File: [local] / OpenXM_contrib2 / asir2000 / engine / Q.c (download)
Revision 1.3, Mon May 29 08:54:46 2000 UTC (24 years ago) by noro
Branch: MAIN
Changes since 1.2: +19 -1
lines
1. F4
array.c, gr.c (still experimental)
2. Weyl algebra
dist.c, distm.c : product of monomials (weyl_mul*)
Q.c : coefficients of the expansion of D^k*x^l (mkwc, mkwcm)
Note that the current implementation and specification are too ad hoc.
If ctrl("do_weyl",1) is executed, then all monomial operations are
done in Weyl algebra. If the length of the exponent of a monomial
is n=2m, then it is regarded as an element of Q<x1,...,xm,Dx1,...,Dxm>.
If the length is n=2m+1, then it is a regarded as an element of
Q[h]<x1,...,xm,Dx1,...,Dxm>, where h is the homogenization variable.
The order specification is the same as in the commutative case, so
one should use matrix order to realize natural orderings in Weyl
algebra. Negative waits have not yet been supported.
|
/* $OpenXM: OpenXM_contrib2/asir2000/engine/Q.c,v 1.3 2000/05/29 08:54:46 noro Exp $ */
#include "ca.h"
#include "base.h"
#include "inline.h"
void addq(n1,n2,nr)
Q n1,n2,*nr;
{
N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m;
int sgn;
unsigned t,t1;
if ( !n1 )
*nr = n2;
else if ( !n2 )
*nr = n1;
else if ( INT(n1) )
if ( INT(n2) ) {
nm1 = NM(n1); nm2 = NM(n2);
if ( SGN(n1) == SGN(n2) ) {
if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
t1 = BD(nm1)[0]; t = t1+BD(nm2)[0];
if ( t < t1 ) {
m = NALLOC(2); PL(m) = 2; BD(m)[0] = t; BD(m)[1] = 1;
} else
UTON(t,m);
} else
addn(NM(n1),NM(n2),&m);
NTOQ(m,SGN(n1),*nr);
} else {
if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
t1 = BD(nm1)[0]; t = t1-BD(nm2)[0];
if ( !t )
sgn = 0;
else if ( t > t1 ) {
sgn = -1; t = -((int)t); UTON(t,m);
} else {
sgn = 1; UTON(t,m);
}
} else
sgn = subn(NM(n1),NM(n2),&m);
if ( sgn ) {
if ( SGN(n1) == sgn )
NTOQ(m,1,*nr);
else
NTOQ(m,-1,*nr);
} else
*nr = 0;
}
} else {
kmuln(NM(n1),DN(n2),&m);
if ( SGN(n1) == SGN(n2) ) {
sgn = SGN(n1); addn(m,NM(n2),&nm);
} else
sgn = SGN(n1)*subn(m,NM(n2),&nm);
if ( sgn ) {
dn = DN(n2); NDTOQ(nm,dn,sgn,*nr);
} else
*nr = 0;
}
else if ( INT(n2) ) {
kmuln(NM(n2),DN(n1),&m);
if ( SGN(n1) == SGN(n2) ) {
sgn = SGN(n1); addn(m,NM(n1),&nm);
} else
sgn = SGN(n1)*subn(NM(n1),m,&nm);
if ( sgn ) {
dn = DN(n1); NDTOQ(nm,dn,sgn,*nr);
} else
*nr = 0;
} else {
gcdn(DN(n1),DN(n2),&g); divsn(DN(n1),g,&dn1); divsn(DN(n2),g,&dn2);
kmuln(NM(n1),dn2,&nm1); kmuln(NM(n2),dn1,&nm2);
if ( SGN(n1) == SGN(n2) ) {
sgn = SGN(n1); addn(nm1,nm2,&nm3);
} else
sgn = SGN(n1)*subn(nm1,nm2,&nm3);
if ( sgn ) {
gcdn(nm3,g,&g1); divsn(nm3,g1,&nm); divsn(g,g1,&g0);
kmuln(dn1,dn2,&dn0); kmuln(g0,dn0,&dn);
if ( UNIN(dn) )
NTOQ(nm,sgn,*nr);
else
NDTOQ(nm,dn,sgn,*nr);
} else
*nr = 0;
}
}
void subq(n1,n2,nr)
Q n1,n2,*nr;
{
Q m;
if ( !n1 )
if ( !n2 )
*nr = 0;
else {
DUPQ(n2,*nr); SGN(*nr) = -SGN(n2);
}
else if ( !n2 )
*nr = n1;
else if ( n1 == n2 )
*nr = 0;
else {
DUPQ(n2,m); SGN(m) = -SGN(n2); addq(n1,m,nr);
}
}
void mulq(n1,n2,nr)
Q n1,n2,*nr;
{
N nm,nm1,nm2,dn,dn1,dn2,g;
int sgn;
unsigned u,l;
if ( !n1 || !n2 ) *nr = 0;
else if ( INT(n1) )
if ( INT(n2) ) {
nm1 = NM(n1); nm2 = NM(n2);
if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
DM(BD(NM(n1))[0],BD(NM(n2))[0],u,l)
if ( u ) {
nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = l; BD(nm)[1] = u;
} else {
nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = l;
}
} else
kmuln(nm1,nm2,&nm);
sgn = (SGN(n1)==SGN(n2)?1:-1); NTOQ(nm,sgn,*nr);
} else {
gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn);
kmuln(nm1,NM(n2),&nm); sgn = SGN(n1)*SGN(n2);
if ( UNIN(dn) )
NTOQ(nm,sgn,*nr);
else
NDTOQ(nm,dn,sgn,*nr);
}
else if ( INT(n2) ) {
gcdn(NM(n2),DN(n1),&g); divsn(NM(n2),g,&nm2); divsn(DN(n1),g,&dn);
kmuln(nm2,NM(n1),&nm); sgn = SGN(n1)*SGN(n2);
if ( UNIN(dn) )
NTOQ(nm,sgn,*nr);
else
NDTOQ(nm,dn,sgn,*nr);
} else {
gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn2);
gcdn(DN(n1),NM(n2),&g); divsn(DN(n1),g,&dn1); divsn(NM(n2),g,&nm2);
kmuln(nm1,nm2,&nm); kmuln(dn1,dn2,&dn); sgn = SGN(n1) * SGN(n2);
if ( UNIN(dn) )
NTOQ(nm,sgn,*nr);
else
NDTOQ(nm,dn,sgn,*nr);
}
}
void divq(n1,n2,nq)
Q n1,n2,*nq;
{
Q m;
if ( !n2 ) {
error("division by 0");
*nq = 0;
return;
} else if ( !n1 )
*nq = 0;
else if ( n1 == n2 )
*nq = ONE;
else {
invq(n2,&m); mulq(n1,m,nq);
}
}
void invq(n,nr)
Q n,*nr;
{
N nm,dn;
if ( INT(n) )
if ( UNIN(NM(n)) )
if ( SGN(n) > 0 )
*nr = ONE;
else {
nm = ONEN; NTOQ(nm,SGN(n),*nr);
}
else {
nm = ONEN; dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
}
else if ( UNIN(NM(n)) ) {
nm = DN(n); NTOQ(nm,SGN(n),*nr);
} else {
nm = DN(n); dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
}
}
void chsgnq(n,nr)
Q n,*nr;
{
Q t;
if ( !n )
*nr = 0;
else {
DUPQ(n,t); SGN(t) = -SGN(t); *nr = t;
}
}
void pwrq(n1,n,nr)
Q n1,n,*nr;
{
N nm,dn;
int sgn;
if ( !n || UNIQ(n1) )
*nr = ONE;
else if ( !n1 )
*nr = 0;
else if ( !INT(n) ) {
error("can't calculate fractional power."); *nr = 0;
} else if ( MUNIQ(n1) )
*nr = BD(NM(n))[0]%2 ? n1 : ONE;
else if ( PL(NM(n)) > 1 ) {
error("exponent too big."); *nr = 0;
} else {
sgn = ((SGN(n1)>0)||EVENN(NM(n))?1:-1);
pwrn(NM(n1),BD(NM(n))[0],&nm);
if ( INT(n1) ) {
if ( UNIN(nm) )
if ( sgn == 1 )
*nr = ONE;
else
NTOQ(ONEN,sgn,*nr);
else if ( SGN(n) > 0 )
NTOQ(nm,sgn,*nr);
else
NDTOQ(ONEN,nm,sgn,*nr);
} else {
pwrn(DN(n1),BD(NM(n))[0],&dn);
if ( SGN(n) > 0 )
NDTOQ(nm,dn,sgn,*nr);
else if ( UNIN(nm) )
NTOQ(dn,sgn,*nr);
else
NDTOQ(dn,nm,sgn,*nr);
}
}
}
int cmpq(q1,q2)
Q q1,q2;
{
int sgn;
N t,s;
if ( !q1 )
if ( !q2 )
return ( 0 );
else
return ( -SGN(q2) );
else if ( !q2 )
return ( SGN(q1) );
else if ( SGN(q1) != SGN(q2) )
return ( SGN(q1) );
else if ( INT(q1) )
if ( INT(q2) ) {
sgn = cmpn(NM(q1),NM(q2));
if ( !sgn )
return ( 0 );
else
return ( SGN(q1)==sgn?1:-1 );
} else {
kmuln(NM(q1),DN(q2),&t); sgn = cmpn(t,NM(q2));
return ( SGN(q1) * sgn );
}
else if ( INT(q2) ) {
kmuln(NM(q2),DN(q1),&t); sgn = cmpn(NM(q1),t);
return ( SGN(q1) * sgn );
} else {
kmuln(NM(q1),DN(q2),&t); kmuln(NM(q2),DN(q1),&s); sgn = cmpn(t,s);
return ( SGN(q1) * sgn );
}
}
void remq(n,m,nr)
Q n,m;
Q *nr;
{
N q,r,s;
if ( !n ) {
*nr = 0;
return;
}
divn(NM(n),NM(m),&q,&r);
if ( !r )
*nr = 0;
else if ( SGN(n) > 0 )
NTOQ(r,1,*nr);
else {
subn(NM(m),r,&s); NTOQ(s,1,*nr);
}
}
/* t = [nC0 nC1 ... nCn] */
void mkbc(n,t)
int n;
Q *t;
{
int i;
N c,d;
for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1;
kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
}
for ( ; i <= n; i++ )
t[i] = t[n-i];
}
/*
* Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
*
* t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
* where W(k,l,i) = i! * kCi * lCi
*/
void mkwc(k,l,t)
int k,l;
Q *t;
{
int i,n,up,low;
N nm,d,c;
n = MIN(k,l);
for ( t[0] = ONE, i = 1; i <= n; i++ ) {
DM(k-i+1,l-i+1,up,low);
if ( up ) {
nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;
} else {
nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;
}
kmuln(NM(t[i-1]),nm,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
}
}
/* mod m table */
/* XXX : should be optimized */
void mkwcm(k,l,m,t)
int k,l,m;
int *t;
{
int i,n;
Q *s;
n = MIN(k,l);
s = (Q *)ALLOCA((n+1)*sizeof(Q));
mkwc(k,l,s);
for ( i = 0; i <= n; i++ ) {
t[i] = rem(NM(s[i]),m);
}
}
void factorial(n,r)
int n;
Q *r;
{
Q t,iq,s;
unsigned int i,m,m0;
unsigned int c;
if ( !n )
*r = ONE;
else if ( n < 0 )
*r = 0;
else {
for ( i = 1, t = ONE; (int)i <= n; ) {
for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) {
m0 = m; DM(m0,i,c,m)
}
if ( c ) {
m = m0; i--;
}
UTOQ(m,iq); mulq(t,iq,&s); t = s;
}
*r = t;
}
}
void invl(a,mod,ar)
Q a,mod,*ar;
{
Q f1,f2,a1,a2,q,m,s;
N qn,rn;
a1 = ONE; a2 = 0;
DUPQ(a,f1); SGN(f1) = 1; f2 = mod;
while ( 1 ) {
divn(NM(f1),NM(f2),&qn,&rn);
if ( !qn )
q = 0;
else
NTOQ(qn,1,q);
f1 = f2;
if ( !rn )
break;
else
NTOQ(rn,1,f2);
mulq(a2,q,&m); subq(a1,m,&s); a1 = a2;
if ( !s )
a2 = 0;
else
remq(s,mod,&a2);
}
if ( SGN(a) < 0 )
chsgnq(a2,&m);
else
m = a2;
if ( SGN(m) >= 0 )
*ar = m;
else
addq(m,mod,ar);
}
int kara_mag=100;
void kmuln(n1,n2,nr)
N n1,n2,*nr;
{
N n,t,s,m,carry;
int d,d1,d2,len,i,l;
int *r,*r0;
if ( !n1 || !n2 ) {
*nr = 0; return;
}
d1 = PL(n1); d2 = PL(n2);
if ( (d1 < kara_mag) || (d2 < kara_mag) ) {
muln(n1,n2,nr); return;
}
if ( d1 < d2 ) {
n = n1; n1 = n2; n2 = n;
d = d1; d1 = d2; d2 = d;
}
if ( d2 > (d1+1)/2 ) {
kmulnmain(n1,n2,nr); return;
}
d = (d1/d2)+((d1%d2)!=0);
len = (d+1)*d2;
r0 = (int *)ALLOCA(len*sizeof(int));
bzero((char *)r0,len*sizeof(int));
for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {
extractn(n1,i*d2,d2,&m);
if ( m ) {
kmuln(m,n2,&t);
addn(t,carry,&s);
copyn(s,d2,r);
extractn(s,d2,d2,&carry);
} else {
copyn(carry,d2,r);
carry = 0;
}
}
copyn(carry,d2,r);
for ( l = len - 1; !r0[l]; l-- );
l++;
*nr = t = NALLOC(l); PL(t) = l;
bcopy((char *)r0,(char *)BD(t),l*sizeof(int));
}
void extractn(n,index,len,nr)
N n;
int index,len;
N *nr;
{
unsigned int *m;
int l;
N t;
if ( !n ) {
*nr = 0; return;
}
m = BD(n)+index;
if ( (l = PL(n)-index) >= len ) {
for ( l = len - 1; (l >= 0) && !m[l]; l-- );
l++;
}
if ( l <= 0 ) {
*nr = 0; return;
} else {
*nr = t = NALLOC(l); PL(t) = l;
bcopy((char *)m,(char *)BD(t),l*sizeof(int));
}
}
void copyn(n,len,p)
N n;
int len;
int *p;
{
if ( n )
bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));
}
void dupn(n,p)
N n;
N p;
{
if ( n )
bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));
}
void kmulnmain(n1,n2,nr)
N n1,n2,*nr;
{
int d1,d2,h,sgn,sgn1,sgn2,len;
N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2;
extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi);
extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi);
kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo);
len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len;
bzero((char *)BD(t1),len*sizeof(int));
if ( lo )
bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int));
if ( hi )
bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int));
addn(hi,lo,&mid1);
sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2);
kmuln(s1,s2,&mid2); sgn = sgn1*sgn2;
if ( sgn > 0 )
addn(mid1,mid2,&mid);
else
subn(mid1,mid2,&mid);
if ( mid ) {
len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len;
bzero((char *)BD(t2),len*sizeof(int));
bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int));
addn(t1,t2,nr);
} else
*nr = t1;
}