mpr_complex.cc File Reference

#include "misc/auxiliary.h"
#include "reporter/reporter.h"
#include "coeffs/coeffs.h"
#include "coeffs/numbers.h"
#include "coeffs/mpr_complex.h"
#include "coeffs/longrat.h"
#include <cmath>

Go to the source code of this file.

Macros
#define	SR_HDL(A)   ((long)(A))
#define	SR_TO_INT(SR)   (((long)SR) >> 2)
#define	SIGN_PLUS   1
#define	SIGN_SPACE   2
#define	SIGN_EMPTY   4
#define	EXTRABYTES   4
#define	DEFPREC   20

Functions
void	setGMPFloatDigits (size_t digits, size_t rest)
	Set size of mantissa digits - the number of output digits (basis 10) the size of mantissa consists of two parts: the "output" part a and the "rest" part b. More...
gmp_float	operator+ (const gmp_float &a, const gmp_float &b)
gmp_float	operator- (const gmp_float &a, const gmp_float &b)
gmp_float	operator* (const gmp_float &a, const gmp_float &b)
gmp_float	operator/ (const gmp_float &a, const gmp_float &b)
bool	operator== (const gmp_float &a, const gmp_float &b)
bool	operator> (const gmp_float &a, const gmp_float &b)
bool	operator< (const gmp_float &a, const gmp_float &b)
bool	operator>= (const gmp_float &a, const gmp_float &b)
bool	operator<= (const gmp_float &a, const gmp_float &b)
gmp_float	operator- (const gmp_float &a)
gmp_float	abs (const gmp_float &a)
gmp_float	sqrt (const gmp_float &a)
gmp_float	sin (const gmp_float &a)
gmp_float	cos (const gmp_float &a)
gmp_float	log (const gmp_float &a)
gmp_float	hypot (const gmp_float &a, const gmp_float &b)
gmp_float	exp (const gmp_float &a)
gmp_float	max (const gmp_float &a, const gmp_float &b)
gmp_float	numberToFloat (number num, const coeffs src)
gmp_float	numberFieldToFloat (number num, int cf)
char *	nicifyFloatStr (char *in, mp_exp_t exponent, size_t oprec, int *size, int thesign)
char *	floatToStr (const gmp_float &r, const unsigned int oprec)
gmp_complex	operator+ (const gmp_complex &a, const gmp_complex &b)
gmp_complex	operator- (const gmp_complex &a, const gmp_complex &b)
gmp_complex	operator* (const gmp_complex &a, const gmp_complex &b)
gmp_complex	operator/ (const gmp_complex &a, const gmp_complex &b)
gmp_complex	sqrt (const gmp_complex &x)
char *	complexToStr (gmp_complex &c, const unsigned int oprec, const coeffs src)
bool	complexNearZero (gmp_complex *c, int digits)

Variables
VAR size_t	gmp_output_digits = DEFPREC
STATIC_VAR gmp_float *	gmpRel =NULL
STATIC_VAR gmp_float *	diff =NULL

Macro Definition Documentation

DEFPREC

#define DEFPREC   20

Definition at line 41 of file mpr_complex.cc.

EXTRABYTES

#define EXTRABYTES   4

Definition at line 39 of file mpr_complex.cc.

SIGN_EMPTY

#define SIGN_EMPTY   4

Definition at line 37 of file mpr_complex.cc.

SIGN_PLUS

#define SIGN_PLUS   1

Definition at line 35 of file mpr_complex.cc.

SIGN_SPACE

#define SIGN_SPACE   2

Definition at line 36 of file mpr_complex.cc.

SR_HDL

#define SR_HDL

(

A

)

((long)(A))

Definition at line 32 of file mpr_complex.cc.

SR_TO_INT

#define SR_TO_INT

(

SR

)

(((long)SR) >> 2)

Definition at line 33 of file mpr_complex.cc.

Function Documentation

abs()

gmp_float abs

(

const gmp_float &

a

)

Definition at line 321 of file mpr_complex.cc.

{
  gmp_float tmp;
  mpf_abs( *(tmp._mpfp()), *a.mpfp() );
  return tmp;
}

complexNearZero()

bool complexNearZero	(	gmp_complex *	c,
		int	digits
	)

Definition at line 765 of file mpr_complex.cc.

