source: git/libpolys/coeffs/gnumpc.cc @ 560a3d

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