Loading...
Searching...
No Matches
amp.h
Go to the documentation of this file.
#ifndef _AMP_R_H
#define _AMP_R_H
#include <gmp.h>
#include <mpfr.h>
#include <stdexcept>
#include <math.h>
#include <string>
#include <stdio.h>
#include <time.h>
#include <memory.h>
#include <vector>
#include <list>
#include <ap.h>
//#define _AMP_NO_TEMPLATE_CONSTRUCTORS
namespace amp
{
class exception {};
class incorrectPrecision : public exception {};
class overflow : public exception {};
class divisionByZero : public exception {};
class sqrtOfNegativeNumber : public exception {};
class invalidConversion : public exception {};
class invalidString : public exception {};
class internalError : public exception {};
class domainError : public exception {};
typedef unsigned long unsigned32;
typedef signed long signed32;
struct mpfr_record
{
unsigned int refCount;
unsigned int Precision;
mpfr_t value;
mpfr_record *next;
};
typedef mpfr_record* mpfr_record_ptr;
//
// storage for mpfr_t instances
//
class mpfr_storage
{
public:
static mpfr_record* newMpfr(unsigned int Precision);
static void deleteMpfr(mpfr_record* ref);
/*static void clearStorage();*/
static gmp_randstate_t* getRandState();
private:
static mpfr_record_ptr& getList(unsigned int Precision);
};
//
// mpfr_t reference
//
class mpfr_reference
{
public:
mpfr_reference();
mpfr_reference(const mpfr_reference& r);
mpfr_reference& operator= (const mpfr_reference &r);
~mpfr_reference();
void initialize(int Precision);
void free();
mpfr_srcptr getReadPtr() const;
mpfr_ptr getWritePtr();
private:
mpfr_record *ref;
};
//
// ampf template
//
template<unsigned int Precision>
class ampf
{
public:
//
// Destructor
//
~ampf()
{
rval->refCount--;
if( rval->refCount==0 )
mpfr_storage::deleteMpfr(rval);
}
//
// Initializing
//
ampf () { InitializeAsZero(); }
ampf(mpfr_record *v) { rval = v; }
ampf (long double v) { InitializeAsDouble(v); }
ampf (double v) { InitializeAsDouble(v); }
ampf (float v) { InitializeAsDouble(v); }
ampf (signed long v) { InitializeAsSLong(v); }
ampf (unsigned long v) { InitializeAsULong(v); }
ampf (signed int v) { InitializeAsSLong(v); }
ampf (unsigned int v) { InitializeAsULong(v); }
ampf (signed short v) { InitializeAsSLong(v); }
ampf (unsigned short v) { InitializeAsULong(v); }
ampf (signed char v) { InitializeAsSLong(v); }
ampf (unsigned char v) { InitializeAsULong(v); }
//
// initializing from string
// string s must have format "X0.hhhhhhhh@eee" or "X-0.hhhhhhhh@eee"
//
ampf (const std::string &s) { InitializeAsString(s.c_str()); }
ampf (const char *s) { InitializeAsString(s); }
//
// copy constructors
//
ampf(const ampf& r)
{
rval = r.rval;
rval->refCount++;
}
#ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
template<unsigned int Precision2>
ampf(const ampf<Precision2>& r)
{
CheckPrecision();
rval = mpfr_storage::newMpfr(Precision);
mpfr_set(getWritePtr(), r.getReadPtr(), GMP_RNDN);
}
#endif
//
// Assignment constructors
//
ampf& operator= (long double v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (double v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (float v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (signed long v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (unsigned long v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (signed int v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (unsigned int v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (signed short v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (unsigned short v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (signed char v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (unsigned char v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN); return *this; }
ampf& operator= (const char *s) { mpfr_strtofr(getWritePtr(), s, NULL, 0, GMP_RNDN); return *this; }
ampf& operator= (const std::string &s) { mpfr_strtofr(getWritePtr(), s.c_str(), NULL, 0, GMP_RNDN); return *this; }
ampf& operator= (const ampf& r)
{
// TODO: may be copy ref
if( this==&r )
return *this;
if( rval==r.rval )
return *this;
rval->refCount--;
if( rval->refCount==0 )
mpfr_storage::deleteMpfr(rval);
rval = r.rval;
rval->refCount++;
//mpfr_set(getWritePtr(), r.getReadPtr(), GMP_RNDN);
return *this;
}
#ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
template<unsigned int Precision2>
ampf& operator= (const ampf<Precision2>& r)
{
if( (void*)this==(void*)(&r) )
return *this;
mpfr_set(getWritePtr(), r.getReadPtr(), GMP_RNDN);
return *this;
}
#endif
//
// in-place operators
// TODO: optimize
//
template<class T> ampf& operator+=(const T& v){ *this = *this + v; return *this; };
template<class T> ampf& operator-=(const T& v){ *this = *this - v; return *this; };
template<class T> ampf& operator*=(const T& v){ *this = *this * v; return *this; };
template<class T> ampf& operator/=(const T& v){ *this = *this / v; return *this; };
//
// MPFR access
//
mpfr_srcptr getReadPtr() const;
mpfr_ptr getWritePtr();
//
// properties and information
//
bool isFiniteNumber() const;
bool isPositiveNumber() const;
bool isZero() const;
bool isNegativeNumber() const;
const ampf getUlpOf();
//
// conversions
//
double toDouble() const;
std::string toHex() const;
std::string toDec() const;
char * toString() const;
//
// static methods
//
static const ampf getUlpOf(const ampf &x);
static const ampf getUlp();
static const ampf getUlp256();
static const ampf getUlp512();
static const ampf getMaxNumber();
static const ampf getMinNumber();
static const ampf getAlgoPascalEpsilon();
static const ampf getAlgoPascalMaxNumber();
static const ampf getAlgoPascalMinNumber();
static const ampf getRandom();
private:
void CheckPrecision();
void InitializeAsZero();
void InitializeAsSLong(signed long v);
void InitializeAsULong(unsigned long v);
void InitializeAsDouble(long double v);
void InitializeAsString(const char *s);
//mpfr_reference ref;
mpfr_record *rval;
};
/*void ampf<Precision>::CheckPrecision()
{
if( Precision<32 )
throw incorrectPrecision();
}***/
template<unsigned int Precision>
void ampf<Precision>::CheckPrecision()
{
if( Precision<32 )
throw incorrectPrecision();
}
template<unsigned int Precision>
void ampf<Precision>::InitializeAsZero()
{
CheckPrecision();
rval = mpfr_storage::newMpfr(Precision);
mpfr_set_ui(getWritePtr(), 0, GMP_RNDN);
}
template<unsigned int Precision>
void ampf<Precision>::InitializeAsSLong(signed long sv)
{
CheckPrecision();
rval = mpfr_storage::newMpfr(Precision);
mpfr_set_si(getWritePtr(), sv, GMP_RNDN);
}
template<unsigned int Precision>
void ampf<Precision>::InitializeAsULong(unsigned long v)
{
CheckPrecision();
rval = mpfr_storage::newMpfr(Precision);
mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
}
template<unsigned int Precision>
void ampf<Precision>::InitializeAsDouble(long double v)
{
CheckPrecision();
rval = mpfr_storage::newMpfr(Precision);
mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
}
template<unsigned int Precision>
void ampf<Precision>::InitializeAsString(const char *s)
{
CheckPrecision();
rval = mpfr_storage::newMpfr(Precision);
mpfr_strtofr(getWritePtr(), s, NULL, 0, GMP_RNDN);
}
template<unsigned int Precision>
mpfr_srcptr ampf<Precision>::getReadPtr() const
{
// TODO: подумать, нужно ли сделать, чтобы и при getRead, и при
// getWrite создавалась новая instance mpfr_t.
// это может быть нужно для корректной обработки ситуаций вида
// mpfr_чего_то_там( a.getWritePtr(), a.getReadPtr())
// вроде бы нужно, а то если там завязано на side-effects...
return rval->value;
}
template<unsigned int Precision>
mpfr_ptr ampf<Precision>::getWritePtr()
{
if( rval->refCount==1 )
return rval->value;
mpfr_record *newrval = mpfr_storage::newMpfr(Precision);
mpfr_set(newrval->value, rval->value, GMP_RNDN);
rval->refCount--;
rval = newrval;
return rval->value;
}
template<unsigned int Precision>
bool ampf<Precision>::isFiniteNumber() const
{
return mpfr_number_p(getReadPtr())!=0;
}
template<unsigned int Precision>
bool ampf<Precision>::isPositiveNumber() const
{
if( !isFiniteNumber() )
return false;
return mpfr_sgn(getReadPtr())>0;
}
template<unsigned int Precision>
bool ampf<Precision>::isZero() const
{
return mpfr_zero_p(getReadPtr())!=0;
}
template<unsigned int Precision>
bool ampf<Precision>::isNegativeNumber() const
{
if( !isFiniteNumber() )
return false;
return mpfr_sgn(getReadPtr())<0;
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getUlpOf()
{
return getUlpOf(*this);
}
template<unsigned int Precision>
double ampf<Precision>::toDouble() const
{
return mpfr_get_d(getReadPtr(), GMP_RNDN);
}
template<unsigned int Precision>
std::string ampf<Precision>::toHex() const
{
//
// some special cases
//
if( !isFiniteNumber() )
{
std::string r;
mp_exp_t _e;
char *ptr;
ptr = mpfr_get_str(NULL, &_e, 16, 0, getReadPtr(), GMP_RNDN);
r = ptr;
mpfr_free_str(ptr);
return r;
}
//
// general case
//
std::string r;
char buf_e[128];
signed long iexpval;
mp_exp_t expval;
char *ptr;
char *ptr2;
ptr = mpfr_get_str(NULL, &expval, 16, 0, getReadPtr(), GMP_RNDN);
ptr2 = ptr;
iexpval = expval;
if( iexpval!=expval )
throw internalError();
sprintf(buf_e, "%ld", long(iexpval));
if( *ptr=='-' )
{
r = "-";
ptr++;
}
r += "0x0.";
r += ptr;
r += "@";
r += buf_e;
mpfr_free_str(ptr2);
return r;
}
template<unsigned int Precision>
std::string ampf<Precision>::toDec() const
{
// TODO: advanced output formatting (zero, integers)
//
// some special cases
//
if( !isFiniteNumber() )
{
std::string r;
mp_exp_t _e;
char *ptr;
ptr = mpfr_get_str(NULL, &_e, 10, 0, getReadPtr(), GMP_RNDN);
r = ptr;
mpfr_free_str(ptr);
return r;
}
//
// general case
//
std::string r;
char buf_e[128];
signed long iexpval;
mp_exp_t expval;
char *ptr;
char *ptr2;
ptr = mpfr_get_str(NULL, &expval, 10, 0, getReadPtr(), GMP_RNDN);
ptr2 = ptr;
iexpval = expval;
if( iexpval!=expval )
throw internalError();
sprintf(buf_e, "%ld", long(iexpval));
if( *ptr=='-' )
{
r = "-";
ptr++;
}
r += "0.";
r += ptr;
r += "E";
r += buf_e;
mpfr_free_str(ptr2);
return r;
}
template<unsigned int Precision>
char * ampf<Precision>::toString() const
{
char *toString_Block=(char *)omAlloc(256);
//
// some special cases
//
if( !isFiniteNumber() )
{
mp_exp_t _e;
char *ptr;
ptr = mpfr_get_str(NULL, &_e, 10, 0, getReadPtr(), GMP_RNDN);
strcpy(toString_Block, ptr);
mpfr_free_str(ptr);
return toString_Block;
}
//
// general case
//
char buf_e[128];
signed long iexpval;
mp_exp_t expval;
char *ptr;
char *ptr2;
ptr = mpfr_get_str(NULL, &expval, 10, 0, getReadPtr(), GMP_RNDN);
ptr2 = ptr;
iexpval = expval;
if( iexpval!=expval )
throw internalError();
sprintf(buf_e, "%ld", long(iexpval));
if( *ptr=='-' )
{
ptr++;
sprintf(toString_Block,"-0.%sE%s",ptr,buf_e);
}
else
sprintf(toString_Block,"0.%sE%s",ptr,buf_e);
mpfr_free_str(ptr2);
return toString_Block;
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getUlpOf(const ampf<Precision> &x)
{
if( !x.isFiniteNumber() )
return x;
if( x.isZero() )
return x;
ampf<Precision> r(1);
mpfr_nextabove(r.getWritePtr());
mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
mpfr_mul_2si(
r.getWritePtr(),
r.getWritePtr(),
mpfr_get_exp(x.getReadPtr()),
GMP_RNDN);
mpfr_div_2si(
r.getWritePtr(),
r.getWritePtr(),
1,
GMP_RNDN);
return r;
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getUlp()
{
ampf<Precision> r(1);
mpfr_nextabove(r.getWritePtr());
mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
return r;
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getUlp256()
{
ampf<Precision> r(1);
mpfr_nextabove(r.getWritePtr());
mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
mpfr_mul_2si(
r.getWritePtr(),
r.getWritePtr(),
8,
GMP_RNDN);
return r;
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getUlp512()
{
ampf<Precision> r(1);
mpfr_nextabove(r.getWritePtr());
mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
mpfr_mul_2si(
r.getWritePtr(),
r.getWritePtr(),
9,
GMP_RNDN);
return r;
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getMaxNumber()
{
ampf<Precision> r(1);
mpfr_nextbelow(r.getWritePtr());
mpfr_set_exp(r.getWritePtr(),mpfr_get_emax());
return r;
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getMinNumber()
{
ampf<Precision> r(1);
mpfr_set_exp(r.getWritePtr(),mpfr_get_emin());
return r;
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getAlgoPascalEpsilon()
{
return getUlp256();
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getAlgoPascalMaxNumber()
{
ampf<Precision> r(1);
mp_exp_t e1 = mpfr_get_emax();
mp_exp_t e2 = -mpfr_get_emin();
mp_exp_t e = e1>e2 ? e1 : e2;
mpfr_set_exp(r.getWritePtr(), e-5);
return r;
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getAlgoPascalMinNumber()
{
ampf<Precision> r(1);
mp_exp_t e1 = mpfr_get_emax();
mp_exp_t e2 = -mpfr_get_emin();
mp_exp_t e = e1>e2 ? e1 : e2;
mpfr_set_exp(r.getWritePtr(), 2-(e-5));
return r;
}
template<unsigned int Precision>
const ampf<Precision> ampf<Precision>::getRandom()
{
ampf<Precision> r;
while(mpfr_urandomb(r.getWritePtr(), *amp::mpfr_storage::getRandState()));
return r;
}
//
// comparison operators
//
template<unsigned int Precision>
const bool operator==(const ampf<Precision>& op1, const ampf<Precision>& op2)
{
return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())==0;
}
template<unsigned int Precision>
const bool operator!=(const ampf<Precision>& op1, const ampf<Precision>& op2)
{
return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())!=0;
}
template<unsigned int Precision>
const bool operator<(const ampf<Precision>& op1, const ampf<Precision>& op2)
{
return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())<0;
}
template<unsigned int Precision>
const bool operator>(const ampf<Precision>& op1, const ampf<Precision>& op2)
{
return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())>0;
}
template<unsigned int Precision>
const bool operator<=(const ampf<Precision>& op1, const ampf<Precision>& op2)
{
return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())<=0;
}
template<unsigned int Precision>
const bool operator>=(const ampf<Precision>& op1, const ampf<Precision>& op2)
{
return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())>=0;
}
//
// arithmetic operators
//
template<unsigned int Precision>
const ampf<Precision> operator+(const ampf<Precision>& op1)
{
return op1;
}
template<unsigned int Precision>
const ampf<Precision> operator-(const ampf<Precision>& op1)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_neg(v->value, op1.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> operator+(const ampf<Precision>& op1, const ampf<Precision>& op2)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_add(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> operator-(const ampf<Precision>& op1, const ampf<Precision>& op2)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_sub(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> operator*(const ampf<Precision>& op1, const ampf<Precision>& op2)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_mul(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> operator/(const ampf<Precision>& op1, const ampf<Precision>& op2)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_div(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
return v;
}
//
// basic functions
//
template<unsigned int Precision>
const ampf<Precision> sqr(const ampf<Precision> &x)
{
// TODO: optimize temporary for return value
ampf<Precision> res;
mpfr_sqr(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
return res;
}
template<unsigned int Precision>
const int sign(const ampf<Precision> &x)
{
int s = mpfr_sgn(x.getReadPtr());
if( s>0 )
return +1;
if( s<0 )
return -1;
return 0;
}
template<unsigned int Precision>
const ampf<Precision> abs(const ampf<Precision> &x)
{
// TODO: optimize temporary for return value
ampf<Precision> res;
mpfr_abs(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
return res;
}
template<unsigned int Precision>
const ampf<Precision> maximum(const ampf<Precision> &x, const ampf<Precision> &y)
{
// TODO: optimize temporary for return value
ampf<Precision> res;
mpfr_max(res.getWritePtr(), x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
return res;
}
template<unsigned int Precision>
const ampf<Precision> minimum(const ampf<Precision> &x, const ampf<Precision> &y)
{
// TODO: optimize temporary for return value
ampf<Precision> res;
mpfr_min(res.getWritePtr(), x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
return res;
}
template<unsigned int Precision>
const ampf<Precision> sqrt(const ampf<Precision> &x)
{
// TODO: optimize temporary for return value
ampf<Precision> res;
mpfr_sqrt(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
return res;
}
template<unsigned int Precision>
const signed long trunc(const ampf<Precision> &x)
{
ampf<Precision> tmp;
signed long r;
mpfr_trunc(tmp.getWritePtr(), x.getReadPtr());
if( mpfr_integer_p(tmp.getReadPtr())==0 )
throw invalidConversion();
mpfr_clear_erangeflag();
r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
if( mpfr_erangeflag_p()!=0 )
throw invalidConversion();
return r;
}
template<unsigned int Precision>
const ampf<Precision> frac(const ampf<Precision> &x)
{
// TODO: optimize temporary for return value
ampf<Precision> r;
mpfr_frac(r.getWritePtr(), x.getReadPtr(), GMP_RNDN);
return r;
}
template<unsigned int Precision>
const signed long floor(const ampf<Precision> &x)
{
ampf<Precision> tmp;
signed long r;
mpfr_floor(tmp.getWritePtr(), x.getReadPtr());
if( mpfr_integer_p(tmp.getReadPtr())==0 )
throw invalidConversion();
mpfr_clear_erangeflag();
r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
if( mpfr_erangeflag_p()!=0 )
throw invalidConversion();
return r;
}
template<unsigned int Precision>
const signed long ceil(const ampf<Precision> &x)
{
ampf<Precision> tmp;
signed long r;
mpfr_ceil(tmp.getWritePtr(), x.getReadPtr());
if( mpfr_integer_p(tmp.getReadPtr())==0 )
throw invalidConversion();
mpfr_clear_erangeflag();
r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
if( mpfr_erangeflag_p()!=0 )
throw invalidConversion();
return r;
}
template<unsigned int Precision>
const signed long round(const ampf<Precision> &x)
{
ampf<Precision> tmp;
signed long r;
mpfr_round(tmp.getWritePtr(), x.getReadPtr());
if( mpfr_integer_p(tmp.getReadPtr())==0 )
throw invalidConversion();
mpfr_clear_erangeflag();
r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
if( mpfr_erangeflag_p()!=0 )
throw invalidConversion();
return r;
}
template<unsigned int Precision>
const ampf<Precision> frexp2(const ampf<Precision> &x, mp_exp_t *exponent)
{
// TODO: optimize temporary for return value
ampf<Precision> r;
if( !x.isFiniteNumber() )
throw invalidConversion();
if( x.isZero() )
{
*exponent = 0;
r = 0;
return r;
}
r = x;
*exponent = mpfr_get_exp(r.getReadPtr());
mpfr_set_exp(r.getWritePtr(),0);
return r;
}
template<unsigned int Precision>
const ampf<Precision> ldexp2(const ampf<Precision> &x, mp_exp_t exponent)
{
// TODO: optimize temporary for return value
ampf<Precision> r;
mpfr_mul_2si(r.getWritePtr(), x.getReadPtr(), exponent, GMP_RNDN);
return r;
}
//
// different types of arguments
//
#define __AMP_BINARY_OPI(type) \
template<unsigned int Precision> const ampf<Precision> operator+(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \
template<unsigned int Precision> const ampf<Precision> operator+(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \
template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const signed type& op2) { return op1+ampf<Precision>(op2); } \
template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const unsigned type& op2) { return op1+ampf<Precision>(op2); } \
template<unsigned int Precision> const ampf<Precision> operator-(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \
template<unsigned int Precision> const ampf<Precision> operator-(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \
template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const signed type& op2) { return op1-ampf<Precision>(op2); } \
template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const unsigned type& op2) { return op1-ampf<Precision>(op2); } \
template<unsigned int Precision> const ampf<Precision> operator*(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \
template<unsigned int Precision> const ampf<Precision> operator*(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \
template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const signed type& op2) { return op1*ampf<Precision>(op2); } \
template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const unsigned type& op2) { return op1*ampf<Precision>(op2); } \
template<unsigned int Precision> const ampf<Precision> operator/(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \
template<unsigned int Precision> const ampf<Precision> operator/(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \
template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const signed type& op2) { return op1/ampf<Precision>(op2); } \
template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const unsigned type& op2) { return op1/ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator==(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \
template<unsigned int Precision> const bool operator==(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \
template<unsigned int Precision> const bool operator==(const ampf<Precision>& op1, const signed type& op2) { return op1==ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator==(const ampf<Precision>& op1, const unsigned type& op2) { return op1==ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator!=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \
template<unsigned int Precision> const bool operator!=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \
template<unsigned int Precision> const bool operator!=(const ampf<Precision>& op1, const signed type& op2) { return op1!=ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator!=(const ampf<Precision>& op1, const unsigned type& op2) { return op1!=ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator<=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \
template<unsigned int Precision> const bool operator<=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \
template<unsigned int Precision> const bool operator<=(const ampf<Precision>& op1, const signed type& op2) { return op1<=ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator<=(const ampf<Precision>& op1, const unsigned type& op2) { return op1<=ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator>=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \
template<unsigned int Precision> const bool operator>=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \
template<unsigned int Precision> const bool operator>=(const ampf<Precision>& op1, const signed type& op2) { return op1>=ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator>=(const ampf<Precision>& op1, const unsigned type& op2) { return op1>=ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator<(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \
template<unsigned int Precision> const bool operator<(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \
template<unsigned int Precision> const bool operator<(const ampf<Precision>& op1, const signed type& op2) { return op1<ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator<(const ampf<Precision>& op1, const unsigned type& op2) { return op1<ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator>(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \
template<unsigned int Precision> const bool operator>(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \
template<unsigned int Precision> const bool operator>(const ampf<Precision>& op1, const signed type& op2) { return op1>ampf<Precision>(op2); } \
template<unsigned int Precision> const bool operator>(const ampf<Precision>& op1, const unsigned type& op2) { return op1>ampf<Precision>(op2); }
__AMP_BINARY_OPI(char)
__AMP_BINARY_OPI(short)
__AMP_BINARY_OPI(long)
__AMP_BINARY_OPI(int)
#undef __AMP_BINARY_OPI
#define __AMP_BINARY_OPF(type) \
template<unsigned int Precision> const ampf<Precision> operator+(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \
template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const type& op2) { return op1+ampf<Precision>(op2); } \
template<unsigned int Precision> const ampf<Precision> operator-(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \
template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const type& op2) { return op1-ampf<Precision>(op2); } \
template<unsigned int Precision> const ampf<Precision> operator*(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \
template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const type& op2) { return op1*ampf<Precision>(op2); } \
template<unsigned int Precision> const ampf<Precision> operator/(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \
template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const type& op2) { return op1/ampf<Precision>(op2); } \
template<unsigned int Precision> bool operator==(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \
template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const type& op2) { return op1==ampf<Precision>(op2); } \
template<unsigned int Precision> bool operator!=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \
template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const type& op2) { return op1!=ampf<Precision>(op2); } \
template<unsigned int Precision> bool operator<=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \
template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const type& op2) { return op1<=ampf<Precision>(op2); } \
template<unsigned int Precision> bool operator>=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \
template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const type& op2) { return op1>=ampf<Precision>(op2); } \
template<unsigned int Precision> bool operator<(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \
template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const type& op2) { return op1<ampf<Precision>(op2); } \
template<unsigned int Precision> bool operator>(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \
template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const type& op2) { return op1>ampf<Precision>(op2); }
__AMP_BINARY_OPF(float)
__AMP_BINARY_OPF(double)
__AMP_BINARY_OPF(long double)
#undef __AMP_BINARY_OPF
//
// transcendent functions
//
template<unsigned int Precision>
const ampf<Precision> pi()
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_const_pi(v->value, GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> halfpi()
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_const_pi(v->value, GMP_RNDN);
mpfr_mul_2si(v->value, v->value, -1, GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> twopi()
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_const_pi(v->value, GMP_RNDN);
mpfr_mul_2si(v->value, v->value, +1, GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> sin(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_sin(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> cos(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_cos(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> tan(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_tan(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> asin(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_asin(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> acos(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_acos(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> atan(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_atan(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> atan2(const ampf<Precision> &y, const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_atan2(v->value, y.getReadPtr(), x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> log(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_log(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> log2(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_log2(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> log10(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_log10(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> exp(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_exp(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> sinh(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_sinh(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> cosh(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_cosh(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> tanh(const ampf<Precision> &x)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_tanh(v->value, x.getReadPtr(), GMP_RNDN);
return v;
}
template<unsigned int Precision>
const ampf<Precision> pow(const ampf<Precision> &x, const ampf<Precision> &y)
{
mpfr_record *v = mpfr_storage::newMpfr(Precision);
mpfr_pow(v->value, x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
return v;
}
//
// complex ampf
//
template<unsigned int Precision>
class campf
{
public:
campf():x(0),y(0){};
campf(long double v) { x=v; y=0; }
campf(double v) { x=v; y=0; }
campf(float v) { x=v; y=0; }
campf(signed long v) { x=v; y=0; }
campf(unsigned long v) { x=v; y=0; }
campf(signed int v) { x=v; y=0; }
campf(unsigned int v) { x=v; y=0; }
campf(signed short v) { x=v; y=0; }
campf(unsigned short v) { x=v; y=0; }
campf(signed char v) { x=v; y=0; }
campf(unsigned char v) { x=v; y=0; }
campf(const ampf<Precision> &_x):x(_x),y(0){};
campf(const ampf<Precision> &_x, const ampf<Precision> &_y):x(_x),y(_y){};
campf(const campf &z):x(z.x),y(z.y){};
#ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
template<unsigned int Prec2>
campf(const campf<Prec2> &z):x(z.x),y(z.y){};
#endif
campf& operator= (long double v) { x=v; y=0; return *this; }
campf& operator= (double v) { x=v; y=0; return *this; }
campf& operator= (float v) { x=v; y=0; return *this; }
campf& operator= (signed long v) { x=v; y=0; return *this; }
campf& operator= (unsigned long v) { x=v; y=0; return *this; }
campf& operator= (signed int v) { x=v; y=0; return *this; }
campf& operator= (unsigned int v) { x=v; y=0; return *this; }
campf& operator= (signed short v) { x=v; y=0; return *this; }
campf& operator= (unsigned short v) { x=v; y=0; return *this; }
campf& operator= (signed char v) { x=v; y=0; return *this; }
campf& operator= (unsigned char v) { x=v; y=0; return *this; }
campf& operator= (const char *s) { x=s; y=0; return *this; }
campf& operator= (const std::string &s) { x=s; y=0; return *this; }
campf& operator= (const campf& r)
{
x = r.x;
y = r.y;
return *this;
}
#ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
template<unsigned int Precision2>
campf& operator= (const campf<Precision2>& r)
{
x = r.x;
y = r.y;
return *this;
}
#endif
ampf<Precision> x, y;
};
//
// complex operations
//
template<unsigned int Precision>
const bool operator==(const campf<Precision>& lhs, const campf<Precision>& rhs)
{ return lhs.x==rhs.x && lhs.y==rhs.y; }
template<unsigned int Precision>
const bool operator!=(const campf<Precision>& lhs, const campf<Precision>& rhs)
{ return lhs.x!=rhs.x || lhs.y!=rhs.y; }
template<unsigned int Precision>
const campf<Precision> operator+(const campf<Precision>& lhs)
{ return lhs; }
template<unsigned int Precision>
campf<Precision>& operator+=(campf<Precision>& lhs, const campf<Precision>& rhs)
{ lhs.x += rhs.x; lhs.y += rhs.y; return lhs; }
template<unsigned int Precision>
const campf<Precision> operator+(const campf<Precision>& lhs, const campf<Precision>& rhs)
{ campf<Precision> r = lhs; r += rhs; return r; }
template<unsigned int Precision>
const campf<Precision> operator-(const campf<Precision>& lhs)
{ return campf<Precision>(-lhs.x, -lhs.y); }
template<unsigned int Precision>
campf<Precision>& operator-=(campf<Precision>& lhs, const campf<Precision>& rhs)
{ lhs.x -= rhs.x; lhs.y -= rhs.y; return lhs; }
template<unsigned int Precision>
const campf<Precision> operator-(const campf<Precision>& lhs, const campf<Precision>& rhs)
{ campf<Precision> r = lhs; r -= rhs; return r; }
template<unsigned int Precision>
campf<Precision>& operator*=(campf<Precision>& lhs, const campf<Precision>& rhs)
{
ampf<Precision> xx(lhs.x*rhs.x), yy(lhs.y*rhs.y), mm((lhs.x+lhs.y)*(rhs.x+rhs.y));
lhs.x = xx-yy;
lhs.y = mm-xx-yy;
return lhs;
}
template<unsigned int Precision>
const campf<Precision> operator*(const campf<Precision>& lhs, const campf<Precision>& rhs)
{ campf<Precision> r = lhs; r *= rhs; return r; }
template<unsigned int Precision>
const campf<Precision> operator/(const campf<Precision>& lhs, const campf<Precision>& rhs)
{
campf<Precision> result;
ampf<Precision> e;
ampf<Precision> f;
if( abs(rhs.y)<abs(rhs.x) )
{
e = rhs.y/rhs.x;
f = rhs.x+rhs.y*e;
result.x = (lhs.x+lhs.y*e)/f;
result.y = (lhs.y-lhs.x*e)/f;
}
else
{
e = rhs.x/rhs.y;
f = rhs.y+rhs.x*e;
result.x = (lhs.y+lhs.x*e)/f;
result.y = (-lhs.x+lhs.y*e)/f;
}
return result;
}
template<unsigned int Precision>
campf<Precision>& operator/=(campf<Precision>& lhs, const campf<Precision>& rhs)
{
lhs = lhs/rhs;
return lhs;
}
template<unsigned int Precision>
const ampf<Precision> abscomplex(const campf<Precision> &z)
{
ampf<Precision> w, xabs, yabs, v;
xabs = abs(z.x);
yabs = abs(z.y);
w = xabs>yabs ? xabs : yabs;
v = xabs<yabs ? xabs : yabs;
if( v==0 )
return w;
else
{
ampf<Precision> t = v/w;
return w*sqrt(1+sqr(t));
}
}
template<unsigned int Precision>
const campf<Precision> conj(const campf<Precision> &z)
{
return campf<Precision>(z.x, -z.y);
}
template<unsigned int Precision>
const campf<Precision> csqr(const campf<Precision> &z)
{
ampf<Precision> t = z.x*z.y;
return campf<Precision>(sqr(z.x)-sqr(z.y), t+t);
}
//
// different types of arguments
//
#define __AMP_BINARY_OPI(type) \
template<unsigned int Precision> const campf<Precision> operator+ (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \
template<unsigned int Precision> const campf<Precision> operator+ (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \
template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \
template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \
template<unsigned int Precision> const campf<Precision> operator- (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \
template<unsigned int Precision> const campf<Precision> operator- (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \
template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \
template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \
template<unsigned int Precision> const campf<Precision> operator* (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \
template<unsigned int Precision> const campf<Precision> operator* (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \
template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \
template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \
template<unsigned int Precision> const campf<Precision> operator/ (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(ampf<Precision>(op1),ampf<Precision>(0))/op2; } \
template<unsigned int Precision> const campf<Precision> operator/ (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(ampf<Precision>(op1),ampf<Precision>(0))/op2; } \
template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \
template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \
template<unsigned int Precision> bool operator==(const signed type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \
template<unsigned int Precision> bool operator==(const unsigned type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \
template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const signed type& op2) { return op1.x==op2 && op1.y==0; } \
template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const unsigned type& op2) { return op1.x==op2 && op1.y==0; } \
template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const signed type& op2) { return op1.x!=op2 || op1.y!=0; } \
template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const unsigned type& op2) { return op1.x!=op2 || op1.y!=0; } \
template<unsigned int Precision> bool operator!=(const signed type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } \
template<unsigned int Precision> bool operator!=(const unsigned type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; }
__AMP_BINARY_OPI(char)
__AMP_BINARY_OPI(short)
__AMP_BINARY_OPI(long)
__AMP_BINARY_OPI(int)
#undef __AMP_BINARY_OPI
#define __AMP_BINARY_OPF(type) \
template<unsigned int Precision> const campf<Precision> operator+ (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \
template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \
template<unsigned int Precision> const campf<Precision> operator- (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \
template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \
template<unsigned int Precision> const campf<Precision> operator* (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \
template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \
template<unsigned int Precision> const campf<Precision> operator/ (const type& op1, const campf<Precision>& op2) { return campf<Precision>(ampf<Precision>(op1),ampf<Precision>(0))/op2; } \
template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \
template<unsigned int Precision> bool operator==(const type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \
template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const type& op2) { return op1.x==op2 && op1.y==0; } \
template<unsigned int Precision> bool operator!=(const type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } \
template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const type& op2) { return op1.x!=op2 || op1.y!=0; }
__AMP_BINARY_OPF(float)
__AMP_BINARY_OPF(double)
__AMP_BINARY_OPF(long double)
__AMP_BINARY_OPF(ampf<Precision>)
#undef __AMP_BINARY_OPF
//
// Real linear algebra
//
template<unsigned int Precision>
ampf<Precision> vDotProduct(ap::const_raw_vector< ampf<Precision> > v1, ap::const_raw_vector< ampf<Precision> > v2)
{
ap::ap_error::make_assertion(v1.GetLength()==v2.GetLength());
int i, cnt = v1.GetLength();
const ampf<Precision> *p1 = v1.GetData();
const ampf<Precision> *p2 = v2.GetData();
mpfr_record *r = NULL;
mpfr_record *t = NULL;
try
{
r = mpfr_storage::newMpfr(Precision);
t = mpfr_storage::newMpfr(Precision);
mpfr_set_ui(r->value, 0, GMP_RNDN);
for(i=0; i<cnt; i++)
{
mpfr_mul(t->value, p1->getReadPtr(), p2->getReadPtr(), GMP_RNDN);
mpfr_add(r->value, r->value, t->value, GMP_RNDN);
p1 += v1.GetStep();
p2 += v2.GetStep();
}
mpfr_storage::deleteMpfr(t);
return r;
}
catch(...)
{
if( r!=NULL )
mpfr_storage::deleteMpfr(r);
if( t!=NULL )
mpfr_storage::deleteMpfr(t);
throw;
}
}
template<unsigned int Precision>
void vMove(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
{
ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
int i, cnt = vDst.GetLength();
ampf<Precision> *pDst = vDst.GetData();
const ampf<Precision> *pSrc = vSrc.GetData();
if( pDst==pSrc )
return;
for(i=0; i<cnt; i++)
{
*pDst = *pSrc;
pDst += vDst.GetStep();
pSrc += vSrc.GetStep();
}
}
template<unsigned int Precision>
void vMoveNeg(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
{
ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
int i, cnt = vDst.GetLength();
ampf<Precision> *pDst = vDst.GetData();
const ampf<Precision> *pSrc = vSrc.GetData();
for(i=0; i<cnt; i++)
{
*pDst = *pSrc;
mpfr_ptr v = pDst->getWritePtr();
mpfr_neg(v, v, GMP_RNDN);
pDst += vDst.GetStep();
pSrc += vSrc.GetStep();
}
}
template<unsigned int Precision, class T2>
void vMove(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
{
ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
int i, cnt = vDst.GetLength();
ampf<Precision> *pDst = vDst.GetData();
const ampf<Precision> *pSrc = vSrc.GetData();
ampf<Precision> a(alpha);
for(i=0; i<cnt; i++)
{
*pDst = *pSrc;
mpfr_ptr v = pDst->getWritePtr();
mpfr_mul(v, v, a.getReadPtr(), GMP_RNDN);
pDst += vDst.GetStep();
pSrc += vSrc.GetStep();
}
}
template<unsigned int Precision>
void vAdd(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
{
ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
int i, cnt = vDst.GetLength();
ampf<Precision> *pDst = vDst.GetData();
const ampf<Precision> *pSrc = vSrc.GetData();
for(i=0; i<cnt; i++)
{
mpfr_ptr v = pDst->getWritePtr();
mpfr_srcptr vs = pSrc->getReadPtr();
mpfr_add(v, v, vs, GMP_RNDN);
pDst += vDst.GetStep();
pSrc += vSrc.GetStep();
}
}
template<unsigned int Precision, class T2>
void vAdd(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
{
ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
int i, cnt = vDst.GetLength();
ampf<Precision> *pDst = vDst.GetData();
const ampf<Precision> *pSrc = vSrc.GetData();
ampf<Precision> a(alpha), tmp;
for(i=0; i<cnt; i++)
{
mpfr_ptr v = pDst->getWritePtr();
mpfr_srcptr vs = pSrc->getReadPtr();
mpfr_mul(tmp.getWritePtr(), a.getReadPtr(), vs, GMP_RNDN);
mpfr_add(v, v, tmp.getWritePtr(), GMP_RNDN);
pDst += vDst.GetStep();
pSrc += vSrc.GetStep();
}
}
template<unsigned int Precision>
void vSub(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
{
ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
int i, cnt = vDst.GetLength();
ampf<Precision> *pDst = vDst.GetData();
const ampf<Precision> *pSrc = vSrc.GetData();
for(i=0; i<cnt; i++)
{
mpfr_ptr v = pDst->getWritePtr();
mpfr_srcptr vs = pSrc->getReadPtr();
mpfr_sub(v, v, vs, GMP_RNDN);
pDst += vDst.GetStep();
pSrc += vSrc.GetStep();
}
}
template<unsigned int Precision, class T2>
void vSub(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
{
vAdd(vDst, vSrc, -alpha);
}
template<unsigned int Precision, class T2>
void vMul(ap::raw_vector< ampf<Precision> > vDst, T2 alpha)
{
int i, cnt = vDst.GetLength();
ampf<Precision> *pDst = vDst.GetData();
ampf<Precision> a(alpha);
for(i=0; i<cnt; i++)
{
mpfr_ptr v = pDst->getWritePtr();
mpfr_mul(v, a.getReadPtr(), v, GMP_RNDN);
pDst += vDst.GetStep();
}
}
}
#endif
__AMP_BINARY_OPF
#define __AMP_BINARY_OPF(type)
__AMP_BINARY_OPI
#define __AMP_BINARY_OPI(type)
Definition amp.h:848
ap.h
i
int i
Definition cfEzgcd.cc:132
x
Variable x
Definition cfModGcd.cc:4090
f
FILE * f
Definition checklibs.c:9
amp::ampf
Definition amp.h:82
amp::ampf::operator/=
ampf & operator/=(const T &v)
Definition amp.h:186
amp::ampf::ampf
ampf(unsigned int v)
Definition amp.h:106
amp::ampf::operator-=
ampf & operator-=(const T &v)
Definition amp.h:184
amp::ampf::getAlgoPascalEpsilon
static const ampf getAlgoPascalEpsilon()
Definition amp.h:565
amp::ampf::ampf
ampf(signed char v)
Definition amp.h:109
amp::ampf::InitializeAsSLong
void InitializeAsSLong(signed long v)
Definition amp.h:259
amp::ampf::InitializeAsULong
void InitializeAsULong(unsigned long v)
Definition amp.h:267
amp::ampf::getUlp256
static const ampf getUlp256()
Definition amp.h:520
amp::ampf::isNegativeNumber
bool isNegativeNumber() const
Definition amp.h:334
amp::ampf::InitializeAsDouble
void InitializeAsDouble(long double v)
Definition amp.h:275
amp::ampf::ampf
ampf(const char *s)
Definition amp.h:117
amp::ampf::InitializeAsZero
void InitializeAsZero()
Definition amp.h:251
amp::ampf::ampf
ampf()
Definition amp.h:97
amp::ampf::ampf
ampf(signed long v)
Definition amp.h:103
amp::ampf::ampf
ampf(mpfr_record *v)
Definition amp.h:98
amp::ampf::ampf
ampf(signed short v)
Definition amp.h:107
amp::ampf::ampf
ampf(unsigned long v)
Definition amp.h:104
amp::ampf::getWritePtr
mpfr_ptr getWritePtr()
Definition amp.h:302
amp::ampf::toString
char * toString() const
Definition amp.h:445
amp::ampf::CheckPrecision
void CheckPrecision()
Definition amp.h:244
amp::ampf::operator+=
ampf & operator+=(const T &v)
Definition amp.h:183
amp::ampf::getAlgoPascalMaxNumber
static const ampf getAlgoPascalMaxNumber()
Definition amp.h:571
amp::ampf::getRandom
static const ampf getRandom()
Definition amp.h:593
amp::ampf::isPositiveNumber
bool isPositiveNumber() const
Definition amp.h:320
amp::ampf::ampf
ampf(double v)
Definition amp.h:101
amp::ampf::InitializeAsString
void InitializeAsString(const char *s)
Definition amp.h:283
amp::ampf::getReadPtr
mpfr_srcptr getReadPtr() const
Definition amp.h:291
amp::ampf::getUlp512
static const ampf getUlp512()
Definition amp.h:534
amp::ampf::toDouble
double toDouble() const
Definition amp.h:348
amp::ampf::getUlpOf
const ampf getUlpOf()
Definition amp.h:342
amp::ampf::isZero
bool isZero() const
Definition amp.h:328
amp::ampf::~ampf
~ampf()
Definition amp.h:87
amp::ampf::ampf
ampf(float v)
Definition amp.h:102
amp::ampf::ampf
ampf(unsigned char v)
Definition amp.h:110
amp::ampf::operator=
ampf & operator=(long double v)
Definition amp.h:140
amp::ampf::rval
mpfr_record * rval
Definition amp.h:234
amp::ampf::ampf
ampf(long double v)
Definition amp.h:100
amp::ampf::getUlp
static const ampf getUlp()
Definition amp.h:511
amp::ampf::ampf
ampf(const ampf &r)
Definition amp.h:122
amp::ampf::toHex
std::string toHex() const
Definition amp.h:354
amp::ampf::operator*=
ampf & operator*=(const T &v)
Definition amp.h:185
amp::ampf::getMinNumber
static const ampf getMinNumber()
Definition amp.h:557
amp::ampf::isFiniteNumber
bool isFiniteNumber() const
Definition amp.h:314
amp::ampf::ampf
ampf(const ampf< Precision2 > &r)
Definition amp.h:129
amp::ampf::toDec
std::string toDec() const
Definition amp.h:399
amp::ampf::ampf
ampf(const std::string &s)
Definition amp.h:116
amp::ampf::ampf
ampf(unsigned short v)
Definition amp.h:108
amp::ampf::getMaxNumber
static const ampf getMaxNumber()
Definition amp.h:548
amp::ampf::ampf
ampf(signed int v)
Definition amp.h:105
amp::ampf::getAlgoPascalMinNumber
static const ampf getAlgoPascalMinNumber()
Definition amp.h:582
amp::campf
Definition amp.h:1074
amp::campf::operator=
campf & operator=(long double v)
Definition amp.h:1096
amp::campf::campf
campf(signed long v)
Definition amp.h:1080
amp::campf::campf
campf(float v)
Definition amp.h:1079
amp::campf::campf
campf(unsigned int v)
Definition amp.h:1083
amp::campf::campf
campf(unsigned char v)
Definition amp.h:1087
amp::campf::campf
campf(signed int v)
Definition amp.h:1082
amp::campf::campf
campf(unsigned short v)
Definition amp.h:1085
amp::campf::y
ampf< Precision > y
Definition amp.h:1125
amp::campf::campf
campf(const ampf< Precision > &_x)
Definition amp.h:1088
amp::campf::campf
campf(double v)
Definition amp.h:1078
amp::campf::campf
campf(const ampf< Precision > &_x, const ampf< Precision > &_y)
Definition amp.h:1089
amp::campf::campf
campf(signed char v)
Definition amp.h:1086
amp::campf::campf
campf(const campf< Prec2 > &z)
Definition amp.h:1093
amp::campf::x
ampf< Precision > x
Definition amp.h:1125
amp::campf::campf
campf(signed short v)
Definition amp.h:1084
amp::campf::campf
campf()
Definition amp.h:1076
amp::campf::campf
campf(const campf &z)
Definition amp.h:1090
amp::campf::campf
campf(long double v)
Definition amp.h:1077
amp::campf::campf
campf(unsigned long v)
Definition amp.h:1081
amp::divisionByZero
Definition amp.h:23
amp::domainError
Definition amp.h:28
amp::exception
Definition amp.h:20
amp::incorrectPrecision
Definition amp.h:21
amp::internalError
Definition amp.h:27
amp::invalidConversion
Definition amp.h:25
amp::invalidString
Definition amp.h:26
amp::mpfr_reference::initialize
void initialize(int Precision)
Definition amp.cpp:115
amp::mpfr_reference::mpfr_reference
mpfr_reference()
Definition amp.cpp:82
amp::mpfr_reference::ref
mpfr_record * ref
Definition amp.h:74
amp::mpfr_reference::operator=
mpfr_reference & operator=(const mpfr_reference &r)
Definition amp.cpp:94
amp::mpfr_reference::getWritePtr
mpfr_ptr getWritePtr()
Definition amp.cpp:142
amp::mpfr_reference::free
void free()
Definition amp.cpp:123
amp::mpfr_reference::~mpfr_reference
~mpfr_reference()
Definition amp.cpp:109
amp::mpfr_reference::getReadPtr
mpfr_srcptr getReadPtr() const
Definition amp.cpp:134
amp::mpfr_storage
Definition amp.h:47
amp::mpfr_storage::getRandState
static gmp_randstate_t * getRandState()
Definition amp.cpp:37
amp::mpfr_storage::deleteMpfr
static void deleteMpfr(mpfr_record *ref)
Definition amp.cpp:30
amp::mpfr_storage::newMpfr
static mpfr_record * newMpfr(unsigned int Precision)
Definition amp.cpp:11
amp::mpfr_storage::getList
static mpfr_record_ptr & getList(unsigned int Precision)
Definition amp.cpp:63
amp::overflow
Definition amp.h:22
amp::sqrtOfNegativeNumber
Definition amp.h:24
ap::ap_error::make_assertion
static void make_assertion(bool bClause)
Definition ap.h:49
ap::const_raw_vector
Definition ap.h:137
ap::raw_vector
Definition ap.h:168
alpha
Variable alpha
Definition facAbsBiFact.cc:52
result
return result
Definition facAbsBiFact.cc:76
s
const CanonicalForm int s
Definition facAbsFact.cc:51
y
const CanonicalForm int const CFList const Variable & y
Definition facAbsFact.cc:53
res
CanonicalForm res
Definition facAbsFact.cc:60
w
const CanonicalForm & w
Definition facAbsFact.cc:51
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
exponent
#define exponent
Definition interpolation.cc:31
T
STATIC_VAR jList * T
Definition janet.cc:30
pi
#define pi
Definition libparse.cc:1145
amp
Definition amp.h:19
amp::operator-=
campf< Precision > & operator-=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition amp.h:1156
amp::operator*=
campf< Precision > & operator*=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition amp.h:1164
amp::abs
const ampf< Precision > abs(const ampf< Precision > &x)
Definition amp.h:713
amp::cos
const ampf< Precision > cos(const ampf< Precision > &x)
Definition amp.h:958
amp::halfpi
const ampf< Precision > halfpi()
Definition amp.h:932
amp::unsigned32
unsigned long unsigned32
Definition amp.h:30
amp::operator/=
campf< Precision > & operator/=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition amp.h:1200
amp::floor
const signed long floor(const ampf< Precision > &x)
Definition amp.h:773
amp::operator>
const bool operator>(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition amp.h:622
amp::operator-
const ampf< Precision > operator-(const ampf< Precision > &op1)
Definition amp.h:649
amp::vAdd
void vAdd(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition amp.h:1379
amp::ceil
const signed long ceil(const ampf< Precision > &x)
Definition amp.h:788
amp::abscomplex
const ampf< Precision > abscomplex(const campf< Precision > &z)
Definition amp.h:1207
amp::operator/
const ampf< Precision > operator/(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition amp.h:682
amp::mpfr_record::value
mpfr_t value
Definition amp.h:37
amp::cosh
const ampf< Precision > cosh(const ampf< Precision > &x)
Definition amp.h:1046
amp::ldexp2
const ampf< Precision > ldexp2(const ampf< Precision > &x, mp_exp_t exponent)
Definition amp.h:837
amp::mpfr_record::Precision
unsigned int Precision
Definition amp.h:36
amp::maximum
const ampf< Precision > maximum(const ampf< Precision > &x, const ampf< Precision > &y)
Definition amp.h:722
amp::exp
const ampf< Precision > exp(const ampf< Precision > &x)
Definition amp.h:1030
amp::sqrt
const ampf< Precision > sqrt(const ampf< Precision > &x)
Definition amp.h:740
amp::conj
const campf< Precision > conj(const campf< Precision > &z)
Definition amp.h:1225
amp::twopi
const ampf< Precision > twopi()
Definition amp.h:941
amp::operator>=
const bool operator>=(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition amp.h:634
amp::vMoveNeg
void vMoveNeg(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition amp.h:1344
amp::atan2
const ampf< Precision > atan2(const ampf< Precision > &y, const ampf< Precision > &x)
Definition amp.h:998
amp::csqr
const campf< Precision > csqr(const campf< Precision > &z)
Definition amp.h:1231
amp::log2
const ampf< Precision > log2(const ampf< Precision > &x)
Definition amp.h:1014
amp::operator<=
const bool operator<=(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition amp.h:628
amp::operator+
const ampf< Precision > operator+(const ampf< Precision > &op1)
Definition amp.h:643
amp::operator!=
const bool operator!=(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition amp.h:610
amp::vMove
void vMove(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition amp.h:1327
amp::log10
const ampf< Precision > log10(const ampf< Precision > &x)
Definition amp.h:1022
amp::asin
const ampf< Precision > asin(const ampf< Precision > &x)
Definition amp.h:974
amp::operator<
const bool operator<(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition amp.h:616
amp::pow
const ampf< Precision > pow(const ampf< Precision > &x, const ampf< Precision > &y)
Definition amp.h:1062
amp::operator+=
campf< Precision > & operator+=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition amp.h:1144
amp::tanh
const ampf< Precision > tanh(const ampf< Precision > &x)
Definition amp.h:1054
amp::frexp2
const ampf< Precision > frexp2(const ampf< Precision > &x, mp_exp_t *exponent)
Definition amp.h:818
amp::vSub
void vSub(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition amp.h:1415
amp::atan
const ampf< Precision > atan(const ampf< Precision > &x)
Definition amp.h:990
amp::vMul
void vMul(ap::raw_vector< ampf< Precision > > vDst, T2 alpha)
Definition amp.h:1438
amp::sign
const int sign(const ampf< Precision > &x)
Definition amp.h:702
amp::sin
const ampf< Precision > sin(const ampf< Precision > &x)
Definition amp.h:950
amp::mpfr_record_ptr
mpfr_record * mpfr_record_ptr
Definition amp.h:41
amp::sqr
const ampf< Precision > sqr(const ampf< Precision > &x)
Definition amp.h:693
amp::minimum
const ampf< Precision > minimum(const ampf< Precision > &x, const ampf< Precision > &y)
Definition amp.h:731
amp::sinh
const ampf< Precision > sinh(const ampf< Precision > &x)
Definition amp.h:1038
amp::operator*
const ampf< Precision > operator*(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition amp.h:674
amp::trunc
const signed long trunc(const ampf< Precision > &x)
Definition amp.h:749
amp::frac
const ampf< Precision > frac(const ampf< Precision > &x)
Definition amp.h:764
amp::operator==
const bool operator==(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition amp.h:604
amp::acos
const ampf< Precision > acos(const ampf< Precision > &x)
Definition amp.h:982
amp::mpfr_record::next
mpfr_record * next
Definition amp.h:38
amp::signed32
signed long signed32
Definition amp.h:31
amp::tan
const ampf< Precision > tan(const ampf< Precision > &x)
Definition amp.h:966
amp::round
const signed long round(const ampf< Precision > &x)
Definition amp.h:803
amp::log
const ampf< Precision > log(const ampf< Precision > &x)
Definition amp.h:1006
amp::mpfr_record::refCount
unsigned int refCount
Definition amp.h:35
amp::mpfr_record
Definition amp.h:34
omAlloc
#define omAlloc(size)
Definition omAllocDecl.h:210
NULL
#define NULL
Definition omList.c:12