{
  gmp_float eps,epsm;
 
  if ( digits < 1 ) return true;
 
  eps=pow(10.0,(int)digits);
  //Print("eps: %s\n",floatToStr(eps,gmp_output_digits));
  eps=(gmp_float)1.0/eps;
  epsm=-eps;
 
  //Print("eps: %s\n",floatToStr(eps,gmp_output_digits));
 
  if ( c->real().sign() > 0 ) // +
    return (c->real() < eps && (c->imag() < eps && c->imag() > epsm));
  else // -
    return (c->real() > epsm && (c->imag() < eps && c->imag() > epsm));
}

complexToStr()

char * complexToStr	(	gmp_complex &	c,
		const unsigned int	oprec,
		const coeffs	src
	)

Definition at line 704 of file mpr_complex.cc.

{
  const char * complex_parameter = "I";
  int N = 1; // strlen(complex_parameter);
 
  if (nCoeff_is_long_C(src))
  {
    complex_parameter = n_ParameterNames(src)[0];
    N = strlen(complex_parameter);
  }
 
  assume( complex_parameter != NULL && N > 0);
 
  char *out,*in_imag,*in_real;
 
  c.SmallToZero();
  if ( !c.imag().isZero() )
  {
 
    in_real=floatToStr( c.real(), oprec );         // get real part
    in_imag=floatToStr( abs(c.imag()), oprec );    // get imaginary part
 
    if (nCoeff_is_long_C(src))
    {
      int len=(strlen(in_real)+strlen(in_imag)+7+N)*sizeof(char);
      out=(char*)omAlloc(len);
      memset(out,0,len);
      if (  !c.real().isZero() )  // (-23-i*5.43) or (15.1+i*5.3)
        sprintf(out,"(%s%s%s*%s)",in_real,c.imag().sign()>=0?"+":"-",complex_parameter,in_imag);
      else // (-i*43) or (i*34)
      {
        if (c.imag().isOne())
          sprintf(out,"%s", complex_parameter);
        else if (c.imag().isMOne())
          sprintf(out,"-%s", complex_parameter);
        else
          sprintf(out,"(%s%s*%s)",c.imag().sign()>=0?"":"-", complex_parameter,in_imag);
      }
    }
    else
    {
      int len=(strlen(in_real)+strlen(in_imag)+9) * sizeof(char);
      out=(char*)omAlloc( len );
      memset(out,0,len);
      if ( !c.real().isZero() )
        sprintf(out,"(%s%s%s)",in_real,c.imag().sign()>=0?"+I*":"-I*",in_imag);
      else
        sprintf(out,"(%s%s)",c.imag().sign()>=0?"I*":"-I*",in_imag);
    }
    omFree( (void *) in_real );
    omFree( (void *) in_imag );
  }
  else
  {
    out= floatToStr( c.real(), oprec );
  }
 
  return out;
}

cos()

gmp_float cos

(

const gmp_float &

a

)

Definition at line 338 of file mpr_complex.cc.

{
  gmp_float tmp( cos((double)a) );
  return tmp;
}

exp()

gmp_float exp

(

const gmp_float &

a

)

Definition at line 357 of file mpr_complex.cc.

{
  gmp_float tmp( exp((double)a) );
  return tmp;
}

floatToStr()

char * floatToStr	(	const gmp_float &	r,
		const unsigned int	oprec
	)

Definition at line 578 of file mpr_complex.cc.

{
#if 1
  mp_exp_t exponent;
  int size,insize;
  char *nout,*out,*in;
 
  insize= (oprec+2) * sizeof(char) + 10;
  in= (char*)omAlloc( insize );
 
  mpf_get_str(in,&exponent,10,oprec,*(r.mpfp()));
 
  //if ( (exponent > 0)
  //&& (exponent < (int)oprec)
  //&& (strlen(in)-(in[0]=='-'?1:0) == oprec) )
  //{
  //  omFree( (void *) in );
  //  insize= (exponent+oprec+2) * sizeof(char) + 10;
  //  in= (char*)omAlloc( insize );
  //  int newprec= exponent+oprec;
  //  mpf_get_str(in,&exponent,10,newprec,*(r.mpfp()));
  //}
  nout= nicifyFloatStr( in, exponent, oprec, &size, SIGN_EMPTY );
  omFree( (void *) in );
  out= (char*)omAlloc( (strlen(nout)+1) * sizeof(char) );
  strcpy( out, nout );
  omFree( (void *) nout );
 
  return out;
#else
  // for testing purpose...
  char *out= (char*)omAlloc( (1024) * sizeof(char) );
  sprintf(out,"% .10f",(double)r);
  return out;
#endif
}

