File: [local] / OpenXM_contrib2 / asir2000 / engine / pari.c (download)
Revision 1.1.1.1 (vendor branch), Fri Dec 3 07:39:08 1999 UTC (24 years, 9 months ago) by noro
Branch: NORO
CVS Tags: RELEASE_20000124, RELEASE_1_1_2, ASIR2000 Changes since 1.1: +0 -0
lines
Imported asir2000 as OpenXM_contrib2/asir2000.
|
/* $OpenXM: OpenXM_contrib2/asir2000/engine/pari.c,v 1.1.1.1 1999/12/03 07:39:08 noro Exp $ */
#include "ca.h"
#if PARI
#include "base.h"
#include <math.h>
#include "genpari.h"
#if defined(THINK_C)
void patori(GEN,Obj *);
void patori_i(GEN,N *);
void ritopa(Obj,GEN *);
void ritopa_i(N,int,GEN *);
#else
void patori();
void patori_i();
void ritopa();
void ritopa_i();
#endif
extern long prec;
extern int paristack;
void risa_pari_init() {
char buf[BUFSIZ];
int i;
pari_init(paristack,2);
prec = 4;
}
void create_pari_variable(index)
int index;
{
static int max_varn;
int i;
char name[BUFSIZ];
if ( index > max_varn ) {
for ( i = max_varn+1; i <= index; i++ ) {
sprintf(name,"x%d",i);
fetch_named_var(name,0);
}
max_varn = index;
}
}
int get_lg(a)
GEN a;
{
return lg(a);
}
void ritopa(a,rp)
Obj a;
GEN *rp;
{
long ltop;
GEN pnm,z,w;
DCP dc;
int i,j,l,row,col;
VL vl;
V v;
if ( !a ) {
*rp = gzero; return;
}
switch ( OID(a) ) {
case O_N:
switch ( NID(a) ) {
case N_Q:
ltop = avma; ritopa_i(NM((Q)a),SGN((Q)a),&pnm);
if ( INT((Q)a) )
*rp = pnm;
else {
*rp = z = cgetg(3,4); z[1] = (long)pnm;
ritopa_i(DN((Q)a),1,(GEN *)&z[2]);
}
break;
case N_R:
*rp = dbltor(BDY((Real)a)); break;
case N_B:
*rp = gcopy((GEN)BDY(((BF)a))); break;
case N_C:
z = cgetg(3,6);
ritopa((Obj)((C)a)->r,(GEN *)&z[1]); ritopa((Obj)((C)a)->i,(GEN *)&z[2]);
*rp = z;
break;
default:
error("ritopa : not implemented"); break;
}
break;
case O_P:
l = UDEG((P)a); *rp = z = cgetg(l+3,10);
setsigne(z,1);
for ( i = 0, vl = CO, v = VR((P)a); vl->v != v;
vl = NEXT(vl), i++ );
create_pari_variable(i);
setvarn(z,i);
setlgef(z,l+3);
for ( i = l+2; i >= 2; i-- )
z[i] = (long)gzero;
for ( dc = DC((P)a); dc; dc = NEXT(dc) )
ritopa((Obj)COEF(dc),(GEN *)&z[QTOS(DEG(dc))+2]);
break;
case O_VECT:
l = ((VECT)a)->len; z = cgetg(l+1,17);
for ( i = 0; i < l; i++ )
ritopa((Obj)BDY((VECT)a)[i],(GEN *)&z[i+1]);
*rp = z;
break;
case O_MAT:
row = ((MAT)a)->row; col = ((MAT)a)->col; z = cgetg(col+1,19);
for ( j = 0; j < col; j++ ) {
w = cgetg(row+1,18);
for ( i = 0; i < row; i++ )
ritopa((Obj)BDY((MAT)a)[i][j],(GEN *)&w[i+1]);
z[j+1] = (long)w;
}
*rp = z;
break;
default:
error("ritopa : not implemented"); break;
}
}
void patori(a,rp)
GEN a;
Obj *rp;
{
Q q,qnm,qdn;
BF r;
C c;
VECT v;
MAT m;
N n,nm,dn;
DCP dc0,dc;
P t;
int s,i,j,l,row,col;
GEN b;
VL vl;
if ( gcmp0(a) )
*rp = 0;
else {
switch ( typ(a) ) {
case 1: /* integer */
patori_i(a,&n); NTOQ(n,(char)signe(a),q); *rp = (Obj)q;
break;
case 2: /* real */
NEWBF(r,lg(a)); bcopy((char *)a,(char *)BDY(r),lg(a)*sizeof(long));
*rp = (Obj)r;
break;
case 4: /* rational; reduced */
patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);
s = signe(a[1])*signe(a[2]);
if ( UNIN(dn) )
NTOQ(nm,s,q);
else
NDTOQ(nm,dn,s,q);
*rp = (Obj)q;
break;
case 5: /* rational; not necessarily reduced */
patori_i((GEN)a[1],&nm); patori_i((GEN)a[2],&dn);
s = signe(a[1])*signe(a[2]);
NTOQ(nm,s,qnm); NTOQ(dn,1,qdn); divq(qnm,qdn,(Q *)rp);
break;
case 6: /* complex */
if ( gcmp0((GEN)a[2]) )
patori((GEN)a[1],rp);
else {
NEWC(c); patori((GEN)a[1],(Obj *)&c->r); patori((GEN)a[2],(Obj *)&c->i);
*rp = (Obj)c;
}
break;
case 10: /* polynomial */
for ( i = lgef(a)-1, dc0 = 0; i >= 2; i-- )
if ( !gcmp0((GEN)a[i]) ) {
NEXTDC(dc0,dc);
patori((GEN)a[i],(Obj *)&COEF(dc)); STOQ(i-2,DEG(dc));
}
if ( !dc0 )
*rp = 0;
else {
/* assuming that CO has not changed. */
for ( s = varn(a), i = 0, vl = CO; i != s;
i++, vl = NEXT(vl) );
NEXT(dc) = 0; MKP(vl->v,dc0,t); *rp = (Obj)t;
}
break;
case 17: /* row vector */
case 18: /* column vector */
l = lg(a)-1; MKVECT(v,l);
for ( i = 0; i < l; i++ )
patori((GEN)a[i+1],(Obj *)&BDY(v)[i]);
*rp = (Obj)v;
break;
case 19: /* matrix */
col = lg(a)-1; row = lg(a[1])-1; MKMAT(m,row,col);
for ( j = 0; j < col; j++ )
for ( i = 0, b = (GEN)a[j+1]; i < row; i++ )
patori((GEN)b[i+1],(Obj *)&BDY(m)[i][j]);
*rp = (Obj)m;
break;
default:
error("patori : not implemented");
break;
}
}
}
#if defined(LONG_IS_32BIT)
void ritopa_i(a,s,rp)
N a;
int s;
GEN *rp;
{
int j,l;
unsigned int *b;
GEN z;
l = PL(a); b = (unsigned int *)BD(a);
z = cgeti(l+2);
bzero((char *)&z[2],l*sizeof(int));
for ( j = 0; j < l; j++ )
z[l-j+1] = b[j];
s = s>0?1:-1;
setsigne(z,s);
setlgefint(z,lg(z));
*rp = z;
}
void patori_i(g,rp)
GEN g;
N *rp;
{
int j,l;
unsigned int *a,*b;
N z;
l = lgef(g)-2;
a = (unsigned int *)g;
*rp = z = NALLOC(l); PL(z) = l;
for ( j = 0, b = (unsigned int *)BD(z); j < l; j++ )
b[l-j-1] = ((unsigned int *)g)[j+2];
}
#endif
#if defined(LONG_IS_64BIT)
void ritopa_i(a,s,rp)
N a;
int s;
GEN *rp;
{
int j,l,words;
unsigned int *b;
GEN z;
l = PL(a); b = BD(a);
words = (l+1)/2;
z = cgeti(words+2);
bzero((char *)&z[2],words*sizeof(long));
for ( j = 0; j < words; j++ )
z[words-j+1] = ((unsigned long)b[2*j])
|(((unsigned long)(2*j+1<l?b[2*j+1]:0))<<32);
s = s>0?1:-1;
setsigne(z,s);
setlgefint(z,lg(z));
*rp = z;
}
void patori_i(g,rp)
GEN g;
N *rp;
{
int j,l,words;
unsigned long t;
unsigned long *a;
unsigned int *b;
N z;
words = lgef(g)-2;
l = 2*words;
a = (unsigned long *)g;
*rp = z = NALLOC(l); PL(z) = l;
for ( j = 0, b = BD(z); j < words; j++ ) {
t = a[words+1-j];
b[2*j] = t&0xffffffff;
b[2*j+1] = t>>32;
}
PL(z) = b[l-1] ? l : l-1;
}
#endif
void strtobf(s,rp)
char *s;
BF *rp;
{
GEN z;
z = lisexpr(s); patori(z,(Obj *)rp); cgiv(z);
}
#endif