source: git/kernel/mpr_complex.cc @ 599326

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