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

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

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