source: git/coeffs/gnumpfl.cc @ 210852

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