[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.10 and 1.11

version 1.10, 2015/08/06 09:12:29 version 1.11, 2015/08/06 23:41:52
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.9 2015/08/06 00:43:20 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.10 2015/08/06 09:12:29 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 44  Num tobf(Num a,int prec)
Line 46  Num tobf(Num a,int prec)
       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;
Line 54  Num tobf(Num a,int prec)
Line 56  Num tobf(Num a,int prec)
         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;
Line 91  void addbf(Num a,Num b,Num *c)
Line 93  void addbf(Num a,Num b,Num *c)
       mpfr_init2(r,BFPREC(a));        mpfr_init2(r,BFPREC(a));
       if ( INT((Q)b) ) {        if ( INT((Q)b) ) {
         z = ztogz((Q)b);          z = ztogz((Q)b);
         mpfr_add_z(r,((BF)a)->body,z->body,MPFR_RNDN);          mpfr_add_z(r,((BF)a)->body,z->body,mpfr_roundmode);
       } else {        } else {
         q = qtogq((Q)b);          q = qtogq((Q)b);
         mpfr_add_q(r,((BF)a)->body,q->body,MPFR_RNDN);          mpfr_add_q(r,((BF)a)->body,q->body,mpfr_roundmode);
       }        }
       break;        break;
     case N_R:      case N_R:
       /* double precision = 53 */        /* double precision = 53 */
       mpfr_init2(r,MAX(BFPREC(a),53));        mpfr_init2(r,MAX(BFPREC(a),53));
       mpfr_add_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);        mpfr_add_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
       break;        break;
     case N_B:      case N_B:
       mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
       mpfr_add(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_add(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     }      }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
Line 115  void addbf(Num a,Num b,Num *c)
Line 117  void addbf(Num a,Num b,Num *c)
       mpfr_init2(r,BFPREC(b));        mpfr_init2(r,BFPREC(b));
       if ( INT((Q)a) ) {        if ( INT((Q)a) ) {
         z = ztogz((Q)a);          z = ztogz((Q)a);
         mpfr_add_z(r,((BF)b)->body,z->body,MPFR_RNDN);          mpfr_add_z(r,((BF)b)->body,z->body,mpfr_roundmode);
       } else {        } else {
         q = qtogq((Q)a);          q = qtogq((Q)a);
         mpfr_add_q(r,((BF)b)->body,q->body,MPFR_RNDN);          mpfr_add_q(r,((BF)b)->body,q->body,mpfr_roundmode);
       }        }
       break;        break;
     case N_R:      case N_R:
       /* double precision = 53 */        /* double precision = 53 */
       mpfr_init2(r,MAX(BFPREC(b),53));        mpfr_init2(r,MAX(BFPREC(b),53));
       mpfr_add_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN);        mpfr_add_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode);
       break;        break;
     }      }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
Line 151  void subbf(Num a,Num b,Num *c)
Line 153  void subbf(Num a,Num b,Num *c)
       mpfr_init2(r,BFPREC(a));        mpfr_init2(r,BFPREC(a));
       if ( INT((Q)b) ) {        if ( INT((Q)b) ) {
         z = ztogz((Q)b);          z = ztogz((Q)b);
         mpfr_sub_z(r,((BF)a)->body,z->body,MPFR_RNDN);          mpfr_sub_z(r,((BF)a)->body,z->body,mpfr_roundmode);
       } else {        } else {
         q = qtogq((Q)b);          q = qtogq((Q)b);
         mpfr_sub_q(r,((BF)a)->body,q->body,MPFR_RNDN);          mpfr_sub_q(r,((BF)a)->body,q->body,mpfr_roundmode);
       }        }
       break;        break;
     case N_R:      case N_R:
       /* double precision = 53 */        /* double precision = 53 */
       mpfr_init2(r,MAX(BFPREC(a),53));        mpfr_init2(r,MAX(BFPREC(a),53));
       mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);        mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
       break;        break;
     case N_B:      case N_B:
       mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
       mpfr_sub(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_sub(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     }      }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
