source: git/libpolys/coeffs/gnumpc.cc @ 61b2e16

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