source: git/libpolys/coeffs/gnumpc.cc @ ce1f78

jengelh-datetimespielwiese
Last change on this file since ce1f78 was ce1f78, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
fixed short/long output of coeffs depending on ShortOut/CanShortOut options of rings Note: the ugly ShortOut correction p_String0Short/p_String0Long is due to Hans
  • Property mode set to 100644
File size: 12.4 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: computations with GMP complex floating-point numbers
7*
8* ngc == number gnu complex
9*/
10
11#include "config.h"
12
13#include <coeffs/coeffs.h>
14#include <coeffs/numbers.h>
15#include <coeffs/longrat.h>
16#include <coeffs/modulop.h>
17#include <coeffs/gnumpc.h>
18#include <coeffs/gnumpfl.h>
19#include <coeffs/mpr_complex.h>
20#include <reporter/reporter.h>
21#include <omalloc/omalloc.h>
22
23
24#include <coeffs/shortfl.h>
25
26/// Our Type!
27static const n_coeffType ID = n_long_C;
28
29
30#ifdef LDEBUG
31// not yet implemented
32BOOLEAN ngcDBTest(number, const char *, const int, const coeffs r)
33{
34  assume( getCoeffType(r) == ID );
35
36  return TRUE;
37}
38#endif
39
40/*
41number   ngcPar(int i, const coeffs r)
42{
43  assume( getCoeffType(r) == ID );
44  assume(i==1);
45
46  gmp_complex* n= new gmp_complex( (long)0, (long)1 );
47  return (number)n;
48}
49*/
50
51/*2
52* n := i
53*/
54number ngcInit (long i, const coeffs r)
55{
56  assume( getCoeffType(r) == ID );
57
58  gmp_complex* n= new gmp_complex( (long)i, (long)0 );
59
60  return (number)n;
61}
62
63/*2
64* convert number to int
65*/
66int ngcInt(number &i, const coeffs r)
67{
68  assume( getCoeffType(r) == ID );
69
70  return (int)((gmp_complex*)i)->real();
71}
72
73int ngcSize(number n, const coeffs R)
74{
75  int r = (int)((gmp_complex*)n)->real();
76  if (r < 0) r = -r;
77  int i = (int)((gmp_complex*)n)->imag();
78  if (i < 0) i = -i;
79  int oneNorm = r + i;
80  /* basically return the 1-norm of n;
81     only if this happens to be zero although n != 0,
82     return 1;
83     (this code ensures that zero has the size zero) */
84  if ((oneNorm == 0.0) & (ngcIsZero(n,R) == FALSE)) oneNorm = 1;
85  return oneNorm;
86}
87
88/*2
89* delete a
90*/
91void ngcDelete (number * a, const coeffs r)
92{
93  assume( getCoeffType(r) == ID );
94
95  if ( *a != NULL )
96  {
97    delete *(gmp_complex**)a;
98    *a=NULL;
99  }
100}
101
102/*2
103 * copy a to b
104*/
105number ngcCopy(number a, const coeffs r)
106{
107  assume( getCoeffType(r) == ID );
108
109  gmp_complex* b= new gmp_complex( *(gmp_complex*)a );
110  return (number)b;
111}
112
113
114/*2
115* za:= - za
116*/
117number ngcNeg (number a, const coeffs R)
118{
119  assume( getCoeffType(R) == ID );
120
121  gmp_complex* r=(gmp_complex*)a;
122  (*r).neg();
123  return (number)a;
124}
125
126/*
127* 1/a
128*/
129number ngcInvers(number a, const coeffs R)
130{
131  assume( getCoeffType(R) == ID );
132
133  gmp_complex* r = NULL;
134  if (((gmp_complex*)a)->isZero())
135  {
136    WerrorS(nDivBy0);
137  }
138  else
139  {
140    r = new gmp_complex( (gmp_complex)1 / (*(gmp_complex*)a) );
141  }
142  return (number)r;
143}
144
145/*2
146* u:= a + b
147*/
148number ngcAdd (number a, number b, const coeffs R)
149{
150  assume( getCoeffType(R) == ID );
151
152  gmp_complex* r= new gmp_complex( (*(gmp_complex*)a) + (*(gmp_complex*)b) );
153  return (number)r;
154}
155
156/*2
157* u:= a - b
158*/
159number ngcSub (number a, number b, const coeffs R)
160{
161  assume( getCoeffType(R) == ID );
162
163  gmp_complex* r= new gmp_complex( (*(gmp_complex*)a) - (*(gmp_complex*)b) );
164  return (number)r;
165}
166
167/*2
168* u := a * b
169*/
170number ngcMult (number a, number b, const coeffs R)
171{
172  assume( getCoeffType(R) == ID );
173
174  gmp_complex* r= new gmp_complex( (*(gmp_complex*)a) * (*(gmp_complex*)b) );
175  return (number)r;
176}
177
178/*2
179* u := a / b
180*/
181number ngcDiv (number a, number b, const coeffs r)
182{
183  assume( getCoeffType(r) == ID );
184
185  if (((gmp_complex*)b)->isZero())
186  {
187    // a/0 = error
188    WerrorS(nDivBy0);
189    return NULL;
190  }
191  gmp_complex* res = new gmp_complex( (*(gmp_complex*)a) / (*(gmp_complex*)b) );
192  return (number)res;
193}
194
195/*2
196* u:= x ^ exp
197*/
198void ngcPower ( number x, int exp, number * u, const coeffs r)
199{
200  assume( getCoeffType(r) == ID );
201
202  if ( exp == 0 )
203  {
204    gmp_complex* n = new gmp_complex(1);
205    *u=(number)n;
206    return;
207  }
208  else if ( exp == 1 )
209  {
210    n_New(u, r);
211    gmp_complex* n = new gmp_complex();
212    *n= *(gmp_complex*)x;
213    *u=(number)n;
214    return;
215  }
216  else if (exp == 2)
217  {
218    n_New(u, r);
219    gmp_complex* n = new gmp_complex();
220    *n= *(gmp_complex*)x;
221    *u=(number)n;
222    *(gmp_complex*)(*u) *= *(gmp_complex*)n;
223    return;
224  }
225  if ( (exp & 1) == 1 )
226  {
227    ngcPower(x,exp-1,u, r);
228    gmp_complex *n = new gmp_complex();
229    *n=*(gmp_complex*)x;
230    *(gmp_complex*)(*u) *= *(gmp_complex*)n;
231    delete n;
232  }
233  else
234  {
235    number w;
236    n_New(&w, r);
237    ngcPower(x,exp/2,&w, r);
238    ngcPower(w,2,u, r);
239    n_Delete(&w, r);
240  }
241}
242
243BOOLEAN ngcIsZero (number a, const coeffs r)
244{
245  assume( getCoeffType(r) == ID );
246
247  return ( ((gmp_complex*)a)->real().isZero() && ((gmp_complex*)a)->imag().isZero());
248}
249
250number ngcRePart(number a, const coeffs r)
251{
252  assume( getCoeffType(r) == ID );
253
254  gmp_complex* n = new gmp_complex(((gmp_complex*)a)->real());
255  return (number)n;
256}
257
258number ngcImPart(number a, const coeffs r)
259{
260  assume( getCoeffType(r) == ID );
261
262  gmp_complex* n = new gmp_complex(((gmp_complex*)a)->imag());
263  return (number)n;
264}
265
266/*2
267* za >= 0 ?
268*/
269BOOLEAN ngcGreaterZero (number a, const coeffs r)
270{
271  assume( getCoeffType(r) == ID );
272
273  if ( ! ((gmp_complex*)a)->imag().isZero() )
274    return ( abs( *(gmp_complex*)a).sign() >= 0 );
275  else
276    return ( ((gmp_complex*)a)->real().sign() >= 0 );
277}
278
279/*2
280* a > b ?
281*/
282BOOLEAN ngcGreater (number a, number b, const coeffs r)
283{
284  assume( getCoeffType(r) == ID );
285
286  gmp_complex *aa=(gmp_complex*)a;
287  gmp_complex *bb=(gmp_complex*)b;
288  return (*aa) > (*bb);
289}
290
291/*2
292* a = b ?
293*/
294BOOLEAN ngcEqual (number a, number b, const coeffs r)
295{
296  assume( getCoeffType(r) == ID );
297
298  gmp_complex *aa=(gmp_complex*)a;
299  gmp_complex *bb=(gmp_complex*)b;
300  return (*aa) == (*bb);
301}
302
303/*2
304* a == 1 ?
305*/
306BOOLEAN ngcIsOne (number a, const coeffs r)
307{
308  assume( getCoeffType(r) == ID );
309
310  return (((gmp_complex*)a)->real().isOne() && ((gmp_complex*)a)->imag().isZero());
311  //return (((gmp_complex*)a)->real().isOne());
312}
313
314/*2
315* a == -1 ?
316*/
317BOOLEAN ngcIsMOne (number a, const coeffs r)
318{
319  assume( getCoeffType(r) == ID );
320
321  return (((gmp_complex*)a)->real().isMOne() && ((gmp_complex*)a)->imag().isZero());
322  //return (((gmp_complex*)a)->real().isMOne());
323}
324
325/*2
326* extracts the number a from s, returns the rest
327*/
328const char * ngcRead (const char * s, number * a, const coeffs r)
329{
330  assume( getCoeffType(r) == ID );
331  assume( r->complex_parameter != NULL );
332
333  if ((*s >= '0') && (*s <= '9'))
334  {
335    gmp_float *re=NULL;
336    s=ngfRead(s,(number *)&re, r);
337    gmp_complex *aa=new gmp_complex(*re);
338    *a=(number)aa;
339    delete re;
340  }
341  else if (strncmp(s, r->complex_parameter,strlen(r->complex_parameter))==0)
342  {
343    s+=strlen(r->complex_parameter);
344    gmp_complex *aa=new gmp_complex((long)0,(long)1);
345    *a=(number)aa;
346  }
347  else
348  {
349    *a=(number) new gmp_complex((long)1);
350  }
351  return s;
352}
353
354
355
356/*2
357* write a floating point number
358*/
359void ngcWrite (number &a, const coeffs r)
360{
361  assume( getCoeffType(r) == ID );
362
363  extern size_t gmp_output_digits; /// comes from mpr_complex.cc
364
365  if (a==NULL)
366    StringAppendS("0");
367  else
368  {
369    char *out;
370    out= complexToStr(*(gmp_complex*)a, gmp_output_digits, r);
371    StringAppendS(out);
372    //    omFreeSize((void *)out, (strlen(out)+1)* sizeof(char) );
373    omFree( (void *)out );
374  }
375}
376
377BOOLEAN ngcCoeffIsEqual (const coeffs r, n_coeffType n, void * parameter)
378{
379  if (n==ID)
380  {
381    LongComplexInfo* p = (LongComplexInfo *)(parameter);
382   
383    if (
384        (p->float_len == r->float_len) &&
385        (p->float_len2 == r->float_len2)
386       )
387      if (strcmp(p->par_name, r->complex_parameter) == 0)
388        return (TRUE);
389  }
390  return (FALSE);
391}
392
393BOOLEAN ngcInitChar(coeffs n, void* parameter)
394{
395  assume( getCoeffType(n) == ID );
396
397  n->cfKillChar = ndKillChar; /* dummy */
398  n->ch = 0;
399
400  n->cfDelete  = ngcDelete;
401  n->cfNormalize=ndNormalize;
402  n->cfInit   = ngcInit;
403  n->cfInt    = ngcInt;
404  n->cfAdd     = ngcAdd;
405  n->cfSub     = ngcSub;
406  n->cfMult    = ngcMult;
407  n->cfDiv     = ngcDiv;
408  n->cfExactDiv= ngcDiv;
409  n->cfNeg     = ngcNeg;
410  n->cfInvers  = ngcInvers;
411  n->cfCopy   = ngcCopy;
412  n->cfGreater = ngcGreater;
413  n->cfEqual   = ngcEqual;
414  n->cfIsZero  = ngcIsZero;
415  n->cfIsOne   = ngcIsOne;
416  n->cfIsMOne  = ngcIsMOne;
417  n->cfGreaterZero = ngcGreaterZero;
418
419  n->cfWriteLong  = ngcWrite;
420  n->cfWriteShort = ngcWrite;
421 
422  n->cfRead    = ngcRead;
423  n->cfPower   = ngcPower;
424  n->cfSetMap = ngcSetMap;
425  n->cfRePart  = ngcRePart;
426  n->cfImPart  = ngcImPart;
427  n->cfCoeffWrite = ngcCoeffWrite;
428    // cfSize  = ndSize;
429#ifdef LDEBUG
430  n->cfDBTest  = ndDBTest; // not yet implemented: ngcDBTest
431#endif
432
433  n->nCoeffIsEqual = ngcCoeffIsEqual;
434
435  n->cfSetChar=ngcSetChar;
436
437  n->cfInit_bigint=ngcMapQ;
438
439/*
440  //r->cfInitChar=nlInitChar;
441  r->cfKillChar=NULL;
442
443  r->cfMult  = nlMult;
444  r->cfSub   = nlSub;
445  r->cfAdd   = nlAdd;
446  r->cfDiv   = nlDiv;
447  r->cfIntDiv= nlIntDiv;
448  r->cfIntMod= nlIntMod;
449  r->cfExactDiv= nlExactDiv;
450  r->cfInit = nlInit;
451  r->cfSize  = nlSize;
452  r->cfInt  = nlInt;
453#ifdef HAVE_RINGS
454  r->cfDivComp = NULL; // only for ring stuff
455  r->cfIsUnit = NULL; // only for ring stuff
456  r->cfGetUnit = NULL; // only for ring stuff
457  r->cfExtGcd = NULL; // only for ring stuff
458#endif
459  r->cfNeg   = nlNeg;
460  r->cfInvers= nlInvers;
461  r->cfCopy  = nl_Copy;
462  r->cfRePart = nl_Copy;
463  r->cfImPart = ndReturn0;
464  r->cfWriteLong = nlWrite;
465  r->cfRead = nlRead;
466  r->cfNormalize=nlNormalize;
467  r->cfGreater = nlGreater;
468#ifdef HAVE_RINGS
469  r->cfDivBy = NULL; // only for ring stuff
470#endif
471  r->cfEqual = nlEqual;
472  r->cfIsZero = nlIsZero;
473  r->cfIsOne = nlIsOne;
474  r->cfIsMOne = nlIsMOne;
475  r->cfGreaterZero = nlGreaterZero;
476  r->cfPower = nlPower;
477  r->cfGetDenom = nlGetDenom;
478  r->cfGetNumerator = nlGetNumerator;
479  r->cfGcd  = nlGcd;
480  r->cfLcm  = nlLcm;
481  r->cfDelete= nlDelete;
482  r->cfSetMap = nlSetMap;
483  r->cfName = ndName;
484  r->cfInpMult=nlInpMult;
485  r->cfInit_bigint=nlCopyMap;
486#ifdef LDEBUG
487  // debug stuff
488  r->cfDBTest=nlDBTest;
489#endif
490
491  // the variables:
492  r->nNULL = INT_TO_SR(0);
493  r->type = n_Q;
494  r->ch = 0;
495  r->has_simple_Alloc=FALSE;
496  r->has_simple_Inverse=FALSE;
497*/
498   
499  if( parameter != NULL)
500  {
501    LongComplexInfo* p = (LongComplexInfo*)parameter;
502    n->complex_parameter = omStrDup(p->par_name);
503    n->float_len = p->float_len;
504    n->float_len2 = p->float_len2;
505   
506  } else // default values, just for testing!
507  {
508    n->complex_parameter = omStrDup("i");
509    n->float_len = SHORT_REAL_LENGTH;
510    n->float_len2 = SHORT_REAL_LENGTH;     
511  }
512   
513  assume( n->float_len <= n->float_len2 );
514  assume( n->float_len2 >= SHORT_REAL_LENGTH );
515  assume( n->complex_parameter != NULL );
516
517  return FALSE;
518}
519
520void ngcSetChar(const coeffs r)
521{
522  setGMPFloatDigits(r->float_len, r->float_len2);
523}
524
525
526
527number ngcMapQ(number from, const coeffs aRing, const coeffs r)
528{
529  assume( getCoeffType(r) == ID );
530  assume( getCoeffType(aRing) == n_Q );
531
532  if ( from != NULL )
533  {
534    gmp_complex *res=new gmp_complex(numberFieldToFloat(from,QTOF,aRing));
535    return (number)res;
536  }
537  else
538    return NULL;
539}
540
541static number ngcMapLongR(number from, const coeffs aRing, const coeffs r)
542{
543  assume( getCoeffType(r) == ID );
544  assume( getCoeffType(aRing) == n_long_R );
545
546  if ( from != NULL )
547  {
548    gmp_complex *res=new gmp_complex(*((gmp_float *)from));
549    return (number)res;
550  }
551  else
552    return NULL;
553}
554
555static number ngcMapR(number from, const coeffs aRing, const coeffs r)
556{
557  assume( getCoeffType(r) == ID );
558  assume( getCoeffType(aRing) == n_R );
559
560  if ( from != NULL )
561  {
562    gmp_complex *res=new gmp_complex((double)nrFloat(from));
563    return (number)res;
564  }
565  else
566    return NULL;
567}
568
569static number ngcMapP(number from, const coeffs aRing, const coeffs r)
570{
571  assume( getCoeffType(r) == ID );
572  assume( getCoeffType(aRing) ==  n_Zp );
573
574  if ( from != NULL )
575    return ngcInit(npInt(from, aRing), r);
576  else
577    return NULL;
578}
579
580static number ngcCopyMap(number from, const coeffs aRing, const coeffs r)
581{
582  assume( getCoeffType(r) == ID );
583  assume( getCoeffType(aRing) ==  ID );
584
585  gmp_complex* b = NULL;
586
587  if ( from !=  NULL )
588  {
589    b = new gmp_complex( *(gmp_complex*)from );
590  }
591  return (number)b;
592}
593
594nMapFunc ngcSetMap(const coeffs src, const coeffs dst)
595{
596  assume( getCoeffType(dst) == ID );
597
598  if (nCoeff_is_Q(src))
599  {
600    return ngcMapQ;
601  }
602  if (nCoeff_is_long_R(src))
603  {
604    return ngcMapLongR;
605  }
606  if (nCoeff_is_long_C(src))
607  {
608    return ngcCopyMap;
609  }
610  if (nCoeff_is_R(src))
611  {
612    return ngcMapR;
613  }
614  if (nCoeff_is_Zp(src))
615  {
616    return ngcMapP;
617  }
618  return NULL;
619}
620
621void    ngcCoeffWrite  (const coeffs r, BOOLEAN details)
622{
623  Print("//   characteristic : 0 (complex:%d digits, additional %d digits)\n",
624        r->float_len, r->float_len2);  /* long C */
625  Print("//   1 parameter    : %s\n", r->complex_parameter);
626  Print("//   minpoly        : (%s^2+1)\n", r->complex_parameter); 
627}
Note: See TracBrowser for help on using the repository browser.