source: git/libpolys/coeffs/gnumpfl.cc @ 06df101

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