[BACK]Return to bfaux.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2018 / builtin

Annotation of OpenXM_contrib2/asir2018/builtin/bfaux.c, Revision 1.3

1.3     ! noro        1: /* $OpenXM: OpenXM_contrib2/asir2018/builtin/bfaux.c,v 1.2 2018/09/28 08:20:27 noro Exp $ */
1.1       noro        2: #include "ca.h"
                      3: #include "parse.h"
                      4:
                      5: void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround();
                      6: void Pmpfr_ai();
                      7: void Pmpfr_eint(), Pmpfr_erf(), Pmpfr_erfc(), Pmpfr_li2();
                      8: void Pmpfr_zeta();
                      9: void Pmpfr_j0(), Pmpfr_j1();
                     10: void Pmpfr_y0(), Pmpfr_y1();
                     11: void Pmpfr_gamma(), Pmpfr_lngamma(), Pmpfr_digamma();
                     12: void Pmpfr_floor(), Pmpfr_round(), Pmpfr_ceil();
                     13: void Prk_ratmat();
                     14: void mp_sin(),mp_cos(),mp_tan(),mp_asin(),mp_acos(),mp_atan();
                     15: void mp_sinh(),mp_cosh(),mp_tanh(),mp_asinh(),mp_acosh(),mp_atanh();
                     16: void mp_exp(),mp_log(),mp_pow();
                     17: void mp_factorial(),mp_abs();
                     18:
                     19: struct ftab bf_tab[] = {
                     20:   {"eval",Peval,-2},
                     21:   {"setprec",Psetprec,-1},
                     22:   {"setbprec",Psetbprec,-1},
                     23:   {"setround",Psetround,-1},
                     24:   {"todouble",Ptodouble,1},
                     25:   {"mpfr_sin",mp_sin,-2},
                     26:   {"mpfr_cos",mp_cos,-2},
                     27:   {"mpfr_tan",mp_tan,-2},
                     28:   {"mpfr_asin",mp_asin,-2},
                     29:   {"mpfr_acos",mp_acos,-2},
                     30:   {"mpfr_atan",mp_atan,-2},
                     31:   {"mpfr_sinh",mp_sinh,-2},
                     32:   {"mpfr_cosh",mp_cosh,-2},
                     33:   {"mpfr_tanh",mp_tanh,-2},
                     34:   {"mpfr_asinh",mp_asinh,-2},
                     35:   {"mpfr_acosh",mp_acosh,-2},
                     36:   {"mpfr_atanh",mp_atanh,-2},
                     37:   {"mpfr_exp",mp_exp,-2},
                     38:   {"mpfr_log",mp_log,-2},
                     39:   {"mpfr_pow",mp_pow,-3},
                     40:   {"mpfr_ai",Pmpfr_ai,-2},
                     41:   {"mpfr_zeta",Pmpfr_zeta,-2},
                     42:   {"mpfr_j0",Pmpfr_j0,-2},
                     43:   {"mpfr_j1",Pmpfr_j1,-2},
                     44:   {"mpfr_y0",Pmpfr_y0,-2},
                     45:   {"mpfr_y1",Pmpfr_y1,-2},
                     46:   {"mpfr_eint",Pmpfr_eint,-2},
                     47:   {"mpfr_erf",Pmpfr_erf,-2},
                     48:   {"mpfr_erfc",Pmpfr_erfc,-2},
                     49:   {"mpfr_li2",Pmpfr_li2,-2},
                     50:   {"mpfr_gamma",Pmpfr_gamma,-2},
                     51:   {"mpfr_lngamma",Pmpfr_gamma,-2},
                     52:   {"mpfr_digamma",Pmpfr_gamma,-2},
                     53:   {"mpfr_floor",Pmpfr_floor,-2},
                     54:   {"mpfr_ceil",Pmpfr_ceil,-2},
                     55:   {"mpfr_round",Pmpfr_round,-2},
                     56:   {"rk_ratmat",Prk_ratmat,7},
                     57:   {0,0,0},
                     58: };
                     59:
                     60: int mpfr_roundmode = MPFR_RNDN;
                     61:
                     62: void todoublen(Num a,Num *rp)
                     63: {
                     64:   double r,i;
                     65:   Real real,imag;
                     66:
                     67:   if ( !a ) {
                     68:     *rp = 0;
                     69:     return;
                     70:   }
                     71:   switch ( NID(a) ) {
                     72:     case N_R: case N_Q: case N_B:
                     73:       r = ToReal(a);
                     74:       MKReal(r,real);
                     75:       *rp = (Num)real;
                     76:       break;
                     77:     case N_C:
                     78:       r = ToReal(((C)a)->r);
                     79:       i = ToReal(((C)a)->i);
                     80:       MKReal(r,real);
                     81:       MKReal(i,imag);
                     82:       reimtocplx((Num)real,(Num)imag,rp);
                     83:       break;
                     84:     default:
                     85:       *rp = a;
                     86:       break;
                     87:   }
                     88: }
                     89:
                     90: void todoublep(P a,P *rp)
                     91: {
                     92:   DCP dc,dcr,dcr0;
                     93:
                     94:   if ( !a ) *rp = 0;
                     95:   else if ( OID(a) == O_N ) todoublen((Num)a,(Num *)rp);
                     96:   else {
                     97:     for ( dcr0 = 0, dc = DC(a); dc; dc = NEXT(dc) ) {
                     98:       NEXTDC(dcr0,dcr);
                     99:       DEG(dcr) = DEG(dc);
                    100:       todoublep(COEF(dc),&COEF(dcr));
                    101:     }
                    102:     NEXT(dcr) = 0;
                    103:     MKP(VR(a),dcr0,*rp);
                    104:   }
                    105: }
                    106:
                    107: void todoubler(R a,R *rp)
                    108: {
                    109:   R b;
                    110:
                    111:   if ( !a ) *rp = 0;
                    112:   else if ( OID(a) <= O_P ) todoublep((P)a,(P *)rp);
                    113:   else {
                    114:     NEWR(b);
                    115:     todoublep(a->nm,&b->nm);
                    116:     todoublep(a->dn,&b->dn);
                    117:     *rp = b;
                    118:   }
                    119: }
                    120:
                    121: void todouble(Obj a,Obj *b)
                    122: {
                    123:   Obj t;
                    124:   LIST l;
                    125:   V v;
                    126:   int row,col,len;
                    127:   VECT vect;
                    128:   MAT mat;
                    129:   int i,j;
                    130:   NODE n0,n,nd;
                    131:   MP m,mp,mp0;
                    132:   DP d;
                    133:
                    134:   if ( !a ) {
                    135:     *b = 0;
                    136:     return;
                    137:   }
                    138:   switch ( OID(a) ) {
                    139:     case O_N:
                    140:       todoublen((Num)a,(Num *)b);
                    141:       break;
                    142:     case O_P:
                    143:       todoublep((P)a,(P *)b);
                    144:       break;
                    145:     case O_R:
                    146:       todoubler((R)a,(R *)b);
                    147:       break;
                    148:     case O_LIST:
                    149:       n0 = 0;
                    150:       for ( nd = BDY((LIST)a); nd; nd = NEXT(nd) ) {
                    151:         NEXTNODE(n0,n);
                    152:         todouble((Obj)BDY(nd),(Obj *)&BDY(n));
                    153:       }
                    154:       if ( n0 )
                    155:         NEXT(n) = 0;
                    156:       MKLIST(l,n0);
                    157:       *b = (Obj)l;
                    158:       break;
                    159:     case O_VECT:
                    160:       len = ((VECT)a)->len;
                    161:       MKVECT(vect,len);
                    162:       for ( i = 0; i < len; i++ ) {
                    163:         todouble((Obj)BDY((VECT)a)[i],(Obj *)&BDY(vect)[i]);
                    164:       }
                    165:       *b = (Obj)vect;
                    166:       break;
                    167:     case O_MAT:
                    168:       row = ((MAT)a)->row;
                    169:       col = ((MAT)a)->col;
                    170:       MKMAT(mat,row,col);
                    171:       for ( i = 0; i < row; i++ )
                    172:         for ( j = 0; j < col; j++ ) {
                    173:           todouble((Obj)BDY((MAT)a)[i][j],(Obj *)&BDY(mat)[i][j]);
                    174:         }
                    175:       *b = (Obj)mat;
                    176:       break;
                    177:     case O_DP:
                    178:       mp0 = 0;
                    179:       for ( m = BDY((DP)a); m; m = NEXT(m) ) {
                    180:         todouble(C(m),&t);
                    181:         if ( t ) {
                    182:           NEXTMP(mp0,mp);
                    183:           C(mp) = t;
                    184:           mp->dl = m->dl;
                    185:         }
                    186:       }
                    187:       if ( mp0 ) {
                    188:         MKDP(NV((DP)a),mp0,d);
                    189:         d->sugar = ((DP)a)->sugar;
                    190:         *b = (Obj)d;
                    191:       } else
                    192:         *b = 0;
                    193:
                    194:       break;
                    195:     default:
                    196:       error("todouble : invalid argument");
                    197:   }
                    198: }
                    199:
                    200: void Ptodouble(NODE arg,Obj *rp)
                    201: {
                    202:   todouble((Obj)ARG0(arg),rp);
                    203: }
                    204:
                    205: void Peval(NODE arg,Obj *rp)
                    206: {
                    207:   long prec;
                    208:
                    209:   asir_assert(ARG0(arg),O_R,"eval");
                    210:   if ( argc(arg) == 2 ) {
1.2       noro      211:     prec = ZTOS((Q)ARG1(arg))*3.32193;
1.1       noro      212:     if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
                    213:     else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
                    214:   } else
                    215:     prec = 0;
                    216:   evalr(CO,(Obj)ARG0(arg),prec,rp);
                    217: }
                    218:
                    219: /* set/get decimal precision */
                    220:
                    221: void Psetprec(NODE arg,Obj *rp)
                    222: {
                    223:   long p;
                    224:   Z q;
                    225:   long prec,dprec;
                    226:
                    227:   prec = mpfr_get_default_prec();
                    228:   /* decimal precision */
                    229:   dprec = prec*0.30103;
1.2       noro      230:   STOZ(dprec,q); *rp = (Obj)q;
1.1       noro      231:   if ( arg ) {
                    232:     asir_assert(ARG0(arg),O_N,"setprec");
1.2       noro      233:     p = ZTOS((Q)ARG0(arg))*3.32193;
1.1       noro      234:     if ( p > 0 )
                    235:       prec = p;
                    236:   }
                    237:   if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
                    238:   else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
                    239:   mpfr_set_default_prec(prec);
                    240: }
                    241:
                    242: /* set/get bit precision */
                    243:
                    244: void Psetbprec(NODE arg,Obj *rp)
                    245: {
                    246:   long p;
                    247:   Z q;
1.3     ! noro      248:   long prec,dprec;
1.1       noro      249:
                    250:   prec = mpfr_get_default_prec();
1.2       noro      251:   STOZ(prec,q); *rp = (Obj)q;
1.1       noro      252:   if ( arg ) {
                    253:     asir_assert(ARG0(arg),O_N,"setbprec");
1.2       noro      254:     p = ZTOS((Q)ARG0(arg));
1.1       noro      255:     if ( p > 0 )
                    256:       prec = p;
                    257:   }
                    258:   if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
                    259:   else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
                    260:   mpfr_set_default_prec(prec);
                    261: }
                    262:
                    263: void Psetround(NODE arg,Z *rp)
                    264: {
                    265:   int round;
                    266:
1.2       noro      267:   STOZ(mpfr_roundmode,*rp);
1.1       noro      268:   if ( arg ) {
                    269:     asir_assert(ARG0(arg),O_N,"setround");
1.2       noro      270:     round = ZTOS((Q)ARG0(arg));
1.1       noro      271:     switch ( round ) {
                    272:     case 0:
                    273:       mpfr_roundmode = MPFR_RNDN;
                    274:       break;
                    275:     case 1:
                    276:       mpfr_roundmode = MPFR_RNDZ;
                    277:       break;
                    278:     case 2:
                    279:       mpfr_roundmode = MPFR_RNDU;
                    280:       break;
                    281:     case 3:
                    282:       mpfr_roundmode = MPFR_RNDD;
                    283:       break;
                    284:     case 4:
                    285:       mpfr_roundmode = MPFR_RNDA;
                    286:       break;
                    287:     case 5:
                    288:       mpfr_roundmode = MPFR_RNDF;
                    289:       break;
                    290:     case 6:
                    291:       mpfr_roundmode = MPFR_RNDNA;
                    292:       break;
                    293:     default:
                    294:       error("setround : invalid rounding mode");
                    295:       break;
                    296:     }
                    297:   }
                    298: }
                    299:
                    300: Num tobf(Num a,int prec);
                    301:
                    302: void mp_pi(NODE arg,BF *rp)
                    303: {
                    304:   int prec;
                    305:   BF r;
                    306:
1.2       noro      307:   prec = arg ? ZTOS((Q)ARG0(arg)) : 0;
1.1       noro      308:   NEWBF(r);
                    309:   prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
                    310:   mpfr_const_pi(r->body,mpfr_roundmode);
                    311:   if ( !cmpbf((Num)r,0) ) r = 0;
                    312:   *rp = r;
                    313: }
                    314:
                    315: void mp_e(NODE arg,BF *rp)
                    316: {
                    317:   int prec;
                    318:   mpfr_t one;
                    319:   BF r;
                    320:
1.2       noro      321:   prec = arg ? ZTOS((Q)ARG0(arg)) : 0;
1.1       noro      322:   NEWBF(r);
                    323:   prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
                    324:   mpfr_init(one);
                    325:   mpfr_set_ui(one,1,mpfr_roundmode);
                    326:   mpfr_exp(r->body,one,mpfr_roundmode);
                    327:   if ( !cmpbf((Num)r,0) ) r = 0;
                    328:   *rp = r;
                    329: }
                    330:
                    331: void mpfr_or_mpc(NODE arg,int (*mpfr_f)(),int (*mpc_f)(),Num *rp)
                    332: {
                    333:   Num a;
                    334:   int prec;
                    335:   BF r,re,im;
                    336:   C c;
                    337:   mpc_t mpc,a1;
                    338:
1.2       noro      339:   prec = NEXT(arg) ? ZTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
1.1       noro      340:   a = tobf(ARG0(arg),prec);
                    341:   if ( a && NID(a)==N_C ) {
                    342:     mpc_init2(mpc,prec); mpc_init2(a1,prec);
                    343:     re = (BF)((C)a)->r; im = (BF)((C)a)->i;
                    344:     mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode);
                    345:     (*mpc_f)(mpc,a1,mpfr_roundmode);
                    346:     MPFRTOBF(mpc_realref(mpc),re);
                    347:     MPFRTOBF(mpc_imagref(mpc),im);
                    348:     if ( !cmpbf((Num)re,0) ) re = 0;
                    349:     if ( !cmpbf((Num)im,0) ) im = 0;
                    350:     if ( !im )
                    351:       *rp = (Num)re;
                    352:     else {
                    353:       NEWC(c); c->r = (Num)re; c->i = (Num)im;
                    354:       *rp = (Num)c;
                    355:     }
                    356:   } else {
                    357:     NEWBF(r);
                    358:     mpfr_init2(r->body,prec);
                    359:     (*mpfr_f)(r->body,((BF)a)->body,mpfr_roundmode);
                    360:     if ( !cmpbf((Num)r,0) ) r = 0;
                    361:     *rp = (Num)r;
                    362:   }
                    363: }
                    364:
                    365: void mp_sin(NODE arg,Num *rp)
                    366: {
                    367:   mpfr_or_mpc(arg,mpfr_sin,mpc_sin,rp);
                    368: }
                    369:
                    370: void mp_cos(NODE arg,Num *rp)
                    371: {
                    372:   mpfr_or_mpc(arg,mpfr_cos,mpc_cos,rp);
                    373: }
                    374:
                    375: void mp_tan(NODE arg,Num *rp)
                    376: {
                    377:   mpfr_or_mpc(arg,mpfr_tan,mpc_tan,rp);
                    378: }
                    379:
                    380: void mp_asin(NODE arg,Num *rp)
                    381: {
                    382:   mpfr_or_mpc(arg,mpfr_asin,mpc_asin,rp);
                    383: }
                    384:
                    385: void mp_acos(NODE arg,Num *rp)
                    386: {
                    387:   mpfr_or_mpc(arg,mpfr_acos,mpc_acos,rp);
                    388: }
                    389:
                    390: void mp_atan(NODE arg,Num *rp)
                    391: {
                    392:   mpfr_or_mpc(arg,mpfr_atan,mpc_atan,rp);
                    393: }
                    394:
                    395: void mp_sinh(NODE arg,Num *rp)
                    396: {
                    397:   mpfr_or_mpc(arg,mpfr_sinh,mpc_sinh,rp);
                    398: }
                    399:
                    400: void mp_cosh(NODE arg,Num *rp)
                    401: {
                    402:   mpfr_or_mpc(arg,mpfr_cosh,mpc_cosh,rp);
                    403: }
                    404:
                    405: void mp_tanh(NODE arg,Num *rp)
                    406: {
                    407:   mpfr_or_mpc(arg,mpfr_tanh,mpc_tanh,rp);
                    408: }
                    409:
                    410: void mp_asinh(NODE arg,Num *rp)
                    411: {
                    412:   mpfr_or_mpc(arg,mpfr_asinh,mpc_asinh,rp);
                    413: }
                    414:
                    415: void mp_acosh(NODE arg,Num *rp)
                    416: {
                    417:   mpfr_or_mpc(arg,mpfr_acosh,mpc_acosh,rp);
                    418: }
                    419:
                    420: void mp_atanh(NODE arg,Num *rp)
                    421: {
                    422:   mpfr_or_mpc(arg,mpfr_atanh,mpc_atanh,rp);
                    423: }
                    424:
                    425: void mp_exp(NODE arg,Num *rp)
                    426: {
                    427:   mpfr_or_mpc(arg,mpfr_exp,mpc_exp,rp);
                    428: }
                    429:
                    430: void mp_log(NODE arg,Num *rp)
                    431: {
                    432:   mpfr_or_mpc(arg,mpfr_log,mpc_log,rp);
                    433: }
                    434:
                    435: void mp_abs(NODE arg,Num *rp)
                    436: {
                    437:   mpfr_or_mpc(arg,mpfr_abs,mpc_abs,rp);
                    438: }
                    439:
                    440: void Pfac(NODE arg,Num *rp);
                    441:
                    442: void mp_factorial(NODE arg,Num *rp)
                    443: {
                    444:   struct oNODE arg0;
                    445:   Num a,a1;
                    446:
                    447:   a = (Num)ARG0(arg);
                    448:   if ( !a ) *rp = (Num)ONE;
                    449:   else if ( INT(a) ) Pfac(arg,rp);
                    450:   else {
                    451:     addnum(0,a,(Num)ONE,&a1);
                    452:     arg0.body = (pointer)a1;
                    453:     arg0.next = arg->next;
                    454:     Pmpfr_gamma(&arg0,rp);
                    455:   }
                    456: }
                    457:
                    458: void mp_pow(NODE arg,Num *rp)
                    459: {
                    460:   Num a,e;
                    461:   int prec;
                    462:   BF r,re,im;
                    463:   C c;
                    464:   mpc_t mpc,a1,e1;
                    465:
1.2       noro      466:   prec = NEXT(NEXT(arg)) ? ZTOS((Q)ARG2(arg)) : mpfr_get_default_prec();
1.1       noro      467:   a = tobf(ARG0(arg),prec);
                    468:   e = tobf(ARG1(arg),prec);
                    469:   if ( NID(a) == N_C || NID(e) == N_C || MPFR_SIGN(((BF)a)->body) < 0 ) {
                    470:     mpc_init2(mpc,prec); mpc_init2(a1,prec); mpc_init2(e1,prec);
                    471:     if ( NID(a) == N_C ) {
                    472:       re = (BF)((C)a)->r; im = (BF)((C)a)->i;
                    473:       mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode);
                    474:     } else {
                    475:       re = (BF)a;
                    476:       mpc_set_fr(a1,re->body,mpfr_roundmode);
                    477:     }
                    478:     if ( NID(e) == N_C ) {
                    479:       re = (BF)((C)e)->r; im = (BF)((C)e)->i;
                    480:       mpc_set_fr_fr(e1,re->body,im->body,mpfr_roundmode);
                    481:     } else {
                    482:       re = (BF)e;
                    483:       mpc_set_fr(e1,re->body,mpfr_roundmode);
                    484:     }
                    485:     mpc_pow(mpc,a1,e1,mpfr_roundmode);
                    486:     MPFRTOBF(mpc_realref(mpc),re);
                    487:     MPFRTOBF(mpc_imagref(mpc),im);
                    488:     if ( !cmpbf((Num)re,0) ) re = 0;
                    489:     if ( !cmpbf((Num)im,0) ) im = 0;
                    490:     if ( !im )
                    491:       *rp = (Num)re;
                    492:     else {
                    493:       NEWC(c); c->r = (Num)re; c->i = (Num)im;
                    494:       *rp = (Num)c;
                    495:     }
                    496:   } else {
                    497:     NEWBF(r);
                    498:     mpfr_init2(r->body,prec);
                    499:     mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
                    500:     *rp = (Num)r;
                    501:   }
                    502: }
                    503:
                    504: #define SETPREC \
