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