[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.8

version 1.7, 2015/08/04 06:20:45 version 1.8, 2015/08/05 03:46:35
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.7 2015/08/04 06:20:45 noro Exp $
  */   */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
Line 33  Num tobf(Num a,int prec)
Line 33  Num tobf(Num a,int prec)
   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 42  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_RNDN);
       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) ) {
Line 67  Num tobf(Num a,int prec)
Line 67  Num tobf(Num a,int prec)
       break;        break;
     default:      default:
       error("tobf : invalid argument");        error("tobf : invalid argument");
           break;        break;
     }      }
   }    }
 }  }
Line 76  void addbf(Num a,Num b,Num *c)
Line 76  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;
Line 84  void addbf(Num a,Num b,Num *c)
Line 86  void addbf(Num a,Num b,Num *c)
   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
     (*addnumt[MIN(NID(a),NID(b))])(a,b,c);      (*addnumt[MIN(NID(a),NID(b))])(a,b,c);
   else {    else {
     if ( NID(a) != N_B ) a = tobf(a,0);      if ( NID(a) == N_B ) {
     if ( NID(b) != N_B ) b = tobf(b,0);        switch ( NID(b) ) {
     mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        case N_Q:
     mpfr_add(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);          mpfr_init2(r,BFPREC(a));
           if ( INT((Q)b) ) {
             z = ztogz((Q)b);
             mpfr_add_z(r,((BF)a)->body,z->body,MPFR_RNDN);
           } else {
             q = qtogq((Q)b);
             mpfr_add_q(r,((BF)a)->body,q->body,MPFR_RNDN);
           }
           break;
         case N_R:
           /* double precision = 53 */
           mpfr_init2(r,MIN(BFPREC(a),53));
           mpfr_add_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
           break;
         case N_B:
           mpfr_init2(r,MIN(BFPREC(a),BFPREC(b)));
           mpfr_add(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
           break;
         }
       } else { /* 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_RNDN);
           } else {
             q = qtogq((Q)a);
             mpfr_add_q(r,((BF)b)->body,q->body,MPFR_RNDN);
           }
           break;
         case N_R:
           /* double precision = 53 */
           mpfr_init2(r,MIN(BFPREC(b),53));
           mpfr_add_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN);
           break;
         }
       }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    }
Line 95  void addbf(Num a,Num b,Num *c)
Line 134  void addbf(Num a,Num b,Num *c)
   
 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 )
Line 105  void subbf(Num a,Num b,Num *c)
Line 146  void subbf(Num a,Num b,Num *c)
   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
     (*subnumt[MIN(NID(a),NID(b))])(a,b,c);      (*subnumt[MIN(NID(a),NID(b))])(a,b,c);
   else {    else {
     if ( NID(a) != N_B ) a = tobf(a,0);      if ( NID(a) == N_B ) {
     if ( NID(b) != N_B ) b = tobf(b,0);        switch ( NID(b) ) {
     mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        case N_Q:
     mpfr_sub(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);          mpfr_init2(r,BFPREC(a));
           if ( INT((Q)b) ) {
             z = ztogz((Q)b);
             mpfr_sub_z(r,((BF)a)->body,z->body,MPFR_RNDN);
           } else {
             q = qtogq((Q)b);
             mpfr_sub_q(r,((BF)a)->body,q->body,MPFR_RNDN);
           }
           break;
         case N_R:
           /* double precision = 53 */
           mpfr_init2(r,MIN(BFPREC(a),53));
           mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
           break;
         case N_B:
           mpfr_init2(r,MIN(BFPREC(a),BFPREC(b)));
           mpfr_sub(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
           break;
         }
       } else { /* NID(b)==N_B */
         switch ( NID(a) ) {
         case N_Q:
           mpfr_init2(s,BFPREC(b));
           mpfr_init2(r,BFPREC(b));
           if ( INT((Q)a) ) {
             z = ztogz((Q)a);
             mpfr_sub_z(s,((BF)b)->body,z->body,MPFR_RNDN);
           } else {
             q = qtogq((Q)a);
             mpfr_sub_q(s,((BF)b)->body,q->body,MPFR_RNDN);
           }
           mpfr_neg(r,s,MPFR_RNDN);
           break;
         case N_R:
           /* double precision = 53 */
           mpfr_init2(r,MIN(BFPREC(b),53));
           mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN);
           break;
         }
       }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    }
Line 117  void subbf(Num a,Num b,Num *c)
Line 197  void subbf(Num a,Num b,Num *c)
 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;
   
Line 125  void mulbf(Num a,Num b,Num *c)
Line 207  void mulbf(Num a,Num b,Num *c)
   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
     (*mulnumt[MIN(NID(a),NID(b))])(a,b,c);      (*mulnumt[MIN(NID(a),NID(b))])(a,b,c);
   else {    else {
     if ( NID(a) != N_B ) a = tobf(a,0);      if ( NID(a) == N_B ) {
     if ( NID(b) != N_B ) b = tobf(b,0);        switch ( NID(b) ) {
     mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        case N_Q:
     mpfr_mul(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);          mpfr_init2(r,BFPREC(a));
           if ( INT((Q)b) ) {
             z = ztogz((Q)b);
             mpfr_mul_z(r,((BF)a)->body,z->body,MPFR_RNDN);
           } else {
             q = qtogq((Q)b);
             mpfr_mul_q(r,((BF)a)->body,q->body,MPFR_RNDN);
           }
           break;
         case N_R:
           /* double precision = 53 */
           mpfr_init2(r,MIN(BFPREC(a),53));
           mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
           break;
         case N_B:
           mpfr_init2(r,MIN(BFPREC(a),BFPREC(b)));
           mpfr_mul(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
           break;
         }
       } else { /* 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_RNDN);
           } else {
             q = qtogq((Q)a);
             mpfr_mul_q(r,((BF)b)->body,q->body,MPFR_RNDN);
           }
           break;
         case N_R:
           /* double precision = 53 */
           mpfr_init2(r,MIN(BFPREC(b),53));
           mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,MPFR_RNDN);
           break;
         }
       }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    }
Line 136  void mulbf(Num a,Num b,Num *c)
Line 255  void mulbf(Num a,Num b,Num *c)
   
 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 )
