source: git/kernel/shortfl.cc @ 3ad53dd

spielwiese
Last change on this file since 3ad53dd was 85e68dd, checked in by Hans Schönemann <hannes@…>, 16 years ago
*hannes: gcc 4.2 git-svn-id: file:///usr/local/Singular/svn/trunk@10634 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 6.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: shortfl.cc,v 1.6 2008-03-19 17:44:12 Singular Exp $ */
5
6/*
7* ABSTRACT:
8*/
9
10#include <string.h>
11#include "mod2.h"
12#include <mylimits.h>
13#include "structs.h"
14#include "febase.h"
15#include "numbers.h"
16#include "longrat.h"
17#include "mpr_complex.h"
18#include "ring.h"
19#include "shortfl.h"
20
21static float nrEps = 1.0e-3;
22union nf
23{
24  float _f;
25  number _n;
26  nf(float f) {_f = f;}
27  nf(number n) {_n = n;}
28  float F() const {return _f;}
29  number N() const {return _n;}
30};
31
32float nrFloat(number n)
33{
34  return nf(n).F();
35}
36
37BOOLEAN nrGreaterZero (number k)
38{
39  return nf(k).F() >= 0.0;
40}
41
42number nrMult (number a,number b)
43{
44  return nf(nf(a).F() * nf(b).F()).N();
45}
46
47/*2
48* create a number from int
49*/
50number nrInit (int i)
51{
52  float r = (float)i;
53  return nf(nf(r).F()).N();
54}
55
56/*2
57* convert a number to int
58*/
59int nrInt(number &n)
60{
61  int i;
62  float r = nf(n).F();
63  if (((float)INT_MIN <= r) || ((float)MAX_INT_VAL >= r))
64    i = (int)r;
65  else
66    i = 0;
67  return i;
68}
69
70number nrAdd (number a, number b)
71{
72  float x = nf(a).F();
73  float y = nf(b).F();
74  float r = x + y;
75  if (x > 0.0)
76  {
77    if (y < 0.0)
78    {
79      x = r / (x - y);
80      if (x < 0.0)
81        x = -x;
82      if (x < nrEps)
83        r = 0.0;
84    }
85  }
86  else
87  {
88    if (y > 0.0)
89    {
90      x = r / (y - x);
91      if (x < 0.0)
92        x = -x;
93      if (x < nrEps)
94        r = 0.0;
95    }
96  }
97  return nf(r).N();
98}
99
100number nrSub (number a, number b)
101{
102  float x = nf(a).F();
103  float y = nf(b).F();
104  float r = x - y;
105  if (x > 0.0)
106  {
107    if (y > 0.0)
108    {
109      x = r / (x + y);
110      if (x < 0.0)
111        x = -x;
112      if (x < nrEps)
113        r = 0.0;
114    }
115  }
116  else
117  {
118    if (y < 0.0)
119    {
120      x = r / (x + y);
121      if (x < 0.0)
122        x = -x;
123      if (x < nrEps)
124        r = 0.0;
125    }
126  }
127  return nf(r).N();
128}
129
130BOOLEAN nrIsZero (number  a)
131{
132  return (0.0 == nf(a).F());
133}
134
135BOOLEAN nrIsOne (number a)
136{
137  float aa=nf(a).F()-1.0;
138  if (aa<0.0) aa=-aa;
139  return (aa<nrEps);
140}
141
142BOOLEAN nrIsMOne (number a)
143{
144  float aa=nf(a).F()+1.0;
145  if (aa<0.0) aa=-aa;
146  return (aa<nrEps);
147}
148
149number nrDiv (number a,number b)
150{
151  float n = nf(b).F();
152  if (n == 0.0)
153  {
154    WerrorS(nDivBy0);
155    return nf((float)0.0).N();
156  }
157  else
158    return nf(nf(a).F() / n).N();
159}
160
161number  nrInvers (number c)
162{
163  float n = nf(c).F();
164  if (n == 0.0)
165  {
166    WerrorS(nDivBy0);
167    return nf((float)0.0).N();
168  }
169  return nf(1.0 / n).N();
170}
171
172number nrNeg (number c)
173{
174  return nf(-nf(c).F()).N();
175}
176
177BOOLEAN nrGreater (number a,number b)
178{
179  return nf(a).F() > nf(b).F();
180}
181
182BOOLEAN nrEqual (number a,number b)
183{
184  number x = nrSub(a,b);
185  return nf(x).F() == nf((float)0.0).F();
186}
187
188void nrWrite (number &a)
189{
190  StringAppend("%9.3e", nf(a).F());
191}
192
193void nrPower (number a, int i, number * result)
194{
195  if (i==0)
196  {
197    *result = nf(nf(1.0).F()).N();
198    return;
199  }
200  if (i==1)
201  {
202    *result = nf(nf(a).F()).N();
203    return;
204  }
205  nrPower(a,i-1,result);
206  *result = nf(nf(a).F() * nf(*result).F()).N();
207}
208
209static const char* nrEatr(const char *s, float *r)
210{
211  int i;
212
213  if    (*s >= '0' && *s <= '9')
214  {
215    *r = 0.0;
216    do
217    {
218      *r *= 10.0;
219      i = *s++ - '0';
220      *r += (float)i;
221    }
222    while (*s >= '0' && *s <= '9');
223  }
224  else *r = 1.0;
225  return s;
226}
227
228const char * nrRead (const char *s, number *a)
229{
230  const char *t;
231  float z1,z2;
232  float n=1.0;
233
234  s = nrEatr(s, &z1);
235  if (*s == '/')
236  {
237    s++;
238    s = nrEatr(s, &z2);
239    if (z2==0.0)
240      WerrorS(nDivBy0);
241    else
242      z1 /= z2;
243  }
244  else if (*s =='.')
245  {
246    s++;
247    t = s;
248    while (*t >= '0' && *t <= '9')
249    {
250      t++;
251      n *= 10.0;
252    }
253    s = nrEatr(s, &z2);
254    z1 = (z1*n + z2) / n;
255    if (*s=='e')
256    {
257      int e=0; /* exponent */
258      int si=1;/* sign of exponent */
259      s++;
260      if (*s=='+') s++;
261      else if (*s=='-') {s++; si=-1; }
262      while (*s >= '0' && *s <= '9')
263      {
264        e=e*10+(*s)-'0';
265        s++;
266      }
267      if (si==1)
268      {
269        while (e>0) {z1*=10.0; e--; }
270      }
271      else
272      {
273        while (e>0) {z1/=10.0; e--; }
274      }
275    }
276  }
277  *a = nf(z1).N();
278  return s;
279}
280
281/*2
282* the last used charcteristic
283*/
284int nrGetChar()
285{
286  return 0;
287}
288
289#ifdef LDEBUG
290/*2
291* test valid numbers: not implemented yet
292*/
293//BOOLEAN nrDBTest(number a, const char *f, const int l)
294//{
295//  return TRUE;
296//}
297#endif
298
299/* in longrat.h
300typedef MP_INT lint;
301#define SR_INT    1
302#define mpz_size1(A) (ABS((A)->_mp_size))
303*/
304#define SR_HDL(A) ((long)(A))
305#define mpz_isNeg(A) ((A)->_mp_size<0)
306#define mpz_limb_size(A) ((A)->_mp_size)
307#define mpz_limb_d(A) ((A)->_mp_d)
308#define MPZ_DIV(A,B,C) mpz_tdiv_q((A),(B),(C))
309#define IS_INT(A) ((A)->s==3)
310#define IS_IMM(A) (SR_HDL(A)&SR_INT)
311#define GET_NOM(A) &((A)->z)
312#define GET_DENOM(A) &((A)->n)
313#define MPZ_INIT mpz_init
314#define MPZ_CLEAR mpz_clear
315
316number nrMapQ(number from)
317{
318  lint h,*g,*z,*n;
319  int i,j,t,s;
320  float ba,rr,rn,y;
321
322  if (IS_IMM(from))
323    return nf((float)nlInt(from)).N();
324  z=GET_NOM(from);
325  s=0X10000;
326  ba=(float)s;
327  ba*=ba;
328  rr=0.0;
329  i=mpz_size1(z);
330  if(IS_INT(from))
331  {
332    if(i>4)
333    {
334      WerrorS("float overflow");
335      return nf(rr).N();
336    }
337    i--;
338    rr=(float)mpz_limb_d(z)[i];
339    while(i>0)
340    {
341      i--;
342      y=(float)mpz_limb_d(z)[i];
343      rr=rr*ba+y;
344    }
345    if(mpz_isNeg(z))
346      rr=-rr;
347    return nf(rr).N();
348  }
349  n=GET_DENOM(from);
350  j=s=mpz_limb_size(n);
351  if(j>i)
352  {
353    g=n; n=z; z=g;
354    t=j; j=i; i=t;
355  }
356  t=i-j;
357  if(t>4)
358  {
359    if(j==s)
360      WerrorS("float overflow");
361    return nf(rr).N();
362  }
363  if(t>1)
364  {
365    g=&h;
366    MPZ_INIT(g);
367    MPZ_DIV(g,z,n);
368    t=mpz_size1(g);
369    if(t>4)
370    {
371      MPZ_CLEAR(g);
372      if(j==s)
373        WerrorS("float overflow");
374      return nf(rr).N();
375    }
376    t--;
377    rr=(float)mpz_limb_d(g)[t];
378    while(t)
379    {
380      t--;
381      y=(float)mpz_limb_d(g)[t];
382      rr=rr*ba+y;
383    }
384    MPZ_CLEAR(g);
385    if(j!=s)
386      rr=1.0/rr;
387    if(mpz_isNeg(z))
388      rr=-rr;
389    return nf(rr).N();
390  }
391  rn=(float)mpz_limb_d(n)[j-1];
392  rr=(float)mpz_limb_d(z)[i-1];
393  if(j>1)
394  {
395    rn=rn*ba+(float)mpz_limb_d(n)[j-2];
396    rr=rr*ba+(float)mpz_limb_d(z)[i-2];
397    i--;
398  }
399  if(t!=0)
400    rr=rr*ba+(float)mpz_limb_d(z)[i-2];
401  if(j==s)
402    rr=rr/rn;
403  else
404    rr=rn/rr;
405  if(mpz_isNeg(z))
406    rr=-rr;
407  return nf(rr).N();
408}
409
410static number nrMapP(number from)
411{
412  int i = (int)((long)from);
413  float r = (float)i;
414  return nf(r).N();
415}
416
417static number nrMapLongR(number from)
418{
419  float t =(float)mpf_get_d((mpf_srcptr)from);
420  return nf(t).N();
421}
422static number nrMapC(number from)
423{
424  gmp_float h = ((gmp_complex*)from)->real();
425  float t =(float)mpf_get_d((mpf_srcptr)&h);
426  return nf(t).N();
427}
428
429nMapFunc nrSetMap(ring src, ring dst)
430{
431  if (rField_is_Q(src))
432  {
433    return nrMapQ;
434  }
435  if (rField_is_long_R(src))
436  {
437    return nrMapLongR;
438  }
439  if (rField_is_R(src))
440  {
441    return ndCopy;
442  }
443  if(rField_is_Zp(src))
444  {
445    return nrMapP;
446  }
447  if (rField_is_long_C(src))
448  {
449    return nrMapC;
450  }
451  return NULL;
452}
Note: See TracBrowser for help on using the repository browser.