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

Annotation of OpenXM_contrib2/asir2000/engine/bf.c, Revision 1.11

1.2       noro        1: /*
1.11    ! noro        2:  * $OpenXM: OpenXM_contrib2/asir2000/engine/bf.c,v 1.10 2015/08/06 09:12:29 noro Exp $
1.7       noro        3:  */
1.1       noro        4: #include "ca.h"
                      5: #include "base.h"
                      6: #include <math.h>
                      7:
1.11    ! noro        8: extern int mpfr_roundmode;
        !             9:
1.7       noro       10: Num tobf(Num,int);
1.1       noro       11:
1.7       noro       12: #define BFPREC(a) (((BF)(a))->body->_mpfr_prec)
                     13:
                     14: void strtobf(char *s,BF *p)
                     15: {
                     16:   BF r;
                     17:   NEWBF(r);
                     18:   mpfr_init(r->body);
1.11    ! noro       19:   mpfr_set_str(r->body,s,10,mpfr_roundmode);
1.7       noro       20:   *p = r;
                     21: }
                     22:
                     23: double mpfrtodbl(mpfr_t a)
                     24: {
1.11    ! noro       25:   return mpfr_get_d(a,mpfr_roundmode);
1.7       noro       26: }
                     27:
                     28: Num tobf(Num a,int prec)
                     29: {
                     30:   mpfr_t r;
                     31:   mpz_t z;
                     32:   mpq_t q;
                     33:   BF d;
                     34:   N nm,dn;
                     35:   int sgn;
                     36:
                     37:   if ( !a ) {
1.8       noro       38:     prec ? mpfr_init2(r,prec) : mpfr_init(r);
1.7       noro       39:     mpfr_set_zero(r,1);
                     40:     MPFRTOBF(r,d);
                     41:     return (Num)d;
                     42:   } else {
                     43:     switch ( NID(a) ) {
                     44:     case N_B:
                     45:       return a;
                     46:       break;
1.8       noro       47:     case N_R:
                     48:       prec ? mpfr_init2(r,prec) : mpfr_init(r);
1.11    ! noro       49:       mpfr_init_set_d(r,((Real)a)->body,mpfr_roundmode);
1.7       noro       50:       MPFRTOBF(r,d);
                     51:       return (Num)d;
1.8       noro       52:       break;
1.7       noro       53:     case N_Q:
                     54:       nm = NM((Q)a); dn = DN((Q)a); sgn = SGN((Q)a);
                     55:       if ( INT((Q)a) ) {
                     56:         mpz_init(z);
                     57:         mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
                     58:         if ( sgn < 0 ) mpz_neg(z,z);
1.11    ! noro       59:         mpfr_init_set_z(r,z,mpfr_roundmode);
1.7       noro       60:       } else {
                     61:         mpq_init(q);
                     62:         mpz_import(mpq_numref(q),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
                     63:         mpz_import(mpq_denref(q),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn));
                     64:         if ( sgn < 0 ) mpq_neg(q,q);
1.11    ! noro       65:         mpfr_init_set_q(r,q,mpfr_roundmode);
1.7       noro       66:       }
                     67:       MPFRTOBF(r,d);
                     68:       return (Num)d;
                     69:       break;
                     70:     default:
                     71:       error("tobf : invalid argument");
1.8       noro       72:       break;
1.7       noro       73:     }
                     74:   }
                     75: }
                     76:
                     77: void addbf(Num a,Num b,Num *c)
                     78: {
                     79:   mpfr_t r;
                     80:   BF d;
1.8       noro       81:   GZ z;
                     82:   GQ q;
1.7       noro       83:
                     84:   if ( !a )
                     85:     *c = b;
                     86:   else if ( !b )
                     87:     *c = a;
                     88:   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
                     89:     (*addnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9       noro       90:   else if ( NID(a) == N_B ) {
                     91:     switch ( NID(b) ) {
                     92:     case N_Q:
                     93:       mpfr_init2(r,BFPREC(a));
                     94:       if ( INT((Q)b) ) {
                     95:         z = ztogz((Q)b);
1.11    ! noro       96:         mpfr_add_z(r,((BF)a)->body,z->body,mpfr_roundmode);
1.9       noro       97:       } else {
                     98:         q = qtogq((Q)b);
1.11    ! noro       99:         mpfr_add_q(r,((BF)a)->body,q->body,mpfr_roundmode);
1.9       noro      100:       }
                    101:       break;
                    102:     case N_R:
                    103:       /* double precision = 53 */
                    104:       mpfr_init2(r,MAX(BFPREC(a),53));
1.11    ! noro      105:       mpfr_add_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
1.9       noro      106:       break;
                    107:     case N_B:
                    108:       mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
1.11    ! noro      109:       mpfr_add(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      110:       break;
                    111:     }
                    112:     MPFRTOBF(r,d);
                    113:     *c = (Num)d;
                    114:   } else { /* NID(b)==N_B */
                    115:     switch ( NID(a) ) {
                    116:     case N_Q:
                    117:       mpfr_init2(r,BFPREC(b));
                    118:       if ( INT((Q)a) ) {
                    119:         z = ztogz((Q)a);
1.11    ! noro      120:         mpfr_add_z(r,((BF)b)->body,z->body,mpfr_roundmode);
1.9       noro      121:       } else {
                    122:         q = qtogq((Q)a);
1.11    ! noro      123:         mpfr_add_q(r,((BF)b)->body,q->body,mpfr_roundmode);
1.8       noro      124:       }
1.9       noro      125:       break;
                    126:     case N_R:
                    127:       /* double precision = 53 */
                    128:       mpfr_init2(r,MAX(BFPREC(b),53));
1.11    ! noro      129:       mpfr_add_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode);
1.9       noro      130:       break;
1.8       noro      131:     }
1.7       noro      132:     MPFRTOBF(r,d);
                    133:     *c = (Num)d;
                    134:   }
                    135: }
                    136:
                    137: void subbf(Num a,Num b,Num *c)
                    138: {
1.8       noro      139:   mpfr_t r,s;
                    140:   GZ z;
                    141:   GQ q;
1.7       noro      142:   BF d;
                    143:
                    144:   if ( !a )
                    145:     (*chsgnnumt[NID(b)])(b,c);
                    146:   else if ( !b )
                    147:     *c = a;
                    148:   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
                    149:     (*subnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9       noro      150:   else if ( NID(a) == N_B ) {
                    151:     switch ( NID(b) ) {
                    152:     case N_Q:
                    153:       mpfr_init2(r,BFPREC(a));
                    154:       if ( INT((Q)b) ) {
                    155:         z = ztogz((Q)b);
1.11    ! noro      156:         mpfr_sub_z(r,((BF)a)->body,z->body,mpfr_roundmode);
1.9       noro      157:       } else {
                    158:         q = qtogq((Q)b);
1.11    ! noro      159:         mpfr_sub_q(r,((BF)a)->body,q->body,mpfr_roundmode);
1.9       noro      160:       }
                    161:       break;
                    162:     case N_R:
                    163:       /* double precision = 53 */
                    164:       mpfr_init2(r,MAX(BFPREC(a),53));
1.11    ! noro      165:       mpfr_sub_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
1.9       noro      166:       break;
                    167:     case N_B:
                    168:       mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
1.11    ! noro      169:       mpfr_sub(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      170:       break;
                    171:     }
                    172:     MPFRTOBF(r,d);
                    173:     *c = (Num)d;
                    174:   } else { /* NID(b)==N_B */
                    175:     switch ( NID(a) ) {
                    176:     case N_Q:
                    177:       mpfr_init2(r,BFPREC(b));
                    178:       if ( INT((Q)a) ) {
                    179:         z = ztogz((Q)a);
1.11    ! noro      180:         mpfr_sub_z(r,((BF)b)->body,z->body,mpfr_roundmode);
1.9       noro      181:       } else {
                    182:         q = qtogq((Q)a);
1.11    ! noro      183:         mpfr_sub_q(r,((BF)b)->body,q->body,mpfr_roundmode);
1.8       noro      184:       }
1.11    ! noro      185:       mpfr_neg(r,r,mpfr_roundmode);
1.9       noro      186:       break;
                    187:     case N_R:
                    188:       /* double precision = 53 */
                    189:       mpfr_init2(r,MAX(BFPREC(b),53));
1.11    ! noro      190:       mpfr_d_sub(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      191:       break;
1.8       noro      192:     }
1.7       noro      193:     MPFRTOBF(r,d);
                    194:     *c = (Num)d;
                    195:   }
                    196: }
                    197:
                    198: void mulbf(Num a,Num b,Num *c)
                    199: {
                    200:   mpfr_t r;
1.8       noro      201:   GZ z;
                    202:   GQ q;
1.7       noro      203:   BF d;
                    204:   int prec;
                    205:
                    206:   if ( !a || !b )
                    207:     *c = 0;
                    208:   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
                    209:     (*mulnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9       noro      210:   else if ( NID(a) == N_B ) {
                    211:     switch ( NID(b) ) {
                    212:     case N_Q:
                    213:       mpfr_init2(r,BFPREC(a));
                    214:       if ( INT((Q)b) ) {
                    215:         z = ztogz((Q)b);
1.11    ! noro      216:         mpfr_mul_z(r,((BF)a)->body,z->body,mpfr_roundmode);
1.9       noro      217:       } else {
                    218:         q = qtogq((Q)b);
1.11    ! noro      219:         mpfr_mul_q(r,((BF)a)->body,q->body,mpfr_roundmode);
1.9       noro      220:       }
                    221:       break;
                    222:     case N_R:
                    223:       /* double precision = 53 */
                    224:       mpfr_init2(r,MAX(BFPREC(a),53));
1.11    ! noro      225:       mpfr_mul_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
1.9       noro      226:       break;
                    227:     case N_B:
                    228:       mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
1.11    ! noro      229:       mpfr_mul(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      230:       break;
                    231:     }
                    232:     MPFRTOBF(r,d);
                    233:     *c = (Num)d;
                    234:   } else { /* NID(b)==N_B */
                    235:     switch ( NID(a) ) {
                    236:     case N_Q:
                    237:       mpfr_init2(r,BFPREC(b));
                    238:       if ( INT((Q)a) ) {
                    239:         z = ztogz((Q)a);
1.11    ! noro      240:         mpfr_mul_z(r,((BF)b)->body,z->body,mpfr_roundmode);
1.9       noro      241:       } else {
                    242:         q = qtogq((Q)a);
1.11    ! noro      243:         mpfr_mul_q(r,((BF)b)->body,q->body,mpfr_roundmode);
1.8       noro      244:       }
1.9       noro      245:       break;
                    246:     case N_R:
                    247:       /* double precision = 53 */
                    248:       mpfr_init2(r,MAX(BFPREC(b),53));
1.11    ! noro      249:       mpfr_mul_d(r,((BF)b)->body,((Real)a)->body,mpfr_roundmode);
1.9       noro      250:       break;
1.8       noro      251:     }
1.7       noro      252:     MPFRTOBF(r,d);
                    253:     *c = (Num)d;
                    254:   }
                    255: }
                    256:
                    257: void divbf(Num a,Num b,Num *c)
                    258: {
1.8       noro      259:   mpfr_t s,r;
                    260:   GZ z;
                    261:   GQ q;
1.7       noro      262:   BF d;
                    263:
                    264:   if ( !b )
                    265:     error("divbf : division by 0");
                    266:   else if ( !a )
                    267:     *c = 0;
                    268:   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
                    269:     (*divnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9       noro      270:   else if ( NID(a) == N_B ) {
                    271:     switch ( NID(b) ) {
                    272:     case N_Q:
                    273:       mpfr_init2(r,BFPREC(a));
                    274:       if ( INT((Q)b) ) {
                    275:         z = ztogz((Q)b);
1.11    ! noro      276:         mpfr_div_z(r,((BF)a)->body,z->body,mpfr_roundmode);
1.9       noro      277:       } else {
                    278:         q = qtogq((Q)b);
1.11    ! noro      279:         mpfr_div_q(r,((BF)a)->body,q->body,mpfr_roundmode);
1.8       noro      280:       }
1.9       noro      281:       break;
                    282:     case N_R:
                    283:       /* double precision = 53 */
                    284:       mpfr_init2(r,MAX(BFPREC(a),53));
1.11    ! noro      285:       mpfr_div_d(r,((BF)a)->body,((Real)b)->body,mpfr_roundmode);
1.9       noro      286:       break;
                    287:     case N_B:
                    288:       mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
1.11    ! noro      289:       mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      290:       break;
                    291:     }
                    292:     MPFRTOBF(r,d);
                    293:     *c = (Num)d;
                    294:   } else { /* NID(b)==N_B */
                    295:     switch ( NID(a) ) {
                    296:     case N_Q:
                    297:       /* XXX : mpfr_z_div and mpfr_q_div are not implemented */
                    298:       a = tobf(a,BFPREC(b));
                    299:       mpfr_init2(r,BFPREC(b));
1.11    ! noro      300:       mpfr_div(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      301:       break;
                    302:     case N_R:
                    303:       /* double precision = 53 */
                    304:       mpfr_init2(r,MAX(BFPREC(b),53));
1.11    ! noro      305:       mpfr_d_div(r,((Real)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      306:       break;
1.8       noro      307:     }
1.7       noro      308:     MPFRTOBF(r,d);
                    309:     *c = (Num)d;
                    310:   }
                    311: }
                    312:
                    313: void pwrbf(Num a,Num b,Num *c)
                    314: {
1.8       noro      315:   int prec;
1.7       noro      316:   mpfr_t r;
1.8       noro      317:   GZ z;
1.7       noro      318:   BF d;
                    319:
                    320:   if ( !b )
                    321:     *c = (Num)ONE;
                    322:   else if ( !a )
                    323:     *c = 0;
                    324:   else if ( (NID(a) <= N_A) && (NID(b) <= N_A ) )
                    325:     (*pwrnumt[MIN(NID(a),NID(b))])(a,b,c);
1.9       noro      326:   else if ( NID(a) == N_B ) {
                    327:     switch ( NID(b) ) {
                    328:     case N_Q:
                    329:       mpfr_init2(r,BFPREC(a));
                    330:       if ( INT((Q)b) ) {
                    331:         z = ztogz((Q)b);
1.11    ! noro      332:         mpfr_pow_z(r,((BF)a)->body,z->body,mpfr_roundmode);
1.9       noro      333:       } else {
                    334:         b = tobf(b,BFPREC(a));
1.11    ! noro      335:         mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.8       noro      336:       }
1.9       noro      337:       break;
                    338:     case N_R:
                    339:       /* double precision = 53 */
                    340:       prec = MAX(BFPREC(a),53);
                    341:       mpfr_init2(r,prec);
                    342:       b = tobf(b,prec);
1.11    ! noro      343:       mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      344:       break;
                    345:     case N_B:
                    346:       mpfr_init2(r,MAX(BFPREC(a),BFPREC(b)));
1.11    ! noro      347:       mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      348:       break;
                    349:     }
                    350:     MPFRTOBF(r,d);
                    351:     *c = (Num)d;
                    352:   } else { /* NID(b)==N_B */
                    353:     switch ( NID(a) ) {
                    354:     case N_Q:
                    355:       mpfr_init2(r,BFPREC(b));
                    356:       a = tobf(a,BFPREC(b));
1.11    ! noro      357:       mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      358:       break;
                    359:     case N_R:
                    360:       /* double precision = 53 */
                    361:       prec = MAX(BFPREC(a),53);
                    362:       mpfr_init2(r,prec);
                    363:       a = tobf(a,prec);
1.11    ! noro      364:       mpfr_pow(r,((BF)a)->body,((BF)b)->body,mpfr_roundmode);
1.9       noro      365:       break;
1.8       noro      366:     }
1.7       noro      367:     MPFRTOBF(r,d);
                    368:     *c = (Num)d;
                    369:   }
                    370: }
                    371:
                    372: void chsgnbf(Num a,Num *c)
                    373: {
                    374:   mpfr_t r;
                    375:   BF d;
                    376:
                    377:   if ( !a )
                    378:     *c = 0;
                    379:   else if ( NID(a) <= N_A )
                    380:     (*chsgnnumt[NID(a)])(a,c);
                    381:   else {
                    382:     mpfr_init2(r,BFPREC(a));
1.11    ! noro      383:     mpfr_neg(r,((BF)a)->body,mpfr_roundmode);
1.7       noro      384:     MPFRTOBF(r,d);
                    385:     *c = (Num)d;
                    386:   }
                    387: }
                    388:
                    389: int cmpbf(Num a,Num b)
                    390: {
1.9       noro      391:   int ret;
                    392:   GZ z;
                    393:   GQ q;
                    394:
1.7       noro      395:   if ( !a ) {
                    396:     if ( !b || (NID(b)<=N_A) )
                    397:       return (*cmpnumt[NID(b)])(a,b);
                    398:     else
                    399:       return -mpfr_sgn(((BF)a)->body);
                    400:   } else if ( !b ) {
                    401:     if ( !a || (NID(a)<=N_A) )
                    402:       return (*cmpnumt[NID(a)])(a,b);
                    403:     else
                    404:       return mpfr_sgn(((BF)a)->body);
1.9       noro      405:   } else if ( NID(a) == N_B ) {
                    406:     switch ( NID(b) ) {
                    407:     case N_Q:
                    408:       if ( INT((Q)b) ) {
                    409:         z = ztogz((Q)b);
                    410:         ret = mpfr_cmp_z(((BF)a)->body,z->body);
                    411:       } else {
                    412:         q = qtogq((Q)b);
                    413:         ret = mpfr_cmp_q(((BF)a)->body,q->body);
                    414:       }
                    415:       break;
                    416:     case N_R:
                    417:       /* double precision = 53 */
                    418:       ret = mpfr_cmp_d(((BF)a)->body,((Real)b)->body);
                    419:       break;
                    420:     case N_B:
                    421:       ret = mpfr_cmp(((BF)a)->body,((BF)b)->body);
                    422:       break;
                    423:     }
                    424:     return ret;
                    425:   } else { /* NID(b)==N_B */
                    426:     switch ( NID(a) ) {
                    427:     case N_Q:
                    428:       if ( INT((Q)a) ) {
                    429:         z = ztogz((Q)a);
                    430:         ret = mpfr_cmp_z(((BF)b)->body,z->body);
                    431:       } else {
                    432:         q = qtogq((Q)a);
                    433:         ret = mpfr_cmp_q(((BF)b)->body,q->body);
                    434:       }
                    435:       break;
                    436:     case N_R:
                    437:       /* double precision = 53 */
                    438:       ret = mpfr_cmp_d(((BF)b)->body,((Real)a)->body);
                    439:       break;
                    440:     }
                    441:     return -ret;
1.7       noro      442:   }
1.1       noro      443: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>