source: git/libpolys/coeffs/mpr_complex.cc @ 560a3d

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