hypot()

gmp_float hypot	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 348 of file mpr_complex.cc.

{
#if 1
  return ( sqrt( (a*a) + (b*b) ) );
#else
  gmp_float tmp( hypot( (double)a, (double)b ) );
  return tmp;
#endif
}

log()

gmp_float log

(

const gmp_float &

a

)

Definition at line 343 of file mpr_complex.cc.

{
  gmp_float tmp( log((double)a) );
  return tmp;
}

max()

gmp_float max	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 362 of file mpr_complex.cc.

{
  gmp_float tmp;
  a > b ? tmp= a : tmp= b;
  return tmp;
}

nicifyFloatStr()

char * nicifyFloatStr	(	char *	in,
		mp_exp_t	exponent,
		size_t	oprec,
		int *	size,
		int	thesign
	)

Definition at line 485 of file mpr_complex.cc.

{
  char *out;
 
  int sign= (in[0] == '-') ? 1 : 0;
  char csign[2];
 
  switch (thesign)
  {
    case SIGN_PLUS:
      sign ? strcpy(csign,"-") : strcpy(csign,"+");  //+123, -123
      break;
    case SIGN_SPACE:
      sign ? strcpy(csign,"-") : strcpy(csign," ");  // 123, -123
      break;
    case SIGN_EMPTY:
    default:
      sign ? strcpy(csign,"-") : strcpy(csign,"");   //123, -123
      break;
  }
 
  if ( strlen(in) == 0 )
  {
    *size= 2*sizeof(char);
    return omStrDup("0");
  }
 
  if ( ((unsigned int)ABS(exponent) <= oprec)
       /*|| (exponent+sign >= (int)strlen(in))*/ )
  {
    if ( exponent+sign < (int)strlen(in) )
    {
      int eexponent= (exponent >= 0) ? 0 : -exponent;
      int eeexponent= (exponent >= 0) ? exponent : 0;
      *size= (strlen(in)+15+eexponent) * sizeof(char);
      out= (char*)omAlloc(*size);
      memset(out,0,*size);
 
      strcpy(out,csign);
      strncat(out,in+sign,eeexponent);
 
      if (exponent == 0)
        strcat(out,"0.");
      else if ( exponent > 0 )
        strcat(out,".");
      else
      {
        strcat(out,"0.");
        memset(out+strlen(out),'0',eexponent);
      }
      strcat(out,in+sign+eeexponent);
    }
    else if ( exponent+sign > (int)strlen(in) )
    {
      *size= (strlen(in)+exponent+12)*sizeof(char);
      out= (char*)omAlloc(*size);
      memset(out,0,*size);
      sprintf(out,"%s%s",csign,in+sign);
      memset(out+strlen(out),'0',exponent-strlen(in)+sign);
    }
    else
    {
      *size= (strlen(in)+2) * sizeof(char) + 10;
      out= (char*)omAlloc(*size);
      memset(out,0,*size);
      sprintf(out,"%s%s",csign,in+sign);
    }
  }
  else
  {
//      if ( exponent > 0 )
//      {
      int c=1,d=10;
      while ( exponent / d > 0 )
      { // count digits
        d*=10;
        c++;
      }
      *size= (strlen(in)+12+c) * sizeof(char) + 10;
      out= (char*)omAlloc(*size);
      memset(out,0,*size);
      sprintf(out,"%s0.%se%s%d",csign,in+sign,exponent>=0?"+":"",(int)exponent);
//      }
//      else
//      {
//        *size=2;
//        out= (char*)omAlloc(*size);
//        strcpy(out,"0");
//      }
  }
  return out;
}

numberFieldToFloat()

gmp_float numberFieldToFloat	(	number	num,
		int	cf
	)

Definition at line 438 of file mpr_complex.cc.

