[8a150b] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
[f543de1] | 4 | /* $Id: mpr_complex.cc,v 1.5 1999-06-24 07:46:50 wenk 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 | |
---|
[f543de1] | 12 | // WARNING! ALWAYS use AllocL and FreeL when alloc. memory for some char* !! |
---|
| 13 | |
---|
[8a150b] | 14 | #include "mod2.h" |
---|
| 15 | //#ifdef HAVE_MPR |
---|
| 16 | #include "structs.h" |
---|
| 17 | #include "febase.h" |
---|
| 18 | #include "mmemory.h" |
---|
| 19 | #include "numbers.h" |
---|
| 20 | #include "longrat.h" |
---|
| 21 | #include <math.h> |
---|
| 22 | #include "mpr_complex.h" |
---|
| 23 | |
---|
[f543de1] | 24 | //%s |
---|
[8a150b] | 25 | // this was copied form longrat0.cc |
---|
[f543de1] | 26 | // and will be used in numberToFloat. |
---|
| 27 | // Make sure that it is up to date!! |
---|
[8a150b] | 28 | #define SR_HDL(A) ((long)(A)) |
---|
| 29 | #define SR_INT 1 |
---|
| 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 | size_t gmp_output_digits= DEFPREC; |
---|
| 39 | |
---|
[f543de1] | 40 | //-> setGMPFloat* |
---|
[8a150b] | 41 | /** Set size of mantissa to <bytes> bytes and guess the number of |
---|
| 42 | * output digits. |
---|
| 43 | * |
---|
| 44 | * Set internal gmp floating point precision to (bytes+EXTRABYTES)*8 bits and |
---|
| 45 | * set external precision (i.e. any number smaller than this is treated as ZERO) |
---|
| 46 | * to bytes*8 bits. |
---|
| 47 | * The precision in bytes is the size of the mantissa! |
---|
| 48 | * Guesses the number of digits <-> precision |
---|
| 49 | */ |
---|
| 50 | void setGMPFloatPrecBytes( unsigned long int bytes ) |
---|
| 51 | { |
---|
| 52 | unsigned long int bits= bytes * 8; |
---|
| 53 | gmp_float::setPrecision( (bytes+EXTRABYTES)*8 ); |
---|
| 54 | gmp_float::setEqualBits( bits ); |
---|
| 55 | // guess the maximal number of digits for this precision |
---|
| 56 | gmp_output_digits= -10 + (size_t) |
---|
| 57 | ( (((bits>64?bits:64)+2*mp_bits_per_limb-1)/mp_bits_per_limb) |
---|
| 58 | * mp_bits_per_limb * (log(2)/log(10))); |
---|
| 59 | } |
---|
| 60 | |
---|
| 61 | unsigned long int getGMPFloatPrecBytes() |
---|
| 62 | { |
---|
| 63 | return gmp_float::getEqualBits()/8; |
---|
| 64 | } |
---|
| 65 | |
---|
[f543de1] | 66 | // Sets the lenght of the mantissa to <digits> digits |
---|
[8a150b] | 67 | void setGMPFloatDigits( size_t digits ) |
---|
| 68 | { |
---|
| 69 | size_t bits= 1 + (size_t) (digits / (log(2)/log(10))); |
---|
| 70 | bits= bits>64?bits:64; |
---|
[f543de1] | 71 | // gmp_float::setPrecision( bits+EXTRABYTES*8 ); |
---|
| 72 | gmp_float::setPrecision( bits+(bits/2) ); |
---|
[8a150b] | 73 | gmp_float::setEqualBits( bits ); |
---|
| 74 | gmp_output_digits= digits; |
---|
| 75 | } |
---|
| 76 | |
---|
| 77 | size_t getGMPFloatDigits() |
---|
| 78 | { |
---|
| 79 | return gmp_output_digits; |
---|
| 80 | } |
---|
[f543de1] | 81 | //<- |
---|
[8a150b] | 82 | |
---|
| 83 | //-> gmp_float::* |
---|
| 84 | unsigned long int gmp_float::gmp_default_prec_bits= GMP_DEFAULT_PREC_BITS; |
---|
| 85 | unsigned long int gmp_float::gmp_needequal_bits= GMP_NEEDEQUAL_BITS; |
---|
| 86 | |
---|
| 87 | // <gmp_float> = <gmp_float> operator <gmp_float> |
---|
| 88 | gmp_float operator + ( const gmp_float & a, const gmp_float & b ) |
---|
| 89 | { |
---|
| 90 | gmp_float tmp( a ); |
---|
| 91 | tmp += b; |
---|
| 92 | return tmp; |
---|
| 93 | } |
---|
| 94 | gmp_float operator - ( const gmp_float & a, const gmp_float & b ) |
---|
| 95 | { |
---|
| 96 | gmp_float tmp( a ); |
---|
| 97 | tmp -= b; |
---|
| 98 | return tmp; |
---|
| 99 | } |
---|
| 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 | |
---|
| 113 | // <gmp_float> == <gmp_float> ?? up to the first gmp_float::gmp_needequal_bits bits |
---|
| 114 | bool operator == ( const gmp_float & a, const gmp_float & b ) |
---|
| 115 | { |
---|
| 116 | //return mpf_cmp( a.t, b.t ) == 0; |
---|
| 117 | return mpf_eq( a.t , b.t , gmp_float::gmp_needequal_bits ); |
---|
| 118 | } |
---|
| 119 | bool operator > ( const gmp_float & a, const gmp_float & b ) |
---|
| 120 | { |
---|
| 121 | return mpf_cmp( a.t, b.t ) > 0; |
---|
| 122 | } |
---|
| 123 | bool operator < ( const gmp_float & a, const gmp_float & b ) |
---|
| 124 | { |
---|
| 125 | return mpf_cmp( a.t, b.t ) < 0; |
---|
| 126 | } |
---|
| 127 | bool operator >= ( const gmp_float & a, const gmp_float & b ) |
---|
| 128 | { |
---|
| 129 | return mpf_cmp( a.t, b.t ) >= 0; |
---|
| 130 | } |
---|
| 131 | bool operator <= ( const gmp_float & a, const gmp_float & b ) |
---|
| 132 | { |
---|
| 133 | return mpf_cmp( a.t, b.t ) <= 0; |
---|
| 134 | } |
---|
| 135 | |
---|
| 136 | // unary - |
---|
| 137 | gmp_float operator - ( const gmp_float & a ) |
---|
| 138 | { |
---|
| 139 | gmp_float tmp; |
---|
| 140 | mpf_neg( *(tmp.mpfp()), *(a.mpfp()) ); |
---|
| 141 | return tmp; |
---|
| 142 | } |
---|
| 143 | |
---|
| 144 | gmp_float abs( const gmp_float & a ) |
---|
| 145 | { |
---|
| 146 | gmp_float *tmp= new gmp_float(); |
---|
| 147 | mpf_abs( *tmp->mpfp(), *a.mpfp() ); |
---|
| 148 | return *tmp; |
---|
| 149 | } |
---|
| 150 | gmp_float sqrt( const gmp_float & a ) |
---|
| 151 | { |
---|
| 152 | gmp_float *tmp= new gmp_float(); |
---|
| 153 | mpf_sqrt( *tmp->mpfp(), *a.mpfp() ); |
---|
| 154 | return *tmp; |
---|
| 155 | } |
---|
| 156 | gmp_float sin( const gmp_float & a ) |
---|
| 157 | { |
---|
| 158 | gmp_float *tmp= new gmp_float( sin((double)a) ); |
---|
| 159 | return *tmp; |
---|
| 160 | } |
---|
| 161 | gmp_float cos( const gmp_float & a ) |
---|
| 162 | { |
---|
| 163 | gmp_float *tmp= new gmp_float( cos((double)a) ); |
---|
| 164 | return *tmp; |
---|
| 165 | } |
---|
| 166 | gmp_float log( const gmp_float & a ) |
---|
| 167 | { |
---|
| 168 | gmp_float *tmp= new gmp_float( log((double)a) ); |
---|
| 169 | return *tmp; |
---|
| 170 | } |
---|
| 171 | gmp_float hypot( const gmp_float & a, const gmp_float & b ) |
---|
| 172 | { |
---|
| 173 | #if 1 |
---|
| 174 | gmp_float *tmp= new gmp_float(); |
---|
| 175 | *tmp= sqrt( (a*a) + (b*b) ); |
---|
| 176 | return *tmp; |
---|
| 177 | #else |
---|
| 178 | gmp_float *tmp= new gmp_float( hypot( (double)a, (double)b ) ); |
---|
| 179 | return *tmp; |
---|
| 180 | #endif |
---|
| 181 | } |
---|
| 182 | gmp_float exp( const gmp_float & a ) |
---|
| 183 | { |
---|
| 184 | gmp_float *tmp= new gmp_float( exp((double)a) ); |
---|
| 185 | return *tmp; |
---|
| 186 | } |
---|
| 187 | gmp_float max( const gmp_float & a, const gmp_float & b ) |
---|
| 188 | { |
---|
| 189 | gmp_float *tmp= new gmp_float(); |
---|
| 190 | a > b ? *tmp= a : *tmp= b; |
---|
| 191 | return *tmp; |
---|
| 192 | } |
---|
| 193 | // |
---|
[f543de1] | 194 | // number to float, number = Q, R, C |
---|
| 195 | // makes a COPY of num! (Ist das gut?) |
---|
[8a150b] | 196 | // |
---|
| 197 | gmp_float numberToFloat( number num ) |
---|
| 198 | { |
---|
| 199 | gmp_float r; |
---|
| 200 | |
---|
[f543de1] | 201 | if ( rField_is_Q() ) { |
---|
| 202 | if ( num != NULL ) |
---|
[0ae5116] | 203 | { |
---|
[f543de1] | 204 | if (SR_HDL(num) & SR_INT) |
---|
| 205 | { |
---|
| 206 | r= SR_TO_INT(num); |
---|
| 207 | } |
---|
| 208 | else |
---|
| 209 | { |
---|
| 210 | if ( num->s == 0 ) |
---|
| 211 | { |
---|
| 212 | nlNormalize( num ); |
---|
| 213 | } |
---|
| 214 | if (SR_HDL(num) & SR_INT) |
---|
| 215 | { |
---|
| 216 | r= SR_TO_INT(num); |
---|
| 217 | } |
---|
| 218 | else |
---|
| 219 | { |
---|
| 220 | if ( num->s != 3 ) |
---|
| 221 | { |
---|
| 222 | r= &num->z; |
---|
| 223 | r/= (gmp_float)&num->n; |
---|
| 224 | } |
---|
| 225 | else |
---|
| 226 | { |
---|
| 227 | r= &num->z; |
---|
| 228 | } |
---|
| 229 | } |
---|
| 230 | } |
---|
[0ae5116] | 231 | } |
---|
[f543de1] | 232 | else |
---|
[0ae5116] | 233 | { |
---|
[f543de1] | 234 | r= 0.0; |
---|
[8a150b] | 235 | } |
---|
[f543de1] | 236 | } else if (rField_is_long_R() || rField_is_long_C()) { |
---|
| 237 | r= *(gmp_float*)num; |
---|
| 238 | } else if ( rField_is_R() ) { |
---|
| 239 | // Add some code here :-) |
---|
| 240 | Werror("Wrong field!"); |
---|
| 241 | } else { |
---|
| 242 | Werror("Wrong field!"); |
---|
[0ae5116] | 243 | } |
---|
[8a150b] | 244 | |
---|
| 245 | return r; |
---|
| 246 | } |
---|
| 247 | |
---|
| 248 | // Do some strange things with the mantissa string and the exponent |
---|
| 249 | // to get some nice output string. |
---|
| 250 | char *nicifyFloatStr( char * in, mp_exp_t exponent, size_t oprec, int *size, int thesign ) |
---|
| 251 | { |
---|
| 252 | char *out; |
---|
| 253 | |
---|
| 254 | int sign= (in[0] == '-') ? 1 : 0; |
---|
| 255 | char csign[2]; |
---|
| 256 | |
---|
| 257 | switch (thesign) |
---|
| 258 | { |
---|
| 259 | case SIGN_PLUS: |
---|
[f543de1] | 260 | sign ? strcpy(csign,"-") : strcpy(csign,"+"); //+123, -123 |
---|
[8a150b] | 261 | break; |
---|
| 262 | case SIGN_SPACE: |
---|
[f543de1] | 263 | sign ? strcpy(csign,"-") : strcpy(csign," "); // 123, -123 |
---|
[8a150b] | 264 | break; |
---|
| 265 | case SIGN_EMPTY: |
---|
| 266 | default: |
---|
[f543de1] | 267 | sign ? strcpy(csign,"-") : strcpy(csign,""); //123, -123 |
---|
[8a150b] | 268 | break; |
---|
| 269 | } |
---|
| 270 | |
---|
| 271 | if ( strlen(in) == 0 ) |
---|
| 272 | { |
---|
[f543de1] | 273 | *size= 2*sizeof(char)+10; |
---|
| 274 | out= (char*)AllocL( *size ); |
---|
[8a150b] | 275 | strcpy(out,"0"); |
---|
| 276 | return out; |
---|
| 277 | } |
---|
| 278 | |
---|
| 279 | if ( ((unsigned int)abs(exponent) <= oprec) |
---|
| 280 | /*|| (exponent+sign >= (int)strlen(in))*/ ) |
---|
| 281 | { |
---|
| 282 | if ( exponent+sign < (int)strlen(in) ) |
---|
| 283 | { |
---|
| 284 | int eexponent= (exponent >= 0) ? 0 : -exponent; |
---|
| 285 | int eeexponent= (exponent >= 0) ? exponent : 0; |
---|
[f543de1] | 286 | *size= (strlen(in)+5+eexponent) * sizeof(char) + 10; |
---|
| 287 | out= (char*)AllocL(*size); |
---|
| 288 | memset(out,'\0',*size); |
---|
[8a150b] | 289 | |
---|
| 290 | strcpy(out,csign); |
---|
| 291 | strncat(out,in+sign,eeexponent); |
---|
| 292 | |
---|
| 293 | if (exponent == 0) |
---|
| 294 | strcat(out,"0."); |
---|
| 295 | else if ( exponent > 0 ) |
---|
| 296 | strcat(out,"."); |
---|
| 297 | else |
---|
| 298 | { |
---|
| 299 | strcat(out,"0."); |
---|
| 300 | memset(out+strlen(out),'0',eexponent); |
---|
| 301 | } |
---|
| 302 | strcat(out,in+sign+eeexponent); |
---|
| 303 | } |
---|
| 304 | else if ( exponent+sign > (int)strlen(in) ) |
---|
| 305 | { |
---|
[f543de1] | 306 | *size= (strlen(in)+exponent+2)*sizeof(char)+10; |
---|
| 307 | out= (char*)AllocL(*size); |
---|
[8a150b] | 308 | sprintf(out,"%s%s",csign,in+sign); |
---|
| 309 | memset(out+strlen(out),'0',exponent-strlen(in)+sign); |
---|
| 310 | } |
---|
| 311 | else |
---|
| 312 | { |
---|
[f543de1] | 313 | *size= (strlen(in)+2) * sizeof(char) + 10; |
---|
| 314 | out= (char*)AllocL(*size); |
---|
[8a150b] | 315 | sprintf(out,"%s%s",csign,in+sign); |
---|
| 316 | } |
---|
| 317 | } |
---|
| 318 | else |
---|
| 319 | { |
---|
[f543de1] | 320 | // if ( exponent > 0 ) |
---|
| 321 | // { |
---|
[8a150b] | 322 | int c=1,d=10; |
---|
| 323 | while ( exponent / d > 0 ) |
---|
| 324 | { // count digits |
---|
| 325 | d*=10; |
---|
| 326 | c++; |
---|
| 327 | } |
---|
[f543de1] | 328 | *size= (strlen(in)+12+c) * sizeof(char) + 10; |
---|
| 329 | out= (char*)AllocL(*size); |
---|
[8a150b] | 330 | sprintf(out,"%s0.%se%d",csign,in+sign,(unsigned int)exponent); |
---|
[f543de1] | 331 | // } |
---|
| 332 | // else |
---|
| 333 | // { |
---|
| 334 | // *size=2; |
---|
| 335 | // out= (char*)AllocL(*size); |
---|
| 336 | // strcpy(out,"0"); |
---|
| 337 | // } |
---|
[8a150b] | 338 | } |
---|
| 339 | return out; |
---|
| 340 | } |
---|
| 341 | |
---|
| 342 | char *floatToStr( const gmp_float & r, const size_t oprec ) |
---|
| 343 | { |
---|
| 344 | #if 1 |
---|
| 345 | mp_exp_t exponent; |
---|
| 346 | int size,insize; |
---|
| 347 | char *nout,*out,*in; |
---|
| 348 | |
---|
[f543de1] | 349 | insize= (oprec+2) * sizeof(char) + 10; |
---|
| 350 | in= (char*)AllocL( insize ); |
---|
[8a150b] | 351 | |
---|
| 352 | mpf_get_str(in,&exponent,10,oprec,*(r.mpfp())); |
---|
[f543de1] | 353 | |
---|
[0ae5116] | 354 | if ( (exponent > 0) |
---|
| 355 | && (exponent < (int)oprec) |
---|
| 356 | && (strlen(in)-(in[0]=='-'?1:0) == oprec) ) |
---|
| 357 | { |
---|
[f543de1] | 358 | FreeL( (ADDRESS) in ); |
---|
| 359 | insize= (exponent+oprec+2) * sizeof(char) + 10; |
---|
| 360 | in= (char*)AllocL( insize ); |
---|
[8a150b] | 361 | int newprec= exponent+oprec; |
---|
| 362 | mpf_get_str(in,&exponent,10,newprec,*(r.mpfp())); |
---|
| 363 | } |
---|
| 364 | nout= nicifyFloatStr( in, exponent, oprec, &size, SIGN_EMPTY ); |
---|
[f543de1] | 365 | FreeL( (ADDRESS) in ); |
---|
| 366 | out= (char*)AllocL( (strlen(nout)+1) * sizeof(char) ); |
---|
[8a150b] | 367 | strcpy( out, nout ); |
---|
[f543de1] | 368 | FreeL( (ADDRESS) nout ); |
---|
| 369 | |
---|
[8a150b] | 370 | return out; |
---|
| 371 | #else |
---|
[f543de1] | 372 | // for testing purpose... |
---|
| 373 | char *out= (char*)AllocL( (1024) * sizeof(char) ); |
---|
[8a150b] | 374 | sprintf(out,"% .10f",(double)r); |
---|
| 375 | return out; |
---|
| 376 | #endif |
---|
| 377 | } |
---|
| 378 | //<- |
---|
| 379 | |
---|
[f543de1] | 380 | //-> gmp_complex::* |
---|
| 381 | // <gmp_complex> = <gmp_complex> operator <gmp_complex> |
---|
[8a150b] | 382 | // |
---|
[f543de1] | 383 | gmp_complex operator + ( const gmp_complex & a, const gmp_complex & b ) |
---|
[8a150b] | 384 | { |
---|
[f543de1] | 385 | return gmp_complex( a.r + b.r, a.i + b.i ); |
---|
[8a150b] | 386 | } |
---|
[f543de1] | 387 | gmp_complex operator - ( const gmp_complex & a, const gmp_complex & b ) |
---|
[8a150b] | 388 | { |
---|
[f543de1] | 389 | return gmp_complex( a.r - b.r, a.i - b.i ); |
---|
[8a150b] | 390 | } |
---|
[f543de1] | 391 | gmp_complex operator * ( const gmp_complex & a, const gmp_complex & b ) |
---|
[8a150b] | 392 | { |
---|
[f543de1] | 393 | return gmp_complex( a.real() * b.real() - a.imag() * b.imag(), |
---|
[8a150b] | 394 | a.real() * b.imag() + a.imag() * b.real()); |
---|
| 395 | } |
---|
[f543de1] | 396 | gmp_complex operator / ( const gmp_complex & a, const gmp_complex & b ) |
---|
[8a150b] | 397 | { |
---|
[f543de1] | 398 | gmp_float ar = abs(b.real()); |
---|
| 399 | gmp_float ai = abs(b.imag()); |
---|
| 400 | gmp_float nr, ni, t, d; |
---|
[0ae5116] | 401 | if (ar <= ai) |
---|
| 402 | { |
---|
[8a150b] | 403 | t = b.real() / b.imag(); |
---|
[f543de1] | 404 | d = b.imag() * ((gmp_float)1 + t*t); |
---|
[8a150b] | 405 | nr = (a.real() * t + a.imag()) / d; |
---|
| 406 | ni = (a.imag() * t - a.real()) / d; |
---|
[0ae5116] | 407 | } |
---|
| 408 | else |
---|
| 409 | { |
---|
[8a150b] | 410 | t = b.imag() / b.real(); |
---|
[f543de1] | 411 | d = b.real() * ((gmp_float)1 + t*t); |
---|
[8a150b] | 412 | nr = (a.real() + a.imag() * t) / d; |
---|
| 413 | ni = (a.imag() - a.real() * t) / d; |
---|
| 414 | } |
---|
[f543de1] | 415 | return gmp_complex( nr, ni ); |
---|
[8a150b] | 416 | } |
---|
| 417 | |
---|
[f543de1] | 418 | // <gmp_complex> operator <gmp_complex> |
---|
[8a150b] | 419 | // |
---|
[f543de1] | 420 | gmp_complex & gmp_complex::operator += ( const gmp_complex & b ) |
---|
[8a150b] | 421 | { |
---|
| 422 | r+=b.r; |
---|
| 423 | i+=b.i; |
---|
| 424 | return *this; |
---|
| 425 | } |
---|
[f543de1] | 426 | gmp_complex & gmp_complex::operator -= ( const gmp_complex & b ) |
---|
[8a150b] | 427 | { |
---|
| 428 | r-=b.r; |
---|
| 429 | i-=b.i; |
---|
| 430 | return *this; |
---|
| 431 | } |
---|
[f543de1] | 432 | gmp_complex & gmp_complex::operator *= ( const gmp_complex & b ) |
---|
[8a150b] | 433 | { |
---|
[f543de1] | 434 | gmp_float f = r * b.r - i * b.i; |
---|
[8a150b] | 435 | i = r * b.i + i * b.r; |
---|
| 436 | r = f; |
---|
| 437 | return *this; |
---|
| 438 | } |
---|
[f543de1] | 439 | gmp_complex & gmp_complex::operator /= ( const gmp_complex & b ) |
---|
[8a150b] | 440 | { |
---|
[f543de1] | 441 | gmp_float ar = abs(b.r); |
---|
| 442 | gmp_float ai = abs(b.i); |
---|
| 443 | gmp_float nr, ni, t, d; |
---|
[0ae5116] | 444 | if (ar <= ai) |
---|
| 445 | { |
---|
[8a150b] | 446 | t = b.r / b.i; |
---|
[f543de1] | 447 | d = b.i * ((gmp_float)1 + t*t); |
---|
[8a150b] | 448 | nr = (r * t + i) / d; |
---|
| 449 | ni = (i * t - r) / d; |
---|
[0ae5116] | 450 | } |
---|
| 451 | else |
---|
| 452 | { |
---|
[8a150b] | 453 | t = b.i / b.r; |
---|
[f543de1] | 454 | d = b.r * ((gmp_float)1 + t*t); |
---|
[8a150b] | 455 | nr = (r + i * t) / d; |
---|
| 456 | ni = (i - r * t) / d; |
---|
| 457 | } |
---|
| 458 | r = nr; |
---|
| 459 | i = ni; |
---|
| 460 | return *this; |
---|
| 461 | } |
---|
| 462 | |
---|
[f543de1] | 463 | // Returns square root of gmp_complex number |
---|
[8a150b] | 464 | // |
---|
[f543de1] | 465 | gmp_complex sqrt( const gmp_complex & x ) |
---|
[8a150b] | 466 | { |
---|
[f543de1] | 467 | gmp_float r = abs(x); |
---|
| 468 | gmp_float nr, ni; |
---|
| 469 | if (r == (gmp_float) 0.0) |
---|
[0ae5116] | 470 | { |
---|
[8a150b] | 471 | nr = ni = r; |
---|
[0ae5116] | 472 | } |
---|
[f543de1] | 473 | else if ( x.real() > (gmp_float)0) |
---|
[0ae5116] | 474 | { |
---|
[f543de1] | 475 | nr = sqrt((gmp_float)0.5 * (r + x.real())); |
---|
| 476 | ni = x.imag() / nr / (gmp_float)2; |
---|
[0ae5116] | 477 | } |
---|
| 478 | else |
---|
| 479 | { |
---|
[f543de1] | 480 | ni = sqrt((gmp_float)0.5 * (r - x.real())); |
---|
| 481 | if (x.imag() < (gmp_float)0) |
---|
[0ae5116] | 482 | { |
---|
[8a150b] | 483 | ni = - ni; |
---|
| 484 | } |
---|
[f543de1] | 485 | nr = x.imag() / ni / (gmp_float)2; |
---|
[8a150b] | 486 | } |
---|
[f543de1] | 487 | gmp_complex *tmp= new gmp_complex(nr, ni); |
---|
[8a150b] | 488 | return *tmp; |
---|
| 489 | } |
---|
| 490 | |
---|
[f543de1] | 491 | // converts a gmp_complex to a string ( <real part> + I * <imaginary part> ) |
---|
[8a150b] | 492 | // |
---|
[f543de1] | 493 | char *complexToStr( const gmp_complex & c, const size_t oprec ) |
---|
[8a150b] | 494 | { |
---|
| 495 | char *out,*in_imag,*in_real; |
---|
| 496 | |
---|
[0ae5116] | 497 | if ( !c.imag().isZero() ) |
---|
| 498 | { |
---|
[8a150b] | 499 | in_real=floatToStr( c.real(), oprec ); // get real part |
---|
| 500 | in_imag=floatToStr( abs(c.imag()), oprec ); // get imaginary part |
---|
[f543de1] | 501 | |
---|
| 502 | if (rField_is_long_C()) { |
---|
| 503 | out=(char*)AllocL((strlen(in_real)+strlen(in_imag)+5+strlen(currRing->parameter[0]))*sizeof(char)); |
---|
| 504 | sprintf(out,"%s%s%s*%s",in_real,c.imag().sign()>=0?"+":"-",currRing->parameter[0],in_imag); |
---|
| 505 | } else { |
---|
| 506 | out=(char*)AllocL( (strlen(in_real)+strlen(in_imag)+8) * sizeof(char)); |
---|
| 507 | sprintf(out,"%s%s%s",in_real,c.imag().sign()>=0?" + I ":" - I ",in_imag); |
---|
| 508 | } |
---|
| 509 | FreeL( (ADDRESS) in_real ); |
---|
| 510 | FreeL( (ADDRESS) in_imag ); |
---|
[0ae5116] | 511 | } |
---|
| 512 | else |
---|
| 513 | { |
---|
[8a150b] | 514 | out= floatToStr( c.real(), oprec ); |
---|
| 515 | } |
---|
| 516 | return out; |
---|
| 517 | } |
---|
| 518 | //<- |
---|
[f543de1] | 519 | //%e |
---|
[8a150b] | 520 | |
---|
| 521 | //#endif // HAVE_MPR |
---|
| 522 | |
---|
| 523 | // local Variables: *** |
---|
| 524 | // folded-file: t *** |
---|
| 525 | // compile-command-1: "make installg" *** |
---|
| 526 | // compile-command-2: "make install" *** |
---|
| 527 | // End: *** |
---|