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