[35aab3] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
[341696] | 4 | /* $Id$ */ |
---|
[35aab3] | 5 | |
---|
| 6 | /* |
---|
| 7 | * ABSTRACT - multipolynomial resultants - real floating-point numbers using gmp |
---|
| 8 | * and complex numbers based on pairs of real floating-point numbers |
---|
| 9 | * |
---|
| 10 | */ |
---|
| 11 | |
---|
| 12 | // WARNING! ALWAYS use omAlloc and FreeL when alloc. memory for some char* !! |
---|
| 13 | |
---|
[691dba] | 14 | #include "config.h" |
---|
[35aab3] | 15 | //#ifdef HAVE_MPR |
---|
[2d805a] | 16 | #include <coeffs/coeffs.h> |
---|
[31213a4] | 17 | #include <reporter/reporter.h> |
---|
| 18 | #include <omalloc/omalloc.h> |
---|
[2d805a] | 19 | #include <coeffs/numbers.h> |
---|
| 20 | #include <coeffs/longrat.h> |
---|
[35aab3] | 21 | #include <math.h> |
---|
[2d805a] | 22 | #include <coeffs/mpr_complex.h> |
---|
[35aab3] | 23 | |
---|
| 24 | //%s |
---|
| 25 | // this was copied form longrat0.cc |
---|
| 26 | // and will be used in numberToFloat. |
---|
| 27 | // Make sure that it is up to date!! |
---|
| 28 | #define SR_HDL(A) ((long)(A)) |
---|
| 29 | #define SR_TO_INT(SR) (((long)SR) >> 2) |
---|
| 30 | |
---|
| 31 | #define SIGN_PLUS 1 |
---|
| 32 | #define SIGN_SPACE 2 |
---|
| 33 | #define SIGN_EMPTY 4 |
---|
| 34 | |
---|
| 35 | #define EXTRABYTES 4 |
---|
| 36 | |
---|
| 37 | #define DEFPREC 20 // minimum number of digits (output operations) |
---|
| 38 | size_t gmp_output_digits= DEFPREC; |
---|
| 39 | |
---|
[8fc1b4d] | 40 | static gmp_float *gmpRel=NULL; |
---|
| 41 | static gmp_float *diff=NULL; |
---|
[35aab3] | 42 | |
---|
| 43 | |
---|
| 44 | /** Set size of mantissa |
---|
| 45 | * digits - the number of output digits (basis 10) |
---|
| 46 | * the size of mantissa consists of two parts: |
---|
| 47 | * the "output" part a and the "rest" part b. |
---|
| 48 | * According to the GMP-precision digits is |
---|
| 49 | * recomputed to bits (basis 2). |
---|
| 50 | * Two numbers a, b are equal if |
---|
| 51 | * | a - b | < | a | * 0.1^digits . |
---|
| 52 | * In this case we have a - b = 0 . |
---|
| 53 | * The epsilon e is e=0.1^(digits+rest) with |
---|
| 54 | * 1+e != 1, but 1+0.1*e = 1. |
---|
| 55 | */ |
---|
| 56 | void setGMPFloatDigits( size_t digits, size_t rest ) |
---|
| 57 | { |
---|
| 58 | size_t bits = 1 + (size_t) ((float)digits * 3.5); |
---|
| 59 | size_t rb = 1 + (size_t) ((float)rest * 3.5); |
---|
| 60 | size_t db = bits+rb; |
---|
| 61 | gmp_output_digits= digits; |
---|
| 62 | mpf_set_default_prec( db ); |
---|
[8fc1b4d] | 63 | if (diff!=NULL) delete diff; |
---|
| 64 | diff=new gmp_float(0.0); |
---|
| 65 | mpf_set_prec(*diff->_mpfp(),32); |
---|
| 66 | if (gmpRel!=NULL) delete gmpRel; |
---|
| 67 | gmpRel=new gmp_float(0.0); |
---|
| 68 | mpf_set_prec(*gmpRel->_mpfp(),32); |
---|
| 69 | mpf_set_d(*gmpRel->_mpfp(),0.1); |
---|
| 70 | mpf_pow_ui(*gmpRel->_mpfp(),*gmpRel->_mpfp(),digits); |
---|
[35aab3] | 71 | } |
---|
| 72 | |
---|
| 73 | size_t getGMPFloatDigits() |
---|
| 74 | { |
---|
| 75 | return gmp_output_digits; |
---|
| 76 | } |
---|
| 77 | |
---|
[c8f8f7] | 78 | #if 1 |
---|
[85e68dd] | 79 | void gmp_float::setFromStr(const char * in ) |
---|
[35aab3] | 80 | { |
---|
[d5a95c] | 81 | BOOLEAN neg=false; |
---|
| 82 | if (*in == '-') { in++; neg=TRUE; } |
---|
| 83 | char *s; |
---|
[587cc5] | 84 | if ((s=strchr((char *)in,'E')) !=NULL) |
---|
[d5a95c] | 85 | { |
---|
| 86 | *s='e'; |
---|
| 87 | } |
---|
| 88 | |
---|
[35aab3] | 89 | // gmp doesn't understand number which begin with "." -- it needs 0. |
---|
| 90 | // so, insert the zero |
---|
| 91 | if (*in == '.') |
---|
| 92 | { |
---|
| 93 | int len = strlen(in)+2; |
---|
| 94 | char* c_in = (char*) omAlloc(len); |
---|
| 95 | *c_in = '0'; |
---|
| 96 | strcpy(&(c_in[1]), in); |
---|
| 97 | |
---|
| 98 | mpf_set_str( t, c_in, 10 ); |
---|
| 99 | omFreeSize((void*)c_in, len); |
---|
| 100 | } |
---|
| 101 | else |
---|
| 102 | { |
---|
| 103 | mpf_set_str( t, in, 10 ); |
---|
| 104 | } |
---|
[d5a95c] | 105 | if (neg) mpf_neg( t, t ); |
---|
[35aab3] | 106 | } |
---|
[c8f8f7] | 107 | #else |
---|
| 108 | // problemns with solve_s.tst |
---|
| 109 | void gmp_float::setFromStr(const char * in ) |
---|
| 110 | { |
---|
| 111 | BOOLEAN neg=false; |
---|
| 112 | BOOLEAN E_found=FALSE; |
---|
| 113 | if (*in == '-') { in++; neg=TRUE; } |
---|
| 114 | char *s; |
---|
| 115 | if ((s=strchr(in,'E')) !=NULL) |
---|
| 116 | { |
---|
| 117 | *s='e'; |
---|
| 118 | E_found=TRUE; |
---|
| 119 | } |
---|
| 120 | // gmp doesn't understand number like 1e1, it need 1e+1 |
---|
| 121 | // so, insert the + |
---|
| 122 | if (E_found ||((s=strchr(in,'e')) !=NULL)) |
---|
| 123 | { |
---|
| 124 | if ((*(s+1)!='+') && (*(s+1)!='-')) |
---|
| 125 | { |
---|
| 126 | int len = strlen(in)+3; |
---|
| 127 | char* c_in = (char*) omAlloc(len); |
---|
| 128 | if (*in == '.') |
---|
| 129 | { |
---|
| 130 | *c_in = '0'; |
---|
| 131 | strcpy(&(c_in[1]), in); |
---|
| 132 | } |
---|
| 133 | else |
---|
| 134 | { |
---|
| 135 | strcpy(c_in, in); |
---|
| 136 | } |
---|
| 137 | char * ss=strchr(c_in,'e'); |
---|
| 138 | memmove(ss+2,s+1,strlen(s+1)); |
---|
| 139 | *(ss+1)+'+'; |
---|
| 140 | |
---|
| 141 | mpf_set_str( t, c_in, 10 ); |
---|
| 142 | omFreeSize((void*)c_in, len); |
---|
| 143 | } |
---|
| 144 | } |
---|
| 145 | |
---|
| 146 | // gmp doesn't understand number which begin with "." -- it needs 0. |
---|
| 147 | // so, insert the zero |
---|
| 148 | else if (*in == '.') |
---|
| 149 | { |
---|
| 150 | int len = strlen(in)+2; |
---|
| 151 | char* c_in = (char*) omAlloc(len); |
---|
| 152 | *c_in = '0'; |
---|
| 153 | strcpy(&(c_in[1]), in); |
---|
| 154 | |
---|
| 155 | mpf_set_str( t, c_in, 10 ); |
---|
| 156 | omFreeSize((void*)c_in, len); |
---|
| 157 | } |
---|
| 158 | else |
---|
| 159 | { |
---|
| 160 | mpf_set_str( t, in, 10 ); |
---|
| 161 | } |
---|
| 162 | if (neg) mpf_neg( t, t ); |
---|
| 163 | } |
---|
| 164 | #endif |
---|
[35aab3] | 165 | |
---|
| 166 | |
---|
| 167 | // <gmp_float> = <gmp_float> operator <gmp_float> |
---|
| 168 | gmp_float operator + ( const gmp_float & a, const gmp_float & b ) |
---|
| 169 | { |
---|
| 170 | gmp_float tmp( a ); |
---|
| 171 | tmp += b; |
---|
| 172 | return tmp; |
---|
| 173 | } |
---|
| 174 | gmp_float operator - ( const gmp_float & a, const gmp_float & b ) |
---|
| 175 | { |
---|
| 176 | gmp_float tmp( a ); |
---|
| 177 | tmp -= b; |
---|
| 178 | return tmp; |
---|
| 179 | } |
---|
| 180 | gmp_float operator * ( const gmp_float & a, const gmp_float & b ) |
---|
| 181 | { |
---|
| 182 | gmp_float tmp( a ); |
---|
| 183 | tmp *= b; |
---|
| 184 | return tmp; |
---|
| 185 | } |
---|
| 186 | gmp_float operator / ( const gmp_float & a, const gmp_float & b ) |
---|
| 187 | { |
---|
| 188 | gmp_float tmp( a ); |
---|
| 189 | tmp /= b; |
---|
| 190 | return tmp; |
---|
| 191 | } |
---|
| 192 | |
---|
| 193 | // <gmp_float> operator <gmp_float> |
---|
| 194 | gmp_float & gmp_float::operator += ( const gmp_float & a ) |
---|
| 195 | { |
---|
| 196 | if (mpf_sgn(t) != -(mpf_sgn(a.t))) |
---|
| 197 | { |
---|
| 198 | mpf_add( t, t, a.t); |
---|
| 199 | return *this; |
---|
| 200 | } |
---|
| 201 | if((mpf_sgn(a.t)==0) && (mpf_sgn(t)==0)) |
---|
| 202 | { |
---|
| 203 | mpf_set_d( t, 0.0); |
---|
| 204 | return *this; |
---|
| 205 | } |
---|
| 206 | mpf_add( t, t, a.t ); |
---|
[8fc1b4d] | 207 | mpf_set(diff->t, t); |
---|
| 208 | mpf_set_prec(diff->t, 32); |
---|
| 209 | mpf_div(diff->t, diff->t, a.t); |
---|
| 210 | mpf_abs(diff->t, diff->t); |
---|
| 211 | if(mpf_cmp(diff->t, gmpRel->t) < 0) |
---|
[35aab3] | 212 | mpf_set_d( t, 0.0); |
---|
| 213 | return *this; |
---|
| 214 | } |
---|
| 215 | gmp_float & gmp_float::operator -= ( const gmp_float & a ) |
---|
| 216 | { |
---|
| 217 | if (mpf_sgn(t) != mpf_sgn(a.t)) |
---|
| 218 | { |
---|
| 219 | mpf_sub( t, t, a.t); |
---|
| 220 | return *this; |
---|
| 221 | } |
---|
| 222 | if((mpf_sgn(a.t)==0) && (mpf_sgn(t)==0)) |
---|
| 223 | { |
---|
| 224 | mpf_set_d( t, 0.0); |
---|
| 225 | return *this; |
---|
| 226 | } |
---|
| 227 | mpf_sub( t, t, a.t ); |
---|
[8fc1b4d] | 228 | mpf_set(diff->t, t); |
---|
| 229 | mpf_set_prec(diff->t, 32); |
---|
| 230 | mpf_div(diff->t, diff->t, a.t); |
---|
| 231 | mpf_abs(diff->t, diff->t); |
---|
| 232 | if(mpf_cmp(diff->t, gmpRel->t) < 0) |
---|
[35aab3] | 233 | mpf_set_d( t, 0.0); |
---|
| 234 | return *this; |
---|
| 235 | } |
---|
| 236 | |
---|
| 237 | // <gmp_float> == <gmp_float> ?? |
---|
| 238 | bool operator == ( const gmp_float & a, const gmp_float & b ) |
---|
| 239 | { |
---|
| 240 | if(mpf_sgn(a.t) != mpf_sgn(b.t)) |
---|
| 241 | return false; |
---|
| 242 | if((mpf_sgn(a.t)==0) && (mpf_sgn(b.t)==0)) |
---|
| 243 | return true; |
---|
[8fc1b4d] | 244 | mpf_sub(diff->t, a.t, b.t); |
---|
| 245 | mpf_div(diff->t, diff->t, a.t); |
---|
| 246 | mpf_abs(diff->t, diff->t); |
---|
| 247 | if(mpf_cmp(diff->t, gmpRel->t) < 0) |
---|
[35aab3] | 248 | return true; |
---|
| 249 | else |
---|
| 250 | return false; |
---|
| 251 | } |
---|
| 252 | // t == 0 ? |
---|
[66b967] | 253 | bool gmp_float::isZero() const |
---|
[35aab3] | 254 | { |
---|
| 255 | return (mpf_sgn( t ) == 0); |
---|
| 256 | } |
---|
| 257 | // t == 1 ? |
---|
[66b967] | 258 | bool gmp_float::isOne() const |
---|
[35aab3] | 259 | { |
---|
| 260 | #ifdef VARIANTE_1 |
---|
| 261 | return (mpf_cmp_ui( t , 1 ) == 0); |
---|
| 262 | #else |
---|
| 263 | if (mpf_sgn(t) <= 0) |
---|
| 264 | return false; |
---|
[8fc1b4d] | 265 | mpf_sub_ui(diff->t, t, 1); |
---|
| 266 | mpf_abs(diff->t, diff->t); |
---|
| 267 | if(mpf_cmp(diff->t, gmpRel->t) < 0) |
---|
[35aab3] | 268 | return true; |
---|
| 269 | else |
---|
| 270 | return false; |
---|
| 271 | #endif |
---|
| 272 | } |
---|
| 273 | // t == -1 ? |
---|
[66b967] | 274 | bool gmp_float::isMOne() const |
---|
[35aab3] | 275 | { |
---|
| 276 | #ifdef VARIANTE_1 |
---|
| 277 | return (mpf_cmp_si( t , -1 ) == 0); |
---|
| 278 | #else |
---|
| 279 | if (mpf_sgn(t) >= 0) |
---|
| 280 | return false; |
---|
[8fc1b4d] | 281 | mpf_add_ui(diff->t, t, 1); |
---|
| 282 | mpf_abs(diff->t, diff->t); |
---|
| 283 | if(mpf_cmp(diff->t, gmpRel->t) < 0) |
---|
[35aab3] | 284 | return true; |
---|
| 285 | else |
---|
| 286 | return false; |
---|
| 287 | #endif |
---|
| 288 | } |
---|
| 289 | bool operator > ( const gmp_float & a, const gmp_float & b ) |
---|
| 290 | { |
---|
| 291 | if (a.t == b.t) |
---|
| 292 | return false; |
---|
| 293 | return mpf_cmp( a.t, b.t ) > 0; |
---|
| 294 | } |
---|
| 295 | bool operator < ( const gmp_float & a, const gmp_float & b ) |
---|
| 296 | { |
---|
| 297 | if (a.t == b.t) |
---|
| 298 | return false; |
---|
| 299 | return mpf_cmp( a.t, b.t ) < 0; |
---|
| 300 | } |
---|
| 301 | bool operator >= ( const gmp_float & a, const gmp_float & b ) |
---|
| 302 | { |
---|
| 303 | if (a.t == b.t) |
---|
| 304 | return true; |
---|
| 305 | return mpf_cmp( a.t, b.t ) >= 0; |
---|
| 306 | } |
---|
| 307 | bool operator <= ( const gmp_float & a, const gmp_float & b ) |
---|
| 308 | { |
---|
| 309 | if (a.t == b.t) |
---|
| 310 | return true; |
---|
| 311 | return mpf_cmp( a.t, b.t ) <= 0; |
---|
| 312 | } |
---|
| 313 | |
---|
| 314 | // unary - |
---|
| 315 | gmp_float operator - ( const gmp_float & a ) |
---|
| 316 | { |
---|
| 317 | gmp_float tmp; |
---|
| 318 | mpf_neg( *(tmp._mpfp()), *(a.mpfp()) ); |
---|
| 319 | return tmp; |
---|
| 320 | } |
---|
| 321 | |
---|
| 322 | gmp_float abs( const gmp_float & a ) |
---|
| 323 | { |
---|
| 324 | gmp_float tmp; |
---|
| 325 | mpf_abs( *(tmp._mpfp()), *a.mpfp() ); |
---|
| 326 | return tmp; |
---|
| 327 | } |
---|
| 328 | gmp_float sqrt( const gmp_float & a ) |
---|
| 329 | { |
---|
| 330 | gmp_float tmp; |
---|
| 331 | mpf_sqrt( *(tmp._mpfp()), *a.mpfp() ); |
---|
| 332 | return tmp; |
---|
| 333 | } |
---|
| 334 | gmp_float sin( const gmp_float & a ) |
---|
| 335 | { |
---|
| 336 | gmp_float tmp( sin((double)a) ); |
---|
| 337 | return tmp; |
---|
| 338 | } |
---|
| 339 | gmp_float cos( const gmp_float & a ) |
---|
| 340 | { |
---|
| 341 | gmp_float tmp( cos((double)a) ); |
---|
| 342 | return tmp; |
---|
| 343 | } |
---|
| 344 | gmp_float log( const gmp_float & a ) |
---|
| 345 | { |
---|
| 346 | gmp_float tmp( log((double)a) ); |
---|
| 347 | return tmp; |
---|
| 348 | } |
---|
| 349 | gmp_float hypot( const gmp_float & a, const gmp_float & b ) |
---|
| 350 | { |
---|
| 351 | #if 1 |
---|
| 352 | return ( sqrt( (a*a) + (b*b) ) ); |
---|
| 353 | #else |
---|
| 354 | gmp_float tmp( hypot( (double)a, (double)b ) ); |
---|
| 355 | return tmp; |
---|
| 356 | #endif |
---|
| 357 | } |
---|
| 358 | gmp_float exp( const gmp_float & a ) |
---|
| 359 | { |
---|
| 360 | gmp_float tmp( exp((double)a) ); |
---|
| 361 | return tmp; |
---|
| 362 | } |
---|
| 363 | gmp_float max( const gmp_float & a, const gmp_float & b ) |
---|
| 364 | { |
---|
| 365 | gmp_float tmp; |
---|
| 366 | a > b ? tmp= a : tmp= b; |
---|
| 367 | return tmp; |
---|
| 368 | } |
---|
| 369 | // |
---|
| 370 | // number to float, number = Q, R, C |
---|
| 371 | // makes a COPY of num! (Ist das gut?) |
---|
| 372 | // |
---|
[ac5686] | 373 | gmp_float numberToFloat( number num, const coeffs src) |
---|
[35aab3] | 374 | { |
---|
| 375 | gmp_float r; |
---|
| 376 | |
---|
[1cce47] | 377 | if ( nCoeff_is_Q(src) ) |
---|
[35aab3] | 378 | { |
---|
| 379 | if ( num != NULL ) |
---|
| 380 | { |
---|
| 381 | if (SR_HDL(num) & SR_INT) |
---|
| 382 | { |
---|
| 383 | r= SR_TO_INT(num); |
---|
| 384 | } |
---|
| 385 | else |
---|
| 386 | { |
---|
| 387 | if ( num->s == 0 ) |
---|
| 388 | { |
---|
[ac5686] | 389 | nlNormalize( num, src ); |
---|
[35aab3] | 390 | } |
---|
| 391 | if (SR_HDL(num) & SR_INT) |
---|
| 392 | { |
---|
| 393 | r= SR_TO_INT(num); |
---|
| 394 | } |
---|
| 395 | else |
---|
| 396 | { |
---|
| 397 | if ( num->s != 3 ) |
---|
| 398 | { |
---|
[a604c3] | 399 | r= num->z; |
---|
| 400 | r/= (gmp_float)num->n; |
---|
[35aab3] | 401 | } |
---|
| 402 | else |
---|
| 403 | { |
---|
[a604c3] | 404 | r= num->z; |
---|
[35aab3] | 405 | } |
---|
| 406 | } |
---|
| 407 | } |
---|
| 408 | } |
---|
| 409 | else |
---|
| 410 | { |
---|
| 411 | r= 0.0; |
---|
| 412 | } |
---|
| 413 | } |
---|
[1cce47] | 414 | else if (nCoeff_is_long_R(src) || nCoeff_is_long_C(src)) |
---|
[35aab3] | 415 | { |
---|
| 416 | r= *(gmp_float*)num; |
---|
| 417 | } |
---|
[1cce47] | 418 | else if ( nCoeff_is_R(src) ) |
---|
[35aab3] | 419 | { |
---|
| 420 | // Add some code here :-) |
---|
| 421 | WerrorS("Ground field not implemented!"); |
---|
| 422 | } |
---|
| 423 | else |
---|
| 424 | { |
---|
| 425 | WerrorS("Ground field not implemented!"); |
---|
| 426 | } |
---|
| 427 | |
---|
| 428 | return r; |
---|
| 429 | } |
---|
| 430 | |
---|
[ac5686] | 431 | gmp_float numberFieldToFloat( number num, int k, const coeffs src) |
---|
[35aab3] | 432 | { |
---|
| 433 | gmp_float r; |
---|
| 434 | |
---|
| 435 | switch (k) |
---|
| 436 | { |
---|
| 437 | case QTOF: |
---|
| 438 | if ( num != NULL ) |
---|
| 439 | { |
---|
| 440 | if (SR_HDL(num) & SR_INT) |
---|
| 441 | { |
---|
| 442 | r= SR_TO_INT(num); |
---|
| 443 | } |
---|
| 444 | else |
---|
| 445 | { |
---|
| 446 | if ( num->s == 0 ) |
---|
| 447 | { |
---|
[ac5686] | 448 | nlNormalize( num, src ); |
---|
[35aab3] | 449 | } |
---|
| 450 | if (SR_HDL(num) & SR_INT) |
---|
| 451 | { |
---|
| 452 | r= SR_TO_INT(num); |
---|
| 453 | } |
---|
| 454 | else |
---|
| 455 | { |
---|
| 456 | if ( num->s != 3 ) |
---|
| 457 | { |
---|
[a604c3] | 458 | r= num->z; |
---|
| 459 | r/= (gmp_float)num->n; |
---|
[35aab3] | 460 | } |
---|
| 461 | else |
---|
| 462 | { |
---|
[a604c3] | 463 | r= num->z; |
---|
[35aab3] | 464 | } |
---|
| 465 | } |
---|
| 466 | } |
---|
| 467 | } |
---|
| 468 | else |
---|
| 469 | { |
---|
| 470 | r= 0.0; |
---|
| 471 | } |
---|
| 472 | break; |
---|
| 473 | case RTOF: |
---|
| 474 | r= *(gmp_float*)num; |
---|
| 475 | break; |
---|
| 476 | case CTOF: |
---|
| 477 | WerrorS("Can not map from field C to field R!"); |
---|
| 478 | break; |
---|
| 479 | case ZTOF: |
---|
| 480 | default: |
---|
| 481 | WerrorS("Ground field not implemented!"); |
---|
| 482 | } // switch |
---|
| 483 | |
---|
| 484 | return r; |
---|
| 485 | } |
---|
| 486 | |
---|
| 487 | // Do some strange things with the mantissa string and the exponent |
---|
| 488 | // to get some nice output string. |
---|
| 489 | char *nicifyFloatStr( char * in, mp_exp_t exponent, size_t oprec, int *size, int thesign ) |
---|
| 490 | { |
---|
| 491 | char *out; |
---|
| 492 | |
---|
| 493 | int sign= (in[0] == '-') ? 1 : 0; |
---|
| 494 | char csign[2]; |
---|
| 495 | |
---|
| 496 | switch (thesign) |
---|
| 497 | { |
---|
| 498 | case SIGN_PLUS: |
---|
| 499 | sign ? strcpy(csign,"-") : strcpy(csign,"+"); //+123, -123 |
---|
| 500 | break; |
---|
| 501 | case SIGN_SPACE: |
---|
| 502 | sign ? strcpy(csign,"-") : strcpy(csign," "); // 123, -123 |
---|
| 503 | break; |
---|
| 504 | case SIGN_EMPTY: |
---|
| 505 | default: |
---|
| 506 | sign ? strcpy(csign,"-") : strcpy(csign,""); //123, -123 |
---|
| 507 | break; |
---|
| 508 | } |
---|
| 509 | |
---|
| 510 | if ( strlen(in) == 0 ) |
---|
| 511 | { |
---|
| 512 | *size= 2*sizeof(char); |
---|
| 513 | return omStrDup("0"); |
---|
| 514 | } |
---|
| 515 | |
---|
| 516 | if ( ((unsigned int)ABS(exponent) <= oprec) |
---|
| 517 | /*|| (exponent+sign >= (int)strlen(in))*/ ) |
---|
| 518 | { |
---|
| 519 | if ( exponent+sign < (int)strlen(in) ) |
---|
| 520 | { |
---|
| 521 | int eexponent= (exponent >= 0) ? 0 : -exponent; |
---|
| 522 | int eeexponent= (exponent >= 0) ? exponent : 0; |
---|
| 523 | *size= (strlen(in)+15+eexponent) * sizeof(char); |
---|
| 524 | out= (char*)omAlloc(*size); |
---|
| 525 | memset(out,0,*size); |
---|
| 526 | |
---|
| 527 | strcpy(out,csign); |
---|
| 528 | strncat(out,in+sign,eeexponent); |
---|
| 529 | |
---|
| 530 | if (exponent == 0) |
---|
| 531 | strcat(out,"0."); |
---|
| 532 | else if ( exponent > 0 ) |
---|
| 533 | strcat(out,"."); |
---|
| 534 | else |
---|
| 535 | { |
---|
| 536 | strcat(out,"0."); |
---|
| 537 | memset(out+strlen(out),'0',eexponent); |
---|
| 538 | } |
---|
| 539 | strcat(out,in+sign+eeexponent); |
---|
| 540 | } |
---|
| 541 | else if ( exponent+sign > (int)strlen(in) ) |
---|
| 542 | { |
---|
| 543 | *size= (strlen(in)+exponent+12)*sizeof(char); |
---|
| 544 | out= (char*)omAlloc(*size); |
---|
| 545 | memset(out,0,*size); |
---|
| 546 | sprintf(out,"%s%s",csign,in+sign); |
---|
| 547 | memset(out+strlen(out),'0',exponent-strlen(in)+sign); |
---|
| 548 | } |
---|
| 549 | else |
---|
| 550 | { |
---|
| 551 | *size= (strlen(in)+2) * sizeof(char) + 10; |
---|
| 552 | out= (char*)omAlloc(*size); |
---|
| 553 | memset(out,0,*size); |
---|
| 554 | sprintf(out,"%s%s",csign,in+sign); |
---|
| 555 | } |
---|
| 556 | } |
---|
| 557 | else |
---|
| 558 | { |
---|
| 559 | // if ( exponent > 0 ) |
---|
| 560 | // { |
---|
| 561 | int c=1,d=10; |
---|
| 562 | while ( exponent / d > 0 ) |
---|
| 563 | { // count digits |
---|
| 564 | d*=10; |
---|
| 565 | c++; |
---|
| 566 | } |
---|
| 567 | *size= (strlen(in)+12+c) * sizeof(char) + 10; |
---|
| 568 | out= (char*)omAlloc(*size); |
---|
| 569 | memset(out,0,*size); |
---|
| 570 | sprintf(out,"%s0.%se%s%d",csign,in+sign,exponent>=0?"+":"",(int)exponent); |
---|
| 571 | // } |
---|
| 572 | // else |
---|
| 573 | // { |
---|
| 574 | // *size=2; |
---|
| 575 | // out= (char*)omAlloc(*size); |
---|
| 576 | // strcpy(out,"0"); |
---|
| 577 | // } |
---|
| 578 | } |
---|
| 579 | return out; |
---|
| 580 | } |
---|
| 581 | |
---|
| 582 | char *floatToStr( const gmp_float & r, const unsigned int oprec ) |
---|
| 583 | { |
---|
| 584 | #if 1 |
---|
| 585 | mp_exp_t exponent; |
---|
| 586 | int size,insize; |
---|
| 587 | char *nout,*out,*in; |
---|
| 588 | |
---|
| 589 | insize= (oprec+2) * sizeof(char) + 10; |
---|
| 590 | in= (char*)omAlloc( insize ); |
---|
| 591 | |
---|
| 592 | mpf_get_str(in,&exponent,10,oprec,*(r.mpfp())); |
---|
| 593 | |
---|
| 594 | if ( (exponent > 0) |
---|
| 595 | && (exponent < (int)oprec) |
---|
| 596 | && (strlen(in)-(in[0]=='-'?1:0) == oprec) ) |
---|
| 597 | { |
---|
[7d90aa] | 598 | omFree( (void *) in ); |
---|
[35aab3] | 599 | insize= (exponent+oprec+2) * sizeof(char) + 10; |
---|
| 600 | in= (char*)omAlloc( insize ); |
---|
| 601 | int newprec= exponent+oprec; |
---|
| 602 | mpf_get_str(in,&exponent,10,newprec,*(r.mpfp())); |
---|
| 603 | } |
---|
| 604 | nout= nicifyFloatStr( in, exponent, oprec, &size, SIGN_EMPTY ); |
---|
[7d90aa] | 605 | omFree( (void *) in ); |
---|
[35aab3] | 606 | out= (char*)omAlloc( (strlen(nout)+1) * sizeof(char) ); |
---|
| 607 | strcpy( out, nout ); |
---|
[7d90aa] | 608 | omFree( (void *) nout ); |
---|
[35aab3] | 609 | |
---|
| 610 | return out; |
---|
| 611 | #else |
---|
| 612 | // for testing purpose... |
---|
| 613 | char *out= (char*)omAlloc( (1024) * sizeof(char) ); |
---|
| 614 | sprintf(out,"% .10f",(double)r); |
---|
| 615 | return out; |
---|
| 616 | #endif |
---|
| 617 | } |
---|
| 618 | //<- |
---|
| 619 | |
---|
| 620 | //-> gmp_complex::* |
---|
| 621 | // <gmp_complex> = <gmp_complex> operator <gmp_complex> |
---|
| 622 | // |
---|
| 623 | gmp_complex operator + ( const gmp_complex & a, const gmp_complex & b ) |
---|
| 624 | { |
---|
| 625 | return gmp_complex( a.r + b.r, a.i + b.i ); |
---|
| 626 | } |
---|
| 627 | gmp_complex operator - ( const gmp_complex & a, const gmp_complex & b ) |
---|
| 628 | { |
---|
| 629 | return gmp_complex( a.r - b.r, a.i - b.i ); |
---|
| 630 | } |
---|
| 631 | gmp_complex operator * ( const gmp_complex & a, const gmp_complex & b ) |
---|
| 632 | { |
---|
| 633 | return gmp_complex( a.r * b.r - a.i * b.i, |
---|
| 634 | a.r * b.i + a.i * b.r); |
---|
| 635 | } |
---|
| 636 | gmp_complex operator / ( const gmp_complex & a, const gmp_complex & b ) |
---|
| 637 | { |
---|
| 638 | gmp_float d = b.r*b.r + b.i*b.i; |
---|
| 639 | return gmp_complex( (a.r * b.r + a.i * b.i) / d, |
---|
| 640 | (a.i * b.r - a.r * b.i) / d); |
---|
| 641 | } |
---|
| 642 | |
---|
| 643 | // <gmp_complex> operator <gmp_complex> |
---|
| 644 | // |
---|
| 645 | gmp_complex & gmp_complex::operator += ( const gmp_complex & b ) |
---|
| 646 | { |
---|
| 647 | r+=b.r; |
---|
| 648 | i+=b.i; |
---|
| 649 | return *this; |
---|
| 650 | } |
---|
| 651 | gmp_complex & gmp_complex::operator -= ( const gmp_complex & b ) |
---|
| 652 | { |
---|
| 653 | r-=b.r; |
---|
| 654 | i-=b.i; |
---|
| 655 | return *this; |
---|
| 656 | } |
---|
| 657 | gmp_complex & gmp_complex::operator *= ( const gmp_complex & b ) |
---|
| 658 | { |
---|
| 659 | gmp_float f = r * b.r - i * b.i; |
---|
| 660 | i = r * b.i + i * b.r; |
---|
| 661 | r = f; |
---|
| 662 | return *this; |
---|
| 663 | } |
---|
[95d8df] | 664 | gmp_complex & gmp_complex::neg ( ) |
---|
| 665 | { |
---|
| 666 | i.neg(); |
---|
| 667 | r.neg(); |
---|
| 668 | return *this; |
---|
| 669 | } |
---|
[35aab3] | 670 | gmp_complex & gmp_complex::operator /= ( const gmp_complex & b ) |
---|
| 671 | { |
---|
| 672 | gmp_float d = b.r*b.r + b.i*b.i; |
---|
| 673 | r = (r * b.r + i * b.i) / d; |
---|
| 674 | i = (i * b.r - r * b.i) / d; |
---|
| 675 | return *this; |
---|
| 676 | } |
---|
| 677 | |
---|
| 678 | // Returns square root of gmp_complex number |
---|
| 679 | // |
---|
| 680 | gmp_complex sqrt( const gmp_complex & x ) |
---|
| 681 | { |
---|
| 682 | gmp_float r = abs(x); |
---|
| 683 | gmp_float nr, ni; |
---|
| 684 | if (r == (gmp_float) 0.0) |
---|
| 685 | { |
---|
| 686 | nr = ni = r; |
---|
| 687 | } |
---|
| 688 | else if ( x.real() > (gmp_float)0) |
---|
| 689 | { |
---|
| 690 | nr = sqrt((gmp_float)0.5 * (r + x.real())); |
---|
| 691 | ni = x.imag() / nr / (gmp_float)2; |
---|
| 692 | } |
---|
| 693 | else |
---|
| 694 | { |
---|
| 695 | ni = sqrt((gmp_float)0.5 * (r - x.real())); |
---|
| 696 | if (x.imag() < (gmp_float)0) |
---|
| 697 | { |
---|
| 698 | ni = - ni; |
---|
| 699 | } |
---|
| 700 | nr = x.imag() / ni / (gmp_float)2; |
---|
| 701 | } |
---|
[88d2e5] | 702 | gmp_complex tmp(nr, ni); |
---|
| 703 | return tmp; |
---|
[35aab3] | 704 | } |
---|
| 705 | |
---|
| 706 | // converts a gmp_complex to a string ( <real part> + I * <imaginary part> ) |
---|
| 707 | // |
---|
[ac5686] | 708 | char *complexToStr( gmp_complex & c, const unsigned int oprec, const coeffs src ) |
---|
[35aab3] | 709 | { |
---|
[7fee876] | 710 | const char * complex_parameter = "I"; |
---|
| 711 | int N = 1; // strlen(complex_parameter); |
---|
| 712 | |
---|
| 713 | if (nCoeff_is_long_C(src)) |
---|
| 714 | { |
---|
| 715 | complex_parameter = n_ParameterNames(src)[0]; |
---|
| 716 | N = strlen(complex_parameter); |
---|
| 717 | } |
---|
| 718 | |
---|
| 719 | assume( complex_parameter != NULL && N > 0); |
---|
[4c6e420] | 720 | |
---|
[35aab3] | 721 | char *out,*in_imag,*in_real; |
---|
| 722 | |
---|
| 723 | c.SmallToZero(); |
---|
| 724 | if ( !c.imag().isZero() ) |
---|
| 725 | { |
---|
| 726 | |
---|
| 727 | in_real=floatToStr( c.real(), oprec ); // get real part |
---|
| 728 | in_imag=floatToStr( abs(c.imag()), oprec ); // get imaginary part |
---|
| 729 | |
---|
[1cce47] | 730 | if (nCoeff_is_long_C(src)) |
---|
[35aab3] | 731 | { |
---|
[7fee876] | 732 | int len=(strlen(in_real)+strlen(in_imag)+7+N)*sizeof(char); |
---|
[35aab3] | 733 | out=(char*)omAlloc(len); |
---|
| 734 | memset(out,0,len); |
---|
| 735 | if ( !c.real().isZero() ) // (-23-i*5.43) or (15.1+i*5.3) |
---|
[7fee876] | 736 | sprintf(out,"(%s%s%s*%s)",in_real,c.imag().sign()>=0?"+":"-",complex_parameter,in_imag); |
---|
[35aab3] | 737 | else // (-i*43) or (i*34) |
---|
| 738 | { |
---|
| 739 | if (c.imag().isOne()) |
---|
[7fee876] | 740 | sprintf(out,"%s", complex_parameter); |
---|
[35aab3] | 741 | else if (c.imag().isMOne()) |
---|
[7fee876] | 742 | sprintf(out,"-%s", complex_parameter); |
---|
[35aab3] | 743 | else |
---|
[7fee876] | 744 | sprintf(out,"(%s%s*%s)",c.imag().sign()>=0?"":"-", complex_parameter,in_imag); |
---|
[35aab3] | 745 | } |
---|
| 746 | } |
---|
| 747 | else |
---|
| 748 | { |
---|
| 749 | int len=(strlen(in_real)+strlen(in_imag)+9) * sizeof(char); |
---|
| 750 | out=(char*)omAlloc( len ); |
---|
| 751 | memset(out,0,len); |
---|
| 752 | if ( !c.real().isZero() ) |
---|
| 753 | sprintf(out,"(%s%s%s)",in_real,c.imag().sign()>=0?"+I*":"-I*",in_imag); |
---|
| 754 | else |
---|
| 755 | sprintf(out,"(%s%s)",c.imag().sign()>=0?"I*":"-I*",in_imag); |
---|
| 756 | } |
---|
[7d90aa] | 757 | omFree( (void *) in_real ); |
---|
| 758 | omFree( (void *) in_imag ); |
---|
[35aab3] | 759 | } |
---|
| 760 | else |
---|
| 761 | { |
---|
| 762 | out= floatToStr( c.real(), oprec ); |
---|
| 763 | } |
---|
| 764 | |
---|
| 765 | return out; |
---|
| 766 | } |
---|
| 767 | //<- |
---|
| 768 | |
---|
| 769 | bool complexNearZero( gmp_complex * c, int digits ) |
---|
| 770 | { |
---|
| 771 | gmp_float eps,epsm; |
---|
| 772 | |
---|
| 773 | if ( digits < 1 ) return true; |
---|
| 774 | |
---|
| 775 | eps=pow(10.0,(int)digits); |
---|
| 776 | //Print("eps: %s\n",floatToStr(eps,gmp_output_digits)); |
---|
| 777 | eps=(gmp_float)1.0/eps; |
---|
| 778 | epsm=-eps; |
---|
| 779 | |
---|
| 780 | //Print("eps: %s\n",floatToStr(eps,gmp_output_digits)); |
---|
| 781 | |
---|
| 782 | if ( c->real().sign() > 0 ) // + |
---|
| 783 | return (c->real() < eps && (c->imag() < eps && c->imag() > epsm)); |
---|
| 784 | else // - |
---|
| 785 | return (c->real() > epsm && (c->imag() < eps && c->imag() > epsm)); |
---|
| 786 | } |
---|
| 787 | |
---|
| 788 | void gmp_complex::SmallToZero() |
---|
| 789 | { |
---|
| 790 | gmp_float ar=this->real(); |
---|
| 791 | gmp_float ai=this->imag(); |
---|
| 792 | if (ar.isZero() || ai.isZero()) return; |
---|
| 793 | mpf_abs(*ar._mpfp(), *ar._mpfp()); |
---|
| 794 | mpf_abs(*ai._mpfp(), *ai._mpfp()); |
---|
| 795 | mpf_set_prec(*ar._mpfp(), 32); |
---|
| 796 | mpf_set_prec(*ai._mpfp(), 32); |
---|
| 797 | if (ar > ai) |
---|
| 798 | { |
---|
| 799 | mpf_div(*ai._mpfp(), *ai._mpfp(), *ar._mpfp()); |
---|
[8fc1b4d] | 800 | if (ai < *gmpRel) this->imag(0.0); |
---|
[35aab3] | 801 | } |
---|
| 802 | else |
---|
| 803 | { |
---|
| 804 | mpf_div(*ar._mpfp(), *ar._mpfp(), *ai._mpfp()); |
---|
[8fc1b4d] | 805 | if (ar < *gmpRel) this->real(0.0); |
---|
[35aab3] | 806 | } |
---|
| 807 | } |
---|
| 808 | |
---|
| 809 | //%e |
---|
| 810 | |
---|
| 811 | //#endif // HAVE_MPR |
---|
| 812 | |
---|
| 813 | // local Variables: *** |
---|
| 814 | // folded-file: t *** |
---|
| 815 | // compile-command-1: "make installg" *** |
---|
| 816 | // compile-command-2: "make install" *** |
---|
| 817 | // End: *** |
---|