{
  gmp_float r;
 
  switch (cf)
  {
  case QTOF:
    if ( num != NULL )
    {
      if (SR_HDL(num) & SR_INT)
      {
        r = gmp_float(SR_TO_INT(num));
      }
      else
      {
        if ( num->s != 3 )
        {
          r= gmp_float(num->z);
          r/= gmp_float(num->n);
        }
        else
        {
          r= num->z;
        }
      }
    }
    else
    {
      r= 0.0;
    }
    break;
  case RTOF:
    r= *(gmp_float*)num;
    break;
  case CTOF:
    WerrorS("Can not map from field C to field R!");
    break;
  case ZTOF:
  default:
    WerrorS("Ground field not implemented!");
  } // switch
 
  return r;
}

numberToFloat()

gmp_float numberToFloat	(	number	num,
		const coeffs	src
	)

Definition at line 372 of file mpr_complex.cc.

{
  gmp_float r;
 
  if ( nCoeff_is_Q(src) )
  {
    if ( num != NULL )
    {
      if (SR_HDL(num) & SR_INT)
      {
        //n_Print(num, src);printf("\n");
        int nn = SR_TO_INT(num);
        if((long)nn == SR_TO_INT(num))
            r = SR_TO_INT(num);
        else
            r = gmp_float(SR_TO_INT(num));
        //int dd = 20;
        //gmp_printf("\nr = %.*Ff\n",dd,*r.mpfp());
        //getchar();
      }
      else
      {
        if ( num->s == 0 )
        {
          nlNormalize( num, src ); // FIXME? TODO? // extern void     nlNormalize(number &x, const coeffs r); // FIXME
        }
        if (SR_HDL(num) & SR_INT)
        {
          r= SR_TO_INT(num);
        }
        else
        {
          if ( num->s != 3 )
          {
            r= num->z;
            r/= (gmp_float)num->n;
          }
          else
          {
            r= num->z;
          }
        }
      }
    }
    else
    {
      r= 0.0;
    }
  }
  else if (nCoeff_is_long_R(src) || nCoeff_is_long_C(src))
  {
    r= *(gmp_float*)num;
  }
  else if ( nCoeff_is_R(src) )
  {
    // Add some code here :-)
    WerrorS("Ground field not implemented!");
  }
  else
  {
    WerrorS("Ground field not implemented!");
  }
 
  return r;
}

operator*() [1/2]

gmp_complex operator*	(	const gmp_complex &	a,
		const gmp_complex &	b
	)

Definition at line 627 of file mpr_complex.cc.

{
  return gmp_complex( a.r * b.r - a.i * b.i,
                  a.r * b.i + a.i * b.r);
}

operator*() [2/2]

gmp_float operator*	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 179 of file mpr_complex.cc.

{
  gmp_float tmp( a );
  tmp *= b;
  return tmp;
}

operator+() [1/2]

gmp_complex operator+	(	const gmp_complex &	a,
		const gmp_complex &	b
	)

Definition at line 619 of file mpr_complex.cc.

{
  return gmp_complex( a.r + b.r, a.i + b.i );
}

operator+() [2/2]

gmp_float operator+	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 167 of file mpr_complex.cc.

{
  gmp_float tmp( a );
  tmp += b;
  return tmp;
}

operator-() [1/3]

gmp_complex operator-	(	const gmp_complex &	a,
		const gmp_complex &	b
	)

Definition at line 623 of file mpr_complex.cc.

{
  return gmp_complex( a.r - b.r, a.i - b.i );
}

operator-() [2/3]

gmp_float operator-

(

const gmp_float &

a

)

Definition at line 314 of file mpr_complex.cc.

{
  gmp_float tmp;
  mpf_neg( *(tmp._mpfp()), *(a.mpfp()) );
  return tmp;
}

operator-() [3/3]

gmp_float operator-	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 173 of file mpr_complex.cc.

{
  gmp_float tmp( a );
  tmp -= b;
  return tmp;
}

operator/() [1/2]

gmp_complex operator/	(	const gmp_complex &	a,
		const gmp_complex &	b
	)

Definition at line 632 of file mpr_complex.cc.

{
  gmp_float d = b.r*b.r + b.i*b.i;
  return gmp_complex( (a.r * b.r + a.i * b.i) / d,
                (a.i * b.r - a.r * b.i) / d);
}

operator/() [2/2]

gmp_float operator/	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 185 of file mpr_complex.cc.

{
  gmp_float tmp( a );
  tmp /= b;
  return tmp;
}

operator<()

bool operator<	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 294 of file mpr_complex.cc.

