My Project  UNKNOWN_GIT_VERSION
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, è ïðè
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
const ampf< Precision > halfpi()
Definition: amp.h:933
ampf(unsigned short v)
Definition: amp.h:109
campf(const campf &z)
Definition: amp.h:1091
ampf(signed char v)
Definition: amp.h:110
mpfr_record * rval
Definition: amp.h:235
void InitializeAsZero()
Definition: amp.h:252
static const ampf getUlp256()
Definition: amp.h:521
const CanonicalForm int s
Definition: facAbsFact.cc:55
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:57
const ampf< Precision > operator/(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:683
void CheckPrecision()
Definition: amp.h:245
mpfr_record * next
Definition: amp.h:39
static const ampf getAlgoPascalMinNumber()
Definition: amp.h:583
void InitializeAsULong(unsigned long v)
Definition: amp.h:268
campf(unsigned long v)
Definition: amp.h:1082
void InitializeAsDouble(long double v)
Definition: amp.h:276
static const ampf getRandom()
Definition: amp.h:594
ampf(signed int v)
Definition: amp.h:106
ampf(const ampf &r)
Definition: amp.h:123
mpfr_srcptr getReadPtr() const
Definition: amp.h:292
ampf & operator*=(const T &v)
Definition: amp.h:186
const bool operator==(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:605
static const ampf getAlgoPascalEpsilon()
Definition: amp.h:566
unsigned int Precision
Definition: amp.h:37
ampf(signed long v)
Definition: amp.h:104
ampf(long double v)
Definition: amp.h:101
const ampf< Precision > log2(const ampf< Precision > &x)
Definition: amp.h:1015
const ampf< Precision > acos(const ampf< Precision > &x)
Definition: amp.h:983
mpfr_ptr getWritePtr()
Definition: amp.cpp:142
void vAdd(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition: amp.h:1380
campf(const ampf< Precision > &_x)
Definition: amp.h:1089
signed long signed32
Definition: amp.h:32
ampf(mpfr_record *v)
Definition: amp.h:99
ampf(const char *s)
Definition: amp.h:118
#define string
Definition: libparse.cc:1250
ampf & operator/=(const T &v)
Definition: amp.h:187
campf()
Definition: amp.h:1077
mpfr_srcptr getReadPtr() const
Definition: amp.cpp:134
campf(const ampf< Precision > &_x, const ampf< Precision > &_y)
Definition: amp.h:1090
const ampf< Precision > atan2(const ampf< Precision > &y, const ampf< Precision > &x)
Definition: amp.h:999
static mpfr_record * newMpfr(unsigned int Precision)
Definition: amp.cpp:11
const ampf< Precision > minimum(const ampf< Precision > &x, const ampf< Precision > &y)
Definition: amp.h:732
campf(signed int v)
Definition: amp.h:1083
const ampf< Precision > cos(const ampf< Precision > &x)
Definition: amp.h:959
static const ampf getAlgoPascalMaxNumber()
Definition: amp.h:572
const ampf< Precision > exp(const ampf< Precision > &x)
Definition: amp.h:1031
void initialize(int Precision)
Definition: amp.cpp:115
Variable alpha
Definition: facAbsBiFact.cc:52
const signed long trunc(const ampf< Precision > &x)
Definition: amp.h:750
static void make_assertion(bool bClause)
Definition: ap.h:49
campf< Precision > & operator/=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition: amp.h:1201
__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
const ampf< Precision > operator*(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:675
#define omAlloc(size)
Definition: omAllocDecl.h:210
ampf & operator=(long double v)
Definition: amp.h:141
ampf & operator+=(const T &v)
Definition: amp.h:184
const bool operator>=(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:635
static const ampf getUlp()
Definition: amp.h:512
mpfr_record * mpfr_record_ptr
Definition: amp.h:42
const ampf< Precision > log(const ampf< Precision > &x)
Definition: amp.h:1007
const signed long round(const ampf< Precision > &x)
Definition: amp.h:804
const signed long floor(const ampf< Precision > &x)
Definition: amp.h:774
ampf(signed short v)
Definition: amp.h:108
const ampf< Precision > sinh(const ampf< Precision > &x)
Definition: amp.h:1039
static void deleteMpfr(mpfr_record *ref)
Definition: amp.cpp:30
bool isFiniteNumber() const
Definition: amp.h:315
mpfr_reference & operator=(const mpfr_reference &r)
Definition: amp.cpp:94
campf(long double v)
Definition: amp.h:1078
campf(signed long v)
Definition: amp.h:1081
ampf(float v)
Definition: amp.h:103
CanonicalForm res
Definition: facAbsFact.cc:64
bool isNegativeNumber() const
Definition: amp.h:335
const ampf< Precision > operator+(const ampf< Precision > &op1)
Definition: amp.h:644
~ampf()
Definition: amp.h:88
const bool operator!=(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:611
campf(float v)
Definition: amp.h:1080
ampf(unsigned int v)
Definition: amp.h:107
void vMul(ap::raw_vector< ampf< Precision > > vDst, T2 alpha)
Definition: amp.h:1439
const ampf< Precision > abscomplex(const campf< Precision > &z)
Definition: amp.h:1208
campf(unsigned char v)
Definition: amp.h:1088
std::string toDec() const
Definition: amp.h:400
const signed long ceil(const ampf< Precision > &x)
Definition: amp.h:789
static const ampf getMinNumber()
Definition: amp.h:558
const ampf getUlpOf()
Definition: amp.h:343
bool isPositiveNumber() const
Definition: amp.h:321
ampf< Precision > y
Definition: amp.h:1126
campf(double v)
Definition: amp.h:1079
const ampf< Precision > sqrt(const ampf< Precision > &x)
Definition: amp.h:741
campf(signed char v)
Definition: amp.h:1087
const ampf< Precision > pow(const ampf< Precision > &x, const ampf< Precision > &y)
Definition: amp.h:1063
const campf< Precision > csqr(const campf< Precision > &z)
Definition: amp.h:1232
ampf(double v)
Definition: amp.h:102
const bool operator>(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:623
void T2(ideal h)
Definition: cohomo.cc:2754
campf< Precision > & operator+=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition: amp.h:1145
mpfr_ptr getWritePtr()
Definition: amp.h:303
char * toString() const
Definition: amp.h:446
FILE * f
Definition: checklibs.c:9
void vMoveNeg(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition: amp.h:1345
int i
Definition: cfEzgcd.cc:125
void vSub(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition: amp.h:1416
void InitializeAsString(const char *s)
Definition: amp.h:284
campf< Precision > & operator-=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition: amp.h:1157
ampf(const ampf< Precision2 > &r)
Definition: amp.h:130
campf(const campf< Prec2 > &z)
Definition: amp.h:1094
ampf(const std::string &s)
Definition: amp.h:117
campf & operator=(long double v)
Definition: amp.h:1097
mpfr_t value
Definition: amp.h:38
static mpfr_record_ptr & getList(unsigned int Precision)
Definition: amp.cpp:63
static const ampf getMaxNumber()
Definition: amp.h:549
const bool operator<(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:617
campf< Precision > & operator*=(campf< Precision > &lhs, const campf< Precision > &rhs)
Definition: amp.h:1165
const ampf< Precision > ldexp2(const ampf< Precision > &x, mp_exp_t exponent)
Definition: amp.h:838
const ampf< Precision > tan(const ampf< Precision > &x)
Definition: amp.h:967
#define pi
Definition: libparse.cc:1143
const campf< Precision > conj(const campf< Precision > &z)
Definition: amp.h:1226
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
ampf()
Definition: amp.h:98
#define __AMP_BINARY_OPF(type)
campf(unsigned int v)
Definition: amp.h:1084
const ampf< Precision > twopi()
Definition: amp.h:942
Definition: amp.h:19
#define NULL
Definition: omList.c:10
unsigned int refCount
Definition: amp.h:36
const ampf< Precision > frac(const ampf< Precision > &x)
Definition: amp.h:765
ampf & operator-=(const T &v)
Definition: amp.h:185
const ampf< Precision > operator-(const ampf< Precision > &op1)
Definition: amp.h:650
campf(unsigned short v)
Definition: amp.h:1086
ampf(unsigned char v)
Definition: amp.h:111
const ampf< Precision > asin(const ampf< Precision > &x)
Definition: amp.h:975
const ampf< Precision > tanh(const ampf< Precision > &x)
Definition: amp.h:1055
const CanonicalForm & w
Definition: facAbsFact.cc:55
Variable x
Definition: cfModGcd.cc:4023
unsigned long unsigned32
Definition: amp.h:31
void InitializeAsSLong(signed long v)
Definition: amp.h:260
campf(signed short v)
Definition: amp.h:1085
void vMove(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
Definition: amp.h:1328
const int sign(const ampf< Precision > &x)
Definition: amp.h:703
bool isZero() const
Definition: amp.h:329
int exponent(const CanonicalForm &f, int q)
int exponent ( const CanonicalForm & f, int q )
static gmp_randstate_t * getRandState()
Definition: amp.cpp:37
const ampf< Precision > sqr(const ampf< Precision > &x)
Definition: amp.h:694
const ampf< Precision > frexp2(const ampf< Precision > &x, mp_exp_t *exponent)
Definition: amp.h:819
mpfr_record * ref
Definition: amp.h:75
const ampf< Precision > maximum(const ampf< Precision > &x, const ampf< Precision > &y)
Definition: amp.h:723
Definition: amp.h:82
static jList * T
Definition: janet.cc:31
const bool operator<=(const ampf< Precision > &op1, const ampf< Precision > &op2)
Definition: amp.h:629
double toDouble() const
Definition: amp.h:349
std::string toHex() const
Definition: amp.h:355
const ampf< Precision > sin(const ampf< Precision > &x)
Definition: amp.h:951
const ampf< Precision > log10(const ampf< Precision > &x)
Definition: amp.h:1023
const ampf< Precision > atan(const ampf< Precision > &x)
Definition: amp.h:991
static const ampf getUlp512()
Definition: amp.h:535
return result
Definition: facAbsBiFact.cc:76
const ampf< Precision > abs(const ampf< Precision > &x)
Definition: amp.h:714
ampf(unsigned long v)
Definition: amp.h:105
ampf< Precision > x
Definition: amp.h:1126
const ampf< Precision > cosh(const ampf< Precision > &x)
Definition: amp.h:1047