Line 175  void subbf(Num a,Num b,Num *c)
Line 177  void subbf(Num a,Num b,Num *c)
       mpfr_init2(r,BFPREC(b));        mpfr_init2(r,BFPREC(b));
       if ( INT((Q)a) ) {        if ( INT((Q)a) ) {
         z = ztogz((Q)a);          z = ztogz((Q)a);
         mpfr_sub_z(r,((BF)b)->body,z->body,MPFR_RNDN);          mpfr_sub_z(r,((BF)b)->body,z->body,mpfr_roundmode);
       } else {        } else {
         q = qtogq((Q)a);          q = qtogq((Q)a);
         mpfr_sub_q(r,((BF)b)->body,q->body,MPFR_RNDN);          mpfr_sub_q(r,((BF)b)->body,q->body,mpfr_roundmode);
       }        }
       mpfr_neg(r,r,MPFR_RNDN);        mpfr_neg(r,r,mpfr_roundmode);
       break;        break;
     case N_R:      case N_R:
       /* double precision = 53 */        /* double precision = 53 */
       mpfr_init2(r,MAX(BFPREC(b),53));        mpfr_init2(r,MAX(BFPREC(b),53));
       mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     }      }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
Line 211  void mulbf(Num a,Num b,Num *c)
Line 213  void mulbf(Num a,Num b,Num *c)
       mpfr_init2(r,BFPREC(a));        mpfr_init2(r,BFPREC(a));
       if ( INT((Q)b) ) {        if ( INT((Q)b) ) {
         z = ztogz((Q)b);          z = ztogz((Q)b);
         mpfr_mul_z(r,((BF)a)->body,z->body,MPFR_RNDN);          mpfr_mul_z(r,((BF)a)->body,z->body,mpfr_roundmode);
       } else {        } else {
         q = qtogq((Q)b);          q = qtogq((Q)b);
         mpfr_mul_q(r,((BF)a)->body,q->body,MPFR_RNDN);          mpfr_mul_q(r,((BF)a)->body,q->body,mpfr_roundmode);
       }        }
       break;        break;
     case N_R:      case N_R:
       /* double precision = 53 */        /* double precision = 53 */
       mpfr_init2(r,MAX(BFPREC(a),53));        mpfr_init2(r,MAX(BFPREC(a),53));
       mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);        mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
       break;        break;
     case N_B:      case N_B:
       mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
       mpfr_mul(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_mul(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     }      }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
Line 235  void mulbf(Num a,Num b,Num *c)
Line 237  void mulbf(Num a,Num b,Num *c)
       mpfr_init2(r,BFPREC(b));        mpfr_init2(r,BFPREC(b));
       if ( INT((Q)a) ) {        if ( INT((Q)a) ) {
         z = ztogz((Q)a);          z = ztogz((Q)a);
         mpfr_mul_z(r,((BF)b)->body,z->body,MPFR_RNDN);          mpfr_mul_z(r,((BF)b)->body,z->body,mpfr_roundmode);
       } else {        } else {
         q = qtogq((Q)a);          q = qtogq((Q)a);
         mpfr_mul_q(r,((BF)b)->body,q->body,MPFR_RNDN);          mpfr_mul_q(r,((BF)b)->body,q->body,mpfr_roundmode);
       }        }
       break;        break;
     case N_R:      case N_R:
       /* double precision = 53 */        /* double precision = 53 */
       mpfr_init2(r,MAX(BFPREC(b),53));        mpfr_init2(r,MAX(BFPREC(b),53));
       mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN);        mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode);
       break;        break;
     }      }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
