source: git/kernel/gnumpfl.cc @ 85e68dd

spielwiese
Last change on this file since 85e68dd was 85e68dd, checked in by Hans Schönemann <hannes@…>, 16 years ago
*hannes: gcc 4.2 git-svn-id: file:///usr/local/Singular/svn/trunk@10634 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 7.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: gnumpfl.cc,v 1.8 2008-03-19 17:44:08 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  if ( exp == 1 )
288  {
289    nNew(u);
290    if ( x == NULL )
291    {
292      gmp_float* n = new gmp_float();
293      *u=(number)n;
294    }
295    else
296    {
297      gmp_float* n = new gmp_float();
298      *n= *(gmp_float*)x;
299      *u=(number)n;
300    }
301    return;
302  }
303  ngfPower(x,exp-1,u);
304
305  gmp_float *n=new gmp_float();
306  *n=*(gmp_float*)x;
307  *(gmp_float*)(*u) *= *(gmp_float*)n;
308  delete (gmp_float*)n;
309}
310
311BOOLEAN ngfIsZero (number a)
312{
313  if ( a == NULL ) return TRUE;
314  return ( ((gmp_float*)a)->isZero() );
315}
316
317/*2
318* za >= 0 ?
319*/
320BOOLEAN ngfGreaterZero (number a)
321{
322  if ( a == NULL ) return TRUE;
323  return ( (*(gmp_float*)a) >= (gmp_float)0.0 );
324}
325
326/*2
327* a > b ?
328*/
329BOOLEAN ngfGreater (number a, number b)
330{
331  if ( a==NULL )
332  {
333    return (((gmp_float*)b)->sign() < 0);
334  }
335  if ( b==NULL )
336  {
337    return (((gmp_float*)a)->sign() > 0);
338  }
339  return ( (*(gmp_float*)a) > (*(gmp_float*)b) );
340}
341
342/*2
343* a = b ?
344*/
345BOOLEAN ngfEqual (number a, number b)
346{
347  if ( a == NULL && b == NULL )
348  {
349    return TRUE;
350  }
351  else if ( a == NULL || b == NULL )
352  {
353    return FALSE;
354  }
355  return ( (*(gmp_float*)a) == (*(gmp_float*)b) );
356}
357
358/*2
359* a == 1 ?
360*/
361BOOLEAN ngfIsOne (number a)
362{
363  if ( a == NULL ) return FALSE;
364  return ((gmp_float*)a)->isOne();
365}
366
367/*2
368* a == -1 ?
369*/
370BOOLEAN ngfIsMOne (number a)
371{
372  if ( a == NULL ) return FALSE;
373  return ((gmp_float*)a)->isMOne();
374}
375
376/*2
377* result =gcd(a,b)
378* dummy, returns 1
379*/
380number ngfGcd(number a, number b)
381{
382  gmp_float *result= new gmp_float( 1 );
383  return (number)result;
384}
385
386static char * ngfEatFloatNExp(char * s )
387{
388  char *start= s;
389
390  // eat floats (mantissa) like:
391  //   0.394394993, 102.203003008,  .300303032, pssibly starting with -
392  if (*s == '-') s++;
393  while ((*s >= '0' && *s <= '9')||(*s == '.')) s++;
394
395  // eat the exponent, starts with 'e' followed by '+', '-'
396  // and digits, like:
397  //   e-202, e+393, accept also E7
398  if ( (s != start) && ((*s == 'e')||(*s=='E')))
399  {
400    if (*s=='E') *s='e';
401    s++; // skip 'e'/'E'
402    if ((*s == '+') || (*s == '-')) s++;
403    while ((*s >= '0' && *s <= '9')) s++;
404  }
405
406  return s;
407}
408
409/*2
410* extracts the number a from s, returns the rest
411*/
412const char * ngfRead (const char * start, number * a)
413{
414  char *s= (char *)start;
415
416  //Print("%s\n",s);
417
418  s= ngfEatFloatNExp( s );
419
420  if (*s=='\0')  // 0
421  {
422    if ( *(gmp_float**)a == NULL ) (*(gmp_float**)a)= new gmp_float();
423    (*(gmp_float**)a)->setFromStr(start);
424  }
425  else if (s==start)  // 1
426  {
427    if ( *(gmp_float**)a != NULL )  delete (*(gmp_float**)a);
428    (*(gmp_float**)a)= new gmp_float(1);
429  }
430  else
431  {
432    gmp_float divisor(1.0);
433    char *start2=s;
434    if ( *s == '/' )
435    {
436      s++;
437      s= ngfEatFloatNExp( (char *)s );
438      if (s!= start2+1)
439      {
440        char tmp_c=*s;
441        *s='\0';
442        divisor.setFromStr(start2+1);
443        *s=tmp_c;
444      }
445      else
446      {
447        Werror("wrong long real format: %s",start2);
448      }
449    }
450    char c=*start2;
451    *start2='\0';
452    if ( *(gmp_float**)a == NULL ) (*(gmp_float**)a)= new gmp_float();
453    (*(gmp_float**)a)->setFromStr(start);
454    *start2=c;
455    if (divisor.isZero())
456    {
457      WerrorS(nDivBy0);
458    }
459    else
460      (**(gmp_float**)a) /= divisor;
461  }
462
463  return s;
464}
465
466/*2
467* write a floating point number
468*/
469void ngfWrite (number &a)
470{
471  char *out;
472  if ( a != NULL )
473  {
474    out= floatToStr(*(gmp_float*)a,gmp_output_digits);
475    StringAppend(out);
476    //omFreeSize((ADDRESS)out, (strlen(out)+1)* sizeof(char) );
477    omFree( (ADDRESS)out );
478  }
479  else
480  {
481    StringAppend("0");
482  }
483}
484
485#ifdef LDEBUG
486//BOOLEAN ngfDBTest(number a, const char *f, const int l)
487//{
488//  return TRUE;
489//}
490#endif
491
492// local Variables: ***
493// folded-file: t ***
494// compile-command: "make installg" ***
495// End: ***
Note: See TracBrowser for help on using the repository browser.