My Project
Loading...
Searching...
No Matches
gnumpfl.cc
Go to the documentation of this file.
1/****************************************
2* Computer Algebra System SINGULAR *
3****************************************/
4/*
5* ABSTRACT: computations with GMP floating-point numbers
6*
7* ngf == number gnu floats
8*/
9
10#include "misc/auxiliary.h"
11
12#include "reporter/reporter.h"
13
14#include "coeffs/coeffs.h"
15#include "coeffs/numbers.h"
16#include "coeffs/mpr_complex.h"
17
18#include "coeffs/longrat.h"
19#include "coeffs/shortfl.h"
20#include "coeffs/gnumpfl.h"
21#include "coeffs/gnumpc.h"
22#include "coeffs/modulop.h"
23
24const char * ngfRead (const char *s, number *a, const coeffs r);
25
26union nf
27{
29 number _n;
31 nf(number n) {_n = n;}
32 SI_FLOAT F() const {return _f;}
33 number N() const {return _n;}
34};
35
36/*2
37* n := i
38*/
39static number ngfInit (long i, const coeffs r)
40{
42
43 gmp_float* n= new gmp_float( (double)i );
44 return (number)n;
45}
46
47/*2
48* convert number to int
49*/
50static long ngfInt(number &i, const coeffs r)
51{
53
54 double d=(double)*(gmp_float*)i;
55 if (d<0.0)
56 return (long)(d-0.5);
57 else
58 return (long)(d+0.5);
59}
60
61static BOOLEAN ngfIsZero (number a, const coeffs r)
62{
64
65 return ( ((gmp_float*)a)->isZero() );
66}
67
68static int ngfSize(number n, const coeffs r)
69{
70 long i = ngfInt(n, r);
71 /* basically return the largest integer in n;
72 only if this happens to be zero although n != 0,
73 return 1;
74 (this code ensures that zero has the size zero) */
75 if ((i == 0) && (ngfIsZero(n,r) == FALSE)) i = 1;
76 return ABS(i);
77}
78
79/*2
80* delete a
81*/
82static void ngfDelete (number * a, const coeffs r)
83{
85
86 if ( *a != NULL )
87 {
88 delete *(gmp_float**)a;
89 *a=NULL;
90 }
91}
92
93/*2
94* copy a to b
95*/
96static number ngfCopy(number a, const coeffs r)
97{
99
100 gmp_float* b= new gmp_float( *(gmp_float*)a );
101 return (number)b;
102}
103
104#if 0
105static number ngfCopyMap(number a, const coeffs r1, const coeffs r2)
106{
107 assume( getCoeffType(r1) == n_long_R );
108 assume( getCoeffType(r2) == n_long_R );
109
110 gmp_float* b= NULL;
111 if ( a != NULL )
112 {
113 b= new gmp_float( *(gmp_float*)a );
114 }
115 return (number)b;
116}
117#endif
118
119/*2
120* za:= - za
121*/
122static number ngfNeg (number a, const coeffs r)
123{
125
126 *(gmp_float*)a= -(*(gmp_float*)a);
127 return (number)a;
128}
129
130/*
131* 1/a
132*/
133static number ngfInvers(number a, const coeffs r)
134{
136
137 gmp_float* f= NULL;
138 if (((gmp_float*)a)->isZero() )
139 {
141 f= new gmp_float( 0 );
142 }
143 else
144 {
145 f= new gmp_float( gmp_float(1) / (*(gmp_float*)a) );
146 }
147 return (number)f;
148}
149
150/*2
151* u:= a + b
152*/
153static number ngfAdd (number a, number b, const coeffs R)
154{
156
157 gmp_float* r= new gmp_float( (*(gmp_float*)a) + (*(gmp_float*)b) );
158 return (number)r;
159}
160
161static void ngfInpAdd (number &a, number b, const coeffs R)
162{
164
165 (*(gmp_float*)a) += (*(gmp_float*)b);
166}
167
168/*2
169* u:= a - b
170*/
171static number ngfSub (number a, number b, const coeffs R)
172{
174
175 gmp_float* r= new gmp_float( (*(gmp_float*)a) - (*(gmp_float*)b) );
176 return (number)r;
177}
178
179/*2
180* u := a * b
181*/
182static number ngfMult (number a, number b, const coeffs R)
183{
185
186 gmp_float* r= new gmp_float( (*(gmp_float*)a) * (*(gmp_float*)b) );
187 return (number)r;
188}
189
190static void ngfInpMult (number &a, number b, const coeffs R)
191{
193
194 (*(gmp_float*)a) *= (*(gmp_float*)b);
195}
196
197/*2
198* u := a / b
199*/
200static number ngfDiv (number a, number b, const coeffs r)
201{
203
204 gmp_float* f;
205 if ( ((gmp_float*)b)->isZero() )
206 {
207 // a/0 = error
209 f= new gmp_float( 0 );
210 }
211 else
212 {
213 f= new gmp_float( (*(gmp_float*)a) / (*(gmp_float*)b) );
214 }
215 return (number)f;
216}
217
218/*2
219* u:= x ^ exp
220*/
221static number ngfPower (number x, int exp, const coeffs r)
222{
224
225 if ( exp == 0 )
226 {
227 gmp_float* n = new gmp_float(1);
228 return (number)n;
229 }
230 else if ( ngfIsZero(x, r) ) // 0^e, e>0
231 {
232 return ngfInit(0, r);
233 }
234 else if ( exp == 1 )
235 {
236 return ngfCopy(x,r);
237 }
238 return (number) ( new gmp_float( (*(gmp_float*)x)^exp ) );
239}
240
241/* kept for compatibility reasons, to be deleted */
242static void ngfPower ( number x, int exp, number * u, const coeffs r )
243{
244 *u = ngfPower(x, exp, r);
245}
246
247/*2
248* za > 0 ?
249*/
250static BOOLEAN ngfGreaterZero (number a, const coeffs r)
251{
253
254 return (((gmp_float*)a)->sign() > 0);
255}
256
257/*2
258* a > b ?
259*/
260static BOOLEAN ngfGreater (number a, number b, const coeffs r)
261{
263
264 return ( (*(gmp_float*)a) > (*(gmp_float*)b) );
265}
266
267/*2
268* a = b ?
269*/
270static BOOLEAN ngfEqual (number a, number b, const coeffs r)
271{
273
274 return ( (*(gmp_float*)a) == (*(gmp_float*)b) );
275}
276
277/*2
278* a == 1 ?
279*/
280static BOOLEAN ngfIsOne (number a, const coeffs r)
281{
283
284 return ((gmp_float*)a)->isOne();
285}
286
287/*2
288* a == -1 ?
289*/
290static BOOLEAN ngfIsMOne (number a, const coeffs r)
291{
293
294 return ((gmp_float*)a)->isMOne();
295}
296
297static char * ngfEatFloatNExp(char * s )
298{
299 char *start= s;
300
301 // eat floats (mantissa) like:
302 // 0.394394993, 102.203003008, .300303032, pssibly starting with -
303 if (*s == '-') s++;
304 while ((*s >= '0' && *s <= '9')||(*s == '.')) s++;
305
306 // eat the exponent, starts with 'e' followed by '+', '-'
307 // and digits, like:
308 // e-202, e+393, accept also E7
309 if ( (s != start) && ((*s == 'e')||(*s=='E')))
310 {
311 if (*s=='E') *s='e';
312 s++; // skip 'e'/'E'
313 if ((*s == '+') || (*s == '-')) s++;
314 while ((*s >= '0' && *s <= '9')) s++;
315 }
316
317 return s;
318}
319
320/*2
321* extracts the number a from s, returns the rest
322*
323* This is also called to print components of complex coefficients.
324* Handle with care!
325*/
326const char * ngfRead (const char * start, number * a, const coeffs r)
327{
329
330 char *s= (char *)start;
331
332 //Print("%s\n",s);
333
334 s= ngfEatFloatNExp( s );
335
336 if (*s=='\0') // 0
337 {
338 if ( *(gmp_float**)a == NULL ) (*(gmp_float**)a)= new gmp_float();
339 (*(gmp_float**)a)->setFromStr(start);
340 }
341 else if (s==start) // 1
342 {
343 if ( *(gmp_float**)a != NULL ) delete (*(gmp_float**)a);
344 (*(gmp_float**)a)= new gmp_float(1);
345 }
346 else
347 {
348 gmp_float divisor(1.0);
349 char *start2=s;
350 if ( *s == '/' )
351 {
352 s++;
353 s= ngfEatFloatNExp( (char *)s );
354 if (s!= start2+1)
355 {
356 char tmp_c=*s;
357 *s='\0';
358 divisor.setFromStr(start2+1);
359 *s=tmp_c;
360 }
361 else
362 {
363 Werror("wrong long real format: %s",start2);
364 }
365 }
366 char c=*start2;
367 *start2='\0';
368 if ( *(gmp_float**)a == NULL ) (*(gmp_float**)a)= new gmp_float();
369 (*(gmp_float**)a)->setFromStr(start);
370 *start2=c;
371 if (divisor.isZero())
372 {
374 }
375 else
376 (**(gmp_float**)a) /= divisor;
377 }
378
379 return s;
380}
381
382/*2
383* write a floating point number
384*/
385static void ngfWrite (number a, const coeffs r)
386{
388
389 char *out;
390 if ( a != NULL )
391 {
392 out= floatToStr(*(gmp_float*)a, r->float_len);
393 StringAppendS(out);
394 //omFreeSize((void *)out, (strlen(out)+1)* sizeof(char) );
395 omFree( (void *)out );
396 }
397 else
398 {
399 StringAppendS("0");
400 }
401}
402
403static BOOLEAN ngfCoeffIsEqual (const coeffs r, n_coeffType n, void * parameter)
404{
405 if (n==n_long_R)
406 {
407 LongComplexInfo* p = (LongComplexInfo *)(parameter);
408 if ((p!=NULL)
409 && (p->float_len == r->float_len)
410 && (p->float_len2 == r->float_len2))
411 return TRUE;
412 }
413 return FALSE;
414}
415
416static void ngfSetChar(const coeffs r)
417{
418 setGMPFloatDigits(r->float_len, r->float_len2);
419}
420
421static char* ngfCoeffName(const coeffs r)
422{
423 STATIC_VAR char ngfCoeffName_buf[30];
424 snprintf(ngfCoeffName_buf,30,"Float(%d,%d)",r->float_len,r->float_len2);
425 return ngfCoeffName_buf;
426}
427
428static number ngfMapQ(number from, const coeffs src, const coeffs dst)
429{
430 assume( getCoeffType(dst) == n_long_R );
431 assume( src->rep == n_rep_gap_rat );
432
434 return (number)res;
435}
436
437static number ngfMapZ(number from, const coeffs aRing, const coeffs r)
438{
440 assume( aRing->rep == n_rep_gmp);
441
442 gmp_float *res=new gmp_float((mpz_ptr)from);
443 return (number)res;
444}
445
446static number ngfMapR(number from, const coeffs src, const coeffs dst)
447{
448 assume( getCoeffType(dst) == n_long_R );
449 assume( getCoeffType(src) == n_R );
450
451 gmp_float *res=new gmp_float((double)nf(from).F());
452 return (number)res;
453}
454
455static number ngfMapP(number from, const coeffs src, const coeffs dst)
456{
457 assume( getCoeffType(dst) == n_long_R );
458 assume( getCoeffType(src) == n_Zp );
459
460 return ngfInit(npInt(from,src), dst); // FIXME? TODO? // extern int npInt (number &n, const coeffs r);
461}
462
463static number ngfMapC(number from, const coeffs src, const coeffs dst)
464{
465 assume( getCoeffType(dst) == n_long_R );
466 assume( getCoeffType(src) == n_long_C );
467
468 gmp_float *res=new gmp_float(((gmp_complex*)from)->real());
469 return (number)res;
470}
471
472static number ngfInitMPZ(mpz_t m, const coeffs)
473{
474 gmp_float *res=new gmp_float(m);
475 return (number)res;
476}
477
478static nMapFunc ngfSetMap(const coeffs src, const coeffs dst)
479{
480 assume( getCoeffType(dst) == n_long_R );
481
482 if (src->rep==n_rep_gap_rat) /*Q, Z*/
483 {
484 return ngfMapQ;
485 }
486 if (src->rep==n_rep_gap_gmp) /*Q, bigint*/
487 {
488 return ngfMapQ;
489 }
490 if (src->rep==n_rep_gmp) /* Z*/
491 {
492 return ngfMapZ;
493 }
494 if ((src->rep==n_rep_gmp_float) && nCoeff_is_long_R(src))
495 {
496 return ndCopyMap; //ngfCopyMap;
497 }
498 if ((src->rep==n_rep_float) && nCoeff_is_R(src))
499 {
500 return ngfMapR;
501 }
502 if ((src->rep==n_rep_gmp_complex) && nCoeff_is_long_C(src))
503 {
504 return ngfMapC;
505 }
506 if ((src->rep==n_rep_int) && nCoeff_is_Zp(src))
507 {
508 return ngfMapP;
509 }
510 return NULL;
511}
512
513BOOLEAN ngfInitChar(coeffs n, void *parameter)
514{
516
517 n->is_field=TRUE;
518 n->is_domain=TRUE;
519 n->rep=n_rep_gmp_float;
520
521 //n->cfKillChar = ndKillChar; /* dummy */
522
523 n->cfSetChar = ngfSetChar;
524 n->ch = 0;
525 n->cfCoeffName=ngfCoeffName;
526
527 n->cfDelete = ngfDelete;
528 //n->cfNormalize=ndNormalize;
529 n->cfInit = ngfInit;
530 n->cfInitMPZ = ngfInitMPZ;
531 n->cfInt = ngfInt;
532 n->cfAdd = ngfAdd;
533 n->cfInpAdd = ngfInpAdd;
534 n->cfSub = ngfSub;
535 n->cfMult = ngfMult;
536 n->cfInpMult = ngfInpMult;
537 n->cfDiv = ngfDiv;
538 n->cfExactDiv= ngfDiv;
539 n->cfInpNeg = ngfNeg;
540 n->cfInvers = ngfInvers;
541 n->cfCopy = ngfCopy;
542 n->cfGreater = ngfGreater;
543 n->cfEqual = ngfEqual;
544 n->cfIsZero = ngfIsZero;
545 n->cfIsOne = ngfIsOne;
546 n->cfIsMOne = ngfIsMOne;
547 n->cfGreaterZero = ngfGreaterZero;
548 n->cfWriteLong = ngfWrite;
549 n->cfRead = ngfRead;
550 n->cfPower = ngfPower;
551 n->cfSetMap = ngfSetMap;
552#ifdef LDEBUG
553 //n->cfDBTest = ndDBTest; // not yet implemented: ngfDBTest
554#endif
555
556 n->nCoeffIsEqual = ngfCoeffIsEqual;
557
558 if( parameter != NULL)
559 {
560 LongComplexInfo* p = (LongComplexInfo*)parameter;
561
562 n->float_len = p->float_len;
563 n->float_len2 = p->float_len2;
564 } else // default values, just for testing!
565 {
566 n->float_len = SHORT_REAL_LENGTH;
567 n->float_len2 = SHORT_REAL_LENGTH;
568 }
569
570 assume( n->float_len2 >= SHORT_REAL_LENGTH );
571
572 assume( n_NumberOfParameters(n) == 0 );
574
575 return FALSE;
576}
All the auxiliary stuff.
static int ABS(int v)
Definition: auxiliary.h:112
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4081
int p
Definition: cfModGcd.cc:4077
CanonicalForm b
Definition: cfModGcd.cc:4102
FILE * f
Definition: checklibs.c:9
gmp_complex numbers based on
Definition: mpr_complex.h:179
void setFromStr(const char *in)
Definition: mpr_complex.cc:78
bool isZero() const
Definition: mpr_complex.cc:252
Coefficient rings, fields and other domains suitable for Singular polynomials.
number ndCopyMap(number a, const coeffs src, const coeffs dst)
Definition: numbers.cc:282
static FORCE_INLINE BOOLEAN nCoeff_is_long_R(const coeffs r)
Definition: coeffs.h:891
n_coeffType
Definition: coeffs.h:27
@ n_R
single prescision (6,6) real numbers
Definition: coeffs.h:31
@ n_long_R
real floating point (GMP) numbers
Definition: coeffs.h:33
@ n_Zp
\F{p < 2^31}
Definition: coeffs.h:29
@ n_long_C
complex floating point (GMP) numbers
Definition: coeffs.h:41
static FORCE_INLINE char const ** n_ParameterNames(const coeffs r)
Returns a (const!) pointer to (const char*) names of parameters.
Definition: coeffs.h:778
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition: coeffs.h:421
static FORCE_INLINE int n_NumberOfParameters(const coeffs r)
Returns the number of parameters.
Definition: coeffs.h:774
static FORCE_INLINE BOOLEAN nCoeff_is_Zp(const coeffs r)
Definition: coeffs.h:800
@ n_rep_gap_rat
(number), see longrat.h
Definition: coeffs.h:111
@ n_rep_gap_gmp
(), see rinteger.h, new impl.
Definition: coeffs.h:112
@ n_rep_float
(float), see shortfl.h
Definition: coeffs.h:116
@ n_rep_int
(int), see modulop.h
Definition: coeffs.h:110
@ n_rep_gmp_float
(gmp_float), see
Definition: coeffs.h:117
@ n_rep_gmp
(mpz_ptr), see rmodulon,h
Definition: coeffs.h:115
@ n_rep_gmp_complex
(gmp_complex), see gnumpc.h
Definition: coeffs.h:118
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:73
static FORCE_INLINE BOOLEAN nCoeff_is_R(const coeffs r)
Definition: coeffs.h:836
static FORCE_INLINE BOOLEAN nCoeff_is_long_C(const coeffs r)
Definition: coeffs.h:894
const CanonicalForm int s
Definition: facAbsFact.cc:51
CanonicalForm res
Definition: facAbsFact.cc:60
bool isZero(const CFArray &A)
checks if entries of A are zero
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define STATIC_VAR
Definition: globaldefs.h:7
static number ngfInit(long i, const coeffs r)
Definition: gnumpfl.cc:39
static number ngfMapC(number from, const coeffs src, const coeffs dst)
Definition: gnumpfl.cc:463
static number ngfCopy(number a, const coeffs r)
Definition: gnumpfl.cc:96
static BOOLEAN ngfGreater(number a, number b, const coeffs r)
Definition: gnumpfl.cc:260
static void ngfSetChar(const coeffs r)
Definition: gnumpfl.cc:416
static void ngfInpAdd(number &a, number b, const coeffs R)
Definition: gnumpfl.cc:161
static number ngfMapZ(number from, const coeffs aRing, const coeffs r)
Definition: gnumpfl.cc:437
static number ngfInvers(number a, const coeffs r)
Definition: gnumpfl.cc:133
static long ngfInt(number &i, const coeffs r)
Definition: gnumpfl.cc:50
static number ngfInitMPZ(mpz_t m, const coeffs)
Definition: gnumpfl.cc:472
static number ngfDiv(number a, number b, const coeffs r)
Definition: gnumpfl.cc:200
static number ngfAdd(number a, number b, const coeffs R)
Definition: gnumpfl.cc:153
static number ngfMapP(number from, const coeffs src, const coeffs dst)
Definition: gnumpfl.cc:455
static BOOLEAN ngfGreaterZero(number a, const coeffs r)
Definition: gnumpfl.cc:250
static BOOLEAN ngfIsMOne(number a, const coeffs r)
Definition: gnumpfl.cc:290
static BOOLEAN ngfIsZero(number a, const coeffs r)
Definition: gnumpfl.cc:61
static void ngfWrite(number a, const coeffs r)
Definition: gnumpfl.cc:385
static number ngfPower(number x, int exp, const coeffs r)
Definition: gnumpfl.cc:221
static BOOLEAN ngfEqual(number a, number b, const coeffs r)
Definition: gnumpfl.cc:270
static int ngfSize(number n, const coeffs r)
Definition: gnumpfl.cc:68
static char * ngfEatFloatNExp(char *s)
Definition: gnumpfl.cc:297
static void ngfDelete(number *a, const coeffs r)
Definition: gnumpfl.cc:82
static number ngfMapQ(number from, const coeffs src, const coeffs dst)
Definition: gnumpfl.cc:428
const char * ngfRead(const char *s, number *a, const coeffs r)
Definition: gnumpfl.cc:326
static BOOLEAN ngfCoeffIsEqual(const coeffs r, n_coeffType n, void *parameter)
Definition: gnumpfl.cc:403
static number ngfNeg(number a, const coeffs r)
Definition: gnumpfl.cc:122
static void ngfInpMult(number &a, number b, const coeffs R)
Definition: gnumpfl.cc:190
static BOOLEAN ngfIsOne(number a, const coeffs r)
Definition: gnumpfl.cc:280
static number ngfMult(number a, number b, const coeffs R)
Definition: gnumpfl.cc:182
static nMapFunc ngfSetMap(const coeffs src, const coeffs dst)
Definition: gnumpfl.cc:478
static char * ngfCoeffName(const coeffs r)
Definition: gnumpfl.cc:421
BOOLEAN ngfInitChar(coeffs n, void *parameter)
Initialize r.
Definition: gnumpfl.cc:513
static number ngfSub(number a, number b, const coeffs R)
Definition: gnumpfl.cc:171
static number ngfMapR(number from, const coeffs src, const coeffs dst)
Definition: gnumpfl.cc:446
#define assume(x)
Definition: mod2.h:389
long npInt(number &n, const coeffs r)
Definition: modulop.cc:85
char * floatToStr(const gmp_float &r, const unsigned int oprec)
Definition: mpr_complex.cc:578
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:357
gmp_float numberFieldToFloat(number num, int cf)
Definition: mpr_complex.cc:438
void setGMPFloatDigits(size_t digits, size_t rest)
Set size of mantissa digits - the number of output digits (basis 10) the size of mantissa consists of...
Definition: mpr_complex.cc:60
#define QTOF
Definition: mpr_complex.h:19
The main handler for Singular numbers which are suitable for Singular polynomials.
const char *const nDivBy0
Definition: numbers.h:88
#define SHORT_REAL_LENGTH
Definition: numbers.h:57
#define omFree(addr)
Definition: omAllocDecl.h:261
#define NULL
Definition: omList.c:12
void StringAppendS(const char *st)
Definition: reporter.cc:107
void Werror(const char *fmt,...)
Definition: reporter.cc:189
static int sign(int x)
Definition: ring.cc:3469
#define SI_FLOAT
Definition: shortfl.h:15
#define R
Definition: sirandom.c:27
Definition: gnumpfl.cc:27
nf(number n)
Definition: gnumpfl.cc:31
nf(SI_FLOAT f)
Definition: gnumpfl.cc:30
SI_FLOAT _f
Definition: gnumpfl.cc:28
number _n
Definition: gnumpfl.cc:29
SI_FLOAT F() const
Definition: gnumpfl.cc:32
number N() const
Definition: gnumpfl.cc:33