Line 271  void divbf(Num a,Num b,Num *c)
Line 273  void divbf(Num a,Num b,Num *c)
       mpfr_init2(r,BFPREC(a));        mpfr_init2(r,BFPREC(a));
       if ( INT((Q)b) ) {        if ( INT((Q)b) ) {
         z = ztogz((Q)b);          z = ztogz((Q)b);
         mpfr_div_z(r,((BF)a)->body,z->body,MPFR_RNDN);          mpfr_div_z(r,((BF)a)->body,z->body,mpfr_roundmode);
       } else {        } else {
         q = qtogq((Q)b);          q = qtogq((Q)b);
         mpfr_div_q(r,((BF)a)->body,q->body,MPFR_RNDN);          mpfr_div_q(r,((BF)a)->body,q->body,mpfr_roundmode);
       }        }
       break;        break;
     case N_R:      case N_R:
       /* double precision = 53 */        /* double precision = 53 */
       mpfr_init2(r,MAX(BFPREC(a),53));        mpfr_init2(r,MAX(BFPREC(a),53));
       mpfr_div_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);        mpfr_div_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
       break;        break;
     case N_B:      case N_B:
       mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
       mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     }      }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
Line 295  void divbf(Num a,Num b,Num *c)
Line 297  void divbf(Num a,Num b,Num *c)
       /* XXX : mpfr_z_div and mpfr_q_div are not implemented */        /* XXX : mpfr_z_div and mpfr_q_div are not implemented */
       a = tobf(a,BFPREC(b));        a = tobf(a,BFPREC(b));
       mpfr_init2(r,BFPREC(b));        mpfr_init2(r,BFPREC(b));
       mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     case N_R:      case N_R:
       /* double precision = 53 */        /* double precision = 53 */
       mpfr_init2(r,MAX(BFPREC(b),53));        mpfr_init2(r,MAX(BFPREC(b),53));
       mpfr_d_div(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_d_div(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     }      }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
Line 327  void pwrbf(Num a,Num b,Num *c)
Line 329  void pwrbf(Num a,Num b,Num *c)
       mpfr_init2(r,BFPREC(a));        mpfr_init2(r,BFPREC(a));
       if ( INT((Q)b) ) {        if ( INT((Q)b) ) {
         z = ztogz((Q)b);          z = ztogz((Q)b);
         mpfr_pow_z(r,((BF)a)->body,z->body,MPFR_RNDN);          mpfr_pow_z(r,((BF)a)->body,z->body,mpfr_roundmode);
       } else {        } else {
         b = tobf(b,BFPREC(a));          b = tobf(b,BFPREC(a));
         mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);          mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
       }        }
       break;        break;
     case N_R:      case N_R:
Line 338  void pwrbf(Num a,Num b,Num *c)
Line 340  void pwrbf(Num a,Num b,Num *c)
       prec = MAX(BFPREC(a),53);        prec = MAX(BFPREC(a),53);
       mpfr_init2(r,prec);        mpfr_init2(r,prec);
       b = tobf(b,prec);        b = tobf(b,prec);
       mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     case N_B:      case N_B:
       mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
       mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     }      }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
Line 352  void pwrbf(Num a,Num b,Num *c)
Line 354  void pwrbf(Num a,Num b,Num *c)
     case N_Q:      case N_Q:
       mpfr_init2(r,BFPREC(b));        mpfr_init2(r,BFPREC(b));
       a = tobf(a,BFPREC(b));        a = tobf(a,BFPREC(b));
       mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     case N_R:      case N_R:
       /* double precision = 53 */        /* double precision = 53 */
       prec = MAX(BFPREC(a),53);        prec = MAX(BFPREC(a),53);
       mpfr_init2(r,prec);        mpfr_init2(r,prec);
       a = tobf(a,prec);        a = tobf(a,prec);
       mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);        mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
       break;        break;
     }      }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
Line 378  void chsgnbf(Num a,Num *c)
Line 380  void chsgnbf(Num a,Num *c)
     (*chsgnnumt[NID(a)])(a,c);      (*chsgnnumt[NID(a)])(a,c);
   else {    else {
     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;
   }    }
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.11
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>