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