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