![[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|
Return to ratint CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
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$