Changeset f543de1 in git


Ignore:
Timestamp:
Jun 24, 1999, 9:46:51 AM (24 years ago)
Author:
Moritz Wenk <wenk@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
e1eec6fcd177391b1ba1d2a6030f298f404ae5f9
Parents:
6b9e9c662b8c35931dfedb7053c0c6d2597e5f68
Message:
*wenk: complex -> gnu_complex


git-svn-id: file:///usr/local/Singular/svn/trunk@3168 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • Singular/gnumpc.h

    r6b9e9c rf543de1  
    44*  Computer Algebra System SINGULAR     *
    55****************************************/
    6 /* $Id: gnumpc.h,v 1.1 1999-05-11 15:42:37 Singular Exp $ */
     6/* $Id: gnumpc.h,v 1.2 1999-06-24 07:46:49 wenk Exp $ */
    77/*
    88* ABSTRACT: computations with GMP floating-point numbers

  • Singular/gnumpfl.h

    r6b9e9c rf543de1  
    44*  Computer Algebra System SINGULAR     *
    55****************************************/
    6 /* $Id: gnumpfl.h,v 1.2 1999-05-10 15:10:50 Singular Exp $ */
     6/* $Id: gnumpfl.h,v 1.3 1999-06-24 07:46:49 wenk Exp $ */
    77/*
    88* ABSTRACT: computations with GMP floating-point numbers
    … …  
    4040BOOLEAN  ngfSetMap(int c, char ** par, int nop, number minpol);
    4141
    42 
     42void setGMPFloatDigits( size_t digits );
    4343void setGMPFloatPrecBytes( unsigned long int bytes );
    4444#endif

  • Singular/mpr_complex.cc

    r6b9e9c rf543de1  
    22*  Computer Algebra System SINGULAR     *
    33****************************************/
    4 /* $Id: mpr_complex.cc,v 1.4 1999-05-17 11:31:47 Singular Exp $ */
     4/* $Id: mpr_complex.cc,v 1.5 1999-06-24 07:46:50 wenk Exp $ */
    55
    66/*
    … …  
    99*
    1010*/
     11
     12// WARNING! ALWAYS use AllocL and FreeL when alloc. memory for some char* !!
    1113
    1214#include "mod2.h"
    … …  
    2022#include "mpr_complex.h"
    2123
     24//%s
    2225// this was copied form longrat0.cc
     26// and will be used in numberToFloat.
     27// Make sure that it is up to date!!
    2328#define SR_HDL(A) ((long)(A))
    2429#define SR_INT    1
    … …  
    3338size_t gmp_output_digits= DEFPREC;
    3439
    35 const gmp_float  gmpOne= 1;
    36 const gmp_float gmpMOne= -1;
    37 const gmp_float gmpZero= 0;
    38 
    39 
     40//-> setGMPFloat*
    4041/** Set size of mantissa to <bytes> bytes and guess the number of
    4142 * output digits.
    … …  
    6364}
    6465
     66// Sets the lenght of the mantissa to <digits> digits
    6567void setGMPFloatDigits( size_t digits )
    6668{
    6769  size_t bits= 1 + (size_t) (digits / (log(2)/log(10)));
    6870  bits= bits>64?bits:64;
    69   gmp_float::setPrecision( bits+EXTRABYTES*8 );
     71  //  gmp_float::setPrecision( bits+EXTRABYTES*8 );
     72  gmp_float::setPrecision( bits+(bits/2) );
    7073  gmp_float::setEqualBits( bits );
    7174  gmp_output_digits= digits;
    … …  
    7679  return gmp_output_digits;
    7780}
    78 
    79 //------------------------------------- gmp_float ----------------------------------------------
     81//<-
    8082
    8183//-> gmp_float::*
    … …  
    8385unsigned long int gmp_float::gmp_needequal_bits= GMP_NEEDEQUAL_BITS;
    8486
    85 gmp_float::gmp_float( const int v )
    86 {
    87   mpf_init2( t, gmp_default_prec_bits );
    88   mpf_set_si( t, (signed long int) v );
    89 }
    90 gmp_float::gmp_float( const long v )
    91 {
    92   mpf_init2( t, gmp_default_prec_bits );
    93   mpf_set_si( t, (signed long int) v );
    94 }
    95 gmp_float::gmp_float( const mprfloat v )
    96 {
    97   mpf_init2( t, gmp_default_prec_bits );
    98   mpf_set_d( t, (double) v );
    99 }
    100 gmp_float::gmp_float( const mpf_t v )
    101 {
    102   mpf_init2( t, gmp_default_prec_bits );
    103   mpf_set( t, v );
    104 }
    105 gmp_float::gmp_float( const mpz_t v )
    106 {
    107   mpf_init2( t, gmp_default_prec_bits );
    108   mpf_set_z( t, v );
    109 }
    110 gmp_float::gmp_float( const gmp_float & v )
    111 {
    112   //mpf_init2( t, mpf_get_prec( v.t ) );
    113   mpf_init2( t, gmp_default_prec_bits );
    114   mpf_set( t, v.t );
    115 }
    116 
    117 gmp_float::~gmp_float()
    118 {
    119   mpf_clear( t );
    120 }
    121 
    12287// <gmp_float> = <gmp_float> operator <gmp_float>
    12388gmp_float operator + ( const gmp_float & a, const gmp_float & b )
    … …  
    146111}
    147112
    148 // <gmp_float> operator <gmp_float>
    149 gmp_float & gmp_float::operator += ( const gmp_float & a )
    150 {
    151   mpf_add( t, t, a.t );
    152   return *this;
    153 }
    154 gmp_float & gmp_float::operator -= ( const gmp_float & a )
    155 {
    156   mpf_sub( t, t, a.t );
    157   return *this;
    158 }
    159 gmp_float & gmp_float::operator *= ( const gmp_float & a )
    160 {
    161   mpf_mul( t, t, a.t );
    162   return *this;
    163 }
    164 gmp_float & gmp_float::operator /= ( const gmp_float & a )
    165 {
    166   mpf_div( t, t, a.t );
    167   return *this;
    168 }
    169 
    170113// <gmp_float> == <gmp_float> ?? up to the first gmp_float::gmp_needequal_bits bits
    171114bool operator == ( const gmp_float & a, const gmp_float & b )
    … …  
    199142}
    200143
    201 // <gmp_float> = <*>
    202 gmp_float & gmp_float::operator = ( const gmp_float & a )
    203 {
    204   mpf_set( t, a.t );
    205   return *this;
    206 }
    207 gmp_float & gmp_float::operator = ( const mpz_t & a )
    208 {
    209   mpf_set_z( t, a );
    210   return *this;
    211 }
    212 gmp_float & gmp_float::operator = ( const mprfloat a )
    213 {
    214   mpf_set_d( t, (double) a );
    215   return *this;
    216 }
    217 gmp_float & gmp_float::operator = ( const long a )
    218 {
    219   mpf_set_si( t, (signed long int) a );
    220   return *this;
    221 }
    222 
    223 // cast to double
    224 gmp_float::operator double()
    225 {
    226   return mpf_get_d( t );
    227 }
    228 gmp_float::operator double() const
    229 {
    230   return mpf_get_d( t );
    231 }
    232 
    233 // cast to int
    234 gmp_float::operator int()
    235 {
    236   return (int)mpf_get_d( t );
    237 }
    238 gmp_float::operator int() const
    239 {
    240   return (int)mpf_get_d( t );
    241 }
    242 
    243 // get sign of real number ( -1: t < 0; 0: t==0; 1: t > 0 )
    244 int gmp_float::sign()
    245 {
    246   return mpf_sgn( t );
    247 }
    248 // t == 0 ?
    249 bool gmp_float::isZero()
    250 {
    251 #ifdef  VARIANTE_1
    252   return (mpf_sgn( t ) == 0);
    253 #else
    254   return  mpf_eq( t , gmpZero.t , gmp_float::gmp_needequal_bits );
    255 #endif
    256 }
    257 // t == 1 ?
    258 bool gmp_float::isOne()
    259 {
    260 #ifdef  VARIANTE_1
    261   return (mpf_cmp_ui( t , 1 ) == 0);
    262 #else
    263   return mpf_eq( t , gmpOne.t , gmp_float::gmp_needequal_bits );
    264 #endif
    265 }
    266 // t == -1 ?
    267 bool gmp_float::isMOne()
    268 {
    269 #ifdef VARIANTE_1
    270   return (mpf_cmp_si( t , -1 ) == 0);
    271 #else
    272   return mpf_eq( t , gmpMOne.t , gmp_float::gmp_needequal_bits );
    273 #endif
    274 }
    275 
    276 bool gmp_float::setFromStr( char * in )
    277 {
    278   return ( mpf_set_str( t, in, 10 ) == 0 );
    279 }
    280 
    281 // access pointer
    282 const mpf_t *gmp_float::mpfp() const
    283 {
    284   return &t;
    285 }
    286 
    287144gmp_float abs( const gmp_float & a )
    288145{
    … …  
    335192}
    336193//
    337 // WARNING:
    338 // supports only Q to float !!!
     194// number to float, number = Q, R, C
     195// makes a COPY of num! (Ist das gut?)
    339196//
    340197gmp_float numberToFloat( number num )
    … …  
    342199  gmp_float r;
    343200
    344   //Print("numberToFloat: ");nPrint(num);
    345 
    346   if ( num != NULL )
    347   {
    348     if (SR_HDL(num) & SR_INT)
    349     {
    350       r= SR_TO_INT(num);
    351     }
     201  if ( rField_is_Q() ) {
     202    if ( num != NULL )
     203      {
     204        if (SR_HDL(num) & SR_INT)
     205          {
     206            r= SR_TO_INT(num);
     207          }
     208        else
     209          {
     210            if ( num->s == 0 )
     211              {
     212                nlNormalize( num );
     213              }
     214            if (SR_HDL(num) & SR_INT)
     215              {
     216                r= SR_TO_INT(num);
     217              }
     218            else
     219              {
     220                if ( num->s != 3 )
     221                  {
     222                    r= &num->z;
     223                    r/= (gmp_float)&num->n;
     224                  }
     225                else
     226                  {
     227                    r= &num->z;
     228                  }
     229              }
     230          }
     231      }
    352232    else
    353     {
    354       if ( num->s == 0 )
    355233      {
    356         nlNormalize( num );
     234        r= 0.0;
    357235      }
    358       if (SR_HDL(num) & SR_INT)
    359       {
    360         r= SR_TO_INT(num);
    361       }
    362       else
    363       {
    364         if ( num->s != 3 )
    365         {
    366           r= &num->z;
    367           r/= (mprfloat_g)&num->n;
    368         }
    369         else
    370         {
    371           r= &num->z;
    372         }
    373       }
    374     }
    375   }
    376   else
    377   {
    378     r= 0.0;
    379   }
    380 
    381   //Print(" --> %s ",floatToStr(r,10));PrintLn();
     236  } else if (rField_is_long_R() || rField_is_long_C()) {
     237    r= *(gmp_float*)num;
     238  } else if ( rField_is_R() ) {
     239    // Add some code here :-)
     240    Werror("Wrong field!");
     241  } else {
     242    Werror("Wrong field!");
     243  }
    382244
    383245  return r;
    … …  
    396258  {
    397259    case SIGN_PLUS:
    398       sign ? strcpy(csign,"-") : strcpy(csign,"+");
     260      sign ? strcpy(csign,"-") : strcpy(csign,"+");  //+123, -123
    399261      break;
    400262    case SIGN_SPACE:
    401       sign ? strcpy(csign,"-") : strcpy(csign," ");
     263      sign ? strcpy(csign,"-") : strcpy(csign," ");  // 123, -123
    402264      break;
    403265    case SIGN_EMPTY:
    404266    default:
    405       sign ? strcpy(csign,"-") : strcpy(csign,"");
     267      sign ? strcpy(csign,"-") : strcpy(csign,"");   //123, -123
    406268      break;
    407269  }
    … …  
    409271  if ( strlen(in) == 0 )
    410272  {
    411     out= (char*)Alloc0( 2*sizeof(char) );
    412     *size= 2*sizeof(char);
     273    *size= 2*sizeof(char)+10;
     274    out= (char*)AllocL( *size );
    413275    strcpy(out,"0");
    414276    return out;
    … …  
    418280       /*|| (exponent+sign >= (int)strlen(in))*/ )
    419281  {
    420 #ifdef mprDEBUG_ALL
    421     Print(" no exponent %d %d\n",abs(exponent),oprec);
    422 #endif
    423282    if ( exponent+sign < (int)strlen(in) )
    424283    {
    425284      int eexponent= (exponent >= 0) ? 0 : -exponent;
    426285      int eeexponent= (exponent >= 0) ? exponent : 0;
    427       *size= (strlen(in)+5+eexponent) * sizeof(char);
    428       out= (char*)Alloc0(*size);
     286      *size= (strlen(in)+5+eexponent) * sizeof(char) + 10;
     287      out= (char*)AllocL(*size);
     288      memset(out,'\0',*size);
    429289
    430290      strcpy(out,csign);
    … …  
    444304    else if ( exponent+sign > (int)strlen(in) )
    445305    {
    446       *size= (strlen(in)+exponent+2)*sizeof(char);
    447       out= (char*)Alloc0(*size);
     306      *size= (strlen(in)+exponent+2)*sizeof(char)+10;
     307      out= (char*)AllocL(*size);
    448308      sprintf(out,"%s%s",csign,in+sign);
    449309      memset(out+strlen(out),'0',exponent-strlen(in)+sign);
    … …  
    451311    else
    452312    {
    453       *size= (strlen(in)+2) * sizeof(char);
    454       out= (char*)Alloc0(*size);
     313      *size= (strlen(in)+2) * sizeof(char) + 10;
     314      out= (char*)AllocL(*size);
    455315      sprintf(out,"%s%s",csign,in+sign);
    456316    }
    … …  
    458318  else
    459319  {
    460 #ifdef mprDEBUG_ALL
    461     Print(" exponent %d %d\n",exponent,oprec);
    462 #endif
    463     if ( exponent > 0 )
    464     {
     320//      if ( exponent > 0 )
     321//      {
    465322      int c=1,d=10;
    466323      while ( exponent / d > 0 )
    … …  
    469326        c++;
    470327      }
    471       *size= (strlen(in)+12+c) * sizeof(char);
    472       out= (char*)Alloc0(*size);
     328      *size= (strlen(in)+12+c) * sizeof(char) + 10;
     329      out= (char*)AllocL(*size);
    473330      sprintf(out,"%s0.%se%d",csign,in+sign,(unsigned int)exponent);
    474     }
    475     else
    476     {
    477       *size=2;
    478       out= (char*)Alloc0(*size);
    479       strcpy(out,"0");
    480     }
     331//      }
     332//      else
     333//      {
     334//        *size=2;
     335//        out= (char*)AllocL(*size);
     336//        strcpy(out,"0");
     337//      }
    481338  }
    482339  return out;
    … …  
    490347  char *nout,*out,*in;
    491348
    492   insize= (oprec+2) * sizeof(char);
    493   in= (char*)Alloc0( insize );
     349  insize= (oprec+2) * sizeof(char) + 10;
     350  in= (char*)AllocL( insize );
    494351
    495352  mpf_get_str(in,&exponent,10,oprec,*(r.mpfp()));
     353
    496354  if ( (exponent > 0)
    497355  && (exponent < (int)oprec)
    498356  && (strlen(in)-(in[0]=='-'?1:0) == oprec) )
    499357  {
    500     in= (char*)ReAlloc( in, insize, (exponent+oprec+2) * sizeof(char) );
    501     insize= (exponent+oprec+2) * sizeof(char);
     358    FreeL( (ADDRESS) in );
     359    insize= (exponent+oprec+2) * sizeof(char) + 10;
     360    in= (char*)AllocL( insize );
    502361    int newprec= exponent+oprec;
    503362    mpf_get_str(in,&exponent,10,newprec,*(r.mpfp()));
    504363  }
    505364  nout= nicifyFloatStr( in, exponent, oprec, &size, SIGN_EMPTY );
    506   Free( (ADDRESS) in, insize );
    507   out= (char*)Alloc0( (strlen(nout)+1) * sizeof(char) );
     365  FreeL( (ADDRESS) in );
     366  out= (char*)AllocL( (strlen(nout)+1) * sizeof(char) );
    508367  strcpy( out, nout );
    509   Free( (ADDRESS) nout, size );
     368  FreeL( (ADDRESS) nout );
     369
    510370  return out;
    511371#else
    512   char *out= (char*)Alloc0( (1024) * sizeof(char) );
     372  // for testing purpose...
     373  char *out= (char*)AllocL( (1024) * sizeof(char) );
    513374  sprintf(out,"% .10f",(double)r);
    514375  return out;
    … …  
    517378//<-
    518379
    519 //------------------------------------- complex ------------------------------------------------
    520 
    521 //-> complex::*
    522 // constructors
     380//-> gmp_complex::*
     381// <gmp_complex> = <gmp_complex> operator <gmp_complex>
    523382//
    524 complex::complex( const mprfloat_g re, const mprfloat_g im )
    525 {
    526   r= re;
    527   i= im;
    528 }
    529 complex::complex( const mprfloat re, const mprfloat im )
    530 {
    531   r= re;
    532   i= im;
    533 }
    534 complex::complex( const long re, const long im )
    535 {
    536   r= re;
    537   i= im;
    538 }
    539 complex::complex( const complex & v )
    540 {
    541   r= v.r;
    542   i= v.i;
    543 }
    544 complex::~complex()
    545 {
    546 }
    547 
    548 // <complex> = <complex> operator <complex>
    549 //
    550 complex operator + ( const complex & a, const complex & b )
    551 {
    552   return complex( a.r + b.r, a.i + b.i );
    553 }
    554 complex operator - ( const complex & a, const complex & b )
    555 {
    556   return complex( a.r - b.r, a.i - b.i );
    557 }
    558 complex operator * ( const complex & a, const complex & b )
    559 {
    560   return complex( a.real() * b.real() - a.imag() * b.imag(),
     383gmp_complex operator + ( const gmp_complex & a, const gmp_complex & b )
     384{
     385  return gmp_complex( a.r + b.r, a.i + b.i );
     386}
     387gmp_complex operator - ( const gmp_complex & a, const gmp_complex & b )
     388{
     389  return gmp_complex( a.r - b.r, a.i - b.i );
     390}
     391gmp_complex operator * ( const gmp_complex & a, const gmp_complex & b )
     392{
     393  return gmp_complex( a.real() * b.real() - a.imag() * b.imag(),
    561394                  a.real() * b.imag() + a.imag() * b.real());
    562395}
    563 complex operator / ( const complex & a, const complex & b )
    564 {
    565   mprfloat_g ar = abs(b.real());
    566   mprfloat_g ai = abs(b.imag());
    567   mprfloat_g nr, ni, t, d;
     396gmp_complex operator / ( const gmp_complex & a, const gmp_complex & b )
     397{
     398  gmp_float ar = abs(b.real());
     399  gmp_float ai = abs(b.imag());
     400  gmp_float nr, ni, t, d;
    568401  if (ar <= ai)
    569402  {
    570403    t = b.real() / b.imag();
    571     d = b.imag() * ((mprfloat_g)1 + t*t);
     404    d = b.imag() * ((gmp_float)1 + t*t);
    572405    nr = (a.real() * t + a.imag()) / d;
    573406    ni = (a.imag() * t - a.real()) / d;
    … …  
    576409  {
    577410    t = b.imag() / b.real();
    578     d = b.real() * ((mprfloat_g)1 + t*t);
     411    d = b.real() * ((gmp_float)1 + t*t);
    579412    nr = (a.real() + a.imag() * t) / d;
    580413    ni = (a.imag() - a.real() * t) / d;
    581414  }
    582   return complex( nr, ni );
    583 }
    584 
    585 // <complex> = <complex> operator <mprfloat_g>
     415  return gmp_complex( nr, ni );
     416}
     417
     418// <gmp_complex> operator <gmp_complex>
    586419//
    587 complex operator + ( const complex & a, const mprfloat_g b_d )
    588 {
    589   return complex( a.r + b_d, a.i );
    590 }
    591 complex operator - ( const complex & a, const mprfloat_g b_d )
    592 {
    593   return complex( a.r - b_d, a.i );
    594 }
    595 complex operator * ( const complex & a, const mprfloat_g b_d )
    596 {
    597   return complex( a.r * b_d, a.i * b_d );
    598 }
    599 complex operator / ( const complex & a, const mprfloat_g b_d )
    600 {
    601   return complex( a.r / b_d, a.i / b_d );
    602 }
    603 
    604 // <complex> operator <complex>
    605 //
    606 complex & complex::operator += ( const complex & b )
     420gmp_complex & gmp_complex::operator += ( const gmp_complex & b )
    607421{
    608422  r+=b.r;
    … …  
    610424  return *this;
    611425}
    612 complex & complex::operator -= ( const complex & b )
     426gmp_complex & gmp_complex::operator -= ( const gmp_complex & b )
    613427{
    614428  r-=b.r;
    … …  
    616430  return *this;
    617431}
    618 complex & complex::operator *= ( const complex & b )
    619 {
    620   mprfloat_g f = r * b.r - i * b.i;
     432gmp_complex & gmp_complex::operator *= ( const gmp_complex & b )
     433{
     434  gmp_float f = r * b.r - i * b.i;
    621435  i = r * b.i + i * b.r;
    622436  r = f;
    623437  return *this;
    624438}
    625 complex & complex::operator /= ( const complex & b )
    626 {
    627   mprfloat_g ar = abs(b.r);
    628   mprfloat_g ai = abs(b.i);
    629   mprfloat_g nr, ni, t, d;
     439gmp_complex & gmp_complex::operator /= ( const gmp_complex & b )
     440{
     441  gmp_float ar = abs(b.r);
     442  gmp_float ai = abs(b.i);
     443  gmp_float nr, ni, t, d;
    630444  if (ar <= ai)
    631445  {
    632446    t = b.r / b.i;
    633     d = b.i * ((mprfloat_g)1 + t*t);
     447    d = b.i * ((gmp_float)1 + t*t);
    634448    nr = (r * t + i) / d;
    635449    ni = (i * t - r) / d;
    … …  
    638452  {
    639453    t = b.i / b.r;
    640     d = b.r * ((mprfloat_g)1 + t*t);
     454    d = b.r * ((gmp_float)1 + t*t);
    641455    nr = (r + i * t) / d;
    642456    ni = (i - r * t) / d;
    … …  
    647461}
    648462
    649 // <complex> == <complex> ?
    650 bool operator == ( const complex & a, const complex & b )
    651 {
    652   return ( b.real() == a.real() ) && ( b.imag() == a.imag() );
    653 }
    654 
    655 // <complex> = <complex>
    656 complex & complex::operator = ( const complex & a )
    657 {
    658   r= a.r;
    659   i= a.i;
    660   return *this;
    661 }
    662 
    663 // Returns absolute value of a complex number
     463// Returns square root of gmp_complex number
    664464//
    665 mprfloat_g abs( const complex & c )
    666 {
    667   return hypot(c.real(),c.imag());
    668 }
    669 
    670 // Returns square root of complex number
    671 //
    672 complex sqrt( const complex & x )
    673 {
    674   mprfloat_g r = abs(x);
    675   mprfloat_g nr, ni;
    676   if (r == (mprfloat_g) 0.0)
     465gmp_complex sqrt( const gmp_complex & x )
     466{
     467  gmp_float r = abs(x);
     468  gmp_float nr, ni;
     469  if (r == (gmp_float) 0.0)
    677470  {
    678471    nr = ni = r;
    679472  }
    680   else if ( x.real() > (mprfloat_g)0)
    681   {
    682     nr = sqrt((mprfloat_g)0.5 * (r + x.real()));
    683     ni = x.imag() / nr / (mprfloat_g)2;
     473  else if ( x.real() > (gmp_float)0)
     474  {
     475    nr = sqrt((gmp_float)0.5 * (r + x.real()));
     476    ni = x.imag() / nr / (gmp_float)2;
    684477  }
    685478  else
    686479  {
    687     ni = sqrt((mprfloat_g)0.5 * (r - x.real()));
    688     if (x.imag() < (mprfloat_g)0)
     480    ni = sqrt((gmp_float)0.5 * (r - x.real()));
     481    if (x.imag() < (gmp_float)0)
    689482    {
    690483      ni = - ni;
    691484    }
    692     nr = x.imag() / ni / (mprfloat_g)2;
    693   }
    694   complex *tmp= new complex(nr, ni);
    695   return *tmp;
    696 }
    697 
    698 // converts a number to a complex
     485    nr = x.imag() / ni / (gmp_float)2;
     486  }
     487  gmp_complex *tmp= new gmp_complex(nr, ni);
     488  return *tmp;
     489}
     490
     491// converts a gmp_complex to a string ( <real part> + I * <imaginary part> )
    699492//
    700 char *complexToStr( const complex & c, const size_t oprec )
     493char *complexToStr( const gmp_complex & c, const size_t oprec )
    701494{
    702495  char *out,*in_imag,*in_real;
    … …  
    706499    in_real=floatToStr( c.real(), oprec );         // get real part
    707500    in_imag=floatToStr( abs(c.imag()), oprec );    // get imaginary part
    708     out= (char*)Alloc0( (strlen(in_real)+strlen(in_imag)+6) * sizeof(char) );
    709     sprintf(out,"%s%s%s",in_real,c.imag().sign() >= 0 ? " + i ":" - i ",in_imag);
    710     Free( in_real, (strlen(in_real)+1)*sizeof(char) );
    711     Free( in_imag, (strlen(in_imag)+1)*sizeof(char) );
     501   
     502    if (rField_is_long_C()) {
     503      out=(char*)AllocL((strlen(in_real)+strlen(in_imag)+5+strlen(currRing->parameter[0]))*sizeof(char));
     504      sprintf(out,"%s%s%s*%s",in_real,c.imag().sign()>=0?"+":"-",currRing->parameter[0],in_imag);
     505    } else {
     506      out=(char*)AllocL( (strlen(in_real)+strlen(in_imag)+8) * sizeof(char));
     507      sprintf(out,"%s%s%s",in_real,c.imag().sign()>=0?" + I ":" - I ",in_imag);
     508    }
     509    FreeL( (ADDRESS) in_real );
     510    FreeL( (ADDRESS) in_imag );
    712511  }
    713512  else
    … …  
    718517}
    719518//<-
     519//%e
    720520
    721521//#endif // HAVE_MPR

  • Singular/mpr_complex.h

    r6b9e9c rf543de1  
    44*  Computer Algebra System SINGULAR     *
    55****************************************/
    6 /* $Id: mpr_complex.h,v 1.3 1999-06-15 16:39:00 Singular Exp $ */
     6/* $Id: mpr_complex.h,v 1.4 1999-06-24 07:46:51 wenk Exp $ */
    77
    88/*
    … …  
    1818}
    1919#include "numbers.h"
     20#include "ring.h"
    2021#include "mpr_global.h"
    2122
    … …  
    2829void setGMPFloatPrecBytes( unsigned long int bytes );
    2930unsigned long int getGMPFloatPrecBytes();
     31
    3032void setGMPFloatDigits( size_t digits );
    3133size_t getGMPFloatDigits();
    … …  
    3840{
    3941public:
    40   gmp_float( const mprfloat v = 0.0 );
    41   gmp_float( const int v );
    42   gmp_float( const long v );
    43   gmp_float( const mpf_t v );
    44   gmp_float( const mpz_t v );
    45   gmp_float( const gmp_float & v );
    46 
    47   ~gmp_float();
     42  gmp_float( const int v = 0 )
     43  {
     44    mpf_init2( t, gmp_default_prec_bits );
     45    mpf_set_si( t, (signed long int) v );
     46  }
     47  gmp_float( const long v )
     48  {
     49    mpf_init2( t, gmp_default_prec_bits );
     50    mpf_set_si( t, (signed long int) v );
     51  }
     52  gmp_float( const mprfloat v ) // double
     53  {
     54    mpf_init2( t, gmp_default_prec_bits );
     55    mpf_set_d( t, (double) v );
     56  }
     57  gmp_float( const mpf_t v )
     58  {
     59    mpf_init2( t, gmp_default_prec_bits );
     60    mpf_set( t, v );
     61  }
     62  gmp_float( const mpz_t v ) // gnu mp Z
     63  {
     64    mpf_init2( t, gmp_default_prec_bits );
     65    mpf_set_z( t, v );
     66  }
     67  gmp_float( const gmp_float & v ) // copy constructor
     68  {
     69    //mpf_init2( t, mpf_get_prec( v.t ) );
     70    mpf_init2( t, gmp_default_prec_bits );
     71    mpf_set( t, v.t );
     72  }
     73
     74  ~gmp_float()
     75  {
     76    mpf_clear( t );
     77  }
    4878
    4979  friend gmp_float operator + ( const gmp_float & a, const gmp_float & b );
    … …  
    5282  friend gmp_float operator / ( const gmp_float & a, const gmp_float & b );
    5383
    54   gmp_float & operator += ( const gmp_float & a );
    55   gmp_float & operator -= ( const gmp_float & a );
    56   gmp_float & operator *= ( const gmp_float & a );
    57   gmp_float & operator /= ( const gmp_float & a );
     84  inline gmp_float & operator += ( const gmp_float & a );
     85  inline gmp_float & operator -= ( const gmp_float & a );
     86  inline gmp_float & operator *= ( const gmp_float & a );
     87  inline gmp_float & operator /= ( const gmp_float & a );
    5888
    5989  friend bool operator == ( const gmp_float & a, const gmp_float & b );
    … …  
    70100  gmp_float & operator = ( const long a );
    71101
    72   int sign();    // t>0:+1, t==0:0, t<0:-1
    73   bool isZero();  // t == 0 ?
    74   bool isOne();   // t == 1 ?
    75   bool isMOne();  // t == -1 ?
    76 
    77   bool setFromStr( char * in );
     102  inline int sign();    // t>0:+1, t==0:0, t<0:-1
     103  inline bool isZero();  // t == 0 ?
     104  inline bool isOne();   // t == 1 ?
     105  inline bool isMOne();  // t == -1 ?
     106
     107  inline bool setFromStr( char * in );
    78108 
    79109  // access
    80110  inline const mpf_t *mpfp() const;
    81111
    82   operator double();
    83   operator double() const;
    84 
    85   operator int();
    86   operator int() const;
     112  inline operator double();
     113  inline operator double() const;
     114
     115  inline operator int();
     116  inline operator int() const;
    87117
    88118public:
    … …  
    104134};
    105135
     136static const gmp_float  gmpOne= 1;
     137static const gmp_float gmpMOne= -1;
     138static const gmp_float gmpZero= 0;
     139
     140// <gmp_float> operator <gmp_float>
     141inline gmp_float & gmp_float::operator += ( const gmp_float & a )
     142{
     143  mpf_add( t, t, a.t );
     144  return *this;
     145}
     146inline gmp_float & gmp_float::operator -= ( const gmp_float & a )
     147{
     148  mpf_sub( t, t, a.t );
     149  return *this;
     150}
     151inline gmp_float & gmp_float::operator *= ( const gmp_float & a )
     152{
     153  mpf_mul( t, t, a.t );
     154  return *this;
     155}
     156inline gmp_float & gmp_float::operator /= ( const gmp_float & a )
     157{
     158  mpf_div( t, t, a.t );
     159  return *this;
     160}
     161
     162// <gmp_float> = <*>
     163inline gmp_float & gmp_float::operator = ( const gmp_float & a )
     164{
     165  mpf_set( t, a.t );
     166  return *this;
     167}
     168inline gmp_float & gmp_float::operator = ( const mpz_t & a )
     169{
     170  mpf_set_z( t, a );
     171  return *this;
     172}
     173inline gmp_float & gmp_float::operator = ( const mprfloat a )
     174{
     175  mpf_set_d( t, (double) a );
     176  return *this;
     177}
     178inline gmp_float & gmp_float::operator = ( const long a )
     179{
     180  mpf_set_si( t, (signed long int) a );
     181  return *this;
     182}
     183
     184// cast to double
     185inline gmp_float::operator double()
     186{
     187  return mpf_get_d( t );
     188}
     189inline gmp_float::operator double() const
     190{
     191  return mpf_get_d( t );
     192}
     193
     194// cast to int
     195inline gmp_float::operator int()
     196{
     197  return (int)mpf_get_d( t );
     198}
     199inline gmp_float::operator int() const
     200{
     201  return (int)mpf_get_d( t );
     202}
     203
     204// get sign of real number ( -1: t < 0; 0: t==0; 1: t > 0 )
     205inline int gmp_float::sign()
     206{
     207  return mpf_sgn( t );
     208}
     209// t == 0 ?
     210inline bool gmp_float::isZero()
     211{
     212#ifdef  VARIANTE_1
     213  return (mpf_sgn( t ) == 0);
     214#else
     215  return  mpf_eq( t , gmpZero.t , gmp_float::gmp_needequal_bits );
     216#endif
     217}
     218// t == 1 ?
     219inline bool gmp_float::isOne()
     220{
     221#ifdef  VARIANTE_1
     222  return (mpf_cmp_ui( t , 1 ) == 0);
     223#else
     224  return mpf_eq( t , gmpOne.t , gmp_float::gmp_needequal_bits );
     225#endif
     226}
     227// t == -1 ?
     228inline bool gmp_float::isMOne()
     229{
     230#ifdef VARIANTE_1
     231  return (mpf_cmp_si( t , -1 ) == 0);
     232#else
     233  return mpf_eq( t , gmpMOne.t , gmp_float::gmp_needequal_bits );
     234#endif
     235}
     236
     237inline bool gmp_float::setFromStr( char * in )
     238{
     239  return ( mpf_set_str( t, in, 10 ) == 0 );
     240}
     241
     242// access pointer
     243inline const mpf_t *gmp_float::mpfp() const
     244{
     245  return &t;
     246}
     247
    106248// built-in functions of GMP
    107249gmp_float abs( const gmp_float & );
    … …  
    118260
    119261gmp_float numberToFloat( number num );
    120 char *floatToStr( const gmp_float & r, const size_t oprec );
    121 
    122 typedef gmp_float mprfloat_g;
     262char *floatToStr( const gmp_float & r, const unsigned int oprec );
    123263//<-
    124264
    125 //-> class complex
     265//-> class gmp_complex
    126266/**
    127  * @short complex numbers based on
     267 * @short gmp_complex numbers based on
    128268 */
    129 class complex
     269class gmp_complex
    130270{
    131271private:
    132   mprfloat_g r, i;
     272  gmp_float r, i;
    133273
    134274public:
    135   complex( const mprfloat_g re= 0.0, const mprfloat_g im= 0.0 );
    136   complex( const mprfloat re, const mprfloat im = 0.0 );
    137   complex( const long re, const long im );
    138   complex( const complex & v );
    139 
    140   ~complex();
    141 
    142   friend complex operator + ( const complex & a, const complex & b );
    143   friend complex operator - ( const complex & a, const complex & b );
    144   friend complex operator * ( const complex & a, const complex & b );
    145   friend complex operator / ( const complex & a, const complex & b );
    146 
    147   friend complex operator + ( const complex & a, const mprfloat_g b_d );
    148   friend complex operator - ( const complex & a, const mprfloat_g b_d );
    149   friend complex operator * ( const complex & a, const mprfloat_g b_d );
    150   friend complex operator / ( const complex & a, const mprfloat_g b_d );
    151 
    152   complex & operator += ( const complex & a );
    153   complex & operator -= ( const complex & a );
    154   complex & operator *= ( const complex & a );
    155   complex & operator /= ( const complex & a );
    156 
    157   friend bool operator == ( const complex & a, const complex & b );
    158 
    159   complex & operator = ( const complex & a );
     275  gmp_complex( const gmp_float re= 0.0, const gmp_float im= 0.0 )
     276  {
     277    r= re;
     278    i= im;
     279  }
     280  gmp_complex( const mprfloat re, const mprfloat im = 0.0 )
     281  {
     282    r= re;
     283    i= im;
     284  }
     285  gmp_complex( const long re, const long im )
     286  {
     287    r= re;
     288    i= im;
     289  }
     290  gmp_complex( const gmp_complex & v )
     291  {
     292    r= v.r;
     293    i= v.i;
     294  }
     295  ~gmp_complex() {}
     296
     297  friend gmp_complex operator + ( const gmp_complex & a, const gmp_complex & b );
     298  friend gmp_complex operator - ( const gmp_complex & a, const gmp_complex & b );
     299  friend gmp_complex operator * ( const gmp_complex & a, const gmp_complex & b );
     300  friend gmp_complex operator / ( const gmp_complex & a, const gmp_complex & b );
     301
     302  // gmp_complex <operator> real
     303  inline friend gmp_complex operator + ( const gmp_complex & a, const gmp_float b_d ); 
     304  inline friend gmp_complex operator - ( const gmp_complex & a, const gmp_float b_d );
     305  inline friend gmp_complex operator * ( const gmp_complex & a, const gmp_float b_d );
     306  inline friend gmp_complex operator / ( const gmp_complex & a, const gmp_float b_d );
     307
     308  gmp_complex & operator += ( const gmp_complex & a );
     309  gmp_complex & operator -= ( const gmp_complex & a );
     310  gmp_complex & operator *= ( const gmp_complex & a );
     311  gmp_complex & operator /= ( const gmp_complex & a );
     312
     313  inline friend bool operator == ( const gmp_complex & a, const gmp_complex & b );
     314
     315  inline gmp_complex & operator = ( const gmp_complex & a );
    160316
    161317  // access to real and imaginary part
    162   inline mprfloat_g real() const { return r; }
    163   inline mprfloat_g imag() const { return i; }
    164 
    165   inline void real( mprfloat_g val ) { r = val; }
    166   inline void imag( mprfloat_g val ) { i = val; }
     318  inline gmp_float real() const { return r; }
     319  inline gmp_float imag() const { return i; }
     320
     321  inline void real( gmp_float val ) { r = val; }
     322  inline void imag( gmp_float val ) { i = val; }
    167323};
    168324
    169 mprfloat_g abs( const complex & c );
    170 complex sqrt( const complex & x );
    171 
    172 inline complex numberToComplex( number num )
    173 {
    174   return complex( numberToFloat(num) );
    175 }
    176 char *complexToStr( const complex & c, const  size_t oprec );
     325// <gmp_complex> = <gmp_complex> operator <gmp_float>
     326//
     327inline gmp_complex operator + ( const gmp_complex & a, const gmp_float b_d )
     328{
     329  return gmp_complex( a.r + b_d, a.i );
     330}
     331inline gmp_complex operator - ( const gmp_complex & a, const gmp_float b_d )
     332{
     333  return gmp_complex( a.r - b_d, a.i );
     334}
     335inline gmp_complex operator * ( const gmp_complex & a, const gmp_float b_d )
     336{
     337  return gmp_complex( a.r * b_d, a.i * b_d );
     338}
     339inline gmp_complex operator / ( const gmp_complex & a, const gmp_float b_d )
     340{
     341  return gmp_complex( a.r / b_d, a.i / b_d );
     342}
     343
     344// <gmp_complex> == <gmp_complex> ?
     345inline bool operator == ( const gmp_complex & a, const gmp_complex & b )
     346{
     347  return ( b.real() == a.real() ) && ( b.imag() == a.imag() );
     348}
     349
     350// <gmp_complex> = <gmp_complex>
     351inline gmp_complex & gmp_complex::operator = ( const gmp_complex & a )
     352{
     353  r= a.r;
     354  i= a.i;
     355  return *this;
     356}
     357
     358// Returns absolute value of a gmp_complex number
     359//
     360inline gmp_float abs( const gmp_complex & c )
     361{
     362  return hypot(c.real(),c.imag());
     363}
     364
     365gmp_complex sqrt( const gmp_complex & x );
     366
     367inline gmp_complex numberToGmp_Complex( number num )
     368{
     369  if (rField_is_long_C()) {
     370    return *(gmp_complex*)num;
     371  } else {
     372    return gmp_complex( numberToFloat(num) );
     373  }
     374}
     375
     376char *complexToStr( const gmp_complex & c, const  unsigned int oprec );
    177377//<-
    178378
    … …  
    184384// compile-command-2: "make install" ***
    185385// End: ***
     386
     387
     388
     389
     390

  • Singular/mpr_global.h

    r6b9e9c rf543de1  
    44*  Computer Algebra System SINGULAR     *
    55****************************************/
    6 /* $Id: mpr_global.h,v 1.1 1999-04-29 11:38:51 Singular Exp $ */
     6/* $Id: mpr_global.h,v 1.2 1999-06-24 07:46:51 wenk Exp $ */
    77
    88/*
    … …  
    1111*/
    1212
    13 // Define WITH_SUBDIV to compute mixed polyhedral subdivision
    14 //#define WITH_SUBDIV
    15 
    16 // <matrix representation>_<determinant algorithm>
    17 #define DENSE_FACTORY    1
    18 #define SPARSE_BAREISS   2
    19 #define SPARSE_GAUSS     3
    20 
    21 #define USE_MATRIX_REP SPARSE_BAREISS
     13// to get detailed timigs, define MPR_TIMING
     14#define MPR_TIMING
    2215
    2316// Set to double or long double. double is recomended.
    24 // Sets the global floating point type used in mpr_roots.cc and mpr_nric.cc.
     17// Sets the global floating point type used in mpr_numeric.cc.
    2518typedef double mprfloat;
    2619
    2720// --------------------------- debugging stuff ----------------------------
     21#if !defined(NDEBUG)
    2822//#define mprDEBUG_ALL
     23#endif
     24
     25#if !defined(NDEBUG) || defined(mprDEBUG_ALL)
    2926#define mprDEBUG_PROT
     27#endif
     28
     29#if !defined(NDEBUG) && !defined(MPR_TIMING)
     30#define MPR_TIMING
     31#endif
     32
    3033#define mprDEBUG_STICKY
    3134
    … …  
    5154
    5255#ifdef mprDEBUG_STICKY
    53 #define mprSTICKYPROT(msg) Print(msg)
    54 //if (BTEST1(OPT_PROT)) Print(msg)
     56// call 'option(prot);' to get status informations
     57#define mprSTICKYPROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
     58#define mprSTICKYPROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
    5559#else
    5660#define mprSTICKYPROT(msg)
     61#define mprSTICKYPROT2(msg,arg)
    5762#endif
    5863
    … …  
    8994// compile-command-2: "make install" ***
    9095// End: ***
     96
     97
     98
Note: See TracChangeset for help on using the changeset viewer.