source: git/kernel/gnumpfl.cc @ ca371d

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