"Fossies" - the Fresh Open Source Software Archive

Member "singular-4.2.1/Singular/svd/libs/amp.h" (9 Jun 2021, 58151 Bytes) of package /linux/misc/singular-4.2.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "amp.h" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 4.2.0p3_vs_4.2.1.

    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;
   38         mpfr_record *next;
   39     };
   40 
   41     typedef mpfr_record* mpfr_record_ptr;
   42 
   43     //
   44     // storage for mpfr_t instances
   45     //
   46     class mpfr_storage
   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     //
   60     class mpfr_reference
   61     {
   62     public:
   63         mpfr_reference();
   64         mpfr_reference(const mpfr_reference& r);
   65         mpfr_reference& operator= (const mpfr_reference &r);
   66         ~mpfr_reference();
   67 
   68         void initialize(int Precision);
   69         void free();
   70 
   71         mpfr_srcptr getReadPtr() const;
   72         mpfr_ptr getWritePtr();
   73     private:
   74         mpfr_record *ref;
   75     };
   76 
   77     //
   78     // ampf template
   79     //
   80     template<unsigned int Precision>
   81     class ampf
   82     {
   83     public:
   84         //
   85         // Destructor
   86         //
   87         ~ampf()
   88         {
   89             rval->refCount--;
   90             if( rval->refCount==0 )
   91                 mpfr_storage::deleteMpfr(rval);
   92         }
   93 
   94         //
   95         // Initializing
   96         //
   97         ampf ()                 { InitializeAsZero(); }
   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>
  129         ampf(const ampf<Precision2>& r)
  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 )
  162                 mpfr_storage::deleteMpfr(rval);
  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>
  170         ampf& operator= (const ampf<Precision2>& r)
  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;
  234         mpfr_record *rval;
  235     };
  236 
  237     /*void ampf<Precision>::CheckPrecision()
  238     {
  239         if( Precision<32 )
  240             throw incorrectPrecision();
  241     }***/
  242 
  243     template<unsigned int Precision>
  244     void ampf<Precision>::CheckPrecision()
  245     {
  246         if( Precision<32 )
  247             throw incorrectPrecision();
  248     }
  249 
  250     template<unsigned int Precision>
  251     void ampf<Precision>::InitializeAsZero()
  252     {
  253         CheckPrecision();
  254         rval = mpfr_storage::newMpfr(Precision);
  255         mpfr_set_ui(getWritePtr(), 0, GMP_RNDN);
  256     }
  257 
  258     template<unsigned int Precision>
  259     void ampf<Precision>::InitializeAsSLong(signed long sv)
  260     {
  261         CheckPrecision();
  262         rval = mpfr_storage::newMpfr(Precision);
  263         mpfr_set_si(getWritePtr(), sv, GMP_RNDN);
  264     }
  265 
  266     template<unsigned int Precision>
  267     void ampf<Precision>::InitializeAsULong(unsigned long v)
  268     {
  269         CheckPrecision();
  270         rval = mpfr_storage::newMpfr(Precision);
  271         mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
  272     }
  273 
  274     template<unsigned int Precision>
  275     void ampf<Precision>::InitializeAsDouble(long double v)
  276     {
  277         CheckPrecision();
  278         rval = mpfr_storage::newMpfr(Precision);
  279         mpfr_set_ld(getWritePtr(), v, GMP_RNDN);
  280     }
  281 
  282     template<unsigned int Precision>
  283     void ampf<Precision>::InitializeAsString(const char *s)
  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>
  302     mpfr_ptr ampf<Precision>::getWritePtr()
  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>
  314     bool ampf<Precision>::isFiniteNumber() const
  315     {
  316         return mpfr_number_p(getReadPtr())!=0;
  317     }
  318 
  319     template<unsigned int Precision>
  320     bool ampf<Precision>::isPositiveNumber() const
  321     {
  322         if( !isFiniteNumber() )
  323             return false;
  324         return mpfr_sgn(getReadPtr())>0;
  325     }
  326 
  327     template<unsigned int Precision>
  328     bool ampf<Precision>::isZero() const
  329     {
  330         return mpfr_zero_p(getReadPtr())!=0;
  331     }
  332 
  333     template<unsigned int Precision>
  334     bool ampf<Precision>::isNegativeNumber() const
  335     {
  336         if( !isFiniteNumber() )
  337             return false;
  338         return mpfr_sgn(getReadPtr())<0;
  339     }
  340 
  341     template<unsigned int Precision>
  342     const ampf<Precision> ampf<Precision>::getUlpOf()
  343     {
  344         return getUlpOf(*this);
  345     }
  346 
  347     template<unsigned int Precision>
  348     double ampf<Precision>::toDouble() const
  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>
  445     char * ampf<Precision>::toString() const
  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>
  488     const ampf<Precision> ampf<Precision>::getUlpOf(const ampf<Precision> &x)
  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>
  511     const ampf<Precision> ampf<Precision>::getUlp()
  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>
  520     const ampf<Precision> ampf<Precision>::getUlp256()
  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>
  534     const ampf<Precision> ampf<Precision>::getUlp512()
  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>
  548     const ampf<Precision> ampf<Precision>::getMaxNumber()
  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>
  557     const ampf<Precision> ampf<Precision>::getMinNumber()
  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>
  565     const ampf<Precision> ampf<Precision>::getAlgoPascalEpsilon()
  566     {
  567         return getUlp256();
  568     }
  569 
  570     template<unsigned int Precision>
  571     const ampf<Precision> ampf<Precision>::getAlgoPascalMaxNumber()
  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>
  582     const ampf<Precision> ampf<Precision>::getAlgoPascalMinNumber()
  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>
  593     const ampf<Precision> ampf<Precision>::getRandom()
  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>
  643     const ampf<Precision> operator+(const ampf<Precision>& op1)
  644     {
  645         return op1;
  646     }
  647 
  648     template<unsigned int Precision>
  649     const ampf<Precision> operator-(const ampf<Precision>& op1)
  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>
  657     const ampf<Precision> operator+(const ampf<Precision>& op1, const ampf<Precision>& op2)
  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>
  665     const ampf<Precision> operator-(const ampf<Precision>& op1, const ampf<Precision>& op2)
  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>
  674     const ampf<Precision> operator*(const ampf<Precision>& op1, const ampf<Precision>& op2)
  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>
  682     const ampf<Precision> operator/(const ampf<Precision>& op1, const ampf<Precision>& op2)
  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>
  693     const ampf<Precision> sqr(const ampf<Precision> &x)
  694     {
  695         // TODO: optimize temporary for return value
  696         ampf<Precision> res;
  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>
  713     const ampf<Precision> abs(const ampf<Precision> &x)
  714     {
  715         // TODO: optimize temporary for return value
  716         ampf<Precision> res;
  717         mpfr_abs(res.getWritePtr(), x.getReadPtr(), GMP_RNDN);
  718         return res;
  719     }
  720 
  721     template<unsigned int Precision>
  722     const ampf<Precision> maximum(const ampf<Precision> &x, const ampf<Precision> &y)
  723     {
  724         // TODO: optimize temporary for return value
  725         ampf<Precision> res;
  726         mpfr_max(res.getWritePtr(), x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
  727         return res;
  728     }
  729 
  730     template<unsigned int Precision>
  731     const ampf<Precision> minimum(const ampf<Precision> &x, const ampf<Precision> &y)
  732     {
  733         // TODO: optimize temporary for return value
  734         ampf<Precision> res;
  735         mpfr_min(res.getWritePtr(), x.getReadPtr(), y.getReadPtr(), GMP_RNDN);
  736         return res;
  737     }
  738 
  739     template<unsigned int Precision>
  740     const ampf<Precision> sqrt(const ampf<Precision> &x)
  741     {
  742         // TODO: optimize temporary for return value
  743         ampf<Precision> res;
  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>
  764     const ampf<Precision> frac(const ampf<Precision> &x)
  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>
  818     const ampf<Precision> frexp2(const ampf<Precision> &x, mp_exp_t *exponent)
  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>
  837     const ampf<Precision> ldexp2(const ampf<Precision> &x, mp_exp_t exponent)
  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); }
  889     __AMP_BINARY_OPI(char)
  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>
  932     const ampf<Precision> halfpi()
  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>
  941     const ampf<Precision> twopi()
  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>
  950     const ampf<Precision> sin(const ampf<Precision> &x)
  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>
  958     const ampf<Precision> cos(const ampf<Precision> &x)
  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>
  966     const ampf<Precision> tan(const ampf<Precision> &x)
  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>
  974     const ampf<Precision> asin(const ampf<Precision> &x)
  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>
  982     const ampf<Precision> acos(const ampf<Precision> &x)
  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>
  990     const ampf<Precision> atan(const ampf<Precision> &x)
  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>
  998     const ampf<Precision> atan2(const ampf<Precision> &y, const ampf<Precision> &x)
  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>
 1006     const ampf<Precision> log(const ampf<Precision> &x)
 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>
 1014     const ampf<Precision> log2(const ampf<Precision> &x)
 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>
 1022     const ampf<Precision> log10(const ampf<Precision> &x)
 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>
 1030     const ampf<Precision> exp(const ampf<Precision> &x)
 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>
 1038     const ampf<Precision> sinh(const ampf<Precision> &x)
 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>
 1046     const ampf<Precision> cosh(const ampf<Precision> &x)
 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>
 1054     const ampf<Precision> tanh(const ampf<Precision> &x)
 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>
 1062     const ampf<Precision> pow(const ampf<Precision> &x, const ampf<Precision> &y)
 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; }
 1109         campf& operator= (const campf& r)
 1110         {
 1111             x = r.x;
 1112             y = r.y;
 1113             return *this;
 1114         }
 1115 #ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
 1116         template<unsigned int Precision2>
 1117         campf& operator= (const campf<Precision2>& r)
 1118         {
 1119             x = r.x;
 1120             y = r.y;
 1121             return *this;
 1122         }
 1123 #endif
 1124 
 1125         ampf<Precision> x, y;
 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>
 1140     const campf<Precision> operator+(const campf<Precision>& lhs)
 1141     { return lhs; }
 1142 
 1143     template<unsigned int Precision>
 1144     campf<Precision>& operator+=(campf<Precision>& lhs, const campf<Precision>& rhs)
 1145     { lhs.x += rhs.x; lhs.y += rhs.y; return lhs; }
 1146 
 1147     template<unsigned int Precision>
 1148     const campf<Precision> operator+(const campf<Precision>& lhs, const campf<Precision>& rhs)
 1149     { campf<Precision> r = lhs; r += rhs; return r; }
 1150 
 1151     template<unsigned int Precision>
 1152     const campf<Precision> operator-(const campf<Precision>& lhs)
 1153     { return campf<Precision>(-lhs.x, -lhs.y); }
 1154 
 1155     template<unsigned int Precision>
 1156     campf<Precision>& operator-=(campf<Precision>& lhs, const campf<Precision>& rhs)
 1157     { lhs.x -= rhs.x; lhs.y -= rhs.y; return lhs; }
 1158 
 1159     template<unsigned int Precision>
 1160     const campf<Precision> operator-(const campf<Precision>& lhs, const campf<Precision>& rhs)
 1161     { campf<Precision> r = lhs; r -= rhs; return r; }
 1162 
 1163     template<unsigned int Precision>
 1164     campf<Precision>& operator*=(campf<Precision>& lhs, const campf<Precision>& rhs)
 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>
 1173     const campf<Precision> operator*(const campf<Precision>& lhs, const campf<Precision>& rhs)
 1174     { campf<Precision> r = lhs; r *= rhs; return r; }
 1175 
 1176     template<unsigned int Precision>
 1177     const campf<Precision> operator/(const campf<Precision>& lhs, const campf<Precision>& rhs)
 1178     {
 1179         campf<Precision> result;
 1180         ampf<Precision> e;
 1181         ampf<Precision> f;
 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>
 1200     campf<Precision>& operator/=(campf<Precision>& lhs, const campf<Precision>& rhs)
 1201     {
 1202         lhs = lhs/rhs;
 1203         return lhs;
 1204     }
 1205 
 1206     template<unsigned int Precision>
 1207     const ampf<Precision> abscomplex(const campf<Precision> &z)
 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>
 1225     const campf<Precision> conj(const campf<Precision> &z)
 1226     {
 1227         return campf<Precision>(z.x, -z.y);
 1228     }
 1229 
 1230     template<unsigned int Precision>
 1231     const campf<Precision> csqr(const campf<Precision> &z)
 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             }
 1313             mpfr_storage::deleteMpfr(t);
 1314             return r;
 1315         }
 1316         catch(...)
 1317         {
 1318             if( r!=NULL )
 1319                 mpfr_storage::deleteMpfr(r);
 1320             if( t!=NULL )
 1321                 mpfr_storage::deleteMpfr(t);
 1322             throw;
 1323         }
 1324     }
 1325 
 1326     template<unsigned int Precision>
 1327     void vMove(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
 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>
 1344     void vMoveNeg(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
 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>
 1361     void vMove(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
 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();
 1367         ampf<Precision>       a(alpha);
 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>
 1379     void vAdd(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
 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>
 1396     void vAdd(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
 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>
 1415     void vSub(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
 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>
 1432     void vSub(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
 1433     {
 1434         vAdd(vDst, vSrc, -alpha);
 1435     }
 1436 
 1437     template<unsigned int Precision, class T2>
 1438     void vMul(ap::raw_vector< ampf<Precision> > vDst, T2 alpha)
 1439     {
 1440         int i, cnt = vDst.GetLength();
 1441         ampf<Precision>       *pDst = vDst.GetData();
 1442         ampf<Precision>       a(alpha);
 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