[BACK]Return to bf.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/bf.c between version 1.7 and 1.14

version 1.7, 2015/08/04 06:20:45 version 1.14, 2017/09/04 01:57:53
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.6 2009/07/15 02:19:47 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.13 2015/09/14 05:47:09 noro Exp $
  */   */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
 #include <math.h>  #include <math.h>
   
   extern int mpfr_roundmode;
   
 Num tobf(Num,int);  Num tobf(Num,int);
   
 #define BFPREC(a) (((BF)(a))->body->_mpfr_prec)  #define BFPREC(a) (((BF)(a))->body->_mpfr_prec)
Line 14  void strtobf(char *s,BF *p)
Line 16  void strtobf(char *s,BF *p)
   BF r;    BF r;
   NEWBF(r);    NEWBF(r);
   mpfr_init(r->body);    mpfr_init(r->body);
   mpfr_set_str(r->body,s,10,MPFR_RNDN);    mpfr_set_str(r->body,s,10,mpfr_roundmode);
   *p = r;    *p = r;
 }  }
   
 double mpfrtodbl(mpfr_t a)  double mpfrtodbl(mpfr_t a)
 {  {
   return mpfr_get_d(a,MPFR_RNDN);    return mpfr_get_d(a,mpfr_roundmode);
 }  }
   
 Num tobf(Num a,int prec)  Num tobf(Num a,int prec)
Line 30  Num tobf(Num a,int prec)
Line 32  Num tobf(Num a,int prec)
   mpq_t q;    mpq_t q;
   BF d;    BF d;
   N nm,dn;    N nm,dn;
     C c;
     Num re,im;
   int sgn;    int sgn;
   
   if ( !a ) {    if ( !a ) {
         prec ? mpfr_init2(r,prec) : mpfr_init(r);      prec ? mpfr_init2(r,prec) : mpfr_init(r);
     mpfr_set_zero(r,1);      mpfr_set_zero(r,1);
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     return (Num)d;      return (Num)d;
Line 42  Num tobf(Num a,int prec)
Line 46  Num tobf(Num a,int prec)
     case N_B:      case N_B:
       return a;        return a;
       break;        break;
         case N_R:      case N_R:
           prec ? mpfr_init2(r,prec) : mpfr_init(r);        prec ? mpfr_init2(r,prec) : mpfr_init(r);
           mpfr_init_set_d(r,((Real)a)->body,MPFR_RNDN);        mpfr_init_set_d(r,((Real)a)->body,mpfr_roundmode);
       MPFRTOBF(r,d);        MPFRTOBF(r,d);
       return (Num)d;        return (Num)d;
           break;        break;
     case N_Q:      case N_Q:
       nm = NM((Q)a); dn = DN((Q)a); sgn = SGN((Q)a);        nm = NM((Q)a); dn = DN((Q)a); sgn = SGN((Q)a);
       if ( INT((Q)a) ) {        if ( INT((Q)a) ) {
         mpz_init(z);          mpz_init(z);
         mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));          mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
         if ( sgn < 0 ) mpz_neg(z,z);          if ( sgn < 0 ) mpz_neg(z,z);
         mpfr_init_set_z(r,z,MPFR_RNDN);          mpfr_init_set_z(r,z,mpfr_roundmode);
       } else {        } else {
         mpq_init(q);          mpq_init(q);
         mpz_import(mpq_numref(q),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));          mpz_import(mpq_numref(q),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
         mpz_import(mpq_denref(q),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn));          mpz_import(mpq_denref(q),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn));
         if ( sgn < 0 ) mpq_neg(q,q);          if ( sgn < 0 ) mpq_neg(q,q);
         mpfr_init_set_q(r,q,MPFR_RNDN);          mpfr_init_set_q(r,q,mpfr_roundmode);
       }        }
       MPFRTOBF(r,d);        MPFRTOBF(r,d);
       return (Num)d;        return (Num)d;
       break;        break;
       case N_C:
         re = tobf(((C)a)->r,prec); im = tobf(((C)a)->i,prec);
         NEWC(c); c->r = re; c->i = im;
         return (Num)c;
         break;
     default:      default:
       error("tobf : invalid argument");        error("tobf : invalid argument");
           break;        break;
     }      }
   }    }
 }  }
Line 76  void addbf(Num a,Num b,Num *c)
Line 85  void addbf(Num a,Num b,Num *c)
 {  {
   mpfr_t r;    mpfr_t r;
   BF d;    BF d;
     GZ z;
     GQ q;
   
   if ( !a )    if ( !a )
     *c = b;      *c = b;
   else if ( !b )    else if ( !b )
     *c = a;      *c = a;
   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
     (*addnumt[MIN(NID(a),NID(b))])(a,b,c);      (*addnumt[MAX(NID(a),NID(b))])(a,b,c);
   else {    else if ( NID(a) == N_B ) {
     if ( NID(a) != N_B ) a = tobf(a,0);      switch ( NID(b) ) {
     if ( NID(b) != N_B ) b = tobf(b,0);      case N_Q:
     mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        mpfr_init2(r,BFPREC(a));
     mpfr_add(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           mpfr_add_z(r,((BF)a)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)b);
           mpfr_add_q(r,((BF)a)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(a),53));
         mpfr_add_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
         break;
       case N_B:
         mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
         mpfr_add(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       default:
         goto err;
         break;
       }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    } else if ( NID(b) == N_B ) {
       switch ( NID(a) ) {
       case N_Q:
         mpfr_init2(r,BFPREC(b));
         if ( INT((Q)a) ) {
           z = ztogz((Q)a);
           mpfr_add_z(r,((BF)b)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)a);
           mpfr_add_q(r,((BF)b)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(b),53));
         mpfr_add_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode);
         break;
       default:
         goto err;
         break;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     } else
       goto err;
     if ( !cmpbf(*c,0) ) *c = 0;
     return;
   
   err: error("addbf : invalid argument");
 }  }
   
 void subbf(Num a,Num b,Num *c)  void subbf(Num a,Num b,Num *c)
 {  {
   mpfr_t r;    mpfr_t r,s;
     GZ z;
     GQ q;
   BF d;    BF d;
   
   if ( !a )    if ( !a )
     (*chsgnnumt[NID(b)])(b,c);      (*chsgnnumt[NID(b)])(b,c);
   else if ( !b )    else if ( !b )
     *c = a;      *c = a;
   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
     (*subnumt[MIN(NID(a),NID(b))])(a,b,c);      (*subnumt[MAX(NID(a),NID(b))])(a,b,c);
   else {    else if ( NID(a) == N_B ) {
     if ( NID(a) != N_B ) a = tobf(a,0);      switch ( NID(b) ) {
     if ( NID(b) != N_B ) b = tobf(b,0);      case N_Q:
     mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        mpfr_init2(r,BFPREC(a));
     mpfr_sub(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           mpfr_sub_z(r,((BF)a)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)b);
           mpfr_sub_q(r,((BF)a)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(a),53));
         mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
         break;
       case N_B:
         mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
         mpfr_sub(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       default:
         goto err;
       }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    } else if ( NID(b)==N_B ) {
       switch ( NID(a) ) {
       case N_Q:
         mpfr_init2(r,BFPREC(b));
         if ( INT((Q)a) ) {
           z = ztogz((Q)a);
           mpfr_sub_z(r,((BF)b)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)a);
           mpfr_sub_q(r,((BF)b)->body,q->body,mpfr_roundmode);
         }
         mpfr_neg(r,r,mpfr_roundmode);
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(b),53));
         mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       default:
         goto err;
       }
   
       MPFRTOBF(r,d);
       *c = (Num)d;
     } else
       goto err;
     if ( !cmpbf(*c,0) ) *c = 0;
     return;
   
   err: error("subbf : invalid argument");
 }  }
   
 void mulbf(Num a,Num b,Num *c)  void mulbf(Num a,Num b,Num *c)
 {  {
   mpfr_t r;    mpfr_t r;
     GZ z;
     GQ q;
   BF d;    BF d;
   int prec;    int prec;
   
   if ( !a || !b )    if ( !a || !b )
     *c = 0;      *c = 0;
   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
     (*mulnumt[MIN(NID(a),NID(b))])(a,b,c);      (*mulnumt[MAX(NID(a),NID(b))])(a,b,c);
   else {    else if ( NID(a) == N_B ) {
     if ( NID(a) != N_B ) a = tobf(a,0);      switch ( NID(b) ) {
     if ( NID(b) != N_B ) b = tobf(b,0);      case N_Q:
     mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        mpfr_init2(r,BFPREC(a));
     mpfr_mul(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           mpfr_mul_z(r,((BF)a)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)b);
           mpfr_mul_q(r,((BF)a)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(a),53));
         mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
         break;
       case N_B:
         mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
         mpfr_mul(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       default:
         goto err;
       }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    } else if ( NID(b) == N_B ) {
       switch ( NID(a) ) {
       case N_Q:
         mpfr_init2(r,BFPREC(b));
         if ( INT((Q)a) ) {
           z = ztogz((Q)a);
           mpfr_mul_z(r,((BF)b)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)a);
           mpfr_mul_q(r,((BF)b)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(b),53));
         mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode);
         break;
       default:
         goto err;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     } else
       goto err;
   
     if ( !cmpbf(*c,0) ) *c = 0;
     return;
   
   err: error("mulbf : invalid argument");
 }  }
   
 void divbf(Num a,Num b,Num *c)  void divbf(Num a,Num b,Num *c)
 {  {
   mpfr_t r;    mpfr_t s,r;
     GZ z;
     GQ q;
   BF d;    BF d;
   
   if ( !b )    if ( !b )
     error("divbf : division by 0");      error("divbf : division by 0");
   else if ( !a )    else if ( !a )
     *c = 0;      *c = 0;
   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
     (*divnumt[MIN(NID(a),NID(b))])(a,b,c);      (*divnumt[MAX(NID(a),NID(b))])(a,b,c);
   else {    else if ( NID(a) == N_B ) {
     if ( NID(a) != N_B ) a = tobf(a,0);      switch ( NID(b) ) {
     if ( NID(b) != N_B ) b = tobf(b,0);      case N_Q:
     mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        mpfr_init2(r,BFPREC(a));
     mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           mpfr_div_z(r,((BF)a)->body,z->body,mpfr_roundmode);
         } else {
           q = qtogq((Q)b);
           mpfr_div_q(r,((BF)a)->body,q->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(a),53));
         mpfr_div_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
         break;
       case N_B:
         mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
         mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       default:
         goto err;
       }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    } else if ( NID(b)==N_B ) {
       switch ( NID(a) ) {
       case N_Q:
         /* XXX : mpfr_z_div and mpfr_q_div are not implemented */
         a = tobf(a,BFPREC(b));
         mpfr_init2(r,BFPREC(b));
         mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       case N_R:
         /* double precision = 53 */
         mpfr_init2(r,MAX(BFPREC(b),53));
         mpfr_d_div(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       default:
         goto err;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     } else
       goto err;
   
     if ( !cmpbf(*c,0) ) *c = 0;
     return;
   
   err: error("mulbf : invalid argument");
 }  }
   
 void pwrbf(Num a,Num b,Num *c)  void pwrbf(Num a,Num b,Num *c)
 {  {
     int prec;
   mpfr_t r;    mpfr_t r;
     GZ z;
   BF d;    BF d;
   
   if ( !b )    if ( !b )
     *c = (Num)ONE;      *c = (Num)ONE;
   else if ( !a )    else if ( !a )
     *c = 0;      *c = 0;
   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_R) && (NID(b) <= N_R ) )
     (*pwrnumt[MIN(NID(a),NID(b))])(a,b,c);      (*pwrnumt[MAX(NID(a),NID(b))])(a,b,c);
   else {    else if ( NID(a) == N_B ) {
     if ( NID(a) != N_B ) a = tobf(a,0);      switch ( NID(b) ) {
     if ( NID(b) != N_B ) b = tobf(b,0);      case N_Q:
     mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        mpfr_init2(r,BFPREC(a));
     mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           mpfr_pow_z(r,((BF)a)->body,z->body,mpfr_roundmode);
         } else {
           b = tobf(b,BFPREC(a));
           mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         }
         break;
       case N_R:
         /* double precision = 53 */
         prec = MAX(BFPREC(a),53);
         mpfr_init2(r,prec);
         b = tobf(b,prec);
         mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       case N_B:
         mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
         mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       default:
         goto err;
       }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    } else if ( NID(b)==N_B ) {
       switch ( NID(a) ) {
       case N_Q:
         mpfr_init2(r,BFPREC(b));
         a = tobf(a,BFPREC(b));
         mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       case N_R:
         /* double precision = 53 */
         prec = MAX(BFPREC(a),53);
         mpfr_init2(r,prec);
         a = tobf(a,prec);
         mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
         break;
       default:
         goto err;
       }
       MPFRTOBF(r,d);
       *c = (Num)d;
     } else
       goto err;
   
     if ( !cmpbf(*c,0) ) *c = 0;
     return;
   
   err: error("pwrbf : invalid argument");
 }  }
   
 void chsgnbf(Num a,Num *c)  void chsgnbf(Num a,Num *c)