{
  if (a.t == b.t)
    return false;
  return mpf_cmp( a.t, b.t ) < 0;
}

operator<=()

bool operator<=	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 306 of file mpr_complex.cc.

{
  if (a.t == b.t)
    return true;
  return mpf_cmp( a.t, b.t ) <= 0;
}

operator==()

bool operator==	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 237 of file mpr_complex.cc.

{
  if(mpf_sgn(a.t) != mpf_sgn(b.t))
    return false;
  if((mpf_sgn(a.t)==0) && (mpf_sgn(b.t)==0))
    return true;
  mpf_sub(diff->t, a.t, b.t);
  mpf_div(diff->t, diff->t, a.t);
  mpf_abs(diff->t, diff->t);
  if(mpf_cmp(diff->t, gmpRel->t) < 0)
    return true;
  else
    return false;
}

operator>()

bool operator>	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 288 of file mpr_complex.cc.

{
  if (a.t == b.t)
    return false;
  return mpf_cmp( a.t, b.t ) > 0;
}

operator>=()

bool operator>=	(	const gmp_float &	a,
		const gmp_float &	b
	)

Definition at line 300 of file mpr_complex.cc.

{
  if (a.t == b.t)
    return true;
  return mpf_cmp( a.t, b.t ) >= 0;
}

setGMPFloatDigits()

void setGMPFloatDigits	(	size_t	digits,
		size_t	rest
	)

Set size of mantissa digits - the number of output digits (basis 10) the size of mantissa consists of two parts: the "output" part a and the "rest" part b.

According to the GMP-precision digits is recomputed to bits (basis 2). Two numbers a, b are equal if | a - b | < | a | * 0.1^digits . In this case we have a - b = 0 . The epsilon e is e=0.1^(digits+rest) with 1+e != 1, but 1+0.1*e = 1.

Definition at line 60 of file mpr_complex.cc.

{
  size_t bits = 1 + (size_t) ((float)digits * 3.5);
  size_t rb = 1 + (size_t) ((float)rest * 3.5);
  size_t db = bits+rb;
  gmp_output_digits= digits;
  mpf_set_default_prec( db );
  if (diff!=NULL) delete diff;
  diff=new gmp_float(0.0);
  mpf_set_prec(*diff->_mpfp(),32);
  if (gmpRel!=NULL) delete gmpRel;
  gmpRel=new gmp_float(0.0);
  mpf_set_prec(*gmpRel->_mpfp(),32);
  mpf_set_d(*gmpRel->_mpfp(),0.1);
  mpf_pow_ui(*gmpRel->_mpfp(),*gmpRel->_mpfp(),digits);
}

sin()

gmp_float sin

(

const gmp_float &

a

)

Definition at line 333 of file mpr_complex.cc.

{
  gmp_float tmp( sin((double)a) );
  return tmp;
}

sqrt() [1/2]

gmp_complex sqrt

(

const gmp_complex &

x

)

Definition at line 676 of file mpr_complex.cc.

{
  gmp_float r = abs(x);
  gmp_float nr, ni;
  if (r == (gmp_float) 0.0)
  {
    nr = ni = r;
  }
  else if ( x.real() > (gmp_float)0)
  {
    nr = sqrt((gmp_float)0.5 * (r + x.real()));
    ni = x.imag() / nr / (gmp_float)2;
  }
  else
  {
    ni = sqrt((gmp_float)0.5 * (r - x.real()));
    if (x.imag() < (gmp_float)0)
    {
      ni = - ni;
    }
    nr = x.imag() / ni / (gmp_float)2;
  }
  gmp_complex tmp(nr, ni);
  return tmp;
}

sqrt() [2/2]

gmp_float sqrt

(

const gmp_float &

a

)

Definition at line 327 of file mpr_complex.cc.

{
  gmp_float tmp;
  mpf_sqrt( *(tmp._mpfp()), *a.mpfp() );
  return tmp;
}

Variable Documentation

diff

STATIC_VAR gmp_float* diff =NULL

Definition at line 45 of file mpr_complex.cc.

gmp_output_digits

VAR size_t gmp_output_digits = DEFPREC

Definition at line 42 of file mpr_complex.cc.

gmpRel

STATIC_VAR gmp_float* gmpRel =NULL

Definition at line 44 of file mpr_complex.cc.