1.2       noro      505:  (prec)=NEXT(arg)?ZTOS((Q)ARG1(arg)):0;\
1.1       noro      506:  (prec)*=3.32193;\
                    507:  (a)=tobf(ARG0(arg),prec);\
                    508:  NEWBF(r);\
                    509:  prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
                    510:
                    511:
                    512: void Pmpfr_gamma(NODE arg,BF *rp)
                    513: {
                    514:   Num a;
                    515:   int prec;
                    516:   BF r;
                    517:
                    518:   SETPREC
                    519:   mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
                    520:   *rp = r;
                    521: }
                    522:
                    523: void Pmpfr_lngamma(NODE arg,BF *rp)
                    524: {
                    525:   Num a;
                    526:   int prec;
                    527:   BF r;
                    528:
                    529:   SETPREC
                    530:   mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode);
                    531:   *rp = r;
                    532: }
                    533:
                    534: void Pmpfr_digamma(NODE arg,BF *rp)
                    535: {
                    536:   Num a;
                    537:   int prec;
                    538:   BF r;
                    539:
                    540:   SETPREC
                    541:   mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode);
                    542:   *rp = r;
                    543: }
                    544:
                    545: void Pmpfr_zeta(NODE arg,BF *rp)
                    546: {
                    547:   Num a;
                    548:   int prec;
                    549:   BF r;
                    550:
                    551:   SETPREC
                    552:   mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode);
                    553:   *rp = r;
                    554: }
                    555:
                    556: void Pmpfr_eint(NODE arg,BF *rp)
                    557: {
                    558:   Num a;
                    559:   int prec;
                    560:   BF r;
                    561:
                    562:   SETPREC
                    563:   mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode);
                    564:   *rp = r;
                    565: }
                    566:
                    567: void Pmpfr_erf(NODE arg,BF *rp)
                    568: {
                    569:   Num a;
                    570:   int prec;
                    571:   BF r;
                    572:
                    573:   SETPREC
                    574:   mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode);
                    575:   *rp = r;
                    576: }
                    577:
                    578: void Pmpfr_erfc(NODE arg,BF *rp)
                    579: {
                    580:   Num a;
                    581:   int prec;
                    582:   BF r;
                    583:
                    584:   SETPREC
                    585:   mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode);
                    586:   *rp = r;
                    587: }
                    588:
                    589: void Pmpfr_j0(NODE arg,BF *rp)
                    590: {
                    591:   Num a;
                    592:   int prec;
                    593:   BF r;
                    594:
                    595:   SETPREC
                    596:   mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode);
                    597:   *rp = r;
                    598: }
                    599:
                    600: void Pmpfr_j1(NODE arg,BF *rp)
                    601: {
                    602:   Num a;
                    603:   int prec;
                    604:   BF r;
                    605:
                    606:   SETPREC
                    607:   mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode);
                    608:   *rp = r;
                    609: }
                    610:
                    611: void Pmpfr_y0(NODE arg,BF *rp)
                    612: {
                    613:   Num a;
                    614:   int prec;
                    615:   BF r;
                    616:
                    617:   SETPREC
                    618:   mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode);
                    619:   *rp = r;
                    620: }
                    621:
                    622: void Pmpfr_y1(NODE arg,BF *rp)
                    623: {
                    624:   Num a;
                    625:   int prec;
                    626:   BF r;
                    627:
                    628:   SETPREC
                    629:   mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode);
                    630:   *rp = r;
                    631: }
                    632:
                    633: void Pmpfr_li2(NODE arg,BF *rp)
                    634: {
                    635:   Num a;
                    636:   int prec;
                    637:   BF r;
                    638:
                    639:   SETPREC
                    640:   mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode);
                    641:   *rp = r;
                    642: }
                    643:
                    644: void Pmpfr_ai(NODE arg,BF *rp)
                    645: {
                    646:   Num a;
                    647:   int prec;
                    648:   BF r;
                    649:
                    650:   SETPREC
                    651:   mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode);
                    652:   *rp = r;
                    653: }
                    654:
                    655: void Pmpfr_floor(NODE arg,Z *rp)
                    656: {
                    657:   Num a;
                    658:   long prec;
                    659:   BF r;
                    660:   mpz_t t;
                    661:
                    662:   SETPREC
                    663:   mpfr_floor(r->body,((BF)a)->body);
                    664:   mpz_init(t);
                    665:   mpfr_get_z(t,r->body,mpfr_roundmode);
                    666:   MPZTOZ(t,*rp);
                    667: }
                    668:
                    669: void Pmpfr_ceil(NODE arg,Z *rp)
                    670: {
                    671:   Num a;
                    672:   long prec;
                    673:   BF r;
                    674:   mpz_t t;
                    675:   Z rz;
                    676:
                    677:   SETPREC
                    678:   mpfr_ceil(r->body,((BF)a)->body);
                    679:   mpz_init(t);
                    680:   mpfr_get_z(t,r->body,mpfr_roundmode);
                    681:   MPZTOZ(t,*rp);
                    682: }
                    683:
                    684: void Pmpfr_round(NODE arg,Z *rp)
                    685: {
                    686:   Num a;
                    687:   long prec;
                    688:   BF r;
                    689:   mpz_t t;
                    690:
                    691:   SETPREC
                    692:   mpfr_round(r->body,((BF)a)->body);
                    693:   mpz_init(t);
                    694:   mpfr_get_z(t,r->body,mpfr_roundmode);
                    695:   MPZTOZ(t,*rp);
                    696: }
                    697:
                    698: double **almat_double(int n)
                    699: {
                    700:   int i;
                    701:   double **a;
                    702:
                    703:   a = (double **)MALLOC(n*sizeof(double *));
                    704:   for ( i = 0; i < n; i++ )
                    705:     a[i] = (double *)MALLOC(n*sizeof(double));
                    706:   return a;
                    707: }
                    708:
                    709: /*
                    710:  *  k <- (A(xi)-(sbeta-mn2/xi))f
                    711:  *  A(t) = (num[0]+num[1]t+...+num[d-1]*t^(d-1))/den(t)
                    712:  */
                    713:
                    714: struct jv {
                    715:   int j;
                    716:   double v;
                    717: };
                    718:
                    719: struct smat {
                    720:   int *rlen;
                    721:   struct jv **row;
                    722: };
                    723:
                    724: void eval_pfaffian2(double *k,int n,int d,struct smat *num,P den,double xi,double *f)
                    725: {
                    726:   struct smat ma;
                    727:   struct jv *maj;
                    728:   int i,j,l,s;
                    729:   double t,dn;
                    730:   P r;
                    731:   Real u;
                    732:
                    733:   memset(k,0,n*sizeof(double));
                    734:   for ( i = d-1; i >= 0; i-- ) {
                    735:     ma = num[i];
                    736:     for ( j = 0; j < n; j++ ) {
                    737:       maj = ma.row[j];
                    738:       l = ma.rlen[j];
                    739:       for ( t = 0, s = 0; s < l; s++, maj++ ) t += maj->v*f[maj->j];
                    740:       k[j] = k[j]*xi+t;
                    741:     }
                    742:   }
                    743:   MKReal(xi,u);
                    744:   substp(CO,den,den->v,(P)u,&r); dn = ToReal(r);
                    745:   for ( j = 0; j < n; j++ )
                    746:     k[j] /= dn;
                    747: }
                    748:
                    749: void Prk_ratmat(NODE arg,LIST *rp)
                    750: {
                    751:   VECT mat;
                    752:   P den;
                    753:   int ord;
                    754:   double sbeta,x0,x1,xi,h,mn2,hd;
                    755:   double a2,a3,a4,a5,a6;
                    756:   double b21,b31,b32,b41,b42,b43,b51,b52,b53,b54,b61,b62,b63,b64,b65;
                    757:   double c1,c2,c3,c4,c5,c6,c7;
                    758:   VECT fv;
                    759:   int step,j,i,k,n,d,len,s;
                    760:   struct smat *num;
                    761:   Obj **b;
                    762:   MAT mati;
                    763:   double *f,*w,*k1,*k2,*k3,*k4,*k5,*k6;
                    764:   NODE nd,nd1;
                    765:   Real x,t;
                    766:   LIST l;
                    767:
1.2       noro      768:   ord = ZTOS((Q)ARG0(arg));
1.1       noro      769:   mat = (VECT)ARG1(arg); den = (P)ARG2(arg);
                    770:   x0 = ToReal((Num)ARG3(arg)); x1 = ToReal((Num)ARG4(arg));
1.2       noro      771:   step = ZTOS((Q)ARG5(arg)); fv = (VECT)ARG6(arg);
1.1       noro      772:   h = (x1-x0)/step;
                    773:
                    774:   n = fv->len;
                    775:   d = mat->len;
                    776:   num = (struct smat *)MALLOC(d*sizeof(struct smat));
                    777:   for ( i = 0; i < d; i++ ) {
                    778:     num[i].rlen = (int *)MALLOC_ATOMIC(n*sizeof(int));
                    779:     num[i].row = (struct jv **)MALLOC(n*sizeof(struct jv *));
                    780:     mati = (MAT)mat->body[i];
                    781:     b = (Obj **)mati->body;
                    782:     for ( j = 0; j < n; j++ ) {
                    783:       for ( len = k = 0; k < n; k++ )
                    784:         if ( b[j][k] ) len++;
                    785:       num[i].rlen[j] = len;
                    786:       if ( !len )
                    787:         num[i].row[j] = 0;
                    788:       else {
                    789:         num[i].row[j] = (struct jv *)MALLOC_ATOMIC((len)*sizeof(struct jv));
                    790:         for ( s = k = 0; k < n; k++ )
                    791:           if ( b[j][k] ) {
                    792:              num[i].row[j][s].j = k;
                    793:              num[i].row[j][s].v = ToReal((Num)b[j][k]);
                    794:              s++;
                    795:           }
                    796:       }
                    797:     }
                    798:   }
                    799:   f = (double *)MALLOC_ATOMIC(n*sizeof(double));
                    800:   for ( j = 0; j < n; j++ )
                    801:     f[j] = ToReal((Num)fv->body[j]);
                    802:   w = (double *)MALLOC_ATOMIC(n*sizeof(double));
                    803:   k1 = (double *)MALLOC_ATOMIC(n*sizeof(double));
                    804:   k2 = (double *)MALLOC_ATOMIC(n*sizeof(double));
                    805:   k3 = (double *)MALLOC_ATOMIC(n*sizeof(double));
                    806:   k4 = (double *)MALLOC_ATOMIC(n*sizeof(double));
                    807:   k5 = (double *)MALLOC_ATOMIC(n*sizeof(double));
                    808:   k6 = (double *)MALLOC(n*sizeof(double));
                    809:   nd = 0;
                    810:   switch ( ord ) {
                    811:   case 4:
                    812:     a2 = 1/2.0*h; b21 = 1/2.0*h;
                    813:     a3 = 1/2.0*h; b31 = 0.0;   b32 = 1/2.0*h;
                    814:     a4 = 1.0*h;   b41 = 0.0;   b42 = 0.0;    b43 = 1.0*h;
                    815:     c1 = 1/6.0*h; c2 = 1/3.0*h;     c3 =  1/3.0*h; c4 = 1/6.0*h;
                    816:     for ( i = 0; i < step; i++ ) {
                    817:       if ( !(i%100000) ) fprintf(stderr,"[%d]",i);
                    818:       xi = x0+i*h;
                    819:       eval_pfaffian2(k1,n,d,num,den,xi,f);
                    820:       memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b21*k1[j];
                    821:       eval_pfaffian2(k2,n,d,num,den,xi+a2,w);
                    822:       memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b31*k1[j]+b32*k2[j];
                    823:       eval_pfaffian2(k3,n,d,num,den,xi+a3,w);
                    824:       memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b41*k1[j]+b42*k2[j]+b43*k3[j];
                    825:       eval_pfaffian2(k4,n,d,num,den,xi+a4,w);
                    826:       memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += c1*k1[j]+c2*k2[j]+c3*k3[j]+c4*k4[j];
                    827:       memcpy(f,w,n*sizeof(double));
                    828:       MKReal(f[0],t);
                    829:       MKReal(xi+h,x);
                    830:       nd1 = mknode(2,x,t);
                    831:       MKLIST(l,nd1);
                    832:       MKNODE(nd1,l,nd);
                    833:       nd = nd1;
                    834:       for ( hd = f[0], j = 0; j < n; j++ ) f[j] /= hd;
                    835:     }
                    836:     MKLIST(*rp,nd);
                    837:     break;
                    838:   case 5:
                    839:   default:
                    840:     a2 = 1/4.0*h; b21 = 1/4.0*h;
                    841:     a3 = 1/4.0*h; b31 = 1/8.0*h; b32 = 1/8.0*h;
                    842:     a4 = 1/2.0*h; b41 = 0.0;   b42 = 0.0;    b43 = 1/2.0*h;
                    843:     a5 = 3/4.0*h; b51 = 3/16.0*h;b52 = -3/8.0*h; b53 = 3/8.0*h;   b54 = 9/16.0*h;
                    844:     a6 = 1.0*h;   b61 = -3/7.0*h;b62 = 8/7.0*h;  b63 = 6/7.0*h;   b64 = -12/7.0*h; b65 = 8/7.0*h;
                    845:     c1 = 7/90.0*h; c2 = 0.0;     c3 =  16/45.0*h; c4 = 2/15.0*h;   c5 = 16/45.0*h; c6 = 7/90.0*h;
                    846:     for ( i = 0; i < step; i++ ) {
                    847:       if ( !(i%100000) ) fprintf(stderr,"[%d]",i);
                    848:       xi = x0+i*h;
                    849:       eval_pfaffian2(k1,n,d,num,den,xi,f);
                    850:       memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b21*k1[j];
                    851:       eval_pfaffian2(k2,n,d,num,den,xi+a2,w);
                    852:       memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b31*k1[j]+b32*k2[j];
                    853:       eval_pfaffian2(k3,n,d,num,den,xi+a3,w);
                    854:       memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b41*k1[j]+b42*k2[j]+b43*k3[j];
                    855:       eval_pfaffian2(k4,n,d,num,den,xi+a4,w);
                    856:       memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b51*k1[j]+b52*k2[j]+b53*k3[j]+b54*k4[j];
                    857:       eval_pfaffian2(k5,n,d,num,den,xi+a5,w);
                    858:       memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b61*k1[j]+b62*k2[j]+b63*k3[j]+b64*k4[j]+b65*k5[j];
                    859:       eval_pfaffian2(k6,n,d,num,den,xi+a6,w);
                    860:       memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += c1*k1[j]+c2*k2[j]+c3*k3[j]+c4*k4[j]+c5*k5[j]+c6*k6[j];
                    861:       memcpy(f,w,n*sizeof(double));
                    862:       MKReal(f[0],t);
                    863:       MKReal(xi+h,x);
                    864:       nd1 = mknode(2,x,t);
                    865:       MKLIST(l,nd1);
                    866:       MKNODE(nd1,l,nd);
                    867:       nd = nd1;
                    868:       for ( hd = f[0], j = 0; j < n; j++ ) f[j] /= hd;
                    869:     }
                    870:     MKLIST(*rp,nd);
                    871:     break;
                    872:   }
                    873: }

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