source: git/libpolys/coeffs/shortfl.cc @ 8c2815d

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