source: git/libpolys/coeffs/mpr_complex.cc @ 0acf3e

jengelh-datetimespielwiese
Last change on this file since 0acf3e was ba5e9e, checked in by Oleksandr Motsak <motsak@…>, 9 years ago
Changed configure-scripts to generate individual public config files for each package: resources, libpolys, singular (main) fix: sources should include correct corresponding config headers.
  • Property mode set to 100644
File size: 17.6 KB
Line 
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
13#ifdef HAVE_CONFIG_H
14#include "libpolysconfig.h"
15#endif /* HAVE_CONFIG_H */
16//#ifdef HAVE_MPR
17#include <coeffs/coeffs.h>
18#include <reporter/reporter.h>
19#include <omalloc/omalloc.h>
20#include <coeffs/numbers.h>
21#include <coeffs/longrat.h>
22#include <math.h>
23#include <coeffs/mpr_complex.h>
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)
39size_t gmp_output_digits= DEFPREC;
40
41static gmp_float *gmpRel=NULL;
42static gmp_float *diff=NULL;
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 */
57void 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 );
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);
72}
73
74size_t getGMPFloatDigits()
75{
76  return gmp_output_digits;
77}
78
79#if 1
80void gmp_float::setFromStr(const char * in )
81{
82  BOOLEAN neg=false;
83  if (*in == '-') { in++; neg=TRUE; }
84  char *s;
85  if ((s=strchr((char *)in,'E')) !=NULL)
86  {
87    *s='e';
88  }
89
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  }
106  if (neg)  mpf_neg( t, t );
107}
108#else
109// problemns with solve_s.tst
110void 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
166
167
168// <gmp_float> = <gmp_float> operator <gmp_float>
169gmp_float operator + ( const gmp_float & a, const gmp_float & b )
170{
171  gmp_float tmp( a );
172  tmp += b;
173  return tmp;
174}
175gmp_float operator - ( const gmp_float & a, const gmp_float & b )
176{
177  gmp_float tmp( a );
178  tmp -= b;
179  return tmp;
180}
181gmp_float operator * ( const gmp_float & a, const gmp_float & b )
182{
183  gmp_float tmp( a );
184  tmp *= b;
185  return tmp;
186}
187gmp_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>
195gmp_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 );
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)
213    mpf_set_d( t, 0.0);
214  return *this;
215}
216gmp_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 );
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)
234    mpf_set_d( t, 0.0);
235  return *this;
236}
237
238// <gmp_float> == <gmp_float> ??
239bool 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;
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)
249    return true;
250  else
251    return false;
252}
253// t == 0 ?
254bool gmp_float::isZero() const
255{
256  return (mpf_sgn( t ) == 0);
257}
258// t == 1 ?
259bool gmp_float::isOne() const
260{
261#ifdef  VARIANTE_1
262  return (mpf_cmp_ui( t , 1 ) == 0);
263#else
264  if (mpf_sgn(t) <= 0)
265    return false;
266  mpf_sub_ui(diff->t, t, 1);
267  mpf_abs(diff->t, diff->t);
268  if(mpf_cmp(diff->t, gmpRel->t) < 0)
269    return true;
270  else
271    return false;
272#endif
273}
274// t == -1 ?
275bool gmp_float::isMOne() const
276{
277#ifdef VARIANTE_1
278  return (mpf_cmp_si( t , -1 ) == 0);
279#else
280  if (mpf_sgn(t) >= 0)
281    return false;
282  mpf_add_ui(diff->t, t, 1);
283  mpf_abs(diff->t, diff->t);
284  if(mpf_cmp(diff->t, gmpRel->t) < 0)
285    return true;
286  else
287    return false;
288#endif
289}
290bool 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}
296bool 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}
302bool 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}
308bool 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 -
316gmp_float operator - ( const gmp_float & a )
317{
318  gmp_float tmp;
319  mpf_neg( *(tmp._mpfp()), *(a.mpfp()) );
320  return tmp;
321}
322
323gmp_float abs( const gmp_float & a )
324{
325  gmp_float tmp;
326  mpf_abs( *(tmp._mpfp()), *a.mpfp() );
327  return tmp;
328}
329gmp_float sqrt( const gmp_float & a )
330{
331  gmp_float tmp;
332  mpf_sqrt( *(tmp._mpfp()), *a.mpfp() );
333  return tmp;
334}
335gmp_float sin( const gmp_float & a )
336{
337  gmp_float tmp( sin((double)a) );
338  return tmp;
339}
340gmp_float cos( const gmp_float & a )
341{
342  gmp_float tmp( cos((double)a) );
343  return tmp;
344}
345gmp_float log( const gmp_float & a )
346{
347  gmp_float tmp( log((double)a) );
348  return tmp;
349}
350gmp_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}
359gmp_float exp( const gmp_float & a )
360{
361  gmp_float tmp( exp((double)a) );
362  return tmp;
363}
364gmp_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//
374gmp_float numberToFloat( number num, const coeffs src)
375{
376  gmp_float r;
377
378  if ( nCoeff_is_Q(src) )
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        {
390          nlNormalize( num, src );
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          {
400            r= num->z;
401            r/= (gmp_float)num->n;
402          }
403          else
404          {
405            r= num->z;
406          }
407        }
408      }
409    }
410    else
411    {
412      r= 0.0;
413    }
414  }
415  else if (nCoeff_is_long_R(src) || nCoeff_is_long_C(src))
416  {
417    r= *(gmp_float*)num;
418  }
419  else if ( nCoeff_is_R(src) )
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
432gmp_float numberFieldToFloat( number num, int k, const coeffs src)
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        {
449          nlNormalize( num, src );
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          {
459            r= num->z;
460            r/= (gmp_float)num->n;
461          }
462          else
463          {
464            r= num->z;
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.
490char *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
583char *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  {
599    omFree( (void *) in );
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 );
606  omFree( (void *) in );
607  out= (char*)omAlloc( (strlen(nout)+1) * sizeof(char) );
608  strcpy( out, nout );
609  omFree( (void *) nout );
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//
624gmp_complex operator + ( const gmp_complex & a, const gmp_complex & b )
625{
626  return gmp_complex( a.r + b.r, a.i + b.i );
627}
628gmp_complex operator - ( const gmp_complex & a, const gmp_complex & b )
629{
630  return gmp_complex( a.r - b.r, a.i - b.i );
631}
632gmp_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}
637gmp_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//
646gmp_complex & gmp_complex::operator += ( const gmp_complex & b )
647{
648  r+=b.r;
649  i+=b.i;
650  return *this;
651}
652gmp_complex & gmp_complex::operator -= ( const gmp_complex & b )
653{
654  r-=b.r;
655  i-=b.i;
656  return *this;
657}
658gmp_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}
665gmp_complex & gmp_complex::neg ( )
666{
667  i.neg();
668  r.neg();
669  return *this;
670}
671gmp_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//
681gmp_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  }
703  gmp_complex tmp(nr, ni);
704  return tmp;
705}
706
707// converts a gmp_complex to a string ( <real part> + I * <imaginary part> )
708//
709char *complexToStr( gmp_complex & c, const unsigned int oprec, const coeffs src )
710{
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);
721 
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
731    if (nCoeff_is_long_C(src))
732    {
733      int len=(strlen(in_real)+strlen(in_imag)+7+N)*sizeof(char);
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)
737        sprintf(out,"(%s%s%s*%s)",in_real,c.imag().sign()>=0?"+":"-",complex_parameter,in_imag);
738      else // (-i*43) or (i*34)
739      {
740        if (c.imag().isOne())
741          sprintf(out,"%s", complex_parameter);
742        else if (c.imag().isMOne())
743          sprintf(out,"-%s", complex_parameter);
744        else
745          sprintf(out,"(%s%s*%s)",c.imag().sign()>=0?"":"-", complex_parameter,in_imag);
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    }
758    omFree( (void *) in_real );
759    omFree( (void *) in_imag );
760  }
761  else
762  {
763    out= floatToStr( c.real(), oprec );
764  }
765
766  return out;
767}
768//<-
769
770bool 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
789void 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());
801    if (ai < *gmpRel) this->imag(0.0);
802  }
803  else
804  {
805    mpf_div(*ar._mpfp(), *ar._mpfp(), *ai._mpfp());
806    if (ar < *gmpRel) this->real(0.0);
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: ***
Note: See TracBrowser for help on using the repository browser.