Annotation of OpenXM_contrib2/asir2018/builtin/bfaux.c, Revision 1.2
1.2 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2018/builtin/bfaux.c,v 1.1 2018/09/19 05:45:05 noro Exp $ */
1.1 noro 2: #include "ca.h"
3: #include "parse.h"
4:
5: void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround();
6: void Pmpfr_ai();
7: void Pmpfr_eint(), Pmpfr_erf(), Pmpfr_erfc(), Pmpfr_li2();
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();
13: void Prk_ratmat();
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();
17: void mp_factorial(),mp_abs();
18:
19: struct ftab bf_tab[] = {
20: {"eval",Peval,-2},
21: {"setprec",Psetprec,-1},
22: {"setbprec",Psetbprec,-1},
23: {"setround",Psetround,-1},
24: {"todouble",Ptodouble,1},
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},
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},
48: {"mpfr_erfc",Pmpfr_erfc,-2},
49: {"mpfr_li2",Pmpfr_li2,-2},
50: {"mpfr_gamma",Pmpfr_gamma,-2},
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},
56: {"rk_ratmat",Prk_ratmat,7},
57: {0,0,0},
58: };
59:
60: int mpfr_roundmode = MPFR_RNDN;
61:
62: void todoublen(Num a,Num *rp)
63: {
64: double r,i;
65: Real real,imag;
66:
67: if ( !a ) {
68: *rp = 0;
69: return;
70: }
71: switch ( NID(a) ) {
72: case N_R: case N_Q: case N_B:
73: r = ToReal(a);
74: MKReal(r,real);
75: *rp = (Num)real;
76: break;
77: case N_C:
78: r = ToReal(((C)a)->r);
79: i = ToReal(((C)a)->i);
80: MKReal(r,real);
81: MKReal(i,imag);
82: reimtocplx((Num)real,(Num)imag,rp);
83: break;
84: default:
85: *rp = a;
86: break;
87: }
88: }
89:
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:
205: void Peval(NODE arg,Obj *rp)
206: {
207: long prec;
208:
209: asir_assert(ARG0(arg),O_R,"eval");
210: if ( argc(arg) == 2 ) {
1.2 ! noro 211: prec = ZTOS((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;
216: evalr(CO,(Obj)ARG0(arg),prec,rp);
217: }
218:
219: /* set/get decimal precision */
220:
221: void Psetprec(NODE arg,Obj *rp)
222: {
223: long p;
224: Z q;
225: long prec,dprec;
226:
227: prec = mpfr_get_default_prec();
228: /* decimal precision */
229: dprec = prec*0.30103;
1.2 ! noro 230: STOZ(dprec,q); *rp = (Obj)q;
1.1 noro 231: if ( arg ) {
232: asir_assert(ARG0(arg),O_N,"setprec");
1.2 ! noro 233: p = ZTOS((Q)ARG0(arg))*3.32193;
1.1 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: }
241:
242: /* set/get bit precision */
243:
244: void Psetbprec(NODE arg,Obj *rp)
245: {
246: long p;
247: Z q;
248: long prec;
249:
250: prec = mpfr_get_default_prec();
1.2 ! noro 251: STOZ(prec,q); *rp = (Obj)q;
1.1 noro 252: if ( arg ) {
253: asir_assert(ARG0(arg),O_N,"setbprec");
1.2 ! noro 254: p = ZTOS((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:
263: void Psetround(NODE arg,Z *rp)
264: {
265: int round;
266:
1.2 ! noro 267: STOZ(mpfr_roundmode,*rp);
1.1 noro 268: if ( arg ) {
269: asir_assert(ARG0(arg),O_N,"setround");
1.2 ! noro 270: round = ZTOS((Q)ARG0(arg));
1.1 noro 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:
300: Num tobf(Num a,int prec);
301:
302: void mp_pi(NODE arg,BF *rp)
303: {
304: int prec;
305: BF r;
306:
1.2 ! noro 307: prec = arg ? ZTOS((Q)ARG0(arg)) : 0;
1.1 noro 308: NEWBF(r);
309: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
310: mpfr_const_pi(r->body,mpfr_roundmode);
311: if ( !cmpbf((Num)r,0) ) r = 0;
312: *rp = r;
313: }
314:
315: void mp_e(NODE arg,BF *rp)
316: {
317: int prec;
318: mpfr_t one;
319: BF r;
320:
1.2 ! noro 321: prec = arg ? ZTOS((Q)ARG0(arg)) : 0;
1.1 noro 322: NEWBF(r);
323: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
324: mpfr_init(one);
325: mpfr_set_ui(one,1,mpfr_roundmode);
326: mpfr_exp(r->body,one,mpfr_roundmode);
327: if ( !cmpbf((Num)r,0) ) r = 0;
328: *rp = r;
329: }
330:
331: void mpfr_or_mpc(NODE arg,int (*mpfr_f)(),int (*mpc_f)(),Num *rp)
332: {
333: Num a;
334: int prec;
335: BF r,re,im;
336: C c;
337: mpc_t mpc,a1;
338:
1.2 ! noro 339: prec = NEXT(arg) ? ZTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
1.1 noro 340: a = tobf(ARG0(arg),prec);
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: }
363: }
364:
365: void mp_sin(NODE arg,Num *rp)
366: {
367: mpfr_or_mpc(arg,mpfr_sin,mpc_sin,rp);
368: }
369:
370: void mp_cos(NODE arg,Num *rp)
371: {
372: mpfr_or_mpc(arg,mpfr_cos,mpc_cos,rp);
373: }
374:
375: void mp_tan(NODE arg,Num *rp)
376: {
377: mpfr_or_mpc(arg,mpfr_tan,mpc_tan,rp);
378: }
379:
380: void mp_asin(NODE arg,Num *rp)
381: {
382: mpfr_or_mpc(arg,mpfr_asin,mpc_asin,rp);
383: }
384:
385: void mp_acos(NODE arg,Num *rp)
386: {
387: mpfr_or_mpc(arg,mpfr_acos,mpc_acos,rp);
388: }
389:
390: void mp_atan(NODE arg,Num *rp)
391: {
392: mpfr_or_mpc(arg,mpfr_atan,mpc_atan,rp);
393: }
394:
395: void mp_sinh(NODE arg,Num *rp)
396: {
397: mpfr_or_mpc(arg,mpfr_sinh,mpc_sinh,rp);
398: }
399:
400: void mp_cosh(NODE arg,Num *rp)
401: {
402: mpfr_or_mpc(arg,mpfr_cosh,mpc_cosh,rp);
403: }
404:
405: void mp_tanh(NODE arg,Num *rp)
406: {
407: mpfr_or_mpc(arg,mpfr_tanh,mpc_tanh,rp);
408: }
409:
410: void mp_asinh(NODE arg,Num *rp)
411: {
412: mpfr_or_mpc(arg,mpfr_asinh,mpc_asinh,rp);
413: }
414:
415: void mp_acosh(NODE arg,Num *rp)
416: {
417: mpfr_or_mpc(arg,mpfr_acosh,mpc_acosh,rp);
418: }
419:
420: void mp_atanh(NODE arg,Num *rp)
421: {
422: mpfr_or_mpc(arg,mpfr_atanh,mpc_atanh,rp);
423: }
424:
425: void mp_exp(NODE arg,Num *rp)
426: {
427: mpfr_or_mpc(arg,mpfr_exp,mpc_exp,rp);
428: }
429:
430: void mp_log(NODE arg,Num *rp)
431: {
432: mpfr_or_mpc(arg,mpfr_log,mpc_log,rp);
433: }
434:
435: void mp_abs(NODE arg,Num *rp)
436: {
437: mpfr_or_mpc(arg,mpfr_abs,mpc_abs,rp);
438: }
439:
440: void Pfac(NODE arg,Num *rp);
441:
442: void mp_factorial(NODE arg,Num *rp)
443: {
444: struct oNODE arg0;
445: Num a,a1;
446:
447: a = (Num)ARG0(arg);
448: if ( !a ) *rp = (Num)ONE;
449: else if ( INT(a) ) Pfac(arg,rp);
450: else {
451: addnum(0,a,(Num)ONE,&a1);
452: arg0.body = (pointer)a1;
453: arg0.next = arg->next;
454: Pmpfr_gamma(&arg0,rp);
455: }
456: }
457:
458: void mp_pow(NODE arg,Num *rp)
459: {
460: Num a,e;
461: int prec;
462: BF r,re,im;
463: C c;
464: mpc_t mpc,a1,e1;
465:
1.2 ! noro 466: prec = NEXT(NEXT(arg)) ? ZTOS((Q)ARG2(arg)) : mpfr_get_default_prec();
1.1 noro 467: a = tobf(ARG0(arg),prec);
468: e = tobf(ARG1(arg),prec);
469: if ( NID(a) == N_C || NID(e) == N_C || MPFR_SIGN(((BF)a)->body) < 0 ) {
470: mpc_init2(mpc,prec); mpc_init2(a1,prec); mpc_init2(e1,prec);
471: if ( NID(a) == N_C ) {
472: re = (BF)((C)a)->r; im = (BF)((C)a)->i;
473: mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode);
474: } else {
475: re = (BF)a;
476: mpc_set_fr(a1,re->body,mpfr_roundmode);
477: }
478: if ( NID(e) == N_C ) {
479: re = (BF)((C)e)->r; im = (BF)((C)e)->i;
480: mpc_set_fr_fr(e1,re->body,im->body,mpfr_roundmode);
481: } else {
482: re = (BF)e;
483: mpc_set_fr(e1,re->body,mpfr_roundmode);
484: }
485: mpc_pow(mpc,a1,e1,mpfr_roundmode);
486: MPFRTOBF(mpc_realref(mpc),re);
487: MPFRTOBF(mpc_imagref(mpc),im);
488: if ( !cmpbf((Num)re,0) ) re = 0;
489: if ( !cmpbf((Num)im,0) ) im = 0;
490: if ( !im )
491: *rp = (Num)re;
492: else {
493: NEWC(c); c->r = (Num)re; c->i = (Num)im;
494: *rp = (Num)c;
495: }
496: } else {
497: NEWBF(r);
498: mpfr_init2(r->body,prec);
499: mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
500: *rp = (Num)r;
501: }
502: }
503:
504: #define SETPREC \
1.2 ! noro 505: (prec)=NEXT(arg)?ZTOS((Q)ARG1(arg)):0;\
1.1 noro 506: (prec)*=3.32193;\
507: (a)=tobf(ARG0(arg),prec);\
508: NEWBF(r);\
509: prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
510:
511:
512: void Pmpfr_gamma(NODE arg,BF *rp)
513: {
514: Num a;
515: int prec;
516: BF r;
517:
518: SETPREC
519: mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
520: *rp = r;
521: }
522:
523: void Pmpfr_lngamma(NODE arg,BF *rp)
524: {
525: Num a;
526: int prec;
527: BF r;
528:
529: SETPREC
530: mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode);
531: *rp = r;
532: }
533:
534: void Pmpfr_digamma(NODE arg,BF *rp)
535: {
536: Num a;
537: int prec;
538: BF r;
539:
540: SETPREC
541: mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode);
542: *rp = r;
543: }
544:
545: void Pmpfr_zeta(NODE arg,BF *rp)
546: {
547: Num a;
548: int prec;
549: BF r;
550:
551: SETPREC
552: mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode);
553: *rp = r;
554: }
555:
556: void Pmpfr_eint(NODE arg,BF *rp)
557: {
558: Num a;
559: int prec;
560: BF r;
561:
562: SETPREC
563: mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode);
564: *rp = r;
565: }
566:
567: void Pmpfr_erf(NODE arg,BF *rp)
568: {
569: Num a;
570: int prec;
571: BF r;
572:
573: SETPREC
574: mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode);
575: *rp = r;
576: }
577:
578: void Pmpfr_erfc(NODE arg,BF *rp)
579: {
580: Num a;
581: int prec;
582: BF r;
583:
584: SETPREC
585: mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode);
586: *rp = r;
587: }
588:
589: void Pmpfr_j0(NODE arg,BF *rp)
590: {
591: Num a;
592: int prec;
593: BF r;
594:
595: SETPREC
596: mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode);
597: *rp = r;
598: }
599:
600: void Pmpfr_j1(NODE arg,BF *rp)
601: {
602: Num a;
603: int prec;
604: BF r;
605:
606: SETPREC
607: mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode);
608: *rp = r;
609: }
610:
611: void Pmpfr_y0(NODE arg,BF *rp)
612: {
613: Num a;
614: int prec;
615: BF r;
616:
617: SETPREC
618: mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode);
619: *rp = r;
620: }
621:
622: void Pmpfr_y1(NODE arg,BF *rp)
623: {
624: Num a;
625: int prec;
626: BF r;
627:
628: SETPREC
629: mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode);
630: *rp = r;
631: }
632:
633: void Pmpfr_li2(NODE arg,BF *rp)
634: {
635: Num a;
636: int prec;
637: BF r;
638:
639: SETPREC
640: mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode);
641: *rp = r;
642: }
643:
644: void Pmpfr_ai(NODE arg,BF *rp)
645: {
646: Num a;
647: int prec;
648: BF r;
649:
650: SETPREC
651: mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode);
652: *rp = r;
653: }
654:
655: void Pmpfr_floor(NODE arg,Z *rp)
656: {
657: Num a;
658: long prec;
659: BF r;
660: mpz_t t;
661:
662: SETPREC
663: mpfr_floor(r->body,((BF)a)->body);
664: mpz_init(t);
665: mpfr_get_z(t,r->body,mpfr_roundmode);
666: MPZTOZ(t,*rp);
667: }
668:
669: void Pmpfr_ceil(NODE arg,Z *rp)
670: {
671: Num a;
672: long prec;
673: BF r;
674: mpz_t t;
675: Z 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: MPZTOZ(t,*rp);
682: }
683:
684: void Pmpfr_round(NODE arg,Z *rp)
685: {
686: Num a;
687: long prec;
688: BF r;
689: mpz_t t;
690:
691: SETPREC
692: mpfr_round(r->body,((BF)a)->body);
693: mpz_init(t);
694: mpfr_get_z(t,r->body,mpfr_roundmode);
695: MPZTOZ(t,*rp);
696: }
697:
698: double **almat_double(int n)
699: {
700: int i;
701: double **a;
702:
703: a = (double **)MALLOC(n*sizeof(double *));
704: for ( i = 0; i < n; i++ )
705: a[i] = (double *)MALLOC(n*sizeof(double));
706: return a;
707: }
708:
709: /*
710: * k <- (A(xi)-(sbeta-mn2/xi))f
711: * A(t) = (num[0]+num[1]t+...+num[d-1]*t^(d-1))/den(t)
712: */
713:
714: struct jv {
715: int j;
716: double v;
717: };
718:
719: struct smat {
720: int *rlen;
721: struct jv **row;
722: };
723:
724: void eval_pfaffian2(double *k,int n,int d,struct smat *num,P den,double xi,double *f)
725: {
726: struct smat ma;
727: struct jv *maj;
728: int i,j,l,s;
729: double t,dn;
730: P r;
731: Real u;
732:
733: memset(k,0,n*sizeof(double));
734: for ( i = d-1; i >= 0; i-- ) {
735: ma = num[i];
736: for ( j = 0; j < n; j++ ) {
737: maj = ma.row[j];
738: l = ma.rlen[j];
739: for ( t = 0, s = 0; s < l; s++, maj++ ) t += maj->v*f[maj->j];
740: k[j] = k[j]*xi+t;
741: }
742: }
743: MKReal(xi,u);
744: substp(CO,den,den->v,(P)u,&r); dn = ToReal(r);
745: for ( j = 0; j < n; j++ )
746: k[j] /= dn;
747: }
748:
749: void Prk_ratmat(NODE arg,LIST *rp)
750: {
751: VECT mat;
752: P den;
753: int ord;
754: double sbeta,x0,x1,xi,h,mn2,hd;
755: double a2,a3,a4,a5,a6;
756: double b21,b31,b32,b41,b42,b43,b51,b52,b53,b54,b61,b62,b63,b64,b65;
757: double c1,c2,c3,c4,c5,c6,c7;
758: VECT fv;
759: int step,j,i,k,n,d,len,s;
760: struct smat *num;
761: Obj **b;
762: MAT mati;
763: double *f,*w,*k1,*k2,*k3,*k4,*k5,*k6;
764: NODE nd,nd1;
765: Real x,t;
766: LIST l;
767:
1.2 ! noro 768: ord = ZTOS((Q)ARG0(arg));
1.1 noro 769: mat = (VECT)ARG1(arg); den = (P)ARG2(arg);
770: x0 = ToReal((Num)ARG3(arg)); x1 = ToReal((Num)ARG4(arg));
1.2 ! noro 771: step = ZTOS((Q)ARG5(arg)); fv = (VECT)ARG6(arg);
1.1 noro 772: h = (x1-x0)/step;
773:
774: n = fv->len;
775: d = mat->len;
776: num = (struct smat *)MALLOC(d*sizeof(struct smat));
777: for ( i = 0; i < d; i++ ) {
778: num[i].rlen = (int *)MALLOC_ATOMIC(n*sizeof(int));
779: num[i].row = (struct jv **)MALLOC(n*sizeof(struct jv *));
780: mati = (MAT)mat->body[i];
781: b = (Obj **)mati->body;
782: for ( j = 0; j < n; j++ ) {
783: for ( len = k = 0; k < n; k++ )
784: if ( b[j][k] ) len++;
785: num[i].rlen[j] = len;
786: if ( !len )
787: num[i].row[j] = 0;
788: else {
789: num[i].row[j] = (struct jv *)MALLOC_ATOMIC((len)*sizeof(struct jv));
790: for ( s = k = 0; k < n; k++ )
791: if ( b[j][k] ) {
792: num[i].row[j][s].j = k;
793: num[i].row[j][s].v = ToReal((Num)b[j][k]);
794: s++;
795: }
796: }
797: }
798: }
799: f = (double *)MALLOC_ATOMIC(n*sizeof(double));
800: for ( j = 0; j < n; j++ )
801: f[j] = ToReal((Num)fv->body[j]);
802: w = (double *)MALLOC_ATOMIC(n*sizeof(double));
803: k1 = (double *)MALLOC_ATOMIC(n*sizeof(double));
804: k2 = (double *)MALLOC_ATOMIC(n*sizeof(double));
805: k3 = (double *)MALLOC_ATOMIC(n*sizeof(double));
806: k4 = (double *)MALLOC_ATOMIC(n*sizeof(double));
807: k5 = (double *)MALLOC_ATOMIC(n*sizeof(double));
808: k6 = (double *)MALLOC(n*sizeof(double));
809: nd = 0;
810: switch ( ord ) {
811: case 4:
812: a2 = 1/2.0*h; b21 = 1/2.0*h;
813: a3 = 1/2.0*h; b31 = 0.0; b32 = 1/2.0*h;
814: a4 = 1.0*h; b41 = 0.0; b42 = 0.0; b43 = 1.0*h;
815: c1 = 1/6.0*h; c2 = 1/3.0*h; c3 = 1/3.0*h; c4 = 1/6.0*h;
816: for ( i = 0; i < step; i++ ) {
817: if ( !(i%100000) ) fprintf(stderr,"[%d]",i);
818: xi = x0+i*h;
819: eval_pfaffian2(k1,n,d,num,den,xi,f);
820: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b21*k1[j];
821: eval_pfaffian2(k2,n,d,num,den,xi+a2,w);
822: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b31*k1[j]+b32*k2[j];
823: eval_pfaffian2(k3,n,d,num,den,xi+a3,w);
824: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b41*k1[j]+b42*k2[j]+b43*k3[j];
825: eval_pfaffian2(k4,n,d,num,den,xi+a4,w);
826: 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];
827: memcpy(f,w,n*sizeof(double));
828: MKReal(f[0],t);
829: MKReal(xi+h,x);
830: nd1 = mknode(2,x,t);
831: MKLIST(l,nd1);
832: MKNODE(nd1,l,nd);
833: nd = nd1;
834: for ( hd = f[0], j = 0; j < n; j++ ) f[j] /= hd;
835: }
836: MKLIST(*rp,nd);
837: break;
838: case 5:
839: default:
840: a2 = 1/4.0*h; b21 = 1/4.0*h;
841: a3 = 1/4.0*h; b31 = 1/8.0*h; b32 = 1/8.0*h;
842: a4 = 1/2.0*h; b41 = 0.0; b42 = 0.0; b43 = 1/2.0*h;
843: 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;
844: 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;
845: 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;
846: for ( i = 0; i < step; i++ ) {
847: if ( !(i%100000) ) fprintf(stderr,"[%d]",i);
848: xi = x0+i*h;
849: eval_pfaffian2(k1,n,d,num,den,xi,f);
850: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b21*k1[j];
851: eval_pfaffian2(k2,n,d,num,den,xi+a2,w);
852: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b31*k1[j]+b32*k2[j];
853: eval_pfaffian2(k3,n,d,num,den,xi+a3,w);
854: memcpy(w,f,n*sizeof(double)); for ( j = 0; j < n; j++ ) w[j] += b41*k1[j]+b42*k2[j]+b43*k3[j];
855: eval_pfaffian2(k4,n,d,num,den,xi+a4,w);
856: 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];
857: eval_pfaffian2(k5,n,d,num,den,xi+a5,w);
858: 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];
859: eval_pfaffian2(k6,n,d,num,den,xi+a6,w);
860: 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];
861: memcpy(f,w,n*sizeof(double));
862: MKReal(f[0],t);
863: MKReal(xi+h,x);
864: nd1 = mknode(2,x,t);
865: MKLIST(l,nd1);
866: MKNODE(nd1,l,nd);
867: nd = nd1;
868: for ( hd = f[0], j = 0; j < n; j++ ) f[j] /= hd;
869: }
870: MKLIST(*rp,nd);
871: break;
872: }
873: }
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 / asir2018 / builtin