Annotation of OpenXM_contrib2/asir2000/lib/ratint, Revision 1.1
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/lib/ratint,v 1.1.1.1 1999/11/10 08:12:30 noro Exp $ */
! 2: /*
! 3: * rational function integration (Trager's algorithm)
! 4: * load("gr"); load("sp");
! 5: * toplevel : ratint(F,V)
! 6: * returns a complicated list s.t.
! 7: * [<rational part>, [[root*log(poly),defpoly],...].
! 8: */
! 9:
! 10: #define FIRST(a) car(a)
! 11: #define SECOND(a) car(cdr(a))
! 12: #define THIRD(a) car(cdr(cdr(a)))
! 13:
! 14:
! 15: def substv(P,Sl)
! 16: {
! 17: for ( A = P; Sl != []; Sl = cdr(Sl) )
! 18: A = subst(A,FIRST(car(Sl)),SECOND(car(Sl)));
! 19: return A;
! 20: }
! 21:
! 22: def co(X,V,D)
! 23: {
! 24: for ( I = 0; I < D; I++ )
! 25: X = diff(X,V);
! 26: return sdiv(subst(X,V,0),fac(D));
! 27: }
! 28:
! 29: def solve(El,Vl)
! 30: /*
! 31: * El : list of linear forms
! 32: * Vl : list of variable
! 33: */
! 34: {
! 35: N = length(El); M = length(Vl);
! 36: Mat = newmat(N,M+1);
! 37: W = newvect(M+1); Index = newvect(N); Vs = newvect(M);
! 38: for ( I = 0, Tl = Vl; I < M; Tl = cdr(Tl), I++ )
! 39: Vs[I] = car(Tl);
! 40: for ( I = 0, Tl = El; I < N; Tl = cdr(Tl), I++ ) {
! 41: ltov(car(Tl),Vl,W);
! 42: for ( J = 0; J <= M; J++ )
! 43: Mat[I][J] = W[J];
! 44: }
! 45: Tl = solvemain(Mat,Index,N,M); L = car(Tl); D = car(cdr(Tl));
! 46: if ( L < 0 )
! 47: return [];
! 48: for ( I = L - 1, S = []; I >= 0; I-- ) {
! 49: for ( J = Index[I]+1, A = 0; J < M; J++ ) {
! 50: A += Mat[I][J]*Vs[J];
! 51: }
! 52: S = cons([Vs[Index[I]],-red((A+Mat[I][M])/D)],S);
! 53: }
! 54: return S;
! 55: }
! 56:
! 57: def solvemain(Mat,Index,N,M)
! 58: /*
! 59: * Mat : matrix of size Nx(M+1)
! 60: * Index : vector of length N
! 61: */
! 62: {
! 63: for ( J = 0, L = 0, D = 1; J < M; J++ ) {
! 64: for ( I = L; I < N && !Mat[I][J]; I++ );
! 65: if ( I == N )
! 66: continue;
! 67: Index[L] = J;
! 68: for ( K = 0; K <= M; K++ ) {
! 69: T = Mat[I][K]; Mat[I][K] = Mat[L][K]; Mat[L][K] = T;
! 70: }
! 71: for ( I = L + 1, V = Mat[L][J]; I < N; I++ )
! 72: for ( K = J, U = Mat[I][J]; K <= M; K++ )
! 73: Mat[I][K] = sdiv(Mat[I][K]*V-Mat[L][K]*U,D);
! 74: D = V; L++;
! 75: }
! 76: for ( I = L; I < N; I++ )
! 77: for ( J = 0; J <= M; J++ )
! 78: if ( Mat[I][J] )
! 79: return -1;
! 80: for ( I = L - 2, W = newvect(M+1); I >= 0; I-- ) {
! 81: for ( J = 0; J <= M; J++ )
! 82: W[J] = 0;
! 83: for ( G = I + 1; G < L; G++ )
! 84: for ( H = Index[G], U = Mat[I][H]; H <= M; H++ )
! 85: W[H] += Mat[G][H]*U;
! 86: for ( J = Index[I], U = Mat[I][J]; J <= M; J++ )
! 87: Mat[I][J] = sdiv(Mat[I][J]*D-W[J],U);
! 88: }
! 89: return [L,D];
! 90: }
! 91:
! 92: def length(L)
! 93: {
! 94: for ( I = 0; L != []; L = cdr(L), I++ );
! 95: return I;
! 96: }
! 97:
! 98: def ltov(P,VL,W)
! 99: {
! 100: for ( I = 0, L = VL; L != []; L = cdr(L), I++ ) {
! 101: W[I] = co(P,car(L),1); P -= W[I]*car(L);
! 102: }
! 103: W[I] = P;
! 104: }
! 105:
! 106: def makeucp(N,V) {
! 107: for ( UCV = [], I = 0; I <= N; I++ )
! 108: UCV = cons(uc(),UCV);
! 109: for ( L = UCV, I = P = 0; I <= N; I++, L = cdr(L) )
! 110: P += car(L)*V^I;
! 111: return [P,UCV];
! 112: }
! 113:
! 114: def ratint(F,V) {
! 115: L = ratintsep(F,V);
! 116: Rat = FIRST(L);
! 117: if ( !SECOND(L) )
! 118: return L;
! 119: else {
! 120: Pf = ratintpf(SECOND(L),V);
! 121: for ( T = Pf, S = []; T != []; T = cdr(T) )
! 122: S = cons(ratintlog(car(T),V),S);
! 123: return [Rat,S];
! 124: }
! 125: }
! 126:
! 127: def ratintlog(F,V) {
! 128: Nm = nm(F); Dn = dn(F);
! 129: C = uc();
! 130: R = res(V,ptozp(Nm-C*diff(Dn,V)),ptozp(Dn));
! 131: Rc = FIRST(SECOND(fctr(R)));
! 132: if ( deg(Rc,C) == 1 ) {
! 133: VC = -co(Rc,C,0)/co(Rc,C,1);
! 134: A = gcd(Nm-VC*diff(Dn,V),Dn);
! 135: return [VC*log(A),0];
! 136: } else {
! 137: Root = newalg(Rc);
! 138: A = gcda(subst(ptozp(Nm-C*diff(Dn,V)),C,Root),subst(ptozp(Dn),C,Root));
! 139: return [Root*log(A),defpoly(Root)];
! 140: }
! 141: }
! 142:
! 143: def ratintsep(F,V) {
! 144: B = dn(F); A = srem(nm(F),B); P = sdiv(nm(F)-R,B);
! 145: IP = polyint(P,V);
! 146: G = gcd(B,diff(B,x));
! 147: if ( type(G) == 1 )
! 148: return [IP,red(A/B)];
! 149: H = sdiv(B,G);
! 150: N = deg(B,V); M = deg(H,V);
! 151: CL = makeucp(N-M-1,V); DL = makeucp(M-1,V);
! 152: C = car(CL); CV = car(cdr(CL));
! 153: D = car(DL); DV = car(cdr(DL));
! 154: UCV = append(CV,DV);
! 155: S = solveuc(A-(diff(C,V)*H-C*sdiv(H*diff(G,V),G)+D*G),V,UCV);
! 156: C = substv(C,S); D = substv(D,S);
! 157: return [IP+C/G,red(D/H)];
! 158: }
! 159:
! 160: def polyint(P,V) {
! 161: if ( !P )
! 162: return 0;
! 163: if ( type(P) == 1 )
! 164: return P*V;
! 165: for ( I = deg(P,V), T = 0; I >= 0; I-- )
! 166: T += coef(P,I)/(I+1)*V^(I+1);
! 167: return T;
! 168: }
! 169:
! 170: def ratintpf(P,V) {
! 171: NmP = nm(P); DnP = dn(P);
! 172: DnPf = fctr(DnP);
! 173: L = length(DnPf) - 1;
! 174: if ( L == 1 )
! 175: return [P];
! 176:
! 177: Lc = FIRST(car(DnPf)); DnPf = cdr(DnPf);
! 178: NmP = sdiv(NmP,Lc); DnP = sdiv(DnP,Lc);
! 179: Nm = newvect(L); Dn = newvect(L);
! 180: for ( I = 0, F = DnPf; I < L; I++, F = cdr(F) )
! 181: Dn[I] = FIRST(car(F));
! 182:
! 183: for ( I = 0, U = -NmP, Vl = []; I < L; I++ ) {
! 184: CL = makeucp(deg(Dn[I],V)-1,V);
! 185: Nm[I] = FIRST(CL); Vl = append(Vl,SECOND(CL));
! 186: U += sdiv(DnP,Dn[I])*Nm[I];
! 187: }
! 188:
! 189: S = solveuc(U,V,Vl);
! 190: for ( I = 0, F = []; I < L; I++ )
! 191: if ( T = substv(Nm[I],S) )
! 192: F = cons(T/Dn[I],F);
! 193: return F;
! 194: }
! 195:
! 196: def solveuc(P,V,L) {
! 197: for ( N = deg(P,V), E = [], I = N; I >= 0; I-- )
! 198: if ( C = coef(P,I) )
! 199: E = cons(C,E);
! 200: EVG = eqsimp(E,L);
! 201: if ( FIRST(EVG) == [] )
! 202: return THIRD(EVG);
! 203: else {
! 204: S = solve(FIRST(EVG),SECOND(EVG));
! 205: for ( T = S, G = THIRD(EVG); G != []; G = cdr(G) ) {
! 206: VV = car(G);
! 207: T = cons([FIRST(VV),substv(SECOND(VV),S)],T);
! 208: }
! 209: return T;
! 210: }
! 211:
! 212: }
! 213:
! 214: #if 0
! 215: def append(A,B) {
! 216: return A == [] ? B : cons(car(A),append(cdr(A),B));
! 217: }
! 218: #endif
! 219:
! 220: def eqsimp(E,Vs) {
! 221: for ( Got = []; ; ) {
! 222: if ( (VV = searchmonic(E,Vs)) == [] )
! 223: return [E,Vs,Got];
! 224: V = FIRST(VV); Val = SECOND(VV);
! 225: Vs = subtract(Vs,V);
! 226: for ( T = []; E != []; E = cdr(E) )
! 227: if ( S = subst(car(E),V,Val) )
! 228: T = cons(S,T);
! 229: E = T;
! 230: for ( T = [VV]; Got != []; Got = cdr(Got) ) {
! 231: VV1 = car(Got);
! 232: T = cons([FIRST(VV1),subst(SECOND(VV1),V,Val)],T);
! 233: }
! 234: Got = T;
! 235: }
! 236: }
! 237:
! 238: def searchmonic(E,Vs) {
! 239: for ( ; E != []; E = cdr(E) )
! 240: for ( P = car(E), T = Vs; T != []; T = cdr(T) ) {
! 241: V = car(T); C = diff(P,V);
! 242: if ( C == 1 )
! 243: return [V,-(P-V)];
! 244: else if ( C == -1 )
! 245: return [V,P+V];
! 246: }
! 247: return [];
! 248: }
! 249:
! 250: def subtract(S,E) {
! 251: for ( T = []; S != []; S = cdr(S) )
! 252: if ( car(S) == E )
! 253: return append(T,cdr(S));
! 254: else
! 255: T = cons(car(S),T);
! 256: }
! 257: end$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ratint CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / lib