source: git/libpolys/coeffs/shortfl.cc @ e76d8cb

jengelh-datetimespielwiese
Last change on this file since e76d8cb was e76d8cb, checked in by Yue Ren <ren@…>, 10 years ago
fix: nonexistant overflow during numbers->float throwing error
  • Property mode set to 100644
File size: 11.8 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5/*
6* ABSTRACT:
7*/
8#include <iostream>
9#include <coeffs/shortfl.h>
10
11#include <string.h>
12#include <math.h>
13#include <coeffs/coeffs.h>
14#include <coeffs/numbers.h>
15#include <reporter/reporter.h>
16#include <coeffs/numbers.h>
17#include <coeffs/longrat.h>
18#include <coeffs/mpr_complex.h>
19
20#include <misc/mylimits.h>
21
22/// Our Type!
23static const n_coeffType ID = n_R;
24
25static const float nrEps = 1.0e-3;
26
27union nf
28{
29  float _f;
30  number _n;
31
32  nf(float f): _f(f){};
33
34  nf(number n): _n(n){};
35
36  inline float F() const {return _f;}
37  inline number N() const {return _n;}
38};
39
40
41
42
43float nrFloat(number n)
44{
45  return nf(n).F();
46}
47
48
49void    nrCoeffWrite  (const coeffs r, BOOLEAN /*details*/)
50{
51  assume( getCoeffType(r) == ID );
52  PrintS("//   characteristic : 0 (real)\n");  /* R */
53}
54
55
56BOOLEAN nrGreaterZero (number k, const coeffs r)
57{
58  assume( getCoeffType(r) == ID );
59
60  return nf(k).F() >= 0.0;
61}
62
63number nrMult (number a,number b, const coeffs r)
64{
65  assume( getCoeffType(r) == ID );
66
67  return nf(nf(a).F() * nf(b).F()).N();
68}
69
70/*2
71* create a number from int
72*/
73number nrInit (long i, const coeffs r)
74{
75  assume( getCoeffType(r) == ID );
76
77  float f = (float)i;
78  return nf(nf(f).F()).N();
79}
80
81/*2
82* convert a number to int
83*/
84int nrInt(number &n, const coeffs r)
85{
86  assume( getCoeffType(r) == ID );
87
88  int i;
89  float f = nf(n).F();
90  if (((float)(-MAX_INT_VAL-1) <= f) || ((float)MAX_INT_VAL >= f))
91    i = (int)f;
92  else
93    i = 0;
94  return i;
95}
96
97int nrSize(number n, const coeffs)
98{
99  float f = nf(n).F();
100  int i = (int)f;
101  /* basically return the largest integer in n;
102     only if this happens to be zero although n != 0,
103     return 1;
104     (this code ensures that zero has the size zero) */
105  if ((f != 0.0) & (i == 0)) i = 1;
106  return i;
107}
108
109number nrAdd (number a, number b, const coeffs r)
110{
111  assume( getCoeffType(r) == ID );
112
113  float x = nf(a).F();
114  float y = nf(b).F();
115  float f = x + y;
116  if (x > 0.0)
117  {
118    if (y < 0.0)
119    {
120      x = f / (x - y);
121      if (x < 0.0)
122        x = -x;
123      if (x < nrEps)
124        f = 0.0;
125    }
126  }
127  else
128  {
129    if (y > 0.0)
130    {
131      x = f / (y - x);
132      if (x < 0.0)
133        x = -x;
134      if (x < nrEps)
135        f = 0.0;
136    }
137  }
138  return nf(f).N();
139}
140
141number nrSub (number a, number b, const coeffs r)
142{
143  assume( getCoeffType(r) == ID );
144
145  float x = nf(a).F();
146  float y = nf(b).F();
147  float f = x - y;
148  if (x > 0.0)
149  {
150    if (y > 0.0)
151    {
152      x = f / (x + y);
153      if (x < 0.0)
154        x = -x;
155      if (x < nrEps)
156        f = 0.0;
157    }
158  }
159  else
160  {
161    if (y < 0.0)
162    {
163      x = f / (x + y);
164      if (x < 0.0)
165        x = -x;
166      if (x < nrEps)
167        f = 0.0;
168    }
169  }
170  return nf(f).N();
171}
172
173BOOLEAN nrIsZero (number  a, const coeffs r)
174{
175  assume( getCoeffType(r) == ID );
176
177  return (0.0 == nf(a).F());
178}
179
180BOOLEAN nrIsOne (number a, const coeffs r)
181{
182  assume( getCoeffType(r) == ID );
183
184  float aa=nf(a).F()-1.0;
185  if (aa<0.0) aa=-aa;
186  return (aa<nrEps);
187}
188
189BOOLEAN nrIsMOne (number a, const coeffs r)
190{
191  assume( getCoeffType(r) == ID );
192
193  float aa=nf(a).F()+1.0;
194  if (aa<0.0) aa=-aa;
195  return (aa<nrEps);
196}
197
198number nrDiv (number a,number b, const coeffs r)
199{
200  assume( getCoeffType(r) == ID );
201
202  float n = nf(b).F();
203  if (n == 0.0)
204  {
205    WerrorS(nDivBy0);
206    return nf((float)0.0).N();
207  }
208  else
209    return nf(nf(a).F() / n).N();
210}
211
212number  nrInvers (number c, const coeffs r)
213{
214  assume( getCoeffType(r) == ID );
215
216  float n = nf(c).F();
217  if (n == 0.0)
218  {
219    WerrorS(nDivBy0);
220    return nf((float)0.0).N();
221  }
222  return nf(1.0 / n).N();
223}
224
225number nrNeg (number c, const coeffs r)
226{
227  assume( getCoeffType(r) == ID );
228
229  return nf(-nf(c).F()).N();
230}
231
232BOOLEAN nrGreater (number a,number b, const coeffs r)
233{
234  assume( getCoeffType(r) == ID );
235
236  return nf(a).F() > nf(b).F();
237}
238
239BOOLEAN nrEqual (number a,number b, const coeffs r)
240{
241  assume( getCoeffType(r) == ID );
242
243  number x = nrSub(a,b,r);
244  return nf(x).F() == nf((float)0.0).F();
245}
246
247void nrWrite (number &a, const coeffs r)
248{
249  assume( getCoeffType(r) == ID );
250
251  StringAppend("%9.3e", nf(a).F());
252}
253
254void nrPower (number a, int i, number * result, const coeffs r)
255{
256  assume( getCoeffType(r) == ID );
257
258  if (i==0)
259  {
260    *result = nf(nf(1.0).F()).N();
261    return;
262  }
263  if (i==1)
264  {
265    *result = nf(nf(a).F()).N();
266    return;
267  }
268  nrPower(a,i-1,result,r);
269  *result = nf(nf(a).F() * nf(*result).F()).N();
270}
271
272namespace {
273  static const char* nrEatr(const char *s, float *r)
274  {
275    int i;
276
277    if    (*s >= '0' && *s <= '9')
278    {
279      *r = 0.0;
280      do
281      {
282        *r *= 10.0;
283        i = *s++ - '0';
284        *r += (float)i;
285      }
286      while (*s >= '0' && *s <= '9');
287    }
288    else *r = 1.0;
289    return s;
290  }
291}
292
293const char * nrRead (const char *s, number *a, const coeffs r)
294{
295
296  assume( getCoeffType(r) == ID );
297
298  static const char *nIllegalChar="illegal character in number";
299
300  const char *t;
301  const char *start=s;
302  float z1,z2;
303  float n=1.0;
304
305  s = nrEatr(s, &z1);
306  if (*s == '/')
307  {
308    if (s==start) { WerrorS(nIllegalChar);return s; }
309    s++;
310    s = nrEatr(s, &z2);
311    if (z2==0.0)
312      WerrorS(nDivBy0);
313    else
314      z1 /= z2;
315  }
316  else if (*s =='.')
317  {
318    if (s==start) { WerrorS(nIllegalChar);return s; }
319    s++;
320    t = s;
321    while (*t >= '0' && *t <= '9')
322    {
323      t++;
324      n *= 10.0;
325    }
326    s = nrEatr(s, &z2);
327    z1 = (z1*n + z2) / n;
328    if (*s=='e')
329    {
330      int e=0; /* exponent */
331      int si=1;/* sign of exponent */
332      s++;
333      if (*s=='+') s++;
334      else if (*s=='-') {s++; si=-1; }
335      while (*s >= '0' && *s <= '9')
336      {
337        e=e*10+(*s)-'0';
338        s++;
339      }
340      if (si==1)
341      {
342        while (e>0) {z1*=10.0; e--; }
343      }
344      else
345      {
346        while (e>0) {z1/=10.0; e--; }
347      }
348    }
349  }
350  *a = nf(z1).N();
351  return s;
352}
353
354
355// the last used charcteristic
356// int nrGetChar(){ return 0; }
357
358
359#ifdef LDEBUG
360/*2
361* test valid numbers: not implemented yet
362*/
363#pragma GCC diagnostic ignored "-Wunused-parameter"
364BOOLEAN  nrDBTest(number a, const char *f, const int l, const coeffs r)
365{
366  assume( getCoeffType(r) == ID );
367
368  return TRUE;
369}
370#endif
371
372static number nrMapP(number from, const coeffs aRing, const coeffs r)
373{
374  assume( getCoeffType(r) == ID );
375  assume( getCoeffType(aRing) ==  n_Zp );
376
377  int i = (int)((long)from);
378  float f = (float)i;
379  return nf(f).N();
380}
381
382static number nrMapLongR(number from, const coeffs aRing, const coeffs r)
383{
384  assume( getCoeffType(r) == ID );
385  assume( getCoeffType(aRing) == n_long_R );
386
387  float t =(float)mpf_get_d((mpf_srcptr)from);
388  return nf(t).N();
389}
390
391static number nrMapC(number from, const coeffs aRing, const coeffs r)
392{
393  assume( getCoeffType(r) == ID );
394  assume( getCoeffType(aRing) == n_long_C );
395
396  gmp_float h = ((gmp_complex*)from)->real();
397  float t =(float)mpf_get_d((mpf_srcptr)&h);
398  return nf(t).N();
399}
400
401
402number nrMapQ(number from, const coeffs aRing, const coeffs r)
403{
404/* in longrat.h
405#define SR_INT    1
406#define mpz_size1(A) (ABS((A)->_mp_size))
407*/
408#define SR_HDL(A) ((long)(A))
409#define IS_INT(A) ((A)->s==3)
410#define IS_IMM(A) (SR_HDL(A)&SR_INT)
411#define GET_NOM(A) ((A)->z)
412#define GET_DENOM(A) ((A)->n)
413
414  assume( getCoeffType(r) == ID );
415  assume( getCoeffType(aRing) == n_Q );
416
417  if (IS_IMM(from))
418    return nf((float)nlInt(from,NULL /* dummy for nlInt*/)).N();
419
420  /* read out the enumerator */
421  mpz_ptr z=GET_NOM(from);
422  int i = mpz_size1(z);
423  mpf_t e;
424  mpf_init(e);
425  mpf_set_z(e,z);
426
427  /* if number was an integer, we are done*/
428  if(IS_INT(from))
429  {
430    if(i>4)
431    {
432      WerrorS("float overflow");
433      return nf(0.0).N();
434    }
435    double basis;
436    signed long int exp;
437    basis = mpf_get_d_2exp(&exp, e);
438    float f = ldexp(basis,exp);
439    mpf_clear(e);
440    return nf(f).N();
441  }
442
443  /* else read out the denominator */
444  mpz_ptr n = GET_DENOM(from);
445  int j = mpz_size1(n);
446  if(j-i>4)
447  {
448    WerrorS("float overflow");
449    mpf_clear(e);
450    return nf(0.0).N();
451  }
452  mpf_t d;
453  mpf_init(d);
454  mpf_set_z(d,n);
455
456  /* and compute the quotient */
457  mpf_t q;
458  mpf_init(q);
459  mpf_div(q,e,d);
460
461  double basis;
462  signed long int exp;
463  basis = mpf_get_d_2exp(&exp, q);
464  float f = ldexp(basis,exp);
465  mpf_clear(e);
466  mpf_clear(d);
467  mpf_clear(q);
468  return nf(f).N();
469}
470
471// old version:
472// number nrMapQ(number from, const coeffs aRing, const coeffs r)
473// {
474// /* in longrat.h
475// #define SR_INT    1
476// #define mpz_size1(A) (ABS((A)->_mp_size))
477// */
478// #define SR_HDL(A) ((long)(A))
479// #define mpz_isNeg(A) ((A)->_mp_size<0)
480// #define mpz_limb_size(A) ((A)->_mp_size)
481// #define mpz_limb_d(A) ((A)->_mp_d)
482// #define MPZ_DIV(A,B,C) mpz_tdiv_q((A),(B),(C))
483// #define IS_INT(A) ((A)->s==3)
484// #define IS_IMM(A) (SR_HDL(A)&SR_INT)
485// #define GET_NOM(A) ((A)->z)
486// #define GET_DENOM(A) ((A)->n)
487// #define MPZ_INIT mpz_init
488// #define MPZ_CLEAR mpz_clear
489
490//   assume( getCoeffType(r) == ID );
491//   assume( getCoeffType(aRing) == n_Q );
492
493//   mpz_t h;
494//   mpz_ptr g,z,n;
495//   int i,j,t,s;
496//   float ba,rr,rn,y;
497
498//   if (IS_IMM(from))
499//     return nf((float)nlInt(from,NULL /* dummy for nlInt*/)).N();
500//   z=GET_NOM(from);
501//   s=0X10000;
502//   ba=(float)s;
503//   ba*=ba;
504//   rr=0.0;
505//   i=mpz_size1(z);
506//   if(IS_INT(from))
507//   {
508//     if(i>4)
509//     {
510//       WerrorS("float overflow");
511//       return nf(rr).N();
512//     }
513//     i--;
514//     rr=(float)mpz_limb_d(z)[i];
515//     while(i>0)
516//     {
517//       i--;
518//       y=(float)mpz_limb_d(z)[i];
519//       rr=rr*ba+y;
520//     }
521//     if(mpz_isNeg(z))
522//       rr=-rr;
523//     return nf(rr).N();
524//   }
525//   n=GET_DENOM(from);
526//   j=s=mpz_limb_size(n);
527//   if(j>i)
528//   {
529//     g=n; n=z; z=g;
530//     t=j; j=i; i=t;
531//   }
532//   t=i-j;
533//   if(t>4)
534//   {
535//     if(j==s)
536//       WerrorS("float overflow");
537//     return nf(rr).N();
538//   }
539//   if(t>1)
540//   {
541//     g=h;
542//     MPZ_INIT(g);
543//     MPZ_DIV(g,z,n);
544//     t=mpz_size1(g);
545//     if(t>4)
546//     {
547//       MPZ_CLEAR(g);
548//       if(j==s)
549//         WerrorS("float overflow");
550//       return nf(rr).N();
551//     }
552//     t--;
553//     rr=(float)mpz_limb_d(g)[t];
554//     while(t)
555//     {
556//       t--;
557//       y=(float)mpz_limb_d(g)[t];
558//       rr=rr*ba+y;
559//     }
560//     MPZ_CLEAR(g);
561//     if(j!=s)
562//       rr=1.0/rr;
563//     if(mpz_isNeg(z))
564//       rr=-rr;
565//     return nf(rr).N();
566//   }
567//   rn=(float)mpz_limb_d(n)[j-1];
568//   rr=(float)mpz_limb_d(z)[i-1];
569//   if(j>1)
570//   {
571//     rn=rn*ba+(float)mpz_limb_d(n)[j-2];
572//     rr=rr*ba+(float)mpz_limb_d(z)[i-2];
573//     i--;
574//   }
575//   if(t!=0)
576//     rr=rr*ba+(float)mpz_limb_d(z)[i-2];
577//   if(j==s)
578//     rr=rr/rn;
579//   else
580//     rr=rn/rr;
581//   if(mpz_isNeg(z))
582//     rr=-rr;
583//   return nf(rr).N();
584// }
585
586nMapFunc nrSetMap(const coeffs src, const coeffs dst)
587{
588  assume( getCoeffType(dst) == ID );
589
590  if (nCoeff_is_Q(src))
591  {
592    return nrMapQ;
593  }
594  if (nCoeff_is_long_R(src))
595  {
596    return nrMapLongR;
597  }
598  if (nCoeff_is_R(src))
599  {
600    return ndCopyMap;
601  }
602  if(nCoeff_is_Zp(src))
603  {
604    return nrMapP;
605  }
606  if (nCoeff_is_long_C(src))
607  {
608    return nrMapC;
609  }
610  return NULL;
611}
612
613BOOLEAN nrInitChar(coeffs n, void* p)
614{
615  assume( getCoeffType(n) == ID );
616
617  assume( p == NULL );
618
619  n->cfKillChar = ndKillChar; /* dummy */
620  n->ch = 0;
621
622  n->cfInit = nrInit;
623  n->cfInt  = nrInt;
624  n->cfAdd   = nrAdd;
625  n->cfSub   = nrSub;
626  n->cfMult  = nrMult;
627  n->cfDiv   = nrDiv;
628  n->cfExactDiv= nrDiv;
629  n->cfNeg   = nrNeg;
630  n->cfInvers= nrInvers;
631  n->cfCopy  = ndCopy;
632  n->cfGreater = nrGreater;
633  n->cfEqual = nrEqual;
634  n->cfIsZero = nrIsZero;
635  n->cfIsOne = nrIsOne;
636  n->cfIsMOne = nrIsMOne;
637  n->cfGreaterZero = nrGreaterZero;
638  n->cfWriteLong = nrWrite;
639  n->cfRead = nrRead;
640  n->cfPower = nrPower;
641  n->cfSetMap = nrSetMap;
642  n->cfCoeffWrite  = nrCoeffWrite;
643  n->cfInit_bigint = nrMapQ;
644
645    /* nName= ndName; */
646    /*nSize  = ndSize;*/
647#ifdef LDEBUG
648  n->cfDBTest=ndDBTest; // not yet implemented: nrDBTest;
649#endif
650
651  n->nCoeffIsEqual = ndCoeffIsEqual;
652
653  n->float_len = SHORT_REAL_LENGTH;
654  n->float_len2 = SHORT_REAL_LENGTH;
655
656  // TODO: Any variables?
657  return FALSE;
658}
Note: See TracBrowser for help on using the repository browser.