source: git/kernel/gnumpfl.cc @ e9c3b2

spielwiese
Last change on this file since e9c3b2 was b1dfaf, checked in by Frank Seelisch <seelisch@…>, 14 years ago
patch from Kai (checked for problems under Windows OS: no problems) git-svn-id: file:///usr/local/Singular/svn/trunk@13210 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 6.8 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 <kernel/mod2.h>
12#include <kernel/structs.h>
13#include <kernel/febase.h>
14#include <omalloc/omalloc.h>
15#include <kernel/numbers.h>
16#include <kernel/modulop.h>
17#include <kernel/longrat.h>
18
19#include <kernel/gnumpfl.h>
20#include <kernel/mpr_complex.h>
21
22extern size_t gmp_output_digits;
23ring ngfMapRing; // to be used also in gnumpc.cc
24
25static number ngfMapP(number from)
26{
27  return ngfInit(npInt(from,ngfMapRing), currRing);
28}
29number ngfMapQ(number from)
30{
31  gmp_float *res=new gmp_float(numberFieldToFloat(from,QTOF));
32  return (number)res;
33}
34union nf
35{
36  float _f;
37  number _n;
38  nf(float f) {_f = f;}
39  nf(number n) {_n = n;}
40  float F() const {return _f;}
41  number N() const {return _n;}
42};
43static number ngfMapR(number from)
44{
45  gmp_float *res=new gmp_float((double)nf(from).F());
46  return (number)res;
47}
48static number ngfMapC(number from)
49{
50  gmp_float *res=new gmp_float(((gmp_complex*)from)->real());
51  return (number)res;
52}
53
54nMapFunc ngfSetMap(const ring src, const ring dst)
55{
56  if (rField_is_Q(src))
57  {
58    return ngfMapQ;
59  }
60  if (rField_is_long_R(src))
61  {
62    return ngfCopy;
63  }
64  if (rField_is_R(src))
65  {
66    return ngfMapR;
67  }
68  if (rField_is_Zp(src))
69  {
70    ngfMapRing=src;
71    return ngfMapP;
72  }
73  if (rField_is_long_C(src))
74  {
75    return ngfMapC;
76  }
77  return NULL;
78}
79
80/*2
81* n := i
82*/
83number ngfInit (int i, const ring r)
84{
85  gmp_float* n= new gmp_float( (double)i );
86  return (number)n;
87}
88
89/*2
90* convert number to int
91*/
92int ngfInt(number &i, const ring r)
93{
94  double d=(double)*(gmp_float*)i;
95  if (d<0.0)
96    return (int)(d-0.5);
97  else
98    return (int)(d+0.5);
99}
100
101int ngfSize(number n)
102{
103  int i = ngfInt(n, currRing);
104  /* basically return the largest integer in n;
105     only if this happens to be zero although n != 0,
106     return 1;
107     (this code ensures that zero has the size zero) */
108  if ((i == 0) && (ngfIsZero(n) == FALSE)) i = 1;
109  return i;
110}
111
112/*2
113* delete a
114*/
115void ngfDelete (number * a, const ring r)
116{
117  if ( *a != NULL )
118  {
119    delete *(gmp_float**)a;
120    *a=NULL;
121  }
122}
123
124/*2
125* copy a to b
126*/
127number ngfCopy(number a)
128{
129  gmp_float* b= new gmp_float( *(gmp_float*)a );
130  return (number)b;
131}
132
133number ngf_Copy(number a, ring r)
134{
135  gmp_float* b= new gmp_float( *(gmp_float*)a );
136  return (number)b;
137}
138
139/*2
140* za:= - za
141*/
142number ngfNeg (number a)
143{
144  *(gmp_float*)a= -(*(gmp_float*)a);
145  return (number)a;
146}
147
148/*
149* 1/a
150*/
151number ngfInvers(number a)
152{
153  gmp_float* r= NULL;
154  if (((gmp_float*)a)->isZero() )
155  {
156    WerrorS(nDivBy0);
157  }
158  else
159  {
160    r= new gmp_float( (gmp_float)1 / (*(gmp_float*)a) );
161  }
162  return (number)r;
163}
164
165/*2
166* u:= a + b
167*/
168number ngfAdd (number a, number b)
169{
170  gmp_float* r= new gmp_float( (*(gmp_float*)a) + (*(gmp_float*)b) );
171  return (number)r;
172}
173
174/*2
175* u:= a - b
176*/
177number ngfSub (number a, number b)
178{
179  gmp_float* r= new gmp_float( (*(gmp_float*)a) - (*(gmp_float*)b) );
180  return (number)r;
181}
182
183/*2
184* u := a * b
185*/
186number ngfMult (number a, number b)
187{
188  gmp_float* r= new gmp_float( (*(gmp_float*)a) * (*(gmp_float*)b) );
189  return (number)r;
190}
191
192/*2
193* u := a / b
194*/
195number ngfDiv (number a, number b)
196{
197  if ( ((gmp_float*)b)->isZero() )
198  {
199    // a/0 = error
200    WerrorS(nDivBy0);
201    return NULL;
202  }
203  gmp_float* r= new gmp_float( (*(gmp_float*)a) / (*(gmp_float*)b) );
204  return (number)r;
205}
206
207/*2
208* u:= x ^ exp
209*/
210void ngfPower ( number x, int exp, number * u )
211{
212  if ( exp == 0 )
213  {
214    gmp_float* n = new gmp_float(1);
215    *u=(number)n;
216    return;
217  }
218  else if ( ngfIsZero(x) ) // 0^e, e>0
219  {
220    *u=ngfInit(0, currRing);
221    return;
222  }
223  else if ( exp == 1 )
224  {
225    nNew(u);
226    gmp_float* n = new gmp_float();
227    *n= *(gmp_float*)x;
228    *u=(number)n;
229    return;
230  }
231  ngfPower(x,exp-1,u);
232
233  gmp_float *n=new gmp_float();
234  *n=*(gmp_float*)x;
235  *(gmp_float*)(*u) *= *(gmp_float*)n;
236  delete (gmp_float*)n;
237}
238
239BOOLEAN ngfIsZero (number a)
240{
241  return ( ((gmp_float*)a)->isZero() );
242}
243
244/*2
245* za >= 0 ?
246*/
247BOOLEAN ngfGreaterZero (number a)
248{
249  return ( (*(gmp_float*)a) >= (gmp_float)0.0 );
250}
251
252/*2
253* a > b ?
254*/
255BOOLEAN ngfGreater (number a, number b)
256{
257  return ( (*(gmp_float*)a) > (*(gmp_float*)b) );
258}
259
260/*2
261* a = b ?
262*/
263BOOLEAN ngfEqual (number a, number b)
264{
265  return ( (*(gmp_float*)a) == (*(gmp_float*)b) );
266}
267
268/*2
269* a == 1 ?
270*/
271BOOLEAN ngfIsOne (number a)
272{
273  return ((gmp_float*)a)->isOne();
274}
275
276/*2
277* a == -1 ?
278*/
279BOOLEAN ngfIsMOne (number a)
280{
281  return ((gmp_float*)a)->isMOne();
282}
283
284/*2
285* result =gcd(a,b)
286* dummy, returns 1
287*/
288number ngfGcd(number a, number b)
289{
290  gmp_float *result= new gmp_float( 1 );
291  return (number)result;
292}
293
294static char * ngfEatFloatNExp(char * s )
295{
296  char *start= s;
297
298  // eat floats (mantissa) like:
299  //   0.394394993, 102.203003008,  .300303032, pssibly starting with -
300  if (*s == '-') s++;
301  while ((*s >= '0' && *s <= '9')||(*s == '.')) s++;
302
303  // eat the exponent, starts with 'e' followed by '+', '-'
304  // and digits, like:
305  //   e-202, e+393, accept also E7
306  if ( (s != start) && ((*s == 'e')||(*s=='E')))
307  {
308    if (*s=='E') *s='e';
309    s++; // skip 'e'/'E'
310    if ((*s == '+') || (*s == '-')) s++;
311    while ((*s >= '0' && *s <= '9')) s++;
312  }
313
314  return s;
315}
316
317/*2
318* extracts the number a from s, returns the rest
319*/
320const char * ngfRead (const char * start, number * a)
321{
322  char *s= (char *)start;
323
324  //Print("%s\n",s);
325
326  s= ngfEatFloatNExp( s );
327
328  if (*s=='\0')  // 0
329  {
330    if ( *(gmp_float**)a == NULL ) (*(gmp_float**)a)= new gmp_float();
331    (*(gmp_float**)a)->setFromStr(start);
332  }
333  else if (s==start)  // 1
334  {
335    if ( *(gmp_float**)a != NULL )  delete (*(gmp_float**)a);
336    (*(gmp_float**)a)= new gmp_float(1);
337  }
338  else
339  {
340    gmp_float divisor(1.0);
341    char *start2=s;
342    if ( *s == '/' )
343    {
344      s++;
345      s= ngfEatFloatNExp( (char *)s );
346      if (s!= start2+1)
347      {
348        char tmp_c=*s;
349        *s='\0';
350        divisor.setFromStr(start2+1);
351        *s=tmp_c;
352      }
353      else
354      {
355        Werror("wrong long real format: %s",start2);
356      }
357    }
358    char c=*start2;
359    *start2='\0';
360    if ( *(gmp_float**)a == NULL ) (*(gmp_float**)a)= new gmp_float();
361    (*(gmp_float**)a)->setFromStr(start);
362    *start2=c;
363    if (divisor.isZero())
364    {
365      WerrorS(nDivBy0);
366    }
367    else
368      (**(gmp_float**)a) /= divisor;
369  }
370
371  return s;
372}
373
374/*2
375* write a floating point number
376*/
377void ngfWrite (number &a, const ring r)
378{
379  char *out;
380  if ( a != NULL )
381  {
382    out= floatToStr(*(gmp_float*)a,gmp_output_digits);
383    StringAppendS(out);
384    //omFreeSize((ADDRESS)out, (strlen(out)+1)* sizeof(char) );
385    omFree( (ADDRESS)out );
386  }
387  else
388  {
389    StringAppendS("0");
390  }
391}
392
393#ifdef LDEBUG
394//BOOLEAN ngfDBTest(number a, const char *f, const int l)
395//{
396//  return TRUE;
397//}
398#endif
399
400// local Variables: ***
401// folded-file: t ***
402// compile-command: "make installg" ***
403// End: ***
Note: See TracBrowser for help on using the repository browser.