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

spielwiese
Last change on this file since 560a3d was cffd3e, checked in by Hans Schoenemann <hannes@…>, 12 years ago
fix: ringlist with real/complex
  • Property mode set to 100644
File size: 10.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: computations with GMP floating-point numbers
6*
7* ngf == number gnu floats
8*/
9
10#include "config.h"
11#include <coeffs/coeffs.h>
12#include <omalloc/omalloc.h>
13#include <reporter/reporter.h>
14#include <coeffs/numbers.h>
15#include <coeffs/modulop.h>
16#include <coeffs/longrat.h>
17#include <coeffs/shortfl.h>
18
19#include <coeffs/gnumpfl.h>
20#include <coeffs/mpr_complex.h>
21
22extern size_t gmp_output_digits;
23//ring ngfMapRing; // to be used also in gnumpc.cc
24
25/// Our Type!
26static const n_coeffType ID = n_long_R;
27
28union nf
29{
30  float _f;
31  number _n;
32  nf(float f) {_f = f;}
33  nf(number n) {_n = n;}
34  float F() const {return _f;}
35  number N() const {return _n;}
36};
37
38/*2
39* n := i
40*/
41number ngfInit (long i, const coeffs r)
42{
43  assume( getCoeffType(r) == ID );
44 
45  gmp_float* n= new gmp_float( (double)i );
46  return (number)n;
47}
48
49/*2
50* convert number to int
51*/
52int ngfInt(number &i, const coeffs r)
53{
54  assume( getCoeffType(r) == ID );
55 
56  double d=(double)*(gmp_float*)i;
57  if (d<0.0)
58    return (int)(d-0.5);
59  else
60    return (int)(d+0.5);
61}
62
63int ngfSize(number n, const coeffs r)
64{
65  int i = ngfInt(n, r);
66  /* basically return the largest integer in n;
67     only if this happens to be zero although n != 0,
68     return 1;
69     (this code ensures that zero has the size zero) */
70  if ((i == 0) && (ngfIsZero(n,r) == FALSE)) i = 1;
71  return i;
72}
73
74/*2
75* delete a
76*/
77void ngfDelete (number * a, const coeffs r)
78{
79  assume( getCoeffType(r) == ID );
80 
81  if ( *a != NULL )
82  {
83    delete *(gmp_float**)a;
84    *a=NULL;
85  }
86}
87
88/*2
89* copy a to b
90*/
91number ngfCopy(number a, const coeffs r)
92{
93  assume( getCoeffType(r) == ID );
94 
95  gmp_float* b= new gmp_float( *(gmp_float*)a );
96  return (number)b;
97}
98
99#if 0
100static number ngfCopyMap(number a, const coeffs r1, const coeffs r2)
101{
102  assume( getCoeffType(r1) == ID );
103  assume( getCoeffType(r2) == ID );
104 
105  gmp_float* b= NULL;
106  if ( a !=  NULL )
107  {
108    b= new gmp_float( *(gmp_float*)a );
109  }
110  return (number)b;
111}
112#endif
113
114/*2
115* za:= - za
116*/
117number ngfNeg (number a, const coeffs r)
118{
119  assume( getCoeffType(r) == ID );
120 
121  *(gmp_float*)a= -(*(gmp_float*)a);
122  return (number)a;
123}
124
125/*
126* 1/a
127*/
128number ngfInvers(number a, const coeffs r)
129{
130  assume( getCoeffType(r) == ID );
131 
132  gmp_float* f= NULL;
133  if (((gmp_float*)a)->isZero() )
134  {
135    WerrorS(nDivBy0);
136  }
137  else
138  {
139    f= new gmp_float( gmp_float(1) / (*(gmp_float*)a) );
140  }
141  return (number)f;
142}
143
144/*2
145* u:= a + b
146*/
147number ngfAdd (number a, number b, const coeffs R)
148{
149  assume( getCoeffType(R) == ID );
150 
151  gmp_float* r= new gmp_float( (*(gmp_float*)a) + (*(gmp_float*)b) );
152  return (number)r;
153}
154
155/*2
156* u:= a - b
157*/
158number ngfSub (number a, number b, const coeffs R)
159{
160  assume( getCoeffType(R) == ID );
161 
162  gmp_float* r= new gmp_float( (*(gmp_float*)a) - (*(gmp_float*)b) );
163  return (number)r;
164}
165
166/*2
167* u := a * b
168*/
169number ngfMult (number a, number b, const coeffs R)
170{
171  assume( getCoeffType(R) == ID );
172 
173  gmp_float* r= new gmp_float( (*(gmp_float*)a) * (*(gmp_float*)b) );
174  return (number)r;
175}
176
177/*2
178* u := a / b
179*/
180number ngfDiv (number a, number b, const coeffs r)
181{
182  assume( getCoeffType(r) == ID );
183 
184  if ( ((gmp_float*)b)->isZero() )
185  {
186    // a/0 = error
187    WerrorS(nDivBy0);
188    return NULL;
189  }
190  gmp_float* f= new gmp_float( (*(gmp_float*)a) / (*(gmp_float*)b) );
191  return (number)f;
192}
193
194/*2
195* u:= x ^ exp
196*/
197number ngfPower (number x, int exp, const coeffs r)
198{
199  assume( getCoeffType(r) == ID );
200 
201  if ( exp == 0 )
202  {
203    gmp_float* n = new gmp_float(1);
204    return (number)n;
205  }
206  else if ( ngfIsZero(x, r) ) // 0^e, e>0
207  {
208    return ngfInit(0, r);
209  }
210  else if ( exp == 1 )
211  {
212    return ngfCopy(x,r);
213  }
214  return (number) ( new gmp_float( (*(gmp_float*)x)^exp ) );
215}
216
217/* kept for compatibility reasons, to be deleted */
218void ngfPower ( number x, int exp, number * u, const coeffs r )
219{
220  *u = ngfPower(x, exp, r);
221} 
222
223BOOLEAN ngfIsZero (number a, const coeffs r)
224{
225  assume( getCoeffType(r) == ID );
226 
227  return ( ((gmp_float*)a)->isZero() );
228}
229
230/*2
231* za > 0 ?
232*/
233BOOLEAN ngfGreaterZero (number a, const coeffs r)
234{
235  assume( getCoeffType(r) == ID );
236 
237  return (((gmp_float*)a)->sign() > 0);
238}
239
240/*2
241* a > b ?
242*/
243BOOLEAN ngfGreater (number a, number b, const coeffs r)
244{
245  assume( getCoeffType(r) == ID );
246 
247  return ( (*(gmp_float*)a) > (*(gmp_float*)b) );
248}
249
250/*2
251* a = b ?
252*/
253BOOLEAN ngfEqual (number a, number b, const coeffs r)
254{
255  assume( getCoeffType(r) == ID );
256 
257  return ( (*(gmp_float*)a) == (*(gmp_float*)b) );
258}
259
260/*2
261* a == 1 ?
262*/
263BOOLEAN ngfIsOne (number a, const coeffs r)
264{
265  assume( getCoeffType(r) == ID );
266 
267  return ((gmp_float*)a)->isOne();
268}
269
270/*2
271* a == -1 ?
272*/
273BOOLEAN ngfIsMOne (number a, const coeffs r)
274{
275  assume( getCoeffType(r) == ID );
276 
277  return ((gmp_float*)a)->isMOne();
278}
279
280static char * ngfEatFloatNExp(char * s )
281{
282  char *start= s;
283
284  // eat floats (mantissa) like:
285  //   0.394394993, 102.203003008,  .300303032, pssibly starting with -
286  if (*s == '-') s++;
287  while ((*s >= '0' && *s <= '9')||(*s == '.')) s++;
288
289  // eat the exponent, starts with 'e' followed by '+', '-'
290  // and digits, like:
291  //   e-202, e+393, accept also E7
292  if ( (s != start) && ((*s == 'e')||(*s=='E')))
293  {
294    if (*s=='E') *s='e';
295    s++; // skip 'e'/'E'
296    if ((*s == '+') || (*s == '-')) s++;
297    while ((*s >= '0' && *s <= '9')) s++;
298  }
299
300  return s;
301}
302
303/*2
304* extracts the number a from s, returns the rest
305*
306* This is also called to print components of complex coefficients.
307* Handle with care!
308*/
309const char * ngfRead (const char * start, number * a, const coeffs r)
310{
311  assume( getCoeffType(r) == ID or getCoeffType(r) == n_long_C);
312 
313  char *s= (char *)start;
314
315  //Print("%s\n",s);
316
317  s= ngfEatFloatNExp( s );
318
319  if (*s=='\0')  // 0
320  {
321    if ( *(gmp_float**)a == NULL ) (*(gmp_float**)a)= new gmp_float();
322    (*(gmp_float**)a)->setFromStr(start);
323  }
324  else if (s==start)  // 1
325  {
326    if ( *(gmp_float**)a != NULL )  delete (*(gmp_float**)a);
327    (*(gmp_float**)a)= new gmp_float(1);
328  }
329  else
330  {
331    gmp_float divisor(1.0);
332    char *start2=s;
333    if ( *s == '/' )
334    {
335      s++;
336      s= ngfEatFloatNExp( (char *)s );
337      if (s!= start2+1)
338      {
339        char tmp_c=*s;
340        *s='\0';
341        divisor.setFromStr(start2+1);
342        *s=tmp_c;
343      }
344      else
345      {
346        Werror("wrong long real format: %s",start2);
347      }
348    }
349    char c=*start2;
350    *start2='\0';
351    if ( *(gmp_float**)a == NULL ) (*(gmp_float**)a)= new gmp_float();
352    (*(gmp_float**)a)->setFromStr(start);
353    *start2=c;
354    if (divisor.isZero())
355    {
356      WerrorS(nDivBy0);
357    }
358    else
359      (**(gmp_float**)a) /= divisor;
360  }
361
362  return s;
363}
364
365/*2
366* write a floating point number
367*/
368void ngfWrite (number &a, const coeffs r)
369{
370  assume( getCoeffType(r) == ID );
371 
372  extern size_t gmp_output_digits;
373  char *out;
374  if ( a != NULL )
375  {
376    out= floatToStr(*(gmp_float*)a, gmp_output_digits);
377    StringAppendS(out);
378    //omFreeSize((void *)out, (strlen(out)+1)* sizeof(char) );
379    omFree( (void *)out );
380  }
381  else
382  {
383    StringAppendS("0");
384  }
385}
386
387BOOLEAN ngfCoeffIsEqual (const coeffs r, n_coeffType n, void * parameter)
388{
389  if (n==ID) {
390    LongComplexInfo* p = (LongComplexInfo *)(parameter);
391    if (p->float_len == r->float_len && p->float_len2 == r->float_len2)
392      return TRUE;
393  }
394  return FALSE;
395}
396
397void ngfSetChar(const coeffs r)
398{
399  setGMPFloatDigits(r->float_len, r->float_len2);
400}
401
402BOOLEAN ngfInitChar(coeffs n, void *parameter)
403{
404  assume( getCoeffType(n) == ID );
405
406  n->cfKillChar = ndKillChar; /* dummy */
407
408  n->cfSetChar = ngfSetChar;
409  n->ch = 0;
410 
411  n->cfDelete  = ngfDelete;
412  n->cfNormalize=ndNormalize;
413  n->cfInit   = ngfInit;
414  n->cfInt    = ngfInt;
415  n->cfAdd     = ngfAdd;
416  n->cfSub     = ngfSub;
417  n->cfMult    = ngfMult;
418  n->cfDiv     = ngfDiv;
419  n->cfExactDiv= ngfDiv;
420  n->cfNeg     = ngfNeg;
421  n->cfInvers  = ngfInvers;
422  n->cfCopy   = ngfCopy;
423  n->cfGreater = ngfGreater;
424  n->cfEqual   = ngfEqual;
425  n->cfIsZero  = ngfIsZero;
426  n->cfIsOne   = ngfIsOne;
427  n->cfIsMOne  = ngfIsMOne;
428  n->cfGreaterZero = ngfGreaterZero;
429  n->cfWriteLong  = ngfWrite;
430  n->cfRead    = ngfRead;
431  n->cfPower   = ngfPower;
432  n->cfSetMap = ngfSetMap;
433  n->cfCoeffWrite = ngfCoeffWrite;
434  n->cfInit_bigint = ngfMapQ;
435#ifdef LDEBUG
436  n->cfDBTest  = ndDBTest; // not yet implemented: ngfDBTest
437#endif
438
439  n->nCoeffIsEqual = ngfCoeffIsEqual;
440
441  if( parameter != NULL)
442  {
443    LongComplexInfo* p = (LongComplexInfo*)parameter;
444     
445    n->float_len = p->float_len;
446    n->float_len2 = p->float_len2;
447  } else // default values, just for testing!
448  {
449    n->float_len = SHORT_REAL_LENGTH;
450    n->float_len2 = SHORT_REAL_LENGTH;
451  }
452   
453  assume( n->float_len <= n->float_len2 );
454  assume( n->float_len2 >= SHORT_REAL_LENGTH );
455 
456  assume( n_NumberOfParameters(n) == 0 );
457  assume( n_ParameterNames(n) == NULL );
458
459  return FALSE;
460}
461
462number ngfMapQ(number from, const coeffs src, const coeffs dst)
463{
464  assume( getCoeffType(dst) == ID );
465  assume( getCoeffType(src) == n_Q );
466 
467  gmp_float *res=new gmp_float(numberFieldToFloat(from,QTOF,dst));
468  return (number)res;
469}
470
471static number ngfMapR(number from, const coeffs src, const coeffs dst)
472{
473  assume( getCoeffType(dst) == ID );
474  assume( getCoeffType(src) == n_R );
475 
476  gmp_float *res=new gmp_float((double)nf(from).F());
477  return (number)res;
478}
479
480static number ngfMapP(number from, const coeffs src, const coeffs dst)
481{
482  assume( getCoeffType(dst) == ID );
483  assume( getCoeffType(src) ==  n_Zp );
484 
485  return ngfInit(npInt(from,src), dst);
486}
487
488static number ngfMapC(number from, const coeffs src, const coeffs dst)
489{
490  assume( getCoeffType(dst) == ID );
491  assume( getCoeffType(src) ==  n_long_C );
492 
493  gmp_float *res=new gmp_float(((gmp_complex*)from)->real());
494  return (number)res;
495}
496
497nMapFunc ngfSetMap(const coeffs src, const coeffs dst)
498{
499  assume( getCoeffType(dst) == ID );
500 
501  if (nCoeff_is_Q(src))
502  {
503    return ngfMapQ;
504  }
505  if (nCoeff_is_long_R(src))
506  {
507    return ndCopyMap; //ngfCopyMap;
508  }
509  if (nCoeff_is_R(src))
510  {
511    return ngfMapR;
512  }
513  if (nCoeff_is_long_C(src))
514  {
515    return ngfMapC;
516  }
517  if (nCoeff_is_Zp(src))
518  {
519    return ngfMapP;
520  }
521  return NULL;
522}
523
524
525void    ngfCoeffWrite  (const coeffs r, BOOLEAN /*details*/)
526{
527  Print("//   characteristic : 0 (real:%d digits, additional %d digits)\n",
528               r->float_len,r->float_len2);  /* long R */
529}
Note: See TracBrowser for help on using the repository browser.