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