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