Annotation of OpenXM_contrib2/asir2000/builtin/parif.c, Revision 1.18
1.3 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.4 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.18 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/parif.c,v 1.17 2011/08/10 04:51:57 saito Exp $
1.3 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "parse.h"
52:
1.14 ohara 53: #if defined(PARI)
1.1 noro 54: #include "genpari.h"
1.16 ohara 55: // PARI_VERSION(2,2,8) == 131592
56: #if PARI_VERSION_CODE >= 131592
57: #define mppgcd(a,b) (gcdii((a),(b)))
58: #endif
1.1 noro 59:
1.16 ohara 60: long get_pariprec();
61: void set_pariprec(long p);
1.1 noro 62:
63: void patori(GEN,Obj *);
64: void patori_i(GEN,N *);
65: void ritopa(Obj,GEN *);
66: void ritopa_i(N,int,GEN *);
67:
1.11 noro 68: void Ptodouble();
1.1 noro 69: void Peval(),Psetprec(),p_pi(),p_e(),p_mul(),p_gcd();
1.5 noro 70: void asir_cgiv(GEN);
1.1 noro 71:
1.7 saito 72: #if defined(INTERVAL) || 1
73: void Psetprecword();
74: #endif
75:
1.18 ! noro 76: #if 0
1.1 noro 77: struct ftab pari_tab[] = {
1.7 saito 78: {"eval",Peval,-2},
79: {"setprec",Psetprec,-1},
1.11 noro 80: {"todouble",Ptodouble,1},
1.7 saito 81: #if defined(INTERVAL) || 1
82: {"setprecword",Psetprecword,-1},
83: #endif
84: {0,0,0},
1.1 noro 85: };
1.18 ! noro 86: #endif
1.1 noro 87:
1.16 ohara 88: #define MKPREC(a,i,b) (argc(a)==(i)?mkprec(QTOS((Q)(b))):get_pariprec())
1.1 noro 89:
90: #define CALLPARI1(f,a,p,r)\
1.5 noro 91: ritopa((Obj)a,&_pt1_); _pt2_ = f(_pt1_,p); patori(_pt2_,r); asir_cgiv(_pt2_); asir_cgiv(_pt1_)
1.1 noro 92: #define CALLPARI2(f,a,b,p,r)\
1.5 noro 93: ritopa((Obj)a,&_pt1_); ritopa((Obj)b,&_pt2_); _pt3_ = f(_pt1_,_pt2_,p); patori(_pt3_,r); asir_cgiv(_pt3_); asir_cgiv(_pt2_); asir_cgiv(_pt1_)
1.1 noro 94:
95: #define PARIF1P(f,pf)\
96: void f(NODE,Obj *);\
97: void f(ar,rp) NODE ar; Obj *rp;\
98: { GEN _pt1_,_pt2_; CALLPARI1(pf,ARG0(ar),MKPREC(ar,2,ARG1(ar)),rp); }
99: #define PARIF2P(f,pf)\
100: void f(NODE,Obj *);\
101: void f(ar,rp) NODE ar; Obj *rp;\
102: { GEN _pt1_,_pt2_,_pt3_; CALLPARI2(pf,ARG0(ar),ARG1(ar),MKPREC(ar,3,ARG2(ar)),rp); }
103:
1.7 saito 104: #if defined(INTERVAL)
1.17 saito 105: #define PREC_CONV (0.103810253/(BYTES_IN_LONG/4)) /* log(10)/(SL*log(2)) */
1.15 ohara 106: #elif SIZEOF_LONG == 4
1.1 noro 107: #define PREC_CONV 0.103810253
1.15 ohara 108: #elif SIZEOF_LONG == 8
1.1 noro 109: #define PREC_CONV 0.051905126
110: #endif
111:
1.5 noro 112: /* XXX : we should be more careful when we free PARI pointers. */
113:
114: void asir_cgiv(ptr)
115: GEN ptr;
116: {
117: if ( ptr != gzero && ptr != gun
118: && ptr != gdeux && ptr != ghalf
119: && ptr != polvar && ptr != gi )
120: cgiv(ptr);
121: }
122:
1.1 noro 123: mkprec(p)
124: int p;
125: {
1.10 noro 126: if ( p <= 0 )
127: p = 1;
128: return (int)(p*PREC_CONV+3);
1.11 noro 129: }
130:
131: void Ptodouble(arg,rp)
132: NODE arg;
133: Num *rp;
134: {
135: double r,i;
136: Real real,imag;
137: Num num;
138:
139: asir_assert(ARG0(arg),O_N,"todouble");
140: num = (Num)ARG0(arg);
141: if ( !num ) {
142: *rp = 0;
143: return;
144: }
145: switch ( NID(num) ) {
146: case N_R: case N_Q: case N_B:
147: r = ToReal(num);
148: MKReal(r,real);
149: *rp = (Num)real;
150: break;
151: case N_C:
152: r = ToReal(((C)num)->r);
153: i = ToReal(((C)num)->i);
154: MKReal(r,real);
155: MKReal(i,imag);
156: reimtocplx((Num)real,(Num)imag,rp);
157: break;
158: default:
159: *rp = num;
160: break;
161: }
1.1 noro 162: }
163:
164: void Peval(arg,rp)
165: NODE arg;
166: Obj *rp;
167: {
168: asir_assert(ARG0(arg),O_R,"eval");
169: evalr(CO,(Obj)ARG0(arg),argc(arg)==2?QTOS((Q)ARG1(arg)):0,rp);
170: }
171:
172: void Psetprec(arg,rp)
173: NODE arg;
174: Obj *rp;
175: {
176: int p;
177: Q q;
1.16 ohara 178: long prec = get_pariprec();
1.1 noro 179:
1.7 saito 180: #if defined(INTERVAL) || 1
181: p = (int)((prec-2)/PREC_CONV); STOQ(p,q); *rp = (Obj)q;
182: if ( arg ) {
183: asir_assert(ARG0(arg),O_N,"setprec");
184: p = QTOS((Q)ARG0(arg));
185: if ( p > 0 )
186: prec = (long)(p*PREC_CONV+3);
187: }
188: #else
1.1 noro 189: p = (int)((prec-3)/PREC_CONV); STOQ(p,q); *rp = (Obj)q;
190: if ( arg ) {
191: asir_assert(ARG0(arg),O_N,"setprec");
192: prec = mkprec(QTOS((Q)ARG0(arg)));
193: }
1.7 saito 194: #endif
1.16 ohara 195: set_pariprec(prec);
1.1 noro 196: }
1.7 saito 197:
198: #if defined(INTERVAL) || 1
199: void Psetprecword(arg,rp)
200: NODE arg;
201: Obj *rp;
202: {
203: int p;
204: Q q;
1.16 ohara 205: long prec = get_pariprec();
1.7 saito 206:
207: p = (int)((prec-2)); STOQ(p,q); *rp = (Obj)q;
208: if ( arg ) {
209: asir_assert(ARG0(arg),O_N,"setprecword");
210: p = QTOS((Q)ARG0(arg));
211: if ( p > 0 ) {
212: prec = p + 2;
213: }
214: }
1.16 ohara 215: set_pariprec(prec);
1.7 saito 216: }
217: #endif
1.1 noro 218:
219: void p_pi(arg,rp)
220: NODE arg;
221: Obj *rp;
222: {
223: GEN x;
224:
225: x = mppi(MKPREC(arg,1,ARG0(arg)));
1.5 noro 226: patori(x,rp); asir_cgiv(x);
1.1 noro 227: }
228:
229: void p_e(arg,rp)
230: NODE arg;
231: Obj *rp;
232: {
233: GEN x;
234:
1.5 noro 235: x = gexp(gun,MKPREC(arg,1,ARG0(arg))); patori(x,rp); asir_cgiv(x);
1.1 noro 236: }
237:
238: void p_mul(a,b,r)
239: Obj a,b,*r;
240: {
241: GEN p1,p2,p3;
242:
243: ritopa((Obj)a,&p1); ritopa((Obj)b,&p2);
244: p3 = mulii(p1,p2);
1.5 noro 245: patori(p3,r); asir_cgiv(p3); asir_cgiv(p2); asir_cgiv(p1);
1.1 noro 246: }
247:
248: void p_gcd(a,b,r)
249: N a,b,*r;
250: {
251: GEN p1,p2,p3;
252:
253: ritopa_i(a,1,&p1); ritopa_i(b,1,&p2);
254: p3 = mppgcd(p1,p2);
1.5 noro 255: patori_i(p3,r); asir_cgiv(p3); asir_cgiv(p2); asir_cgiv(p1);
1.1 noro 256: }
257:
258: PARIF1P(p_sin,gsin) PARIF1P(p_cos,gcos) PARIF1P(p_tan,gtan)
259: PARIF1P(p_asin,gasin) PARIF1P(p_acos,gacos) PARIF1P(p_atan,gatan)
260: PARIF1P(p_sinh,gsh) PARIF1P(p_cosh,gch) PARIF1P(p_tanh,gth)
261: PARIF1P(p_asinh,gash) PARIF1P(p_acosh,gach) PARIF1P(p_atanh,gath)
262: PARIF1P(p_exp,gexp) PARIF1P(p_log,glog)
263: PARIF1P(p_dilog,dilog) PARIF1P(p_erf,gerfc)
264: PARIF1P(p_eigen,eigen) PARIF1P(p_roots,roots)
265:
266: PARIF2P(p_pow,gpui)
267:
268: pointer evalparif(f,arg)
269: FUNC f;
270: NODE arg;
271: {
272: GEN a,v;
273: long ltop,lbot;
274: pointer r;
1.9 noro 275: int ac,opt,intarg,ret;
1.1 noro 276: char buf[BUFSIZ];
1.9 noro 277: Q q;
1.8 noro 278: GEN (*dmy)();
1.1 noro 279:
280: if ( !f->f.binf ) {
281: sprintf(buf,"pari : %s undefined.",f->name);
282: error(buf);
1.10 noro 283: /* NOTREACHED */
284: return 0;
1.1 noro 285: }
286: switch ( f->type ) {
1.9 noro 287: case 0: /* in/out : integer */
288: ac = argc(arg);
289: if ( ac > 2 ) {
290: fprintf(stderr,"argument mismatch in %s()\n",NAME(f));
291: error("");
1.10 noro 292: /* NOTREACHED */
293: return 0;
1.9 noro 294: }
295: intarg = !ac ? 0 : QTOS((Q)ARG0(arg));
296: dmy = (GEN (*)())f->f.binf;
1.10 noro 297: ret = (int)(*dmy)(intarg);
1.9 noro 298: STOQ(ret,q);
299: return (pointer)q;
300:
1.1 noro 301: case 1:
302: ac = argc(arg);
303: if ( !ac || ( ac > 2 ) ) {
304: fprintf(stderr,"argument mismatch in %s()\n",NAME(f));
305: error("");
1.10 noro 306: /* NOTREACHED */
307: return 0;
1.1 noro 308: }
309: ltop = avma;
310: ritopa((Obj)ARG0(arg),&a);
311: dmy = (GEN (*)())f->f.binf;
312: v = (*dmy)(a,MKPREC(arg,2,ARG1(arg)));
313: lbot = avma;
314: patori(v,(Obj *)&r); gerepile(ltop,lbot,0);
315: return r;
1.8 noro 316:
317: case 2:
318: ac = argc(arg);
319: if ( !ac || ( ac > 2 ) ) {
320: fprintf(stderr,"argument mismatch in %s()\n",NAME(f));
321: error("");
1.10 noro 322: /* NOTREACHED */
323: return 0;
1.8 noro 324: }
325: if ( ac == 1 )
326: opt = 0;
327: else
328: opt = QTOS((Q)ARG1(arg));
329: ltop = avma;
330: ritopa((Obj)ARG0(arg),&a);
331: dmy = (GEN (*)())f->f.binf;
332: v = (*dmy)(a,opt);
333: lbot = avma;
334: patori(v,(Obj *)&r); gerepile(ltop,lbot,0);
335: return r;
336:
1.1 noro 337: default:
338: error("evalparif : not implemented yet.");
1.10 noro 339: /* NOTREACHED */
340: return 0;
1.1 noro 341: }
342: }
343:
344: struct pariftab {
345: char *name;
346: GEN (*f)();
347: int type;
348: };
349:
1.8 noro 350: /*
351: * type = 1 => argc = 1, second arg = precision
352: * type = 2 => argc = 1, second arg = optional (long int)
353: *
354: */
355:
1.1 noro 356: struct pariftab pariftab[] = {
1.10 noro 357: {"allocatemem",(GEN(*)())allocatemoremem,0},
1.1 noro 358: {"abs",(GEN (*)())gabs,1},
359: {"adj",adj,1},
360: {"arg",garg,1},
361: {"bigomega",gbigomega,1},
362: {"binary",binaire,1},
363: {"ceil",gceil,1},
364: {"centerlift",centerlift,1},
365: {"cf",gcf,1},
366: {"classno",classno,1},
367: {"classno2",classno2,1},
368: {"conj",gconj,1},
369: {"content",content,1},
370: {"denom",denom,1},
371: {"det",det,1},
372: {"det2",det2,1},
373: {"dilog",dilog,1},
374: {"disc",discsr,1},
375: {"discf",discf,1},
376: {"divisors",divisors,1},
377: {"eigen",eigen,1},
378: {"eintg1",eint1,1},
379: {"erfc",gerfc,1},
380: {"eta",eta,1},
381: {"floor",gfloor,1},
382: {"frac",gfrac,1},
383: {"galois",galois,1},
384: {"galoisconj",galoisconj,1},
385: {"gamh",ggamd,1},
386: {"gamma",ggamma,1},
387: {"hclassno",classno3,1},
388: {"hermite",hnf,1},
389: {"hess",hess,1},
390: {"imag",gimag,1},
391: {"image",image,1},
392: {"image2",image2,1},
393: {"indexrank",indexrank,1},
394: {"indsort",indexsort,1},
395: {"initalg",initalg,1},
396: {"isfund",gisfundamental,1},
1.13 noro 397: {"isprime",gisprime,2},
1.1 noro 398: {"ispsp",gispsp,1},
399: {"isqrt",racine,1},
400: {"issqfree",gissquarefree,1},
401: {"issquare",gcarreparfait,1},
402: {"jacobi",jacobi,1},
403: {"jell",jell,1},
404: {"ker",ker,1},
405: {"keri",keri,1},
406: {"kerint",kerint,1},
407: {"kerintg1",kerint1,1},
408: {"length",(GEN(*)())glength,1},
409: {"lexsort",lexsort,1},
410: {"lift",lift,1},
411: {"lindep",lindep,1},
412: {"lll",lll,1},
413: {"lllgen",lllgen,1},
414: {"lllgram",lllgram,1},
415: {"lllgramgen",lllgramgen,1},
416: {"lllgramint",lllgramint,1},
417: {"lllgramkerim",lllgramkerim,1},
418: {"lllgramkerimgen",lllgramkerimgen,1},
419: {"lllint",lllint,1},
420: {"lllkerim",lllkerim,1},
421: {"lllkerimgen",lllkerimgen,1},
422: {"lngamma",glngamma,1},
423: {"logagm",glogagm,1},
424: {"mat",gtomat,1},
425: {"matrixqz2",matrixqz2,1},
426: {"matrixqz3",matrixqz3,1},
427: {"matsize",matsize,1},
428: {"modreverse",polymodrecip,1},
429: {"mu",gmu,1},
430: {"nextprime",nextprime,1},
431: {"norm",gnorm,1},
432: {"norml2",gnorml2,1},
433: {"numdiv",numbdiv,1},
434: {"numer",numer,1},
435: {"omega",gomega,1},
436: {"order",order,1},
437: {"ordred",ordred,1},
438: {"phi",phi,1},
439: {"pnqn",pnqn,1},
440: {"polred",polred,1},
441: {"polred2",polred2,1},
442: {"primroot",gener,1},
443: {"psi",gpsi,1},
444: {"quadgen",quadgen ,1},
445: {"quadpoly",quadpoly ,1},
446: {"real",greal,1},
447: {"recip",polrecip ,1},
448: {"redreal",redreal ,1},
449: {"regula",regula ,1},
450: {"reorder",reorder ,1},
451: {"reverse",recip ,1},
452: {"rhoreal",rhoreal ,1},
453: {"roots",roots,1},
454: {"round",ground,1},
455: {"sigma",sumdiv,1},
456: {"signat",signat,1},
457: {"simplify",simplify,1},
458: {"smalldiscf",smalldiscf,1},
459: {"smallfact",smallfact,1},
460: {"smallpolred",smallpolred,1},
461: {"smallpolred2",smallpolred2,1},
462: {"smith",smith,1},
463: {"smith2",smith2,1},
464: {"sort",sort,1},
465: {"sqr",gsqr,1},
466: {"sqred",sqred,1},
467: {"sqrt",gsqrt,1},
468: {"supplement",suppl,1},
469: {"trace",gtrace,1},
470: {"trans",gtrans,1},
471: {"trunc",gtrunc,1},
472: {"unit",fundunit,1},
473: {"vec",gtovec,1},
474: {"wf",wf,1},
475: {"wf2",wf2,1},
476: {"zeta",gzeta,1},
1.9 noro 477: {"factor",factor,1},
1.8 noro 478: {"factorint",factorint,2},
1.1 noro 479: {0,0,0},
480: };
481:
482: void parif_init() {
483: int i;
484:
485: for ( i = 0, parif = 0; pariftab[i].name; i++ )
486: appendparif(&parif,pariftab[i].name, (int (*)())pariftab[i].f,pariftab[i].type);
487: }
488: #else /* PARI */
489:
490: struct ftab pari_tab[] = {
491: {0,0,0},
492: };
493:
494: void parif_init() {}
495:
1.2 noro 496: pointer evalparif(f,arg)
497: FUNC f;
498: NODE arg;
499: {
1.1 noro 500: error("evalparif : PARI is not combined.");
501: }
502: #endif /*PARI */
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to parif.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin