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