![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pari.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to pari.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.2 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.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 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.11 ! ohara 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/pari.c,v 1.10 2009/09/10 01:51:54 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
1.9 ohara 51:
1.6 ohara 52: #if defined(PARI)
1.1 noro 53: #include "base.h"
54: #include <math.h>
55: #include "genpari.h"
56:
57: void patori(GEN,Obj *);
58: void patori_i(GEN,N *);
59: void ritopa(Obj,GEN *);
60: void ritopa_i(N,int,GEN *);
61:
1.9 ohara 62: // PARI_VERSION(2,2,12) == 131596
63: #if PARI_VERSION_CODE >= 131596
64: #define prec precreal
65: #endif
1.1 noro 66: extern long prec;
67: extern int paristack;
68:
1.9 ohara 69: long get_pariprec() {
70: return prec;
71: }
72: void set_pariprec(long p) {
73: prec = p;
74: }
75:
1.1 noro 76: void risa_pari_init() {
77: pari_init(paristack,2);
1.9 ohara 78: set_pariprec(4);
1.1 noro 79: }
80:
81: void create_pari_variable(index)
82: int index;
83: {
84: static int max_varn;
85: int i;
86: char name[BUFSIZ];
87:
88: if ( index > max_varn ) {
89: for ( i = max_varn+1; i <= index; i++ ) {
90: sprintf(name,"x%d",i);
1.7 saito 91: #if (PARI_VERSION_CODE < 131594)
1.1 noro 92: fetch_named_var(name,0);
1.7 saito 93: #else
94: fetch_named_var(name);
95: #endif
1.1 noro 96: }
97: max_varn = index;
98: }
99: }
100:
101: int get_lg(a)
102: GEN a;
103: {
104: return lg(a);
105: }
106:
1.9 ohara 107: void gpui_ri(Obj a, Obj e, Obj *c)
108: {
109: GEN pa,pe,z;
1.10 noro 110: long ltop,lbot;
1.9 ohara 111:
112: ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;
113: z = gerepile(ltop,lbot,gpui(pa,pe,get_pariprec()));
114: patori(z,c); cgiv(z);
115: }
116:
1.1 noro 117: void ritopa(a,rp)
118: Obj a;
119: GEN *rp;
120: {
121: long ltop;
1.11 ! ohara 122: GEN pnm,z,w,u;
1.1 noro 123: DCP dc;
124: int i,j,l,row,col;
125: VL vl;
126: V v;
127:
128: if ( !a ) {
129: *rp = gzero; return;
130: }
131: switch ( OID(a) ) {
132: case O_N:
133: switch ( NID(a) ) {
134: case N_Q:
135: ltop = avma; ritopa_i(NM((Q)a),SGN((Q)a),&pnm);
136: if ( INT((Q)a) )
137: *rp = pnm;
138: else {
139: *rp = z = cgetg(3,4); z[1] = (long)pnm;
1.11 ! ohara 140: ritopa_i(DN((Q)a),1,&u); z[2] = u;
1.1 noro 141: }
142: break;
143: case N_R:
144: *rp = dbltor(BDY((Real)a)); break;
145: case N_B:
146: *rp = gcopy((GEN)BDY(((BF)a))); break;
147: case N_C:
148: z = cgetg(3,6);
1.11 ! ohara 149: ritopa((Obj)((C)a)->r,&u); z[1] = u;
! 150: ritopa((Obj)((C)a)->i,&u); z[2] = u;
1.1 noro 151: *rp = z;
152: break;
153: default:
154: error("ritopa : not implemented"); break;
155: }
156: break;
157: case O_P:
158: l = UDEG((P)a); *rp = z = cgetg(l+3,10);
159: setsigne(z,1);
160: for ( i = 0, vl = CO, v = VR((P)a); vl->v != v;
161: vl = NEXT(vl), i++ );
162: create_pari_variable(i);
163: setvarn(z,i);
164: setlgef(z,l+3);
165: for ( i = l+2; i >= 2; i-- )
166: z[i] = (long)gzero;
1.11 ! ohara 167: for ( dc = DC((P)a); dc; dc = NEXT(dc) ) {
! 168: ritopa((Obj)COEF(dc),&u); z[QTOS(DEG(dc))+2] = u;
! 169: }
1.1 noro 170: break;
171: case O_VECT:
172: l = ((VECT)a)->len; z = cgetg(l+1,17);
1.11 ! ohara 173: for ( i = 0; i < l; i++ ) {
! 174: ritopa((Obj)BDY((VECT)a)[i],&u); z[i+1] = u;
! 175: }
1.1 noro 176: *rp = z;
177: break;
178: case O_MAT:
179: row = ((MAT)a)->row; col = ((MAT)a)->col; z = cgetg(col+1,19);
180: for ( j = 0; j < col; j++ ) {
181: w = cgetg(row+1,18);
1.11 ! ohara 182: for ( i = 0; i < row; i++ ) {
! 183: ritopa((Obj)BDY((MAT)a)[i][j],&u); w[i+1] = u;
! 184: }
1.1 noro 185: z[j+1] = (long)w;
186: }
187: *rp = z;
188: break;
189: default:
190: error("ritopa : not implemented"); break;
191: }
192: }
193:
194: void patori(a,rp)
195: GEN a;
196: Obj *rp;
197: {
198: Q q,qnm,qdn;
199: BF r;
200: C c;
201: VECT v;
202: MAT m;
203: N n,nm,dn;
204: DCP dc0,dc;
205: P t;
206: int s,i,j,l,row,col;
207: GEN b;
208: VL vl;
209:
210: if ( gcmp0(a) )
211: *rp = 0;
212: else {
213: switch ( typ(a) ) {
214: case 1: /* integer */
215: patori_i(a,&n); NTOQ(n,(char)signe(a),q); *rp = (Obj)q;
216: break;
217: case 2: /* real */
218: NEWBF(r,lg(a)); bcopy((char *)a,(char *)BDY(r),lg(a)*sizeof(long));
219: *rp = (Obj)r;
220: break;
221: case 4: /* rational; reduced */
222: patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);
223: s = signe(a[1])*signe(a[2]);
224: if ( UNIN(dn) )
225: NTOQ(nm,s,q);
226: else
227: NDTOQ(nm,dn,s,q);
228: *rp = (Obj)q;
229: break;
230: case 5: /* rational; not necessarily reduced */
231: patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);
232: s = signe(a[1])*signe(a[2]);
233: NTOQ(nm,s,qnm); NTOQ(dn,1,qdn); divq(qnm,qdn,(Q *)rp);
234: break;
235: case 6: /* complex */
236: if ( gcmp0((GEN)a[2]) )
237: patori((GEN)a[1],rp);
238: else {
239: NEWC(c); patori((GEN)a[1],(Obj *)&c->r); patori((GEN)a[2],(Obj *)&c->i);
240: *rp = (Obj)c;
241: }
242: break;
243: case 10: /* polynomial */
244: for ( i = lgef(a)-1, dc0 = 0; i >= 2; i-- )
245: if ( !gcmp0((GEN)a[i]) ) {
246: NEXTDC(dc0,dc);
247: patori((GEN)a[i],(Obj *)&COEF(dc)); STOQ(i-2,DEG(dc));
248: }
249: if ( !dc0 )
250: *rp = 0;
251: else {
252: /* assuming that CO has not changed. */
253: for ( s = varn(a), i = 0, vl = CO; i != s;
254: i++, vl = NEXT(vl) );
255: NEXT(dc) = 0; MKP(vl->v,dc0,t); *rp = (Obj)t;
256: }
257: break;
258: case 17: /* row vector */
259: case 18: /* column vector */
260: l = lg(a)-1; MKVECT(v,l);
261: for ( i = 0; i < l; i++ )
262: patori((GEN)a[i+1],(Obj *)&BDY(v)[i]);
263: *rp = (Obj)v;
264: break;
265: case 19: /* matrix */
266: col = lg(a)-1; row = lg(a[1])-1; MKMAT(m,row,col);
267: for ( j = 0; j < col; j++ )
268: for ( i = 0, b = (GEN)a[j+1]; i < row; i++ )
269: patori((GEN)b[i+1],(Obj *)&BDY(m)[i][j]);
270: *rp = (Obj)m;
271: break;
272: default:
273: error("patori : not implemented");
274: break;
275: }
276: }
277: }
278:
1.8 ohara 279: #if SIZEOF_LONG == 4
1.1 noro 280: void ritopa_i(a,s,rp)
281: N a;
282: int s;
283: GEN *rp;
284: {
285: int j,l;
286: unsigned int *b;
287: GEN z;
288:
289: l = PL(a); b = (unsigned int *)BD(a);
290: z = cgeti(l+2);
291: bzero((char *)&z[2],l*sizeof(int));
292: for ( j = 0; j < l; j++ )
293: z[l-j+1] = b[j];
294: s = s>0?1:-1;
295: setsigne(z,s);
296: setlgefint(z,lg(z));
297: *rp = z;
298: }
299:
300: void patori_i(g,rp)
301: GEN g;
302: N *rp;
303: {
304: int j,l;
305: unsigned int *a,*b;
306: N z;
307:
308: l = lgef(g)-2;
309: a = (unsigned int *)g;
310: *rp = z = NALLOC(l); PL(z) = l;
311: for ( j = 0, b = (unsigned int *)BD(z); j < l; j++ )
312: b[l-j-1] = ((unsigned int *)g)[j+2];
313: }
314:
1.8 ohara 315: #elif SIZEOF_LONG == 8
1.1 noro 316: void ritopa_i(a,s,rp)
317: N a;
318: int s;
319: GEN *rp;
320: {
321: int j,l,words;
322: unsigned int *b;
323: GEN z;
324:
325: l = PL(a); b = BD(a);
326: words = (l+1)/2;
327: z = cgeti(words+2);
328: bzero((char *)&z[2],words*sizeof(long));
329: for ( j = 0; j < words; j++ )
330: z[words-j+1] = ((unsigned long)b[2*j])
331: |(((unsigned long)(2*j+1<l?b[2*j+1]:0))<<32);
332: s = s>0?1:-1;
333: setsigne(z,s);
334: setlgefint(z,lg(z));
335: *rp = z;
336: }
337:
338: void patori_i(g,rp)
339: GEN g;
340: N *rp;
341: {
342: int j,l,words;
343: unsigned long t;
344: unsigned long *a;
345: unsigned int *b;
346: N z;
347:
348: words = lgef(g)-2;
349: l = 2*words;
350: a = (unsigned long *)g;
351: *rp = z = NALLOC(l); PL(z) = l;
352: for ( j = 0, b = BD(z); j < words; j++ ) {
353: t = a[words+1-j];
354: b[2*j] = t&0xffffffff;
355: b[2*j+1] = t>>32;
356: }
357: PL(z) = b[l-1] ? l : l-1;
358: }
359: #endif
360:
361: void strtobf(s,rp)
362: char *s;
363: BF *rp;
364: {
365: GEN z;
366:
367: z = lisexpr(s); patori(z,(Obj *)rp); cgiv(z);
368: }
369: #endif