source: git/libpolys/coeffs/gnumpc.cc @ 014b65

fieker-DuValspielwiese
Last change on this file since 014b65 was 014b65, checked in by Mohamed Barakat <mohamed.barakat@…>, 14 years ago
- moved misc,reporter,resources,coeffs,polys -> (new) libpolys (Hans agreed) - migrated to automake in coeffs, misc status: everything builds (except polys) todo: . migrate resources and reporter to automake . create autoconf macros for omalloc, factory, and libpolys
  • Property mode set to 100644
File size: 11.2 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 <coeffs/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 
330  if ((*s >= '0') && (*s <= '9'))
331  {
332    gmp_float *re=NULL;
333    s=ngfRead(s,(number *)&re, r);
334    gmp_complex *aa=new gmp_complex(*re);
335    *a=(number)aa;
336    delete re;
337  }
338  else if (strncmp(s, r->parameter[0],strlen(r->parameter[0]))==0)
339  {
340    s+=strlen(r->parameter[0]);
341    gmp_complex *aa=new gmp_complex((long)0,(long)1);
342    *a=(number)aa;
343  }
344  else
345  {
346    *a=(number) new gmp_complex((long)1);
347  }
348  return s;
349}
350
351
352
353/*2
354* write a floating point number
355*/
356void ngcWrite (number &a, const coeffs r)
357{
358  assume( getCoeffType(r) == ID );
359
360  extern size_t gmp_output_digits; /// comes from mpr_complex.cc
361
362  if (a==NULL)
363    StringAppendS("0");
364  else
365  {
366    char *out;
367    out= complexToStr(*(gmp_complex*)a, gmp_output_digits, r);
368    StringAppendS(out);
369    //    omFreeSize((void *)out, (strlen(out)+1)* sizeof(char) );
370    omFree( (void *)out );
371  }
372}
373
374
375
376/// test, whether r is an instance of nInitCoeffs(n, parameter)
377static BOOLEAN ngcCoeffsEqual(const coeffs r, n_coeffType n, void*)
378{
379  assume( getCoeffType(r) == ID );
380 
381  return (n == ID);
382}
383
384BOOLEAN ngcInitChar(coeffs n, void*)
385{
386  assume( getCoeffType(n) == ID );
387
388  n->cfDelete  = ngcDelete;
389  n->cfNormalize=ndNormalize;
390  n->cfInit   = ngcInit;
391  n->cfInt    = ngcInt;
392  n->cfAdd     = ngcAdd;
393  n->cfSub     = ngcSub;
394  n->cfMult    = ngcMult;
395  n->cfDiv     = ngcDiv;
396  n->cfExactDiv= ngcDiv;
397  n->cfNeg     = ngcNeg;
398  n->cfInvers  = ngcInvers;
399  n->cfCopy   = ngcCopy;
400  n->cfGreater = ngcGreater;
401  n->cfEqual   = ngcEqual;
402  n->cfIsZero  = ngcIsZero;
403  n->cfIsOne   = ngcIsOne;
404  n->cfIsMOne  = ngcIsMOne;
405  n->cfGreaterZero = ngcGreaterZero;
406  n->cfWrite  = ngcWrite;
407  n->cfRead    = ngcRead;
408  n->cfPower   = ngcPower;
409  n->cfSetMap = ngcSetMap;
410  n->cfPar     = ngcPar;
411  n->cfRePart  = ngcRePart;
412  n->cfImPart  = ngcImPart;
413    // cfSize  = ndSize;
414#ifdef LDEBUG
415  n->cfDBTest  = ndDBTest; // not yet implemented: ngcDBTest
416#endif
417
418  n->nCoeffIsEqual = ngcCoeffsEqual;
419
420
421 
422/* 
423  //r->cfInitChar=nlInitChar;
424  r->cfKillChar=NULL;
425  r->cfSetChar=NULL;
426  r->nCoeffIsEqual=nlCoeffsEqual;
427
428  r->cfMult  = nlMult;
429  r->cfSub   = nlSub;
430  r->cfAdd   = nlAdd;
431  r->cfDiv   = nlDiv;
432  r->cfIntDiv= nlIntDiv;
433  r->cfIntMod= nlIntMod;
434  r->cfExactDiv= nlExactDiv;
435  r->cfInit = nlInit;
436  r->cfPar = ndPar;
437  r->cfParDeg = ndParDeg;
438  r->cfSize  = nlSize;
439  r->cfInt  = nlInt;
440#ifdef HAVE_RINGS
441  r->cfDivComp = NULL; // only for ring stuff
442  r->cfIsUnit = NULL; // only for ring stuff
443  r->cfGetUnit = NULL; // only for ring stuff
444  r->cfExtGcd = NULL; // only for ring stuff
445#endif
446  r->cfNeg   = nlNeg;
447  r->cfInvers= nlInvers;
448  r->cfCopy  = nl_Copy;
449  r->cfRePart = nl_Copy;
450  r->cfImPart = ndReturn0;
451  r->cfWrite = nlWrite;
452  r->cfRead = nlRead;
453  r->cfNormalize=nlNormalize;
454  r->cfGreater = nlGreater;
455#ifdef HAVE_RINGS
456  r->cfDivBy = NULL; // only for ring stuff
457#endif
458  r->cfEqual = nlEqual;
459  r->cfIsZero = nlIsZero;
460  r->cfIsOne = nlIsOne;
461  r->cfIsMOne = nlIsMOne;
462  r->cfGreaterZero = nlGreaterZero;
463  r->cfPower = nlPower;
464  r->cfGetDenom = nlGetDenom;
465  r->cfGetNumerator = nlGetNumerator;
466  r->cfGcd  = nlGcd;
467  r->cfLcm  = nlLcm;
468  r->cfDelete= nlDelete;
469  r->cfSetMap = nlSetMap;
470  r->cfName = ndName;
471  r->cfInpMult=nlInpMult;
472  r->cfInit_bigint=nlCopyMap;
473#ifdef LDEBUG
474  // debug stuff
475  r->cfDBTest=nlDBTest;
476#endif
477
478  // the variables:
479  r->nNULL = INT_TO_SR(0);
480  r->type = n_Q;
481  r->ch = 0;
482  r->has_simple_Alloc=FALSE;
483  r->has_simple_Inverse=FALSE;
484*/
485
486/// TODO: Any variables?
487  return FALSE;
488}
489
490
491
492
493
494number ngcMapQ(number from, const coeffs aRing, const coeffs r)
495{
496  assume( getCoeffType(r) == ID );
497  assume( getCoeffType(aRing) == n_Q );
498
499  if ( from != NULL )
500  {
501    gmp_complex *res=new gmp_complex(numberFieldToFloat(from,QTOF,aRing));
502    return (number)res;
503  }
504  else
505    return NULL;
506}
507
508static number ngcMapLongR(number from, const coeffs aRing, const coeffs r)
509{
510  assume( getCoeffType(r) == ID );
511  assume( getCoeffType(aRing) == n_long_R );
512
513  if ( from != NULL )
514  {
515    gmp_complex *res=new gmp_complex(*((gmp_float *)from));
516    return (number)res;
517  }
518  else
519    return NULL;
520}
521
522static number ngcMapR(number from, const coeffs aRing, const coeffs r)
523{
524  assume( getCoeffType(r) == ID );
525  assume( getCoeffType(aRing) == n_R );
526
527  if ( from != NULL )
528  {
529    gmp_complex *res=new gmp_complex((double)nrFloat(from));
530    return (number)res;
531  }
532  else
533    return NULL;
534}
535
536static number ngcMapP(number from, const coeffs aRing, const coeffs r)
537{
538  assume( getCoeffType(r) == ID );
539  assume( getCoeffType(aRing) ==  n_Zp );
540
541  if ( from != NULL )
542    return ngcInit(npInt(from, aRing), r);
543  else
544    return NULL;
545}
546
547static number ngcCopyMap(number from, const coeffs aRing, const coeffs r)
548{
549  assume( getCoeffType(r) == ID );
550  assume( getCoeffType(aRing) ==  ID );
551
552  gmp_complex* b = NULL;
553
554  if ( from !=  NULL )
555  { 
556    b = new gmp_complex( *(gmp_complex*)from );
557  }
558  return (number)b; 
559}
560
561nMapFunc ngcSetMap(const coeffs src, const coeffs dst)
562{
563  assume( getCoeffType(dst) == ID );
564
565  if (nCoeff_is_Q(src))
566  {
567    return ngcMapQ;
568  }
569  if (nCoeff_is_long_R(src))
570  {
571    return ngcMapLongR;
572  }
573  if (nCoeff_is_long_C(src))
574  {
575    return ngcCopyMap;
576  }
577  if (nCoeff_is_R(src))
578  {
579    return ngcMapR;
580  }
581  if (nCoeff_is_Zp(src))
582  {
583    return ngcMapP;
584  }
585  return NULL;
586}
587
Note: See TracBrowser for help on using the repository browser.