My Project
Loading...
Searching...
No Matches
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
18namespace amp
19{
20 class exception {};
21 class incorrectPrecision : public exception {};
22 class overflow : public exception {};
23 class divisionByZero : 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
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 //
99
100 ampf (long double v) { InitializeAsDouble(v); }
101 ampf (double 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 {
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; }
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>
354 std::string ampf<Precision>::toHex() const
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>
399 std::string ampf<Precision>::toDec() const
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 {
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
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
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
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)
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 {
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)
#define __AMP_BINARY_OPI(type)
Definition: amp.h:1240
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4082
FILE * f
Definition: checklibs.c:9
Definition: amp.h:82
ampf & operator/=(const T &v)
Definition: amp.h:186
ampf(unsigned int v)
Definition: amp.h:106
ampf & operator-=(const T &v)
Definition: amp.h:184
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
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
void CheckPrecision()
Definition: amp.h:244
ampf & operator+=(const T &v)
Definition: amp.h:183
static const ampf getAlgoPascalMaxNumber()
Definition: amp.h:571
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
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
ampf & operator=(long double v)
Definition: amp.h:140
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
ampf & operator*=(const T &v)
Definition: amp.h:185
static const ampf getMinNumber()
Definition: amp.h:557
bool isFiniteNumber() const
Definition: amp.h:314
ampf(const ampf< Precision2 > &r)
Definition: amp.h:129
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 & operator=(long double v)
Definition: amp.h:1096
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
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
#define exponent
STATIC_VAR jList * T
Definition: janet.cc:30
#define pi
Definition: libparse.cc:1145
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
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