source: git/libpolys/coeffs/gnumpc.cc @ 2544e7

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