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