source: git/libpolys/coeffs/gnumpc.cc @ 45cc512

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