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