Line 183  void chsgnbf(Num a,Num *c)
Line 434  void chsgnbf(Num a,Num *c)
   
   if ( !a )    if ( !a )
     *c = 0;      *c = 0;
   else if ( NID(a) <= N_A )    else if ( NID(a) <= N_R )
     (*chsgnnumt[NID(a)])(a,c);      (*chsgnnumt[NID(a)])(a,c);
   else {    else if ( NID(a) == N_B ) {
     mpfr_init2(r,BFPREC(a));      mpfr_init2(r,BFPREC(a));
     mpfr_neg(r,((BF)a)->body,MPFR_RNDN);      mpfr_neg(r,((BF)a)->body,mpfr_roundmode);
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    } else
       error("chsgnbf : invalid argument");
 }  }
   
 int cmpbf(Num a,Num b)  int cmpbf(Num a,Num b)
 {  {
     int ret;
     GZ z;
     GQ q;
   
   if ( !a ) {    if ( !a ) {
     if ( !b || (NID(b)<=N_A) )      if ( !b ) return 0;
       else if ( NID(b)<=N_R )
       return (*cmpnumt[NID(b)])(a,b);        return (*cmpnumt[NID(b)])(a,b);
       else if ( NID(b)==N_B )
         return -mpfr_sgn(((BF)b)->body);
     else      else
       return -mpfr_sgn(((BF)a)->body);        goto err;
   } else if ( !b ) {    } else if ( !b ) {
     if ( !a || (NID(a)<=N_A) )      if ( NID(a)<=N_R )
       return (*cmpnumt[NID(a)])(a,b);        return (*cmpnumt[NID(a)])(a,b);
     else      else if ( NID(a)==N_B )
       return mpfr_sgn(((BF)a)->body);        return mpfr_sgn(((BF)a)->body);
   } else {      else
     if ( NID(a) != N_B ) a = tobf(a,0);        goto err;
     if ( NID(b) != N_B ) b = tobf(b,0);    } else if ( NID(a) <= N_R && NID(b) <= N_R )
     return mpfr_cmp(((BF)a)->body,((BF)b)->body);      return (*cmpnumt[MAX(NID(a),NID(b))])(a,b);
     else if ( NID(a) == N_B ) {
       switch ( NID(b) ) {
       case N_Q:
         if ( INT((Q)b) ) {
           z = ztogz((Q)b);
           ret = mpfr_cmp_z(((BF)a)->body,z->body);
         } else {
           q = qtogq((Q)b);
           ret = mpfr_cmp_q(((BF)a)->body,q->body);
         }
         break;
       case N_R:
         /* double precision = 53 */
         ret = mpfr_cmp_d(((BF)a)->body,((Real)b)->body);
         break;
       case N_B:
         ret = mpfr_cmp(((BF)a)->body,((BF)b)->body);
         break;
       default:
         goto err;
       }
       return ret;
     } else if ( NID(b)==N_B ) {
       switch ( NID(a) ) {
       case N_Q:
         if ( INT((Q)a) ) {
           z = ztogz((Q)a);
           ret = mpfr_cmp_z(((BF)b)->body,z->body);
         } else {
           q = qtogq((Q)a);
           ret = mpfr_cmp_q(((BF)b)->body,q->body);
         }
         break;
       case N_R:
         /* double precision = 53 */
         ret = mpfr_cmp_d(((BF)b)->body,((Real)a)->body);
         break;
       default:
         goto err;
       }
       return -ret;
   }    }
   err: error("cmpbf : cannot compare");
 }  }
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.14
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>