source: git/libpolys/coeffs/shortfl.cc @ 20e3062

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