Line 146  void divbf(Num a,Num b,Num *c)
Line 267  void divbf(Num a,Num b,Num *c)
   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
     (*divnumt[MIN(NID(a),NID(b))])(a,b,c);      (*divnumt[MIN(NID(a),NID(b))])(a,b,c);
   else {    else {
     if ( NID(a) != N_B ) a = tobf(a,0);      if ( NID(a) == N_B ) {
     if ( NID(b) != N_B ) b = tobf(b,0);        switch ( NID(b) ) {
     mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        case N_Q:
     mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);          mpfr_init2(r,BFPREC(a));
           if ( INT((Q)b) ) {
             z = ztogz((Q)b);
             mpfr_div_z(r,((BF)a)->body,z->body,MPFR_RNDN);
           } else {
             q = qtogq((Q)b);
             mpfr_div_q(r,((BF)a)->body,q->body,MPFR_RNDN);
           }
           break;
         case N_R:
           /* double precision = 53 */
           mpfr_init2(r,MIN(BFPREC(a),53));
           mpfr_div_d(r,((BF)a)->body,((Real)b)->body,MPFR_RNDN);
           break;
         case N_B:
           mpfr_init2(r,MIN(BFPREC(a),BFPREC(b)));
           mpfr_div(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
           break;
         }
       } else { /* 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_RNDN);
           break;
         case N_R:
           /* double precision = 53 */
           mpfr_init2(r,MIN(BFPREC(b),53));
           mpfr_d_div(r,((Real)a)->body,((BF)b)->body,MPFR_RNDN);
           break;
         }
       }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    }
Line 157  void divbf(Num a,Num b,Num *c)
Line 311  void divbf(Num a,Num b,Num *c)
   
 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 )
Line 167  void pwrbf(Num a,Num b,Num *c)
Line 323  void pwrbf(Num a,Num b,Num *c)
   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )    else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
     (*pwrnumt[MIN(NID(a),NID(b))])(a,b,c);      (*pwrnumt[MIN(NID(a),NID(b))])(a,b,c);
   else {    else {
     if ( NID(a) != N_B ) a = tobf(a,0);      if ( NID(a) == N_B ) {
     if ( NID(b) != N_B ) b = tobf(b,0);        switch ( NID(b) ) {
     mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));        case N_Q:
     mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);          mpfr_init2(r,BFPREC(a));
           if ( INT((Q)b) ) {
             z = ztogz((Q)b);
             mpfr_pow_z(r,((BF)a)->body,z->body,MPFR_RNDN);
           } else {
             b = tobf(b,BFPREC(a));
             mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
           }
           break;
         case N_R:
           /* double precision = 53 */
           prec = MIN(BFPREC(a),53);
           mpfr_init2(r,prec);
           b = tobf(b,prec);
           mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
           break;
         case N_B:
           mpfr_init2(r,MIN(BFPREC(a),BFPREC(b)));
           mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
           break;
         }
       } else { /* 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_RNDN);
           break;
         case N_R:
           /* double precision = 53 */
           prec = MIN(BFPREC(a),53);
           mpfr_init2(r,prec);
           a = tobf(a,prec);
           mpfr_pow(r,((BF)a)->body,((BF)b)->body,MPFR_RNDN);
           break;
         }
       }
     MPFRTOBF(r,d);      MPFRTOBF(r,d);
     *c = (Num)d;      *c = (Num)d;
   }    }
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>