source: git/kernel/shortfl.cc @ 6b9532

spielwiese
Last change on this file since 6b9532 was 35aab3, checked in by Hans Schönemann <hannes@…>, 20 years ago
This commit was generated by cvs2svn to compensate for changes in r6879, which included commits to RCS files with non-trunk default branches. git-svn-id: file:///usr/local/Singular/svn/trunk@6880 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 6.8 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: shortfl.cc,v 1.1.1.1 2003-10-06 12:16:03 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("a/0");
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("1/0");
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
209char* nrEatr(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
228char * nrRead (char *s, number *a)
229{
230  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    z1 /= z2;
240  }
241  else if (*s =='.')
242  {
243    s++;
244    t = s;
245    while (*t >= '0' && *t <= '9')
246    {
247      t++;
248      n *= 10.0;
249    }
250    s = nrEatr(s, &z2);
251    z1 = (z1*n + z2) / n;
252    if (*s=='e')
253    {
254      int e=0; /* exponent */
255      int si=1;/* sign of exponent */
256      s++;
257      if (*s=='+') s++;
258      else if (*s=='-') {s++; si=-1; }
259      while (*s >= '0' && *s <= '9')
260      {
261        e=e*10+(*s)-'0';
262        s++;
263      }
264      if (si==1)
265      {
266        while (e>0) {z1*=10.0; e--; }
267      }
268      else
269      {
270        while (e>0) {z1/=10.0; e--; }
271      }
272    }
273  }
274  *a = nf(z1).N();
275  return s;
276}
277
278/*2
279* the last used charcteristic
280*/
281int nrGetChar()
282{
283  return 0;
284}
285
286#ifdef LDEBUG
287/*2
288* test valid numbers: not implemented yet
289*/
290BOOLEAN nrDBTest(number a, char *f, int l)
291{
292  return TRUE;
293}
294#endif
295
296/* in longrat.h
297typedef MP_INT lint;
298#define SR_INT    1
299#define mpz_size1(A) (ABS((A)->_mp_size))
300*/
301#define SR_HDL(A) ((long)(A))
302#define mpz_isNeg(A) ((A)->_mp_size<0)
303#define mpz_limb_size(A) ((A)->_mp_size)
304#define mpz_limb_d(A) ((A)->_mp_d)
305#define MPZ_DIV(A,B,C) mpz_tdiv_q((A),(B),(C))
306#define IS_INT(A) ((A)->s==3)
307#define IS_IMM(A) (SR_HDL(A)&SR_INT)
308#define GET_NOM(A) &((A)->z)
309#define GET_DENOM(A) &((A)->n)
310#define MPZ_INIT mpz_init
311#define MPZ_CLEAR mpz_clear
312
313static number nrMapQ(number from)
314{
315  lint h,*g,*z,*n;
316  int i,j,t,s;
317  float ba,rr,rn,y;
318
319  if (IS_IMM(from))
320    return nf((float)nlInt(from)).N();
321  z=GET_NOM(from);
322  s=0X10000;
323  ba=(float)s;
324  ba*=ba;
325  rr=0.0;
326  i=mpz_size1(z);
327  if(IS_INT(from))
328  {
329    if(i>4)
330    {
331      WerrorS("float overflow");
332      return nf(rr).N();
333    }
334    i--;
335    rr=(float)mpz_limb_d(z)[i];
336    while(i>0)
337    {
338      i--;
339      y=(float)mpz_limb_d(z)[i];
340      rr=rr*ba+y;
341    }
342    if(mpz_isNeg(z))
343      rr=-rr;
344    return nf(rr).N();
345  }
346  n=GET_DENOM(from);
347  j=s=mpz_limb_size(n);
348  if(j>i)
349  {
350    g=n; n=z; z=g;
351    t=j; j=i; i=t;
352  }
353  t=i-j;
354  if(t>4)
355  {
356    if(j==s)
357      WerrorS("float overflow");
358    return nf(rr).N();
359  }
360  if(t>1)
361  {
362    g=&h;
363    MPZ_INIT(g);
364    MPZ_DIV(g,z,n);
365    t=mpz_size1(g);
366    if(t>4)
367    {
368      MPZ_CLEAR(g);
369      if(j==s)
370        WerrorS("float overflow");
371      return nf(rr).N();
372    }
373    t--;
374    rr=(float)mpz_limb_d(g)[t];
375    while(t)
376    {
377      t--;
378      y=(float)mpz_limb_d(g)[t];
379      rr=rr*ba+y;
380    }
381    MPZ_CLEAR(g);
382    if(j!=s)
383      rr=1.0/rr;
384    if(mpz_isNeg(z))
385      rr=-rr;
386    return nf(rr).N();
387  }
388  rn=(float)mpz_limb_d(n)[j-1];
389  rr=(float)mpz_limb_d(z)[i-1];
390  if(j>1)
391  {
392    rn=rn*ba+(float)mpz_limb_d(n)[j-2];
393    rr=rr*ba+(float)mpz_limb_d(z)[i-2];
394    i--;
395  }
396  if(t!=0)
397    rr=rr*ba+(float)mpz_limb_d(z)[i-2];
398  if(j==s)
399    rr=rr/rn;
400  else
401    rr=rn/rr;
402  if(mpz_isNeg(z))
403    rr=-rr;
404  return nf(rr).N();
405}
406
407static number nrMapP(number from)
408{
409  int i = (int)from;
410  float r = (float)i;
411  return nf(r).N();
412}
413
414static number nrMapLongR(number from)
415{
416  float t =(float)mpf_get_d((mpf_srcptr)from);
417  return nf(t).N();
418}
419static number nrMapC(number from)
420{
421  gmp_float h = ((gmp_complex*)from)->real();
422  float t =(float)mpf_get_d((mpf_srcptr)&h);
423  return nf(t).N();
424}
425
426nMapFunc nrSetMap(ring src, ring dst)
427{
428  if (rField_is_Q(src))
429  {
430    return nrMapQ;
431  }
432  if (rField_is_long_R(src))
433  {
434    return nrMapLongR;
435  }
436  if (rField_is_R(src))
437  {
438    return ndCopy;
439  }
440  if(rField_is_Zp(src))
441  {
442    return nrMapP;
443  }
444  if (rField_is_long_C(src))
445  {
446    return nrMapC;
447  }
448  return NULL;
449}
Note: See TracBrowser for help on using the repository browser.