source: git/libpolys/coeffs/gnumpc.cc

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