source: git/kernel/mpr_complex.cc @ 3a67ea7

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