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