Annotation of OpenXM_contrib2/asir2000/builtin/bfaux.c, Revision 1.16
1.16 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.15 2017/08/31 04:21:48 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[] = {
20: {"eval",Peval,-2},
21: {"setprec",Psetprec,-1},
1.3 noro 22: {"setbprec",Psetbprec,-1},
1.4 noro 23: {"setround",Psetround,-1},
1.1 noro 24: {"todouble",Ptodouble,1},
1.14 noro 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},
1.8 noro 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},
1.9 noro 48: {"mpfr_erfc",Pmpfr_erfc,-2},
1.8 noro 49: {"mpfr_li2",Pmpfr_li2,-2},
1.5 noro 50: {"mpfr_gamma",Pmpfr_gamma,-2},
1.8 noro 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},
1.12 noro 56: {"rk_ratmat",Prk_ratmat,7},
1.1 noro 57: {0,0,0},
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: {
64: double r,i;
65: Real real,imag;
66:
1.15 noro 67: if ( !a ) {
1.1 noro 68: *rp = 0;
69: return;
70: }
1.15 noro 71: switch ( NID(a) ) {
1.1 noro 72: case N_R: case N_Q: case N_B:
1.15 noro 73: r = ToReal(a);
1.1 noro 74: MKReal(r,real);
75: *rp = (Num)real;
76: break;
77: case N_C:
1.15 noro 78: r = ToReal(((C)a)->r);
79: i = ToReal(((C)a)->i);
1.1 noro 80: MKReal(r,real);
81: MKReal(i,imag);
82: reimtocplx((Num)real,(Num)imag,rp);
83: break;
84: default:
1.15 noro 85: *rp = a;
1.1 noro 86: break;
87: }
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: {
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:
1.4 noro 205: void Peval(NODE arg,Obj *rp)
1.1 noro 206: {
207: int prec;
208:
209: asir_assert(ARG0(arg),O_R,"eval");
210: if ( argc(arg) == 2 ) {
1.3 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.2 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: {
223: int p;
224: Q q;
1.3 noro 225: int prec,dprec;
226:
227: prec = mpfr_get_default_prec();
228: /* decimal precision */
229: dprec = prec*0.30103;
230: STOQ(dprec,q); *rp = (Obj)q;
231: if ( arg ) {
232: asir_assert(ARG0(arg),O_N,"setprec");
1.11 ohara 233: p = QTOS((Q)ARG0(arg))*3.32193;
1.3 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: }
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: {
246: int p;
247: Q q;
248: int prec;
249:
250: prec = mpfr_get_default_prec();
251: STOQ(prec,q); *rp = (Obj)q;
1.1 noro 252: if ( arg ) {
1.3 noro 253: asir_assert(ARG0(arg),O_N,"setbprec");
1.11 ohara 254: p = QTOS((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:
1.4 noro 263: void Psetround(NODE arg,Q *rp)
264: {
265: int round;
266:
267: STOQ(mpfr_roundmode,*rp);
268: if ( arg ) {
269: asir_assert(ARG0(arg),O_N,"setround");
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.1 noro 305: BF r;
306:
307: prec = arg ? QTOS((Q)ARG0(arg)) : 0;
308: NEWBF(r);
309: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4 noro 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.1 noro 318: mpfr_t one;
319: BF r;
320:
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);
1.4 noro 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: {
333: Num a;
1.10 noro 334: int prec;
335: BF r,re,im;
336: C c;
337: mpc_t mpc,a1;
1.1 noro 338:
1.10 noro 339: prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
1.1 noro 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 {
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: }
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: {
458: Num a,e;
1.10 noro 459: int prec;
460: BF r,re,im;
461: C c;
462: mpc_t mpc,a1,e1;
1.1 noro 463:
1.10 noro 464: prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : mpfr_get_default_prec();
1.1 noro 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 {
495: NEWBF(r);
496: mpfr_init2(r->body,prec);
497: mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
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: {
512: Num a;
513: int prec;
514: BF r;
515:
1.8 noro 516: SETPREC
1.5 noro 517: mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
518: *rp = r;
519: }
1.6 takayama 520:
1.8 noro 521: void Pmpfr_lngamma(NODE arg,BF *rp)
522: {
523: Num a;
524: int prec;
525: BF r;
526:
527: SETPREC
528: mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode);
529: *rp = r;
530: }
531:
532: void Pmpfr_digamma(NODE arg,BF *rp)
533: {
534: Num a;
535: int prec;
536: BF r;
537:
538: SETPREC
539: mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode);
540: *rp = r;
541: }
542:
543: void Pmpfr_zeta(NODE arg,BF *rp)
544: {
545: Num a;
546: int prec;
547: BF r;
548:
549: SETPREC
550: mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode);
551: *rp = r;
552: }
553:
554: void Pmpfr_eint(NODE arg,BF *rp)
555: {
556: Num a;
557: int prec;
558: BF r;
559:
560: SETPREC
561: mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode);
562: *rp = r;
563: }
564:
565: void Pmpfr_erf(NODE arg,BF *rp)
566: {
567: Num a;
568: int prec;
569: BF r;
570:
571: SETPREC
572: mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode);
573: *rp = r;
574: }
575:
1.9 noro 576: void Pmpfr_erfc(NODE arg,BF *rp)
577: {
578: Num a;
579: int prec;
580: BF r;
581:
582: SETPREC
583: mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode);
584: *rp = r;
585: }
586:
1.8 noro 587: void Pmpfr_j0(NODE arg,BF *rp)
588: {
589: Num a;
590: int prec;
591: BF r;
592:
593: SETPREC
594: mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode);
595: *rp = r;
596: }
597:
598: void Pmpfr_j1(NODE arg,BF *rp)
599: {
600: Num a;
601: int prec;
602: BF r;
603:
604: SETPREC
605: mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode);
606: *rp = r;
607: }
608:
609: void Pmpfr_y0(NODE arg,BF *rp)
610: {
611: Num a;
612: int prec;
613: BF r;
614:
615: SETPREC
616: mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode);
617: *rp = r;
618: }
619:
620: void Pmpfr_y1(NODE arg,BF *rp)
621: {
622: Num a;
623: int prec;
624: BF r;
625:
626: SETPREC
627: mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode);
628: *rp = r;
629: }
630:
631: void Pmpfr_li2(NODE arg,BF *rp)
632: {
633: Num a;
634: int prec;
635: BF r;
636:
637: SETPREC
638: mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode);
639: *rp = r;
640: }
641:
642: void Pmpfr_ai(NODE arg,BF *rp)
643: {
644: Num a;
645: int prec;
646: BF r;
647:
648: SETPREC
649: mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode);
650: *rp = r;
651: }
652:
1.7 takayama 653: void Pmpfr_floor(NODE arg,Q *rp)
1.6 takayama 654: {
655: Num a;
656: int prec;
657: BF r;
1.7 takayama 658: mpz_t t;
659: GZ rz;
1.6 takayama 660:
1.8 noro 661: SETPREC
1.6 takayama 662: mpfr_floor(r->body,((BF)a)->body);
1.7 takayama 663: mpz_init(t);
664: mpfr_get_z(t,r->body,mpfr_roundmode);
665: MPZTOGZ(t,rz);
666: *rp = gztoz(rz);
667: }
668:
1.8 noro 669: void Pmpfr_ceil(NODE arg,Q *rp)
670: {
671: Num a;
672: int prec;
673: BF r;
674: mpz_t t;
675: GZ 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: MPZTOGZ(t,rz);
682: *rp = gztoz(rz);
683: }
684:
1.7 takayama 685: void Pmpfr_round(NODE arg,Q *rp)
686: {
687: Num a;
688: int prec;
689: BF r;
690: mpz_t t;
691: GZ rz;
692:
1.8 noro 693: SETPREC
1.7 takayama 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