source: git/kernel/shortfl.cc @ 493225

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