source: git/libpolys/coeffs/gnumpc.cc @ 73a9ffb

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