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