Annotation of OpenXM_contrib2/asir2018/engine/pari.c, Revision 1.1
1.1 ! 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@sec.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$
! 49: */
! 50: #include "ca.h"
! 51:
! 52: #if defined(PARI)
! 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:
! 62: // PARI_VERSION(2,2,12) == 131596
! 63: #if PARI_VERSION_CODE >= 131596
! 64: #define prec precreal
! 65: #endif
! 66: extern long prec;
! 67: extern int paristack;
! 68:
! 69: long get_pariprec() {
! 70: return prec;
! 71: }
! 72: void set_pariprec(long p) {
! 73: prec = p;
! 74: }
! 75:
! 76: void risa_pari_init() {
! 77: pari_init(paristack,2);
! 78: set_pariprec(4);
! 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);
! 91: #if (PARI_VERSION_CODE < 131594)
! 92: fetch_named_var(name,0);
! 93: #else
! 94: fetch_named_var(name);
! 95: #endif
! 96: }
! 97: max_varn = index;
! 98: }
! 99: }
! 100:
! 101: int get_lg(a)
! 102: GEN a;
! 103: {
! 104: return lg(a);
! 105: }
! 106:
! 107: void gpui_ri(Obj a, Obj e, Obj *c)
! 108: {
! 109: GEN pa,pe,z;
! 110: long ltop,lbot;
! 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:
! 117: void ritopa(a,rp)
! 118: Obj a;
! 119: GEN *rp;
! 120: {
! 121: long ltop;
! 122: GEN pnm,z,w,u;
! 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;
! 140: ritopa_i(DN((Q)a),1,&u); z[2] = u;
! 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);
! 149: ritopa((Obj)((C)a)->r,&u); z[1] = u;
! 150: ritopa((Obj)((C)a)->i,&u); z[2] = u;
! 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;
! 167: for ( dc = DC((P)a); dc; dc = NEXT(dc) ) {
! 168: ritopa((Obj)COEF(dc),&u); z[QTOS(DEG(dc))+2] = u;
! 169: }
! 170: break;
! 171: case O_VECT:
! 172: l = ((VECT)a)->len; z = cgetg(l+1,17);
! 173: for ( i = 0; i < l; i++ ) {
! 174: ritopa((Obj)BDY((VECT)a)[i],&u); z[i+1] = u;
! 175: }
! 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);
! 182: for ( i = 0; i < row; i++ ) {
! 183: ritopa((Obj)BDY((MAT)a)[i][j],&u); w[i+1] = u;
! 184: }
! 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,ty;
! 207: GEN b;
! 208: VL vl;
! 209:
! 210: if ( gcmp0(a) )
! 211: *rp = 0;
! 212: else {
! 213: switch ( ty = 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 3: /* integermod */
! 222: MKVECT(v,2);
! 223: patori((GEN)a[1],(Obj *)&BDY(v)[0]);
! 224: patori((GEN)a[2],(Obj *)&BDY(v)[1]);
! 225: *rp = (Obj)v;
! 226: break;
! 227: case 4: /* rational; reduced */
! 228: patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);
! 229: s = signe(a[1])*signe(a[2]);
! 230: if ( UNIN(dn) )
! 231: NTOQ(nm,s,q);
! 232: else
! 233: NDTOQ(nm,dn,s,q);
! 234: *rp = (Obj)q;
! 235: break;
! 236: case 5: /* rational; not necessarily reduced */
! 237: patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);
! 238: s = signe(a[1])*signe(a[2]);
! 239: NTOQ(nm,s,qnm); NTOQ(dn,1,qdn); divq(qnm,qdn,(Q *)rp);
! 240: break;
! 241: case 6: /* complex */
! 242: if ( gcmp0((GEN)a[2]) )
! 243: patori((GEN)a[1],rp);
! 244: else {
! 245: NEWC(c); patori((GEN)a[1],(Obj *)&c->r); patori((GEN)a[2],(Obj *)&c->i);
! 246: *rp = (Obj)c;
! 247: }
! 248: break;
! 249: case 10: /* polynomial */
! 250: for ( i = lgef(a)-1, dc0 = 0; i >= 2; i-- )
! 251: if ( !gcmp0((GEN)a[i]) ) {
! 252: NEXTDC(dc0,dc);
! 253: patori((GEN)a[i],(Obj *)&COEF(dc)); STOQ(i-2,DEG(dc));
! 254: }
! 255: if ( !dc0 )
! 256: *rp = 0;
! 257: else {
! 258: /* assuming that CO has not changed. */
! 259: for ( s = varn(a), i = 0, vl = CO; i != s;
! 260: i++, vl = NEXT(vl) );
! 261: NEXT(dc) = 0; MKP(vl->v,dc0,t); *rp = (Obj)t;
! 262: }
! 263: break;
! 264: case 17: /* row vector */
! 265: case 18: /* column vector */
! 266: l = lg(a)-1; MKVECT(v,l);
! 267: for ( i = 0; i < l; i++ )
! 268: patori((GEN)a[i+1],(Obj *)&BDY(v)[i]);
! 269: *rp = (Obj)v;
! 270: break;
! 271: case 19: /* matrix */
! 272: col = lg(a)-1; row = lg(a[1])-1; MKMAT(m,row,col);
! 273: for ( j = 0; j < col; j++ )
! 274: for ( i = 0, b = (GEN)a[j+1]; i < row; i++ )
! 275: patori((GEN)b[i+1],(Obj *)&BDY(m)[i][j]);
! 276: *rp = (Obj)m;
! 277: break;
! 278: default:
! 279: error("patori : not implemented");
! 280: break;
! 281: }
! 282: }
! 283: }
! 284:
! 285: #if SIZEOF_LONG == 4
! 286: void ritopa_i(a,s,rp)
! 287: N a;
! 288: int s;
! 289: GEN *rp;
! 290: {
! 291: int j,l;
! 292: unsigned int *b;
! 293: GEN z;
! 294:
! 295: l = PL(a); b = (unsigned int *)BD(a);
! 296: z = cgeti(l+2);
! 297: bzero((char *)&z[2],l*sizeof(int));
! 298: for ( j = 0; j < l; j++ )
! 299: z[l-j+1] = b[j];
! 300: s = s>0?1:-1;
! 301: setsigne(z,s);
! 302: setlgefint(z,lg(z));
! 303: *rp = z;
! 304: }
! 305:
! 306: void patori_i(g,rp)
! 307: GEN g;
! 308: N *rp;
! 309: {
! 310: int j,l;
! 311: unsigned int *a,*b;
! 312: N z;
! 313:
! 314: l = lgef(g)-2;
! 315: a = (unsigned int *)g;
! 316: *rp = z = NALLOC(l); PL(z) = l;
! 317: for ( j = 0, b = (unsigned int *)BD(z); j < l; j++ )
! 318: b[l-j-1] = ((unsigned int *)g)[j+2];
! 319: }
! 320:
! 321: #elif SIZEOF_LONG == 8
! 322: void ritopa_i(a,s,rp)
! 323: N a;
! 324: int s;
! 325: GEN *rp;
! 326: {
! 327: int j,l,words;
! 328: unsigned int *b;
! 329: GEN z;
! 330:
! 331: l = PL(a); b = BD(a);
! 332: words = (l+1)/2;
! 333: z = cgeti(words+2);
! 334: bzero((char *)&z[2],words*sizeof(long));
! 335: for ( j = 0; j < words; j++ )
! 336: z[words-j+1] = ((unsigned long)b[2*j])
! 337: |(((unsigned long)(2*j+1<l?b[2*j+1]:0))<<32);
! 338: s = s>0?1:-1;
! 339: setsigne(z,s);
! 340: setlgefint(z,lg(z));
! 341: *rp = z;
! 342: }
! 343:
! 344: void patori_i(g,rp)
! 345: GEN g;
! 346: N *rp;
! 347: {
! 348: int j,l,words;
! 349: unsigned long t;
! 350: unsigned long *a;
! 351: unsigned int *b;
! 352: N z;
! 353:
! 354: words = lgef(g)-2;
! 355: l = 2*words;
! 356: a = (unsigned long *)g;
! 357: *rp = z = NALLOC(l); PL(z) = l;
! 358: for ( j = 0, b = BD(z); j < words; j++ ) {
! 359: t = a[words+1-j];
! 360: b[2*j] = t&0xffffffff;
! 361: b[2*j+1] = t>>32;
! 362: }
! 363: PL(z) = b[l-1] ? l : l-1;
! 364: }
! 365: #endif
! 366:
! 367: void strtobf(s,rp)
! 368: char *s;
! 369: BF *rp;
! 370: {
! 371: GEN z;
! 372:
! 373: z = lisexpr(s); patori(z,(Obj *)rp); cgiv(z);
! 374: }
! 375: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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 / asir2018 / engine