![[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