My Project  debian-1:4.1.1-p2+ds-4build2
amp.h
Go to the documentation of this file.
1 #ifndef _AMP_R_H
2 #define _AMP_R_H
3 
4 #include "omalloc/omalloc.h"
5 #include <gmp.h>
6 #include <mpfr.h>
7 #include <stdexcept>
8 #include <math.h>
9 #include <string>
10 #include <stdio.h>
11 #include <time.h>
12 #include <memory.h>
13 #include <vector>
14 #include <list>
15 #include <ap.h>
16 
17 //#define _AMP_NO_TEMPLATE_CONSTRUCTORS
18 
19 namespace amp
20 {
21  class exception {};
22  class incorrectPrecision : public exception {};
23  class overflow : public exception {};
24  class divisionByZero : public exception {};
25  class sqrtOfNegativeNumber : public exception {};
26  class invalidConversion : public exception {};
27  class invalidString : public exception {};
28  class internalError : public exception {};
29  class domainError : public exception {};
30 
31  typedef unsigned long unsigned32;
32  typedef signed long signed32;
33 
34  struct mpfr_record
35  {
36  unsigned int refCount;
37  unsigned int Precision;
38  mpfr_t value;
40  };
41 
43 
44  //
45  // storage for mpfr_t instances
46  //
48  {
49  public:
50  static mpfr_record* newMpfr(unsigned int Precision);
51  static void deleteMpfr(mpfr_record* ref);
52  /*static void clearStorage();*/
53  static gmp_randstate_t* getRandState();
54  private:
55  static mpfr_record_ptr& getList(unsigned int Precision);
56  };
57 
58  //
59  // mpfr_t reference
60  //
62  {
63  public:
68 
69  void initialize(int Precision);
70  void free();
71 
72  mpfr_srcptr getReadPtr() const;
73  mpfr_ptr getWritePtr();
74  private:
76  };
77 
78  //
79  // ampf template
80  //
81  template<unsigned int Precision>
82  class ampf
83  {
84  public:
85  //
86  // Destructor
87  //
89  {
90  rval->refCount--;
91  if( rval->refCount==0 )
93  }
94 
95  //
96  // Initializing
97  //
99  ampf(mpfr_record *v) { rval = v; }
100 
101  ampf (long double v) { InitializeAsDouble(v); }
102  ampf (double v) { InitializeAsDouble(v); }
103  ampf (float v) { InitializeAsDouble(v); }
104  ampf (signed long v) { InitializeAsSLong(v); }
105  ampf (unsigned long v) { InitializeAsULong(v); }
106  ampf (signed int v) { InitializeAsSLong(v); }
107  ampf (unsigned int v) { InitializeAsULong(v); }
108  ampf (signed short v) { InitializeAsSLong(v); }
109  ampf (unsigned short v) { InitializeAsULong(v); }
110  ampf (signed char v) { InitializeAsSLong(v); }
111  ampf (unsigned char v) { InitializeAsULong(v); }
112 
113  //
114  // initializing from string
115  // string s must have format "X0.hhhhhhhh@eee" or "X-0.hhhhhhhh@eee"
116  //
117  ampf (const std::string &s) { InitializeAsString(s.c_str()); }
118  ampf (const char *s) { InitializeAsString(s); }
119 
120  //
121  // copy constructors
122  //
123  ampf(const ampf& r)
124  {
125  rval = r.rval;
126  rval->refCount++;
127  }
128 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
129  template<unsigned int Precision2>
131  {
132  CheckPrecision();
133  rval = mpfr_storage::newMpfr(Precision);
134  mpfr_set(getWritePtr(), r.getReadPtr(), GMP_RNDN);
135  }
136 #endif
137 
138  //
139  // Assignment constructors
140  //
141  ampf& operator= (long double v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN); return *this; }
142  ampf& operator= (double v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN); return *this; }
143  ampf& operator= (float v) { mpfr_set_ld(getWritePtr(), v, GMP_RNDN); return *this; }
144  ampf& operator= (signed long v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN); return *this; }
145  ampf& operator= (unsigned long v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN); return *this; }
146  ampf& operator= (signed int v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN); return *this; }
147  ampf& operator= (unsigned int v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN); return *this; }
148  ampf& operator= (signed short v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN); return *this; }
149  ampf& operator= (unsigned short v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN); return *this; }
150  ampf& operator= (signed char v) { mpfr_set_si(getWritePtr(), v, GMP_RNDN); return *this; }
151  ampf& operator= (unsigned char v) { mpfr_set_ui(getWritePtr(), v, GMP_RNDN); return *this; }
152  ampf& operator= (const char *s) { mpfr_strtofr(getWritePtr(), s, NULL, 0, GMP_RNDN); return *this; }
153  ampf& operator= (const std::string &s) { mpfr_strtofr(getWritePtr(), s.c_str(), NULL, 0, GMP_RNDN); return *this; }
154  ampf& operator= (const ampf& r)
155  {
156  // TODO: may be copy ref
157  if( this==&r )
158  return *this;
159  if( rval==r.rval )
160  return *this;
161  rval->refCount--;
162  if( rval->refCount==0 )
164  rval = r.rval;
165  rval->refCount++;
166  //mpfr_set(getWritePtr(), r.getReadPtr(), GMP_RNDN);
167  return *this;
168  }
169 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
170  template<unsigned int Precision2>
172  {
173  if( (void*)this==(void*)(&r) )
174  return *this;
175  mpfr_set(getWritePtr(), r.getReadPtr(), GMP_RNDN);
176  return *this;
177  }
178 #endif
179 
180  //
181  // in-place operators
182  // TODO: optimize
183  //
184  template<class T> ampf& operator+=(const T& v){ *this = *this + v; return *this; };
185  template<class T> ampf& operator-=(const T& v){ *this = *this - v; return *this; };
186  template<class T> ampf& operator*=(const T& v){ *this = *this * v; return *this; };
187  template<class T> ampf& operator/=(const T& v){ *this = *this / v; return *this; };
188 
189  //
190  // MPFR access
191  //
192  mpfr_srcptr getReadPtr() const;
193  mpfr_ptr getWritePtr();
194 
195  //
196  // properties and information
197  //
198  bool isFiniteNumber() const;
199  bool isPositiveNumber() const;
200  bool isZero() const;
201  bool isNegativeNumber() const;
202  const ampf getUlpOf();
203 
204  //
205  // conversions
206  //
207  double toDouble() const;
208  std::string toHex() const;
209  std::string toDec() const;
210  char * toString() const;
211 
212 
213  //
214  // static methods
215  //
216  static const ampf getUlpOf(const ampf &x);
217  static const ampf getUlp();
218  static const ampf getUlp256();
219  static const ampf getUlp512();
220  static const ampf getMaxNumber();
221  static const ampf getMinNumber();
222  static const ampf getAlgoPascalEpsilon();
223  static const ampf getAlgoPascalMaxNumber();
224  static const ampf getAlgoPascalMinNumber();
225  static const ampf getRandom();
226  private:
227  void CheckPrecision();
228  void InitializeAsZero();
229  void InitializeAsSLong(signed long v);
230  void InitializeAsULong(unsigned long v);
231  void InitializeAsDouble(long double v);
232  void InitializeAsString(const char *s);
233 
234  //mpfr_reference ref;
236  };
237 
238  /*void ampf<Precision>::CheckPrecision()
239  {
240  if( Precision<32 )
241  throw incorrectPrecision();
242  }***/
243 
244  template<unsigned int Precision>
246  {
247  if( Precision<32 )
248  throw incorrectPrecision();
249  }
250 
251  template<unsigned int Precision>
253  {
254  CheckPrecision();
255  rval = mpfr_storage::newMpfr(Precision);
256  mpfr_set_ui(getWritePtr(), 0, GMP_RNDN);
257  }
258 
259  template<unsigned int Precision>
261  {
262  CheckPrecision();
263  rval = mpfr_storage::newMpfr(Precision);
264  mpfr_set_si(getWritePtr(), sv, GMP_RNDN);
265  }
266 
267  template<unsigned int Precision>
269  {
270  CheckPrecision();
271  rval = mpfr_storage::newMpfr(Precision);
272  mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
273  }
274 
275  template<unsigned int Precision>
277  {
278  CheckPrecision();
279  rval = mpfr_storage::newMpfr(Precision);
280  mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
281  }
282 
283  template<unsigned int Precision>
285  {
286  CheckPrecision();
287  rval = mpfr_storage::newMpfr(Precision);
288  mpfr_strtofr(getWritePtr(), s, NULL, 0, GMP_RNDN);
289  }
290 
291  template<unsigned int Precision>
292  mpfr_srcptr ampf<Precision>::getReadPtr() const
293  {
294  // 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
295  // getWrite ñîçäàâàëàñü íîâàÿ instance mpfr_t.
296  // ýòî ìîæåò áûòü íóæíî äëÿ êîððåêòíîé îáðàáîòêè ñèòóàöèé âèäà
297  // mpfr_÷åãî_òî_òàì( a.getWritePtr(), a.getReadPtr())
298  // âðîäå áû íóæíî, à òî åñëè òàì çàâÿçàíî íà side-effects...
299  return rval->value;
300  }
301 
302  template<unsigned int Precision>
304  {
305  if( rval->refCount==1 )
306  return rval->value;
307  mpfr_record *newrval = mpfr_storage::newMpfr(Precision);
308  mpfr_set(newrval->value, rval->value, GMP_RNDN);
309  rval->refCount--;
310  rval = newrval;
311  return rval->value;
312  }
313 
314  template<unsigned int Precision>
316  {
317  return mpfr_number_p(getReadPtr())!=0;
318  }
319 
320  template<unsigned int Precision>
322  {
323  if( !isFiniteNumber() )
324  return false;
325  return mpfr_sgn(getReadPtr())>0;
326  }
327 
328  template<unsigned int Precision>
330  {
331  return mpfr_zero_p(getReadPtr())!=0;
332  }
333 
334  template<unsigned int Precision>
336  {
337  if( !isFiniteNumber() )
338  return false;
339  return mpfr_sgn(getReadPtr())<0;
340  }
341 
342  template<unsigned int Precision>
344  {
345  return getUlpOf(*this);
346  }
347 
348  template<unsigned int Precision>
350  {
351  return mpfr_get_d(getReadPtr(), GMP_RNDN);
352  }
353 
354  template<unsigned int Precision>
356  {
357  //
358  // some special cases
359  //
360  if( !isFiniteNumber() )
361  {
362  std::string r;
363  mp_exp_t _e;
364  char *ptr;
365  ptr = mpfr_get_str(NULL, &_e, 16, 0, getReadPtr(), GMP_RNDN);
366  r = ptr;
367  mpfr_free_str(ptr);
368  return r;
369  }
370 
371  //
372  // general case
373  //
374  std::string r;
375  char buf_e[128];
376  signed long iexpval;
377  mp_exp_t expval;
378  char *ptr;
379  char *ptr2;
380  ptr = mpfr_get_str(NULL, &expval, 16, 0, getReadPtr(), GMP_RNDN);
381  ptr2 = ptr;
382  iexpval = expval;
383  if( iexpval!=expval )
384  throw internalError();
385  sprintf(buf_e, "%ld", long(iexpval));
386  if( *ptr=='-' )
387  {
388  r = "-";
389  ptr++;
390  }
391  r += "0x0.";
392  r += ptr;
393  r += "@";
394  r += buf_e;
395  mpfr_free_str(ptr2);
396  return r;
397  }
398 
399  template<unsigned int Precision>
401  {
402  // TODO: advanced output formatting (zero, integers)
403 
404  //
405  // some special cases
406  //
407  if( !isFiniteNumber() )
408  {
409  std::string r;
410  mp_exp_t _e;
411  char *ptr;
412  ptr = mpfr_get_str(NULL, &_e, 10, 0, getReadPtr(), GMP_RNDN);
413  r = ptr;
414  mpfr_free_str(ptr);
415  return r;
416  }
417 
418  //
419  // general case
420  //
421  std::string r;
422  char buf_e[128];
423  signed long iexpval;
424  mp_exp_t expval;
425  char *ptr;
426  char *ptr2;
427  ptr = mpfr_get_str(NULL, &expval, 10, 0, getReadPtr(), GMP_RNDN);
428  ptr2 = ptr;
429  iexpval = expval;
430  if( iexpval!=expval )
431  throw internalError();
432  sprintf(buf_e, "%ld", long(iexpval));
433  if( *ptr=='-' )
434  {
435  r = "-";
436  ptr++;
437  }
438  r += "0.";
439  r += ptr;
440  r += "E";
441  r += buf_e;
442  mpfr_free_str(ptr2);
443  return r;
444  }
445  template<unsigned int Precision>
447  {
448  char *toString_Block=(char *)omAlloc(256);
449  //
450  // some special cases
451  //
452  if( !isFiniteNumber() )
453  {
454  mp_exp_t _e;
455  char *ptr;
456  ptr = mpfr_get_str(NULL, &_e, 10, 0, getReadPtr(), GMP_RNDN);
457  strcpy(toString_Block, ptr);
458  mpfr_free_str(ptr);
459  return toString_Block;
460  }
461 
462  //
463  // general case
464  //
465 
466  char buf_e[128];
467  signed long iexpval;
468  mp_exp_t expval;
469  char *ptr;
470  char *ptr2;
471  ptr = mpfr_get_str(NULL, &expval, 10, 0, getReadPtr(), GMP_RNDN);
472  ptr2 = ptr;
473  iexpval = expval;
474  if( iexpval!=expval )
475  throw internalError();
476  sprintf(buf_e, "%ld", long(iexpval));
477  if( *ptr=='-' )
478  {
479  ptr++;
480  sprintf(toString_Block,"-0.%sE%s",ptr,buf_e);
481  }
482  else
483  sprintf(toString_Block,"0.%sE%s",ptr,buf_e);
484  mpfr_free_str(ptr2);
485  return toString_Block;
486  }
487 
488  template<unsigned int Precision>
490  {
491  if( !x.isFiniteNumber() )
492  return x;
493  if( x.isZero() )
494  return x;
495  ampf<Precision> r(1);
496  mpfr_nextabove(r.getWritePtr());
497  mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
498  mpfr_mul_2si(
499  r.getWritePtr(),
500  r.getWritePtr(),
501  mpfr_get_exp(x.getReadPtr()),
502  GMP_RNDN);
503  mpfr_div_2si(
504  r.getWritePtr(),
505  r.getWritePtr(),
506  1,
507  GMP_RNDN);
508  return r;
509  }
510 
511  template<unsigned int Precision>
513  {
514  ampf<Precision> r(1);
515  mpfr_nextabove(r.getWritePtr());
516  mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
517  return r;
518  }
519 
520  template<unsigned int Precision>
522  {
523  ampf<Precision> r(1);
524  mpfr_nextabove(r.getWritePtr());
525  mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
526  mpfr_mul_2si(
527  r.getWritePtr(),
528  r.getWritePtr(),
529  8,
530  GMP_RNDN);
531  return r;
532  }
533 
534  template<unsigned int Precision>
536  {
537  ampf<Precision> r(1);
538  mpfr_nextabove(r.getWritePtr());
539  mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
540  mpfr_mul_2si(
541  r.getWritePtr(),
542  r.getWritePtr(),
543  9,
544  GMP_RNDN);
545  return r;
546  }
547 
548  template<unsigned int Precision>
550  {
551  ampf<Precision> r(1);
552  mpfr_nextbelow(r.getWritePtr());
553  mpfr_set_exp(r.getWritePtr(),mpfr_get_emax());
554  return r;
555  }
556 
557  template<unsigned int Precision>
559  {
560  ampf<Precision> r(1);
561  mpfr_set_exp(r.getWritePtr(),mpfr_get_emin());
562  return r;
563  }
564 
565  template<unsigned int Precision>
567  {
568  return getUlp256();
569  }
570 
571  template<unsigned int Precision>
573  {
574  ampf<Precision> r(1);
575  mp_exp_t e1 = mpfr_get_emax();
576  mp_exp_t e2 = -mpfr_get_emin();
577  mp_exp_t e = e1>e2 ? e1 : e2;
578  mpfr_set_exp(r.getWritePtr(), e-5);
579  return r;
580  }
581 
582  template<unsigned int Precision>
584  {
585  ampf<Precision> r(1);
586  mp_exp_t e1 = mpfr_get_emax();
587  mp_exp_t e2 = -mpfr_get_emin();
588  mp_exp_t e = e1>e2 ? e1 : e2;
589  mpfr_set_exp(r.getWritePtr(), 2-(e-5));
590  return r;
591  }
592 
593  template<unsigned int Precision>
595  {
596  ampf<Precision> r;
597  while(mpfr_urandomb(r.getWritePtr(), *amp::mpfr_storage::getRandState()));
598  return r;
599  }
600 
601  //
602  // comparison operators
603  //
604  template<unsigned int Precision>
605  const bool operator==(const ampf<Precision>& op1, const ampf<Precision>& op2)
606  {
607  return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())==0;
608  }
609 
610  template<unsigned int Precision>
611  const bool operator!=(const ampf<Precision>& op1, const ampf<Precision>& op2)
612  {
613  return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())!=0;
614  }
615 
616  template<unsigned int Precision>
617  const bool operator<(const ampf<Precision>& op1, const ampf<Precision>& op2)
618  {
619  return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())<0;
620  }
621 
622  template<unsigned int Precision>
623  const bool operator>(const ampf<Precision>& op1, const ampf<Precision>& op2)
624  {
625  return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())>0;
626  }
627 
628  template<unsigned int Precision>
629  const bool operator<=(const ampf<Precision>& op1, const ampf<Precision>& op2)
630  {
631  return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())<=0;
632  }
633 
634  template<unsigned int Precision>
635  const bool operator>=(const ampf<Precision>& op1, const ampf<Precision>& op2)
636  {
637  return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())>=0;
638  }
639 
640  //
641  // arithmetic operators
642  //
643  template<unsigned int Precision>
645  {
646  return op1;
647  }
648 
649  template<unsigned int Precision>
651  {
652  mpfr_record *v = mpfr_storage::newMpfr(Precision);
653  mpfr_neg(v->value, op1.getReadPtr(), GMP_RNDN);
654  return v;
655  }
656 
657  template<unsigned int Precision>
659  {
660  mpfr_record *v = mpfr_storage::newMpfr(Precision);
661  mpfr_add(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
662  return v;
663  }
664 
665  template<unsigned int Precision>
667  {
668  mpfr_record *v = mpfr_storage::newMpfr(Precision);
669  mpfr_sub(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
670  return v;
671  }
672 
673 
674  template<unsigned int Precision>
676  {
677  mpfr_record *v = mpfr_storage::newMpfr(Precision);
678  mpfr_mul(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
679  return v;
680  }
681 
682  template<unsigned int Precision>
684  {
685  mpfr_record *v = mpfr_storage::newMpfr(Precision);
686  mpfr_div(v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
687  return v;
688  }
689 
690  //
691  // basic functions
692  //
693  template<unsigned int Precision>
695  {
696  // TODO: optimize temporary for return value
698  mpfr_sqr(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
699  return res;
700  }
701 
702  template<unsigned int Precision>
703  const int sign(const ampf<Precision> &x)
704  {
705  int s = mpfr_sgn(x.getReadPtr());
706  if( s>0 )
707  return +1;
708  if( s<0 )
709  return -1;
710  return 0;
711  }
712 
713  template<unsigned int Precision>
715  {
716  // TODO: optimize temporary for return value
718  mpfr_abs(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
719  return res;
720  }
721 
722  template<unsigned int Precision>
724  {
725  // TODO: optimize temporary for return value
727  mpfr_max(res.getWritePtr(), x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
728  return res;
729  }
730 
731  template<unsigned int Precision>
733  {
734  // TODO: optimize temporary for return value
736  mpfr_min(res.getWritePtr(), x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
737  return res;
738  }
739 
740  template<unsigned int Precision>
742  {
743  // TODO: optimize temporary for return value
745  mpfr_sqrt(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
746  return res;
747  }
748 
749  template<unsigned int Precision>
750  const signed long trunc(const ampf<Precision> &x)
751  {
752  ampf<Precision> tmp;
753  signed long r;
754  mpfr_trunc(tmp.getWritePtr(), x.getReadPtr());
755  if( mpfr_integer_p(tmp.getReadPtr())==0 )
756  throw invalidConversion();
757  mpfr_clear_erangeflag();
758  r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
759  if( mpfr_erangeflag_p()!=0 )
760  throw invalidConversion();
761  return r;
762  }
763 
764  template<unsigned int Precision>
766  {
767  // TODO: optimize temporary for return value
768  ampf<Precision> r;
769  mpfr_frac(r.getWritePtr(), x.getReadPtr(), GMP_RNDN);
770  return r;
771  }
772 
773  template<unsigned int Precision>
774  const signed long floor(const ampf<Precision> &x)
775  {
776  ampf<Precision> tmp;
777  signed long r;
778  mpfr_floor(tmp.getWritePtr(), x.getReadPtr());
779  if( mpfr_integer_p(tmp.getReadPtr())==0 )
780  throw invalidConversion();
781  mpfr_clear_erangeflag();
782  r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
783  if( mpfr_erangeflag_p()!=0 )
784  throw invalidConversion();
785  return r;
786  }
787 
788  template<unsigned int Precision>
789  const signed long ceil(const ampf<Precision> &x)
790  {
791  ampf<Precision> tmp;
792  signed long r;
793  mpfr_ceil(tmp.getWritePtr(), x.getReadPtr());
794  if( mpfr_integer_p(tmp.getReadPtr())==0 )
795  throw invalidConversion();
796  mpfr_clear_erangeflag();
797  r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
798  if( mpfr_erangeflag_p()!=0 )
799  throw invalidConversion();
800  return r;
801  }
802 
803  template<unsigned int Precision>
804  const signed long round(const ampf<Precision> &x)
805  {
806  ampf<Precision> tmp;
807  signed long r;
808  mpfr_round(tmp.getWritePtr(), x.getReadPtr());
809  if( mpfr_integer_p(tmp.getReadPtr())==0 )
810  throw invalidConversion();
811  mpfr_clear_erangeflag();
812  r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
813  if( mpfr_erangeflag_p()!=0 )
814  throw invalidConversion();
815  return r;
816  }
817 
818  template<unsigned int Precision>
820  {
821  // TODO: optimize temporary for return value
822  ampf<Precision> r;
823  if( !x.isFiniteNumber() )
824  throw invalidConversion();
825  if( x.isZero() )
826  {
827  *exponent = 0;
828  r = 0;
829  return r;
830  }
831  r = x;
832  *exponent = mpfr_get_exp(r.getReadPtr());
833  mpfr_set_exp(r.getWritePtr(),0);
834  return r;
835  }
836 
837  template<unsigned int Precision>
839  {
840  // TODO: optimize temporary for return value
841  ampf<Precision> r;
842  mpfr_mul_2si(r.getWritePtr(), x.getReadPtr(), exponent, GMP_RNDN);
843  return r;
844  }
845 
846  //
847  // different types of arguments
848  //
849  #define __AMP_BINARY_OPI(type) \
850  template<unsigned int Precision> const ampf<Precision> operator+(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \
851  template<unsigned int Precision> const ampf<Precision> operator+(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \
852  template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const signed type& op2) { return op1+ampf<Precision>(op2); } \
853  template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const unsigned type& op2) { return op1+ampf<Precision>(op2); } \
854  template<unsigned int Precision> const ampf<Precision> operator-(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \
855  template<unsigned int Precision> const ampf<Precision> operator-(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \
856  template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const signed type& op2) { return op1-ampf<Precision>(op2); } \
857  template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const unsigned type& op2) { return op1-ampf<Precision>(op2); } \
858  template<unsigned int Precision> const ampf<Precision> operator*(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \
859  template<unsigned int Precision> const ampf<Precision> operator*(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \
860  template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const signed type& op2) { return op1*ampf<Precision>(op2); } \
861  template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const unsigned type& op2) { return op1*ampf<Precision>(op2); } \
862  template<unsigned int Precision> const ampf<Precision> operator/(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \
863  template<unsigned int Precision> const ampf<Precision> operator/(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \
864  template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const signed type& op2) { return op1/ampf<Precision>(op2); } \
865  template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const unsigned type& op2) { return op1/ampf<Precision>(op2); } \
866  template<unsigned int Precision> const bool operator==(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \
867  template<unsigned int Precision> const bool operator==(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \
868  template<unsigned int Precision> const bool operator==(const ampf<Precision>& op1, const signed type& op2) { return op1==ampf<Precision>(op2); } \
869  template<unsigned int Precision> const bool operator==(const ampf<Precision>& op1, const unsigned type& op2) { return op1==ampf<Precision>(op2); } \
870  template<unsigned int Precision> const bool operator!=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \
871  template<unsigned int Precision> const bool operator!=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \
872  template<unsigned int Precision> const bool operator!=(const ampf<Precision>& op1, const signed type& op2) { return op1!=ampf<Precision>(op2); } \
873  template<unsigned int Precision> const bool operator!=(const ampf<Precision>& op1, const unsigned type& op2) { return op1!=ampf<Precision>(op2); } \
874  template<unsigned int Precision> const bool operator<=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \
875  template<unsigned int Precision> const bool operator<=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \
876  template<unsigned int Precision> const bool operator<=(const ampf<Precision>& op1, const signed type& op2) { return op1<=ampf<Precision>(op2); } \
877  template<unsigned int Precision> const bool operator<=(const ampf<Precision>& op1, const unsigned type& op2) { return op1<=ampf<Precision>(op2); } \
878  template<unsigned int Precision> const bool operator>=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \
879  template<unsigned int Precision> const bool operator>=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \
880  template<unsigned int Precision> const bool operator>=(const ampf<Precision>& op1, const signed type& op2) { return op1>=ampf<Precision>(op2); } \
881  template<unsigned int Precision> const bool operator>=(const ampf<Precision>& op1, const unsigned type& op2) { return op1>=ampf<Precision>(op2); } \
882  template<unsigned int Precision> const bool operator<(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \
883  template<unsigned int Precision> const bool operator<(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \
884  template<unsigned int Precision> const bool operator<(const ampf<Precision>& op1, const signed type& op2) { return op1<ampf<Precision>(op2); } \
885  template<unsigned int Precision> const bool operator<(const ampf<Precision>& op1, const unsigned type& op2) { return op1<ampf<Precision>(op2); } \
886  template<unsigned int Precision> const bool operator>(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \
887  template<unsigned int Precision> const bool operator>(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \
888  template<unsigned int Precision> const bool operator>(const ampf<Precision>& op1, const signed type& op2) { return op1>ampf<Precision>(op2); } \
889  template<unsigned int Precision> const bool operator>(const ampf<Precision>& op1, const unsigned type& op2) { return op1>ampf<Precision>(op2); }
891  __AMP_BINARY_OPI(short)
892  __AMP_BINARY_OPI(long)
893  __AMP_BINARY_OPI(int)
894  #undef __AMP_BINARY_OPI
895  #define __AMP_BINARY_OPF(type) \
896  template<unsigned int Precision> const ampf<Precision> operator+(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \
897  template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const type& op2) { return op1+ampf<Precision>(op2); } \
898  template<unsigned int Precision> const ampf<Precision> operator-(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \
899  template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const type& op2) { return op1-ampf<Precision>(op2); } \
900  template<unsigned int Precision> const ampf<Precision> operator*(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \
901  template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const type& op2) { return op1*ampf<Precision>(op2); } \
902  template<unsigned int Precision> const ampf<Precision> operator/(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \
903  template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const type& op2) { return op1/ampf<Precision>(op2); } \
904  template<unsigned int Precision> bool operator==(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \
905  template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const type& op2) { return op1==ampf<Precision>(op2); } \
906  template<unsigned int Precision> bool operator!=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \
907  template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const type& op2) { return op1!=ampf<Precision>(op2); } \
908  template<unsigned int Precision> bool operator<=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \
909  template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const type& op2) { return op1<=ampf<Precision>(op2); } \
910  template<unsigned int Precision> bool operator>=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \
911  template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const type& op2) { return op1>=ampf<Precision>(op2); } \
912  template<unsigned int Precision> bool operator<(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \
913  template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const type& op2) { return op1<ampf<Precision>(op2); } \
914  template<unsigned int Precision> bool operator>(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \
915  template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const type& op2) { return op1>ampf<Precision>(op2); }
916  __AMP_BINARY_OPF(float)
917  __AMP_BINARY_OPF(double)
918  __AMP_BINARY_OPF(long double)
919  #undef __AMP_BINARY_OPF
920 
921  //
922  // transcendent functions
923  //
924  template<unsigned int Precision>
925  const ampf<Precision> pi()
926  {
927  mpfr_record *v = mpfr_storage::newMpfr(Precision);
928  mpfr_const_pi(v->value, GMP_RNDN);
929  return v;
930  }
931 
932  template<unsigned int Precision>
934  {
935  mpfr_record *v = mpfr_storage::newMpfr(Precision);
936  mpfr_const_pi(v->value, GMP_RNDN);
937  mpfr_mul_2si(v->value, v->value, -1, GMP_RNDN);
938  return v;
939  }
940 
941  template<unsigned int Precision>
943  {
944  mpfr_record *v = mpfr_storage::newMpfr(Precision);
945  mpfr_const_pi(v->value, GMP_RNDN);
946  mpfr_mul_2si(v->value, v->value, +1, GMP_RNDN);
947  return v;
948  }
949 
950  template<unsigned int Precision>
952  {
953  mpfr_record *v = mpfr_storage::newMpfr(Precision);
954  mpfr_sin(v->value, x.getReadPtr(), GMP_RNDN);
955  return v;
956  }
957 
958  template<unsigned int Precision>
960  {
961  mpfr_record *v = mpfr_storage::newMpfr(Precision);
962  mpfr_cos(v->value, x.getReadPtr(), GMP_RNDN);
963  return v;
964  }
965 
966  template<unsigned int Precision>
968  {
969  mpfr_record *v = mpfr_storage::newMpfr(Precision);
970  mpfr_tan(v->value, x.getReadPtr(), GMP_RNDN);
971  return v;
972  }
973 
974  template<unsigned int Precision>
976  {
977  mpfr_record *v = mpfr_storage::newMpfr(Precision);
978  mpfr_asin(v->value, x.getReadPtr(), GMP_RNDN);
979  return v;
980  }
981 
982  template<unsigned int Precision>
984  {
985  mpfr_record *v = mpfr_storage::newMpfr(Precision);
986  mpfr_acos(v->value, x.getReadPtr(), GMP_RNDN);
987  return v;
988  }
989 
990  template<unsigned int Precision>
992  {
993  mpfr_record *v = mpfr_storage::newMpfr(Precision);
994  mpfr_atan(v->value, x.getReadPtr(), GMP_RNDN);
995  return v;
996  }
997 
998  template<unsigned int Precision>
1000  {
1001  mpfr_record *v = mpfr_storage::newMpfr(Precision);
1002  mpfr_atan2(v->value, y.getReadPtr(), x.getReadPtr(), GMP_RNDN);
1003  return v;
1004  }
1005 
1006  template<unsigned int Precision>
1008  {
1009  mpfr_record *v = mpfr_storage::newMpfr(Precision);
1010  mpfr_log(v->value, x.getReadPtr(), GMP_RNDN);
1011  return v;
1012  }
1013 
1014  template<unsigned int Precision>
1016  {
1017  mpfr_record *v = mpfr_storage::newMpfr(Precision);
1018  mpfr_log2(v->value, x.getReadPtr(), GMP_RNDN);
1019  return v;
1020  }
1021 
1022  template<unsigned int Precision>
1024  {
1025  mpfr_record *v = mpfr_storage::newMpfr(Precision);
1026  mpfr_log10(v->value, x.getReadPtr(), GMP_RNDN);
1027  return v;
1028  }
1029 
1030  template<unsigned int Precision>
1032  {
1033  mpfr_record *v = mpfr_storage::newMpfr(Precision);
1034  mpfr_exp(v->value, x.getReadPtr(), GMP_RNDN);
1035  return v;
1036  }
1037 
1038  template<unsigned int Precision>
1040  {
1041  mpfr_record *v = mpfr_storage::newMpfr(Precision);
1042  mpfr_sinh(v->value, x.getReadPtr(), GMP_RNDN);
1043  return v;
1044  }
1045 
1046  template<unsigned int Precision>
1048  {
1049  mpfr_record *v = mpfr_storage::newMpfr(Precision);
1050  mpfr_cosh(v->value, x.getReadPtr(), GMP_RNDN);
1051  return v;
1052  }
1053 
1054  template<unsigned int Precision>
1056  {
1057  mpfr_record *v = mpfr_storage::newMpfr(Precision);
1058  mpfr_tanh(v->value, x.getReadPtr(), GMP_RNDN);
1059  return v;
1060  }
1061 
1062  template<unsigned int Precision>
1064  {
1065  mpfr_record *v = mpfr_storage::newMpfr(Precision);
1066  mpfr_pow(v->value, x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
1067  return v;
1068  }
1069 
1070  //
1071  // complex ampf
1072  //
1073  template<unsigned int Precision>
1074  class campf
1075  {
1076  public:
1077  campf():x(0),y(0){};
1078  campf(long double v) { x=v; y=0; }
1079  campf(double v) { x=v; y=0; }
1080  campf(float v) { x=v; y=0; }
1081  campf(signed long v) { x=v; y=0; }
1082  campf(unsigned long v) { x=v; y=0; }
1083  campf(signed int v) { x=v; y=0; }
1084  campf(unsigned int v) { x=v; y=0; }
1085  campf(signed short v) { x=v; y=0; }
1086  campf(unsigned short v) { x=v; y=0; }
1087  campf(signed char v) { x=v; y=0; }
1088  campf(unsigned char v) { x=v; y=0; }
1089  campf(const ampf<Precision> &_x):x(_x),y(0){};
1090  campf(const ampf<Precision> &_x, const ampf<Precision> &_y):x(_x),y(_y){};
1091  campf(const campf &z):x(z.x),y(z.y){};
1092 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
1093  template<unsigned int Prec2>
1094  campf(const campf<Prec2> &z):x(z.x),y(z.y){};
1095 #endif
1096 
1097  campf& operator= (long double v) { x=v; y=0; return *this; }
1098  campf& operator= (double v) { x=v; y=0; return *this; }
1099  campf& operator= (float v) { x=v; y=0; return *this; }
1100  campf& operator= (signed long v) { x=v; y=0; return *this; }
1101  campf& operator= (unsigned long v) { x=v; y=0; return *this; }
1102  campf& operator= (signed int v) { x=v; y=0; return *this; }
1103  campf& operator= (unsigned int v) { x=v; y=0; return *this; }
1104  campf& operator= (signed short v) { x=v; y=0; return *this; }
1105  campf& operator= (unsigned short v) { x=v; y=0; return *this; }
1106  campf& operator= (signed char v) { x=v; y=0; return *this; }
1107  campf& operator= (unsigned char v) { x=v; y=0; return *this; }
1108  campf& operator= (const char *s) { x=s; y=0; return *this; }
1109  campf& operator= (const std::string &s) { x=s; y=0; return *this; }
1111  {
1112  x = r.x;
1113  y = r.y;
1114  return *this;
1115  }
1116 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
1117  template<unsigned int Precision2>
1119  {
1120  x = r.x;
1121  y = r.y;
1122  return *this;
1123  }
1124 #endif
1125 
1127  };
1128 
1129  //
1130  // complex operations
1131  //
1132  template<unsigned int Precision>
1133  const bool operator==(const campf<Precision>& lhs, const campf<Precision>& rhs)
1134  { return lhs.x==rhs.x && lhs.y==rhs.y; }
1135 
1136  template<unsigned int Precision>
1137  const bool operator!=(const campf<Precision>& lhs, const campf<Precision>& rhs)
1138  { return lhs.x!=rhs.x || lhs.y!=rhs.y; }
1139 
1140  template<unsigned int Precision>
1142  { return lhs; }
1143 
1144  template<unsigned int Precision>
1146  { lhs.x += rhs.x; lhs.y += rhs.y; return lhs; }
1147 
1148  template<unsigned int Precision>
1150  { campf<Precision> r = lhs; r += rhs; return r; }
1151 
1152  template<unsigned int Precision>
1154  { return campf<Precision>(-lhs.x, -lhs.y); }
1155 
1156  template<unsigned int Precision>
1158  { lhs.x -= rhs.x; lhs.y -= rhs.y; return lhs; }
1159 
1160  template<unsigned int Precision>
1162  { campf<Precision> r = lhs; r -= rhs; return r; }
1163 
1164  template<unsigned int Precision>
1166  {
1167  ampf<Precision> xx(lhs.x*rhs.x), yy(lhs.y*rhs.y), mm((lhs.x+lhs.y)*(rhs.x+rhs.y));
1168  lhs.x = xx-yy;
1169  lhs.y = mm-xx-yy;
1170  return lhs;
1171  }
1172 
1173  template<unsigned int Precision>
1175  { campf<Precision> r = lhs; r *= rhs; return r; }
1176 
1177  template<unsigned int Precision>
1179  {
1181  ampf<Precision> e;
1183  if( abs(rhs.y)<abs(rhs.x) )
1184  {
1185  e = rhs.y/rhs.x;
1186  f = rhs.x+rhs.y*e;
1187  result.x = (lhs.x+lhs.y*e)/f;
1188  result.y = (lhs.y-lhs.x*e)/f;
1189  }
1190  else
1191  {
1192  e = rhs.x/rhs.y;
1193  f = rhs.y+rhs.x*e;
1194  result.x = (lhs.y+lhs.x*e)/f;
1195  result.y = (-lhs.x+lhs.y*e)/f;
1196  }
1197  return result;
1198  }
1199 
1200  template<unsigned int Precision>
1202  {
1203  lhs = lhs/rhs;
1204  return lhs;
1205  }
1206 
1207  template<unsigned int Precision>
1209  {
1210  ampf<Precision> w, xabs, yabs, v;
1211 
1212  xabs = abs(z.x);
1213  yabs = abs(z.y);
1214  w = xabs>yabs ? xabs : yabs;
1215  v = xabs<yabs ? xabs : yabs;
1216  if( v==0 )
1217  return w;
1218  else
1219  {
1220  ampf<Precision> t = v/w;
1221  return w*sqrt(1+sqr(t));
1222  }
1223  }
1224 
1225  template<unsigned int Precision>
1227  {
1228  return campf<Precision>(z.x, -z.y);
1229  }
1230 
1231  template<unsigned int Precision>
1233  {
1234  ampf<Precision> t = z.x*z.y;
1235  return campf<Precision>(sqr(z.x)-sqr(z.y), t+t);
1236  }
1237 
1238  //
1239  // different types of arguments
1240  //
1241  #define __AMP_BINARY_OPI(type) \
1242  template<unsigned int Precision> const campf<Precision> operator+ (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \
1243  template<unsigned int Precision> const campf<Precision> operator+ (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \
1244  template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \
1245  template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \
1246  template<unsigned int Precision> const campf<Precision> operator- (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \
1247  template<unsigned int Precision> const campf<Precision> operator- (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \
1248  template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \
1249  template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \
1250  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); } \
1251  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); } \
1252  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); } \
1253  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); } \
1254  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; } \
1255  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; } \
1256  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); } \
1257  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); } \
1258  template<unsigned int Precision> bool operator==(const signed type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \
1259  template<unsigned int Precision> bool operator==(const unsigned type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \
1260  template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const signed type& op2) { return op1.x==op2 && op1.y==0; } \
1261  template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const unsigned type& op2) { return op1.x==op2 && op1.y==0; } \
1262  template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const signed type& op2) { return op1.x!=op2 || op1.y!=0; } \
1263  template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const unsigned type& op2) { return op1.x!=op2 || op1.y!=0; } \
1264  template<unsigned int Precision> bool operator!=(const signed type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } \
1265  template<unsigned int Precision> bool operator!=(const unsigned type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; }
1266  __AMP_BINARY_OPI(char)
1267  __AMP_BINARY_OPI(short)
1268  __AMP_BINARY_OPI(long)
1269  __AMP_BINARY_OPI(int)
1270  #undef __AMP_BINARY_OPI
1271  #define __AMP_BINARY_OPF(type) \
1272  template<unsigned int Precision> const campf<Precision> operator+ (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \
1273  template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \
1274  template<unsigned int Precision> const campf<Precision> operator- (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \
1275  template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \
1276  template<unsigned int Precision> const campf<Precision> operator* (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \
1277  template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \
1278  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; } \
1279  template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \
1280  template<unsigned int Precision> bool operator==(const type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \
1281  template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const type& op2) { return op1.x==op2 && op1.y==0; } \
1282  template<unsigned int Precision> bool operator!=(const type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } \
1283  template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const type& op2) { return op1.x!=op2 || op1.y!=0; }
1284  __AMP_BINARY_OPF(float)
1285  __AMP_BINARY_OPF(double)
1286  __AMP_BINARY_OPF(long double)
1287  __AMP_BINARY_OPF(ampf<Precision>)
1288  #undef __AMP_BINARY_OPF
1289 
1290  //
1291  // Real linear algebra
1292  //
1293  template<unsigned int Precision>
1294  ampf<Precision> vDotProduct(ap::const_raw_vector< ampf<Precision> > v1, ap::const_raw_vector< ampf<Precision> > v2)
1295  {
1296  ap::ap_error::make_assertion(v1.GetLength()==v2.GetLength());
1297  int i, cnt = v1.GetLength();
1298  const ampf<Precision> *p1 = v1.GetData();
1299  const ampf<Precision> *p2 = v2.GetData();
1300  mpfr_record *r = NULL;
1301  mpfr_record *t = NULL;
1302  try
1303  {
1304  r = mpfr_storage::newMpfr(Precision);
1305  t = mpfr_storage::newMpfr(Precision);
1306  mpfr_set_ui(r->value, 0, GMP_RNDN);
1307  for(i=0; i<cnt; i++)
1308  {
1309  mpfr_mul(t->value, p1->getReadPtr(), p2->getReadPtr(), GMP_RNDN);
1310  mpfr_add(r->value, r->value, t->value, GMP_RNDN);
1311  p1 += v1.GetStep();
1312  p2 += v2.GetStep();
1313  }
1315  return r;
1316  }
1317  catch(...)
1318  {
1319  if( r!=NULL )
1321  if( t!=NULL )
1323  throw;
1324  }
1325  }
1326 
1327  template<unsigned int Precision>
1329  {
1330  ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
1331  int i, cnt = vDst.GetLength();
1332  ampf<Precision> *pDst = vDst.GetData();
1333  const ampf<Precision> *pSrc = vSrc.GetData();
1334  if( pDst==pSrc )
1335  return;
1336  for(i=0; i<cnt; i++)
1337  {
1338  *pDst = *pSrc;
1339  pDst += vDst.GetStep();
1340  pSrc += vSrc.GetStep();
1341  }
1342  }
1343 
1344  template<unsigned int Precision>
1346  {
1347  ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
1348  int i, cnt = vDst.GetLength();
1349  ampf<Precision> *pDst = vDst.GetData();
1350  const ampf<Precision> *pSrc = vSrc.GetData();
1351  for(i=0; i<cnt; i++)
1352  {
1353  *pDst = *pSrc;
1354  mpfr_ptr v = pDst->getWritePtr();
1355  mpfr_neg(v, v, GMP_RNDN);
1356  pDst += vDst.GetStep();
1357  pSrc += vSrc.GetStep();
1358  }
1359  }
1360 
1361  template<unsigned int Precision, class T2>
1363  {
1364  ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
1365  int i, cnt = vDst.GetLength();
1366  ampf<Precision> *pDst = vDst.GetData();
1367  const ampf<Precision> *pSrc = vSrc.GetData();
1369  for(i=0; i<cnt; i++)
1370  {
1371  *pDst = *pSrc;
1372  mpfr_ptr v = pDst->getWritePtr();
1373  mpfr_mul(v, v, a.getReadPtr(), GMP_RNDN);
1374  pDst += vDst.GetStep();
1375  pSrc += vSrc.GetStep();
1376  }
1377  }
1378 
1379  template<unsigned int Precision>
1381  {
1382  ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
1383  int i, cnt = vDst.GetLength();
1384  ampf<Precision> *pDst = vDst.GetData();
1385  const ampf<Precision> *pSrc = vSrc.GetData();
1386  for(i=0; i<cnt; i++)
1387  {
1388  mpfr_ptr v = pDst->getWritePtr();
1389  mpfr_srcptr vs = pSrc->getReadPtr();
1390  mpfr_add(v, v, vs, GMP_RNDN);
1391  pDst += vDst.GetStep();
1392  pSrc += vSrc.GetStep();
1393  }
1394  }
1395 
1396  template<unsigned int Precision, class T2>
1398  {
1399  ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
1400  int i, cnt = vDst.GetLength();
1401  ampf<Precision> *pDst = vDst.GetData();
1402  const ampf<Precision> *pSrc = vSrc.GetData();
1403  ampf<Precision> a(alpha), tmp;
1404  for(i=0; i<cnt; i++)
1405  {
1406  mpfr_ptr v = pDst->getWritePtr();
1407  mpfr_srcptr vs = pSrc->getReadPtr();
1408  mpfr_mul(tmp.getWritePtr(), a.getReadPtr(), vs, GMP_RNDN);
1409  mpfr_add(v, v, tmp.getWritePtr(), GMP_RNDN);
1410  pDst += vDst.GetStep();
1411  pSrc += vSrc.GetStep();
1412  }
1413  }
1414 
1415  template<unsigned int Precision>
1417  {
1418  ap::ap_error::make_assertion(vDst.GetLength()==vSrc.GetLength());
1419  int i, cnt = vDst.GetLength();
1420  ampf<Precision> *pDst = vDst.GetData();
1421  const ampf<Precision> *pSrc = vSrc.GetData();
1422  for(i=0; i<cnt; i++)
1423  {
1424  mpfr_ptr v = pDst->getWritePtr();
1425  mpfr_srcptr vs = pSrc->getReadPtr();
1426  mpfr_sub(v, v, vs, GMP_RNDN);
1427  pDst += vDst.GetStep();
1428  pSrc += vSrc.GetStep();
1429  }
1430  }
1431 
1432  template<unsigned int Precision, class T2>
1434  {
1435  vAdd(vDst, vSrc, -alpha);
1436  }
1437 
1438  template<unsigned int Precision, class T2>
1440  {
1441  int i, cnt = vDst.GetLength();
1442  ampf<Precision> *pDst = vDst.GetData();
1444  for(i=0; i<cnt; i++)
1445  {
1446  mpfr_ptr v = pDst->getWritePtr();
1447  mpfr_mul(v, a.getReadPtr(), v, GMP_RNDN);
1448  pDst += vDst.GetStep();
1449  }
1450  }
1451 }
1452 
1453 #endif
amp::cosh
const ampf< Precision > cosh(const ampf< Precision > &x)
Definition: amp.h:1047
amp::campf::campf
campf(const campf< Prec2 > &z)
Definition: amp.h:1094
amp::ampf::isZero
bool isZero() const
Definition: amp.h:329
omalloc.h
exponent
int exponent(const CanonicalForm &f, int q)
int exponent ( const CanonicalForm & f, int q )
Definition: gengftables-conway.cc:92
amp::ampf::ampf
ampf(const char *s)
Definition: amp.h:118
amp::ampf::InitializeAsDouble
void InitializeAsDouble(long double v)
Definition: amp.h:276
f
FILE * f
Definition: checklibs.c:9
amp::operator*=
campf< Precision > & operator*=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition: amp.h:1165
amp::ampf::rval
mpfr_record * rval
Definition: amp.h:235
amp::ampf::toHex
std::string toHex() const
Definition: amp.h:355
amp::operator+
const ampf< Precision > operator+(const ampf< Precision > &op1)
Definition: amp.h:644
amp::vMove
void vMove(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition: amp.h:1328
amp::sign
const int sign(const ampf< Precision > &x)
Definition: amp.h:703
amp::ampf::ampf
ampf(signed long v)
Definition: amp.h:104
amp::ampf::getAlgoPascalEpsilon
static const ampf getAlgoPascalEpsilon()
Definition: amp.h:566
amp::vMul
void vMul(ap::raw_vector< ampf< Precision > > vDst, T2 alpha)
Definition: amp.h:1439
amp::abscomplex
const ampf< Precision > abscomplex(const campf< Precision > &z)
Definition: amp.h:1208
x
Variable x
Definition: cfModGcd.cc:4023
y
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:57
amp::trunc
const signed long trunc(const ampf< Precision > &x)
Definition: amp.h:750
amp::ampf::CheckPrecision
void CheckPrecision()
Definition: amp.h:245
result
return result
Definition: facAbsBiFact.cc:76
amp::frac
const ampf< Precision > frac(const ampf< Precision > &x)
Definition: amp.h:765
amp::mpfr_reference::initialize
void initialize(int Precision)
Definition: amp.cpp:115
amp::abs
const ampf< Precision > abs(const ampf< Precision > &x)
Definition: amp.h:714
amp::exp
const ampf< Precision > exp(const ampf< Precision > &x)
Definition: amp.h:1031
amp::sqrt
const ampf< Precision > sqrt(const ampf< Precision > &x)
Definition: amp.h:741
amp::sinh
const ampf< Precision > sinh(const ampf< Precision > &x)
Definition: amp.h:1039
amp::ampf::ampf
ampf(unsigned short v)
Definition: amp.h:109
amp::ampf::operator/=
ampf & operator/=(const T &v)
Definition: amp.h:187
amp::ampf::getRandom
static const ampf getRandom()
Definition: amp.h:594
amp::campf::campf
campf(signed char v)
Definition: amp.h:1087
amp::ampf::ampf
ampf(signed int v)
Definition: amp.h:106
amp::operator-
const ampf< Precision > operator-(const ampf< Precision > &op1)
Definition: amp.h:650
amp::ampf::ampf
ampf(const ampf &r)
Definition: amp.h:123
amp::divisionByZero
Definition: amp.h:24
amp::campf::operator=
campf & operator=(long double v)
Definition: amp.h:1097
amp::tanh
const ampf< Precision > tanh(const ampf< Precision > &x)
Definition: amp.h:1055
amp::vMoveNeg
void vMoveNeg(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition: amp.h:1345
amp::atan
const ampf< Precision > atan(const ampf< Precision > &x)
Definition: amp.h:991
amp::ampf::getReadPtr
mpfr_srcptr getReadPtr() const
Definition: amp.h:292
amp::ampf
Definition: amp.h:82
amp::mpfr_storage
Definition: amp.h:47
amp::campf::campf
campf(float v)
Definition: amp.h:1080
amp::incorrectPrecision
Definition: amp.h:22
ap::raw_vector
Definition: ap.h:172
amp::tan
const ampf< Precision > tan(const ampf< Precision > &x)
Definition: amp.h:967
amp::floor
const signed long floor(const ampf< Precision > &x)
Definition: amp.h:774
amp::conj
const campf< Precision > conj(const campf< Precision > &z)
Definition: amp.h:1226
amp::ampf::isFiniteNumber
bool isFiniteNumber() const
Definition: amp.h:315
amp::campf::y
ampf< Precision > y
Definition: amp.h:1126
amp::campf::campf
campf(unsigned int v)
Definition: amp.h:1084
amp::mpfr_storage::getRandState
static gmp_randstate_t * getRandState()
Definition: amp.cpp:37
amp::operator/=
campf< Precision > & operator/=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition: amp.h:1201
amp::ampf::operator*=
ampf & operator*=(const T &v)
Definition: amp.h:186
amp::ampf::operator+=
ampf & operator+=(const T &v)
Definition: amp.h:184
amp::ampf::operator=
ampf & operator=(long double v)
Definition: amp.h:141
w
const CanonicalForm & w
Definition: facAbsFact.cc:55
amp::ampf::getUlp
static const ampf getUlp()
Definition: amp.h:512
amp::ampf::ampf
ampf(long double v)
Definition: amp.h:101
amp::operator==
const bool operator==(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:605
amp::ampf::~ampf
~ampf()
Definition: amp.h:88
amp::log2
const ampf< Precision > log2(const ampf< Precision > &x)
Definition: amp.h:1015
amp::ampf::getAlgoPascalMaxNumber
static const ampf getAlgoPascalMaxNumber()
Definition: amp.h:572
amp::campf
Definition: amp.h:1074
amp::ampf::toDec
std::string toDec() const
Definition: amp.h:400
amp::campf::campf
campf(const ampf< Precision > &_x)
Definition: amp.h:1089
amp::ceil
const signed long ceil(const ampf< Precision > &x)
Definition: amp.h:789
amp::ampf::ampf
ampf(mpfr_record *v)
Definition: amp.h:99
amp::operator!=
const bool operator!=(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:611
amp::vAdd
void vAdd(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition: amp.h:1380
pi
#define pi
Definition: libparse.cc:1143
amp::acos
const ampf< Precision > acos(const ampf< Precision > &x)
Definition: amp.h:983
amp::campf::campf
campf(signed int v)
Definition: amp.h:1083
amp::mpfr_reference::mpfr_reference
mpfr_reference()
Definition: amp.cpp:82
amp::mpfr_record::refCount
unsigned int refCount
Definition: amp.h:36
amp::pow
const ampf< Precision > pow(const ampf< Precision > &x, const ampf< Precision > &y)
Definition: amp.h:1063
ap.h
i
int i
Definition: cfEzgcd.cc:125
amp::log
const ampf< Precision > log(const ampf< Precision > &x)
Definition: amp.h:1007
amp::halfpi
const ampf< Precision > halfpi()
Definition: amp.h:933
amp::campf::campf
campf(const ampf< Precision > &_x, const ampf< Precision > &_y)
Definition: amp.h:1090
res
CanonicalForm res
Definition: facAbsFact.cc:64
amp::exception
Definition: amp.h:21
amp::mpfr_record::next
mpfr_record * next
Definition: amp.h:39
amp::ampf::toString
char * toString() const
Definition: amp.h:446
amp::ampf::getWritePtr
mpfr_ptr getWritePtr()
Definition: amp.h:303
amp::ampf::ampf
ampf(unsigned int v)
Definition: amp.h:107
alpha
Variable alpha
Definition: facAbsBiFact.cc:52
amp::ampf::getMinNumber
static const ampf getMinNumber()
Definition: amp.h:558
T
static jList * T
Definition: janet.cc:31
amp::ldexp2
const ampf< Precision > ldexp2(const ampf< Precision > &x, mp_exp_t exponent)
Definition: amp.h:838
amp::invalidConversion
Definition: amp.h:26
amp::internalError
Definition: amp.h:28
amp::mpfr_reference::free
void free()
Definition: amp.cpp:123
amp::sqr
const ampf< Precision > sqr(const ampf< Precision > &x)
Definition: amp.h:694
amp::operator/
const ampf< Precision > operator/(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:683
T2
void T2(ideal h)
Definition: cohomo.cc:2754
amp::campf::campf
campf(const campf &z)
Definition: amp.h:1091
amp::operator<
const bool operator<(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:617
amp::sin
const ampf< Precision > sin(const ampf< Precision > &x)
Definition: amp.h:951
omAlloc
#define omAlloc(size)
Definition: omAllocDecl.h:210
amp::mpfr_reference::ref
mpfr_record * ref
Definition: amp.h:75
amp::mpfr_reference::getWritePtr
mpfr_ptr getWritePtr()
Definition: amp.cpp:142
amp::ampf::ampf
ampf(signed short v)
Definition: amp.h:108
amp::mpfr_record_ptr
mpfr_record * mpfr_record_ptr
Definition: amp.h:42
amp::ampf::ampf
ampf(const ampf< Precision2 > &r)
Definition: amp.h:130
amp::mpfr_reference::~mpfr_reference
~mpfr_reference()
Definition: amp.cpp:109
amp::invalidString
Definition: amp.h:27
amp::twopi
const ampf< Precision > twopi()
Definition: amp.h:942
amp::sqrtOfNegativeNumber
Definition: amp.h:25
amp::mpfr_reference::operator=
mpfr_reference & operator=(const mpfr_reference &r)
Definition: amp.cpp:94
amp::campf::campf
campf(long double v)
Definition: amp.h:1078
amp::operator+=
campf< Precision > & operator+=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition: amp.h:1145
amp::ampf::ampf
ampf(float v)
Definition: amp.h:103
amp::mpfr_reference
Definition: amp.h:61
amp::ampf::ampf
ampf(unsigned char v)
Definition: amp.h:111
ap::ap_error::make_assertion
static void make_assertion(bool bClause)
Definition: ap.h:61
amp::ampf::isNegativeNumber
bool isNegativeNumber() const
Definition: amp.h:335
amp::minimum
const ampf< Precision > minimum(const ampf< Precision > &x, const ampf< Precision > &y)
Definition: amp.h:732
amp::ampf::InitializeAsString
void InitializeAsString(const char *s)
Definition: amp.h:284
amp::ampf::getUlpOf
const ampf getUlpOf()
Definition: amp.h:343
amp::ampf::isPositiveNumber
bool isPositiveNumber() const
Definition: amp.h:321
amp::campf::campf
campf(signed short v)
Definition: amp.h:1085
amp::ampf::InitializeAsSLong
void InitializeAsSLong(signed long v)
Definition: amp.h:260
amp::atan2
const ampf< Precision > atan2(const ampf< Precision > &y, const ampf< Precision > &x)
Definition: amp.h:999
amp::ampf::ampf
ampf(const std::string &s)
Definition: amp.h:117
amp::__AMP_BINARY_OPI
__AMP_BINARY_OPI(char) __AMP_BINARY_OPI(short) __AMP_BINARY_OPI(long) __AMP_BINARY_OPI(int) __AMP_BINARY_OPF(float) __AMP_BINARY_OPF(double) __AMP_BINARY_OPF(long double) template< unsigned int Precision > const ampf< Precision > pi()
Definition: amp.h:890
amp::campf::campf
campf()
Definition: amp.h:1077
amp::ampf::ampf
ampf(double v)
Definition: amp.h:102
amp::overflow
Definition: amp.h:23
amp::ampf::operator-=
ampf & operator-=(const T &v)
Definition: amp.h:185
amp::maximum
const ampf< Precision > maximum(const ampf< Precision > &x, const ampf< Precision > &y)
Definition: amp.h:723
amp::mpfr_record::Precision
unsigned int Precision
Definition: amp.h:37
amp::ampf::getMaxNumber
static const ampf getMaxNumber()
Definition: amp.h:549
amp::operator<=
const bool operator<=(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:629
amp::frexp2
const ampf< Precision > frexp2(const ampf< Precision > &x, mp_exp_t *exponent)
Definition: amp.h:819
amp::log10
const ampf< Precision > log10(const ampf< Precision > &x)
Definition: amp.h:1023
amp::operator>=
const bool operator>=(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:635
amp::asin
const ampf< Precision > asin(const ampf< Precision > &x)
Definition: amp.h:975
NULL
#define NULL
Definition: omList.c:10
string
#define string
Definition: libparse.cc:1250
amp::mpfr_record
Definition: amp.h:34
amp::vSub
void vSub(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition: amp.h:1416
amp::ampf::getUlp512
static const ampf getUlp512()
Definition: amp.h:535
amp::campf::x
ampf< Precision > x
Definition: amp.h:1126
amp::ampf::ampf
ampf(unsigned long v)
Definition: amp.h:105
amp::campf::campf
campf(signed long v)
Definition: amp.h:1081
amp::unsigned32
unsigned long unsigned32
Definition: amp.h:31
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
amp::mpfr_storage::newMpfr
static mpfr_record * newMpfr(unsigned int Precision)
Definition: amp.cpp:11
__AMP_BINARY_OPF
#define __AMP_BINARY_OPF(type)
amp::ampf::ampf
ampf()
Definition: amp.h:98
amp::operator-=
campf< Precision > & operator-=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition: amp.h:1157
amp::campf::campf
campf(unsigned long v)
Definition: amp.h:1082
amp::ampf::ampf
ampf(signed char v)
Definition: amp.h:110
amp::domainError
Definition: amp.h:29
s
const CanonicalForm int s
Definition: facAbsFact.cc:55
amp::campf::campf
campf(unsigned short v)
Definition: amp.h:1086
amp::mpfr_reference::getReadPtr
mpfr_srcptr getReadPtr() const
Definition: amp.cpp:134
amp::ampf::getUlp256
static const ampf getUlp256()
Definition: amp.h:521
amp::mpfr_record::value
mpfr_t value
Definition: amp.h:38
amp::mpfr_storage::deleteMpfr
static void deleteMpfr(mpfr_record *ref)
Definition: amp.cpp:30
amp::cos
const ampf< Precision > cos(const ampf< Precision > &x)
Definition: amp.h:959
ap::const_raw_vector
Definition: ap.h:141
amp::operator*
const ampf< Precision > operator*(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:675
amp
Definition: amp.h:19
amp::ampf::getAlgoPascalMinNumber
static const ampf getAlgoPascalMinNumber()
Definition: amp.h:583
amp::signed32
signed long signed32
Definition: amp.h:32
amp::ampf::InitializeAsULong
void InitializeAsULong(unsigned long v)
Definition: amp.h:268
amp::campf::campf
campf(unsigned char v)
Definition: amp.h:1088
amp::mpfr_storage::getList
static mpfr_record_ptr & getList(unsigned int Precision)
Definition: amp.cpp:63
amp::ampf::toDouble
double toDouble() const
Definition: amp.h:349
amp::operator>
const bool operator>(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:623
amp::round
const signed long round(const ampf< Precision > &x)
Definition: amp.h:804
amp::campf::campf
campf(double v)
Definition: amp.h:1079
amp::ampf::InitializeAsZero
void InitializeAsZero()
Definition: amp.h:252
amp::csqr
const campf< Precision > csqr(const campf< Precision > &z)
Definition: amp.h:1232