Annotation of OpenXM_contrib2/asir2018/builtin/parif.c, Revision 1.5
1.5 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2018/builtin/parif.c,v 1.4 2020/08/26 06:40:36 noro Exp $ */
1.1 noro 2: #include "ca.h"
3: #include "parse.h"
4: #include "ox.h"
5:
6: Q ox_pari_stream;
7: int ox_pari_stream_initialized = 0;
8: int ox_get_pari_result = 0;
9: P ox_pari_starting_function = 0;
10:
11: typedef void (*mpfr_func)(NODE,Obj *);
12:
13: void Pmpfr_ai();
14: void Pmpfr_eint(), Pmpfr_erf(),Pmpfr_li2();
15: void Pmpfr_zeta();
16: void Pmpfr_j0(), Pmpfr_j1();
17: void Pmpfr_y0(), Pmpfr_y1();
18: void Pmpfr_gamma(), Pmpfr_lngamma(), Pmpfr_digamma();
19: void Pmpfr_floor(), Pmpfr_round(), Pmpfr_ceil();
20:
21: void Pox_shutdown(NODE arg,Q *rp);
22: void Pox_launch_nox(NODE arg,Obj *rp);
1.4 noro 23: void Pox_launch(NODE arg,Obj *rp);
1.1 noro 24: void Pox_push_cmo(NODE arg,Obj *rp);
25: void Pox_pop_cmo(NODE arg,Obj *rp);
26: void Pox_execute_function(NODE arg,Obj *rp);
27:
1.4 noro 28: int debug_pari;
29:
1.1 noro 30: struct mpfr_tab_rec {
31: char *name;
32: mpfr_func func;
33: } mpfr_tab[] = {
34: {"ai",Pmpfr_ai},
35: {"zeta",Pmpfr_zeta},
36: {"j0",Pmpfr_j0},
37: {"j1",Pmpfr_j1},
38: {"y0",Pmpfr_y0},
39: {"y1",Pmpfr_y1},
40: {"eint",Pmpfr_eint},
41: {"erf",Pmpfr_erf},
42: {"li2",Pmpfr_li2},
43: {"gamma",Pmpfr_gamma},
44: {"lngamma",Pmpfr_gamma},
45: {"digamma",Pmpfr_gamma},
46: {"floor",Pmpfr_floor},
47: {"ceil",Pmpfr_ceil},
48: {"round",Pmpfr_round},
49: };
50:
51: mpfr_func mpfr_search(char *name)
52: {
53: int i,n;
54:
55: n = sizeof(mpfr_tab)/sizeof(struct mpfr_tab_rec);
56: for ( i = 0; i < n; i++ )
57: if ( !strcmp(name,mpfr_tab[i].name) )
58: return mpfr_tab[i].func;
59: return 0;
60: }
61:
62: Obj list_to_vect(Obj a)
63: {
64: int len,i;
65: VECT v;
66: NODE nd;
67:
68: if ( !a || OID(a) != O_LIST ) return a;
69: len = length(BDY((LIST)a));
70: MKVECT(v,len);
71: for ( i = 0, nd = BDY((LIST)a); nd; nd = NEXT(nd), i++ )
72: v->body[i] = (pointer)list_to_vect((Obj)BDY(nd));
73: return (Obj)v;
74: }
75:
76: Obj vect_to_mat(VECT v)
77: {
78: MAT m;
79: int len,col,i,j;
80:
81: len = v->len;
82: if ( v->body[0] && OID((Obj)v->body[0]) == O_VECT ) {
83: col = ((VECT)v->body[0])->len;
84: for ( i = 1; i < len; i++ )
85: if ( !v->body[i] || OID((Obj)v->body[i]) != O_VECT
86: || ((VECT)v->body[i])->len != col )
87: break;
88: if ( i == len ) {
89: /* convert to a matrix */
90: MKMAT(m,len,col);
91: for ( i = 0; i < len; i++ )
92: for ( j = 0; j < col; j++ )
93: m->body[i][j] = ((VECT)v->body[i])->body[j];
94: return (Obj)m;
95: }
96: }
97: return (Obj)v;
98: }
99:
100: void reset_ox_pari()
101: {
102: NODE nd;
103: Q r;
104:
105: if ( ox_get_pari_result ) {
106: nd = mknode(1,ox_pari_stream);
107: Pox_shutdown(nd,&r);
108: ox_get_pari_result = 0;
109: ox_pari_stream_initialized = 0;
110: }
111: }
112:
1.5 ! noro 113: void pari_setprec(long n)
! 114: {
! 115: struct oFUNC f;
! 116: Z z;
! 117: NODE arg;
! 118:
! 119: f.name = f.fullname = "pari_setprec";
! 120: f.id = A_PARI;
! 121: f.argc = 1;
! 122: f.f.binf = 0;
! 123: STOZ(n,z);
! 124: arg = mknode(1,z);
! 125: evalparif(&f,arg);
! 126: }
! 127:
! 128: /* decimal precision */
! 129: long mpfrprec = 0;
! 130:
1.1 noro 131: pointer evalparif(FUNC f,NODE arg)
132: {
1.5 ! noro 133: long prec;
! 134: int ac,intarg,opt;
1.1 noro 135: Q q,r,cmd;
136: Z narg;
137: Real sec;
138: NODE nd,oxarg,t,t1,n;
139: STRING name;
140: USINT ui;
141: LIST list;
142: Obj ret,dmy;
143: mpfr_func mpfr_function;
144: V v;
145:
1.3 noro 146: if ( arg && (!ARG0(arg) || (NUM(ARG0(arg)) && NID((Num)ARG0(arg)) != N_C))
1.1 noro 147: && (mpfr_function = mpfr_search(f->name)) ) {
148: (*mpfr_function)(arg,&ret);
149: return (pointer) ret;
150: }
151:
152: if ( !ox_pari_stream_initialized ) {
153: if ( ox_pari_starting_function && OID(ox_pari_starting_function) == O_P ) {
154: v = VR(ox_pari_starting_function);
155: if ( (long)v->attr != V_SR ) {
156: error("pari : no handler.");
157: }
158: MKNODE(nd,0,0);
159: r = (Q)bevalf((FUNC)v->priv,0);
160: }else {
161: #if !defined(VISUAL)
162: MKSTR(name,"ox_pari");
163: nd = mknode(2,NULL,name);
1.4 noro 164: if ( debug_pari )
165: Pox_launch(nd,(Obj *)&r);
166: else
167: Pox_launch_nox(nd,(Obj *)&r);
1.1 noro 168: #else
169: error("Please load names.rr from latest asir-contrib library before using pari functions.");
170: #endif
171: }
172: ox_pari_stream = r;
173: ox_pari_stream_initialized = 1;
174: }
175:
1.5 ! noro 176: prec = mpfr_get_default_prec()*0.30103+1;
! 177: if ( prec != mpfrprec ) {
! 178: mpfrprec = prec;
! 179: pari_setprec(prec);
! 180: }
1.1 noro 181: ac = argc(arg);
182: /* reverse the arg list */
183: for ( n = arg, t = 0; n; n = NEXT(n) ) {
184: MKNODE(t1,BDY(n),t); t = t1;
185: }
186: /* push the reversed arg list */
187: for ( ; t; t = NEXT(t) ) {
188: oxarg = mknode(2,ox_pari_stream,BDY(t));
189: Pox_push_cmo(oxarg,&dmy);
190: }
191: MKSTR(name,f->name);
1.2 noro 192: STOZ(ac,narg);
1.1 noro 193: oxarg = mknode(3,ox_pari_stream,name,narg);
194: Pox_execute_function(oxarg,&dmy);
195: ox_get_pari_result = 1;
196: #if defined(VISUAL) || defined(__MINGW32__)
197: #define SM_popCMO 262
1.2 noro 198: STOZ(SM_popCMO,cmd);
1.1 noro 199: oxarg = mknode(2,ox_pari_stream,cmd);
200: Pox_push_cmd(oxarg,&dmy);
201: nd = mknode(1,ox_pari_stream);
202: MKLIST(list,nd);
203: MKReal(1.0/8,sec);
204: oxarg = mknode(2,list,sec);
205: ret=0;
206: do {
207: check_intr();
208: Pox_select(oxarg,&list);
209: oxarg = mknode(1,list);
210: Plength(oxarg,&ret);
211: }while (!ret);
212: oxarg = mknode(1,ox_pari_stream);
213: Pox_get(oxarg,&ret);
214: #else
215: oxarg = mknode(1,ox_pari_stream);
216: Pox_pop_cmo(oxarg,&ret);
217: #endif
218: ox_get_pari_result = 0;
219: if ( ret && OID(ret) == O_ERR ) {
220: char buf[BUFSIZ];
221: soutput_init(buf);
222: sprintexpr(CO,((ERR)ret)->body);
223: error(buf);
224: }
225: if ( ret && OID(ret) == O_LIST ) {
226: ret = list_to_vect(ret);
227: ret = vect_to_mat((VECT)ret);
228: }
229: return ret;
230: }
231:
232: struct pariftab {
233: char *name;
234: int dmy;
235: int type;
236: };
237:
238: /*
239: * type = 1 => argc = 1, second arg = precision
240: * type = 2 => argc = 1, second arg = (long int)0
241: *
242: */
243: /*
244: {"abs",0,1},
245: {"adj",0,1},
246: */
247:
248: struct pariftab pariftab[] = {
249: {"arg",0,1},
250: {"bigomega",0,1},
251: {"binary",0,1},
252: {"ceil",0,1},
253: {"centerlift",0,1},
254: {"cf",0,1},
255: {"classno",0,1},
256: {"classno2",0,1},
257: {"conj",0,1},
258: {"content",0,1},
259: {"denom",0,1},
260: {"det",0,1},
261: {"det2",0,1},
262: {"dilog",0,1},
263: {"disc",0,1},
264: {"discf",0,1},
265: {"divisors",0,1},
266: {"eigen",0,1},
267: {"eintg1",0,1},
268: {"erfc",0,1},
269: {"eta",0,1},
270: {"floor",0,1},
271: {"frac",0,1},
272: {"galois",0,1},
273: {"galoisconj",0,1},
274: {"gamh",0,1},
275: {"gamma",0,1},
276: {"hclassno",0,1},
277: {"hermite",0,1},
278: {"hess",0,1},
279: {"imag",0,1},
280: {"image",0,1},
281: {"image2",0,1},
282: {"indexrank",0,1},
283: {"indsort",0,1},
284: {"initalg",0,1},
285: {"isfund",0,1},
286: {"ispsp",0,1},
287: {"isqrt",0,1},
288: {"issqfree",0,1},
289: {"issquare",0,1},
290: {"jacobi",0,1},
291: {"jell",0,1},
292: {"ker",0,1},
293: {"keri",0,1},
294: {"kerint",0,1},
295: {"kerintg1",0,1},
296: {"length",0,1},
297: {"lexsort",0,1},
298: {"lift",0,1},
299: {"lindep",0,1},
300: {"lll",0,1},
301: {"lllgen",0,1},
302: {"lllgram",0,1},
303: {"lllgramgen",0,1},
304: {"lllgramint",0,1},
305: {"lllgramkerim",0,1},
306: {"lllgramkerimgen",0,1},
307: {"lllint",0,1},
308: {"lllkerim",0,1},
309: {"lllkerimgen",0,1},
310: {"lngamma",0,1},
311: {"logagm",0,1},
312: {"mat",0,1},
313: {"matrixqz2",0,1},
314: {"matrixqz3",0,1},
315: {"matsize",0,1},
316: {"modreverse",0,1},
317: {"mu",0,1},
318: {"nextprime",0,1},
319: {"norm",0,1},
320: {"norml2",0,1},
321: {"numdiv",0,1},
322: {"numer",0,1},
323: {"omega",0,1},
324: {"order",0,1},
325: {"ordred",0,1},
326: {"phi",0,1},
327: {"pnqn",0,1},
328: {"polred",0,1},
329: {"polred2",0,1},
330: {"primroot",0,1},
331: {"psi",0,1},
332: {"quadgen",0,1},
333: {"quadpoly",0,1},
334: {"real",0,1},
335: {"recip",0,1},
336: {"redreal",0,1},
337: {"regula",0,1},
338: {"reorder",0,1},
339: {"reverse",0,1},
340: {"rhoreal",0,1},
341: {"roots",0,1},
342: {"round",0,1},
343: {"sigma",0,1},
344: {"signat",0,1},
345: {"simplify",0,1},
346: {"smalldiscf",0,1},
347: {"smallfact",0,1},
348: {"smallpolred",0,1},
349: {"smallpolred2",0,1},
350: {"smith",0,1},
351: {"smith2",0,1},
352: {"sort",0,1},
353: {"sqr",0,1},
354: {"sqred",0,1},
355: {"sqrt",0,1},
356: {"supplement",0,1},
357: {"trace",0,1},
358: {"trans",0,1},
359: {"trunc",0,1},
360: {"unit",0,1},
361: {"vec",0,1},
362: {"wf",0,1},
363: {"wf2",0,1},
364: {"zeta",0,1},
365: {"factor",0,1},
366:
367: {"allocatemem",0,0},
368:
369: {"factpol",0,2},
370: {"isprime",0,2},
371: {"factorint",0,2},
372: {0,0,0},
373: };
374:
375: void parif_init() {
376: int i;
377:
378: for ( i = 0, parif = 0; pariftab[i].name; i++ )
379: appendparif(&parif,pariftab[i].name, 0,pariftab[i].type);
380: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>