Annotation of OpenXM_contrib2/asir2000/builtin/bfaux.c, Revision 1.17
1.17 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.16 2018/03/28 05:27:22 noro Exp $ */
1.1 noro 2: #include "ca.h"
3: #include "parse.h"
4:
1.4 noro 5: void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround();
1.8 noro 6: void Pmpfr_ai();
1.9 noro 7: void Pmpfr_eint(), Pmpfr_erf(), Pmpfr_erfc(), Pmpfr_li2();
1.8 noro 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();
1.12 noro 13: void Prk_ratmat();
1.14 noro 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();
1.16 noro 17: void mp_factorial(),mp_abs();
1.1 noro 18:
19: struct ftab bf_tab[] = {
1.17 ! noro 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},
1.1 noro 58: };
59:
1.4 noro 60: int mpfr_roundmode = MPFR_RNDN;
61:
1.15 noro 62: void todoublen(Num a,Num *rp)
1.1 noro 63: {
1.17 ! noro 64: double r,i;
! 65: Real real,imag;
1.1 noro 66:
1.17 ! noro 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: }
1.1 noro 88: }
89:
1.15 noro 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: {
1.17 ! noro 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:
1.15 noro 140: todoublen((Num)a,(Num *)b);
141: break;
142: case O_P:
143: todoublep((P)a,(P *)b);
144: break;
145: case O_R:
1.17 ! noro 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: }
1.15 noro 198: }
199:
200: void Ptodouble(NODE arg,Obj *rp)
201: {
202: todouble((Obj)ARG0(arg),rp);
203: }
204:
1.4 noro 205: void Peval(NODE arg,Obj *rp)
1.1 noro 206: {
207: int prec;
208:
1.17 ! noro 209: asir_assert(ARG0(arg),O_R,"eval");
1.1 noro 210: if ( argc(arg) == 2 ) {
1.17 ! noro 211: prec = QTOS((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;
1.17 ! noro 216: evalr(CO,(Obj)ARG0(arg),prec,rp);
1.1 noro 217: }
218:
1.3 noro 219: /* set/get decimal precision */
1.1 noro 220:
221: void Psetprec(NODE arg,Obj *rp)
222: {
1.17 ! noro 223: int p;
! 224: Q q;
! 225: int prec,dprec;
1.3 noro 226:
227: prec = mpfr_get_default_prec();
228: /* decimal precision */
229: dprec = prec*0.30103;
230: STOQ(dprec,q); *rp = (Obj)q;
1.17 ! noro 231: if ( arg ) {
! 232: asir_assert(ARG0(arg),O_N,"setprec");
! 233: p = QTOS((Q)ARG0(arg))*3.32193;
! 234: if ( p > 0 )
! 235: prec = p;
! 236: }
1.3 noro 237: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
238: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
1.17 ! noro 239: mpfr_set_default_prec(prec);
1.3 noro 240: }
1.1 noro 241:
1.3 noro 242: /* set/get bit precision */
1.1 noro 243:
1.3 noro 244: void Psetbprec(NODE arg,Obj *rp)
245: {
1.17 ! noro 246: int p;
! 247: Q q;
! 248: int prec;
1.3 noro 249:
250: prec = mpfr_get_default_prec();
251: STOQ(prec,q); *rp = (Obj)q;
1.17 ! noro 252: if ( arg ) {
! 253: asir_assert(ARG0(arg),O_N,"setbprec");
! 254: p = QTOS((Q)ARG0(arg));
! 255: if ( p > 0 )
! 256: prec = p;
! 257: }
1.1 noro 258: if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
259: else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
1.17 ! noro 260: mpfr_set_default_prec(prec);
1.1 noro 261: }
262:
1.4 noro 263: void Psetround(NODE arg,Q *rp)
264: {
265: int round;
266:
267: STOQ(mpfr_roundmode,*rp);
268: if ( arg ) {
1.17 ! noro 269: asir_assert(ARG0(arg),O_N,"setround");
1.4 noro 270: round = QTOS((Q)ARG0(arg));
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:
1.1 noro 300: Num tobf(Num a,int prec);
301:
302: void mp_pi(NODE arg,BF *rp)
303: {
1.10 noro 304: int prec;
1.17 ! noro 305: BF r;
1.1 noro 306:
1.17 ! noro 307: prec = arg ? QTOS((Q)ARG0(arg)) : 0;
! 308: NEWBF(r);
! 309: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
! 310: mpfr_const_pi(r->body,mpfr_roundmode);
1.10 noro 311: if ( !cmpbf((Num)r,0) ) r = 0;
312: *rp = r;
1.1 noro 313: }
314:
315: void mp_e(NODE arg,BF *rp)
316: {
1.10 noro 317: int prec;
1.17 ! noro 318: mpfr_t one;
! 319: BF r;
1.1 noro 320:
1.17 ! noro 321: prec = arg ? QTOS((Q)ARG0(arg)) : 0;
! 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);
1.10 noro 327: if ( !cmpbf((Num)r,0) ) r = 0;
328: *rp = r;
1.1 noro 329: }
330:
1.10 noro 331: void mpfr_or_mpc(NODE arg,int (*mpfr_f)(),int (*mpc_f)(),Num *rp)
1.1 noro 332: {
1.17 ! noro 333: Num a;
1.10 noro 334: int prec;
1.17 ! noro 335: BF r,re,im;
1.10 noro 336: C c;
337: mpc_t mpc,a1;
1.1 noro 338:
1.17 ! noro 339: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
! 340: a = tobf(ARG0(arg),prec);
1.10 noro 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 {
1.17 ! noro 357: NEWBF(r);
! 358: mpfr_init2(r->body,prec);
! 359: (*mpfr_f)(r->body,((BF)a)->body,mpfr_roundmode);
1.10 noro 360: if ( !cmpbf((Num)r,0) ) r = 0;
361: *rp = (Num)r;
362: }
1.1 noro 363: }
364:
1.10 noro 365: void mp_sin(NODE arg,Num *rp)
1.1 noro 366: {
1.10 noro 367: mpfr_or_mpc(arg,mpfr_sin,mpc_sin,rp);
368: }
1.1 noro 369:
1.10 noro 370: void mp_cos(NODE arg,Num *rp)
371: {
372: mpfr_or_mpc(arg,mpfr_cos,mpc_cos,rp);
1.1 noro 373: }
374:
1.10 noro 375: void mp_tan(NODE arg,Num *rp)
1.1 noro 376: {
1.10 noro 377: mpfr_or_mpc(arg,mpfr_tan,mpc_tan,rp);
378: }
1.1 noro 379:
1.10 noro 380: void mp_asin(NODE arg,Num *rp)
381: {
382: mpfr_or_mpc(arg,mpfr_asin,mpc_asin,rp);
1.1 noro 383: }
384:
1.10 noro 385: void mp_acos(NODE arg,Num *rp)
1.1 noro 386: {
1.10 noro 387: mpfr_or_mpc(arg,mpfr_acos,mpc_acos,rp);
388: }
1.1 noro 389:
1.10 noro 390: void mp_atan(NODE arg,Num *rp)
391: {
392: mpfr_or_mpc(arg,mpfr_atan,mpc_atan,rp);
1.1 noro 393: }
394:
1.10 noro 395: void mp_sinh(NODE arg,Num *rp)
1.1 noro 396: {
1.10 noro 397: mpfr_or_mpc(arg,mpfr_sinh,mpc_sinh,rp);
1.1 noro 398: }
399:
1.10 noro 400: void mp_cosh(NODE arg,Num *rp)
1.1 noro 401: {
1.10 noro 402: mpfr_or_mpc(arg,mpfr_cosh,mpc_cosh,rp);
1.1 noro 403: }
404:
1.10 noro 405: void mp_tanh(NODE arg,Num *rp)
1.1 noro 406: {
1.10 noro 407: mpfr_or_mpc(arg,mpfr_tanh,mpc_tanh,rp);
1.1 noro 408: }
409:
1.10 noro 410: void mp_asinh(NODE arg,Num *rp)
1.1 noro 411: {
1.10 noro 412: mpfr_or_mpc(arg,mpfr_asinh,mpc_asinh,rp);
1.1 noro 413: }
414:
1.10 noro 415: void mp_acosh(NODE arg,Num *rp)
1.1 noro 416: {
1.10 noro 417: mpfr_or_mpc(arg,mpfr_acosh,mpc_acosh,rp);
1.1 noro 418: }
419:
1.10 noro 420: void mp_atanh(NODE arg,Num *rp)
1.1 noro 421: {
1.10 noro 422: mpfr_or_mpc(arg,mpfr_atanh,mpc_atanh,rp);
1.1 noro 423: }
424:
1.10 noro 425: void mp_exp(NODE arg,Num *rp)
1.1 noro 426: {
1.10 noro 427: mpfr_or_mpc(arg,mpfr_exp,mpc_exp,rp);
1.1 noro 428: }
429:
1.10 noro 430: void mp_log(NODE arg,Num *rp)
1.1 noro 431: {
1.10 noro 432: mpfr_or_mpc(arg,mpfr_log,mpc_log,rp);
1.1 noro 433: }
434:
1.16 noro 435: void mp_abs(NODE arg,Num *rp)
436: {
437: mpfr_or_mpc(arg,mpfr_abs,mpc_abs,rp);
438: }
439:
440: void mp_factorial(NODE arg,Num *rp)
441: {
442: struct oNODE arg0;
443: Num a,a1;
444:
445: a = (Num)ARG0(arg);
446: if ( !a ) *rp = (Num)ONE;
447: else if ( INT(a) ) Pfac(arg,rp);
448: else {
449: addnum(0,a,(Num)ONE,&a1);
450: arg0.body = (pointer)a1;
451: arg0.next = arg->next;
452: Pmpfr_gamma(&arg0,rp);
453: }
454: }
455:
1.10 noro 456: void mp_pow(NODE arg,Num *rp)
1.1 noro 457: {
1.17 ! noro 458: Num a,e;
1.10 noro 459: int prec;
1.17 ! noro 460: BF r,re,im;
1.10 noro 461: C c;
462: mpc_t mpc,a1,e1;
1.1 noro 463:
1.17 ! noro 464: prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : mpfr_get_default_prec();
! 465: a = tobf(ARG0(arg),prec);
! 466: e = tobf(ARG1(arg),prec);
1.10 noro 467: if ( NID(a) == N_C || NID(e) == N_C || MPFR_SIGN(((BF)a)->body) < 0 ) {
468: mpc_init2(mpc,prec); mpc_init2(a1,prec); mpc_init2(e1,prec);
469: if ( NID(a) == N_C ) {
470: re = (BF)((C)a)->r; im = (BF)((C)a)->i;
471: mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode);
472: } else {
473: re = (BF)a;
474: mpc_set_fr(a1,re->body,mpfr_roundmode);
475: }
476: if ( NID(e) == N_C ) {
477: re = (BF)((C)e)->r; im = (BF)((C)e)->i;
478: mpc_set_fr_fr(e1,re->body,im->body,mpfr_roundmode);
479: } else {
480: re = (BF)e;
481: mpc_set_fr(e1,re->body,mpfr_roundmode);
482: }
483: mpc_pow(mpc,a1,e1,mpfr_roundmode);
484: MPFRTOBF(mpc_realref(mpc),re);
485: MPFRTOBF(mpc_imagref(mpc),im);
486: if ( !cmpbf((Num)re,0) ) re = 0;
487: if ( !cmpbf((Num)im,0) ) im = 0;
488: if ( !im )
489: *rp = (Num)re;
490: else {
491: NEWC(c); c->r = (Num)re; c->i = (Num)im;
492: *rp = (Num)c;
493: }
494: } else {
1.17 ! noro 495: NEWBF(r);
! 496: mpfr_init2(r->body,prec);
! 497: mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
1.10 noro 498: *rp = (Num)r;
499: }
1.1 noro 500: }
1.5 noro 501:
1.8 noro 502: #define SETPREC \
503: (prec)=NEXT(arg)?QTOS((Q)ARG1(arg)):0;\
504: (prec)*=3.32193;\
505: (a)=tobf(ARG0(arg),prec);\
506: NEWBF(r);\
507: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
508:
509:
1.5 noro 510: void Pmpfr_gamma(NODE arg,BF *rp)
511: {
1.17 ! noro 512: Num a;
1.5 noro 513: int prec;
1.17 ! noro 514: BF r;
1.5 noro 515:
1.8 noro 516: SETPREC
1.17 ! noro 517: mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
1.5 noro 518: *rp = r;
519: }
1.6 takayama 520:
1.8 noro 521: void Pmpfr_lngamma(NODE arg,BF *rp)
522: {
1.17 ! noro 523: Num a;
1.8 noro 524: int prec;
1.17 ! noro 525: BF r;
1.8 noro 526:
527: SETPREC
1.17 ! noro 528: mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 529: *rp = r;
530: }
531:
532: void Pmpfr_digamma(NODE arg,BF *rp)
533: {
1.17 ! noro 534: Num a;
1.8 noro 535: int prec;
1.17 ! noro 536: BF r;
1.8 noro 537:
538: SETPREC
1.17 ! noro 539: mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 540: *rp = r;
541: }
542:
543: void Pmpfr_zeta(NODE arg,BF *rp)
544: {
1.17 ! noro 545: Num a;
1.8 noro 546: int prec;
1.17 ! noro 547: BF r;
1.8 noro 548:
549: SETPREC
1.17 ! noro 550: mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 551: *rp = r;
552: }
553:
554: void Pmpfr_eint(NODE arg,BF *rp)
555: {
1.17 ! noro 556: Num a;
1.8 noro 557: int prec;
1.17 ! noro 558: BF r;
1.8 noro 559:
560: SETPREC
1.17 ! noro 561: mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 562: *rp = r;
563: }
564:
565: void Pmpfr_erf(NODE arg,BF *rp)
566: {
1.17 ! noro 567: Num a;
1.8 noro 568: int prec;
1.17 ! noro 569: BF r;
1.8 noro 570:
571: SETPREC
1.17 ! noro 572: mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 573: *rp = r;
574: }
575:
1.9 noro 576: void Pmpfr_erfc(NODE arg,BF *rp)
577: {
1.17 ! noro 578: Num a;
1.9 noro 579: int prec;
1.17 ! noro 580: BF r;
1.9 noro 581:
582: SETPREC
1.17 ! noro 583: mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode);
1.9 noro 584: *rp = r;
585: }
586:
1.8 noro 587: void Pmpfr_j0(NODE arg,BF *rp)
588: {
1.17 ! noro 589: Num a;
1.8 noro 590: int prec;
1.17 ! noro 591: BF r;
1.8 noro 592:
593: SETPREC
1.17 ! noro 594: mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 595: *rp = r;
596: }
597:
598: void Pmpfr_j1(NODE arg,BF *rp)
599: {
1.17 ! noro 600: Num a;
1.8 noro 601: int prec;
1.17 ! noro 602: BF r;
1.8 noro 603:
604: SETPREC
1.17 ! noro 605: mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 606: *rp = r;
607: }
608:
609: void Pmpfr_y0(NODE arg,BF *rp)
610: {
1.17 ! noro 611: Num a;
1.8 noro 612: int prec;
1.17 ! noro 613: BF r;
1.8 noro 614:
615: SETPREC
1.17 ! noro 616: mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 617: *rp = r;
618: }
619:
620: void Pmpfr_y1(NODE arg,BF *rp)
621: {
1.17 ! noro 622: Num a;
1.8 noro 623: int prec;
1.17 ! noro 624: BF r;
1.8 noro 625:
626: SETPREC
1.17 ! noro 627: mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 628: *rp = r;
629: }
630:
631: void Pmpfr_li2(NODE arg,BF *rp)
632: {
1.17 ! noro 633: Num a;
1.8 noro 634: int prec;
1.17 ! noro 635: BF r;
1.8 noro 636:
637: SETPREC
1.17 ! noro 638: mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 639: *rp = r;
640: }
641:
642: void Pmpfr_ai(NODE arg,BF *rp)
643: {
1.17 ! noro 644: Num a;
1.8 noro 645: int prec;
1.17 ! noro 646: BF r;
1.8 noro 647:
648: SETPREC
1.17 ! noro 649: mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode);
1.8 noro 650: *rp = r;
651: }
652:
1.7 takayama 653: void Pmpfr_floor(NODE arg,Q *rp)
1.6 takayama 654: {
1.17 ! noro 655: Num a;
1.6 takayama 656: int prec;
1.17 ! noro 657: BF r;
! 658: mpz_t t;
! 659: GZ rz;
1.6 takayama 660:
1.8 noro 661: SETPREC
1.17 ! noro 662: mpfr_floor(r->body,((BF)a)->body);
! 663: mpz_init(t);
! 664: mpfr_get_z(t,r->body,mpfr_roundmode);
! 665: MPZTOGZ(t,rz);
! 666: *rp = gztoz(rz);
1.7 takayama 667: }
668:
1.8 noro 669: void Pmpfr_ceil(NODE arg,Q *rp)
670: {
1.17 ! noro 671: Num a;
1.8 noro 672: int prec;
1.17 ! noro 673: BF r;
! 674: mpz_t t;
! 675: GZ rz;
1.8 noro 676:
677: SETPREC
1.17 ! noro 678: mpfr_ceil(r->body,((BF)a)->body);
! 679: mpz_init(t);
! 680: mpfr_get_z(t,r->body,mpfr_roundmode);
! 681: MPZTOGZ(t,rz);
! 682: *rp = gztoz(rz);
1.8 noro 683: }
684:
1.7 takayama 685: void Pmpfr_round(NODE arg,Q *rp)
686: {
1.17 ! noro 687: Num a;
1.7 takayama 688: int prec;
1.17 ! noro 689: BF r;
! 690: mpz_t t;
! 691: GZ rz;
1.7 takayama 692:
1.8 noro 693: SETPREC
1.17 ! noro 694: mpfr_round(r->body,((BF)a)->body);
! 695: mpz_init(t);
! 696: mpfr_get_z(t,r->body,mpfr_roundmode);
! 697: MPZTOGZ(t,rz);
! 698: *rp = gztoz(rz);
1.6 takayama 699: }
1.12 noro 700:
701: double **almat_double(int n)
702: {
703: int i;
704: double **a;
705:
706: a = (double **)MALLOC(n*sizeof(double *));
707: for ( i = 0; i < n; i++ )
708: a[i] = (double *)MALLOC(n*sizeof(double));
709: return a;
710: }
711:
712: /*
713: * k <- (A(xi)-(sbeta-mn2/xi))f
714: * A(t) = (num[0]+num[1]t+...+num[d-1]*t^(d-1))/den(t)
715: */
716:
717: struct jv {
718: int j;
719: double v;
720: };
721:
722: struct smat {
723: int *rlen;
724: struct jv **row;
725: };
726:
727: void eval_pfaffian2(double *k,int n,int d,struct smat *num,P den,double xi,double *f)
728: {
729: struct smat ma;
730: struct jv *maj;
731: int i,j,l,s;
732: double t,dn;
733: P r;
734: Real u;
735:
736: memset(k,0,n*sizeof(double));
737: for ( i = d-1; i >= 0; i-- ) {
738: ma = num[i];
739: for ( j = 0; j < n; j++ ) {
740: maj = ma.row[j];
741: l = ma.rlen[j];
742: for ( t = 0, s = 0; s < l; s++, maj++ ) t += maj->v*f[maj->j];
743: k[j] = k[j]*xi+t;
744: }
745: }
746: MKReal(xi,u);
747: substp(CO,den,den->v,(P)u,&r); dn = ToReal(r);
748: for ( j = 0; j < n; j++ )
749: k[j] /= dn;
750: }
751:
752: void Prk_ratmat(NODE arg,LIST *rp)
753: {
754: VECT mat;
755: P den;
756: int ord;
757: double sbeta,x0,x1,xi,h,mn2,hd;
758: double a2,a3,a4,a5,a6;
759: double b21,b31,b32,b41,b42,b43,b51,b52,b53,b54,b61,b62,b63,b64,b65;
760: double c1,c2,c3,c4,c5,c6,c7;
761: VECT fv;
762: int step,j,i,k,n,d,len,s;
763: struct smat *num;
764: Obj **b;
765: MAT mati;
766: double *f,*w,*k1,*k2,*k3,*k4,*k5,*k6;
767: NODE nd,nd1;
768: Real x,t;
769: LIST l;
770:
771: ord = QTOS((Q)ARG0(arg));
772: mat = (VECT)ARG1(arg); den = (P)ARG2(arg);
773: x0 = ToReal((Num)ARG3(arg)); x1 = ToReal((Num)ARG4(arg));
774: step = QTOS((Q)ARG5(arg)); fv = (VECT)ARG6(arg);
775: h = (x1-x0)/step;
776:
777: n = fv->len;
778: d = mat->len;
1.13 noro 779: num = (struct smat *)MALLOC(d*sizeof(struct smat));
1.12 noro 780: for ( i = 0; i < d; i++ ) {
1.13 noro 781: num[i].rlen = (int *)MALLOC_ATOMIC(n*sizeof(int));
1.12 noro 782: num[i].row = (struct jv **)MALLOC(n*sizeof(struct jv *));
783: mati = (MAT)mat->body[i];
784: b = (Obj **)mati->body;
785: for ( j = 0; j < n; j++ ) {
786: for ( len = k = 0; k < n; k++ )
787: if ( b[j][k] ) len++;
788: num[i].rlen[j] = len;
1.13 noro 789: if ( !len )
790: num[i].row[j] = 0;
791: else {
792: num[i].row[j] = (struct jv *)MALLOC_ATOMIC((len)*sizeof(struct jv));
793: for ( s = k = 0; k < n; k++ )
794: if ( b[j][k] ) {
795: num[i].row[j][s].j = k;
796: num[i].row[j][s].v = ToReal((Num)b[j][k]);
797: s++;
798: }
799: }
1.12 noro 800: }
801: }
1.13 noro 802: f = (double *)MALLOC_ATOMIC(n*sizeof(double));
1.12 noro 803: for ( j = 0; j < n; j++ )
804: f[j] = ToReal((Num)fv->body[j]);
1.13 noro 805: w = (double *)MALLOC_ATOMIC(n*sizeof(double));
806: k1 = (double *)MALLOC_ATOMIC(n*sizeof(double));
807: k2 = (double *)MALLOC_ATOMIC(n*sizeof(double));
808: k3 = (double *)MALLOC_ATOMIC(n*sizeof(double));
809: k4 = (double *)MALLOC_ATOMIC(n*sizeof(double));
810: k5 = (double *)MALLOC_ATOMIC(n*sizeof(double));
1.12 noro 811: k6 = (double *)MALLOC(n*sizeof(double));
812: nd = 0;
813: switch ( ord ) {
814: case 4:
815: a2 = 1/2.0*h; b21 = 1/2.0*h;
816: a3 = 1/2.0*h; b31 = 0.0; b32 = 1/2.0*h;
817: a4 = 1.0*h; b41 = 0.0; b42 = 0.0; b43 = 1.0*h;
818: c1 = 1/6.0*h; c2 = 1/3.0*h; c3 = 1/3.0*h; c4 = 1/6.0*h;
819: for ( i = 0; i < step; i++ ) {
820: if ( !(i%100000) ) fprintf(stderr,"[%d]",i);
821: xi = x0+i*h;
822: eval_pfaffian2(k1,n,d,num,den,xi,f);
823: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b21*k1[j];
824: eval_pfaffian2(k2,n,d,num,den,xi+a2,w);
825: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b31*k1[j]+b32*k2[j];
826: eval_pfaffian2(k3,n,d,num,den,xi+a3,w);
827: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b41*k1[j]+b42*k2[j]+b43*k3[j];
828: eval_pfaffian2(k4,n,d,num,den,xi+a4,w);
829: 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];
830: memcpy(f,w,n*sizeof(double));
831: MKReal(f[0],t);
832: MKReal(xi+h,x);
833: nd1 = mknode(2,x,t);
834: MKLIST(l,nd1);
835: MKNODE(nd1,l,nd);
836: nd = nd1;
837: for ( hd = f[0], j = 0; j < n; j++ ) f[j] /= hd;
838: }
839: MKLIST(*rp,nd);
840: break;
841: case 5:
842: default:
843: a2 = 1/4.0*h; b21 = 1/4.0*h;
844: a3 = 1/4.0*h; b31 = 1/8.0*h; b32 = 1/8.0*h;
845: a4 = 1/2.0*h; b41 = 0.0; b42 = 0.0; b43 = 1/2.0*h;
846: 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;
847: 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;
848: 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;
849: for ( i = 0; i < step; i++ ) {
850: if ( !(i%100000) ) fprintf(stderr,"[%d]",i);
851: xi = x0+i*h;
852: eval_pfaffian2(k1,n,d,num,den,xi,f);
853: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b21*k1[j];
854: eval_pfaffian2(k2,n,d,num,den,xi+a2,w);
855: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b31*k1[j]+b32*k2[j];
856: eval_pfaffian2(k3,n,d,num,den,xi+a3,w);
857: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b41*k1[j]+b42*k2[j]+b43*k3[j];
858: eval_pfaffian2(k4,n,d,num,den,xi+a4,w);
859: 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];
860: eval_pfaffian2(k5,n,d,num,den,xi+a5,w);
861: 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];
862: eval_pfaffian2(k6,n,d,num,den,xi+a6,w);
863: 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];
864: memcpy(f,w,n*sizeof(double));
865: MKReal(f[0],t);
866: MKReal(xi+h,x);
867: nd1 = mknode(2,x,t);
868: MKLIST(l,nd1);
869: MKNODE(nd1,l,nd);
870: nd = nd1;
871: for ( hd = f[0], j = 0; j < n; j++ ) f[j] /= hd;
872: }
873: MKLIST(*rp,nd);
874: break;
875: }
876: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bfaux.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin