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