source: git/libpolys/coeffs/mpr_complex.cc @ 06df101

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