![[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.2 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/lib/ratint,v 1.1.1.1 1999/12/03 07:39:11 noro Exp $ */
1.1 noro 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 ltov(P,VL,W)
93: {
94: for ( I = 0, L = VL; L != []; L = cdr(L), I++ ) {
95: W[I] = co(P,car(L),1); P -= W[I]*car(L);
96: }
97: W[I] = P;
98: }
99:
100: def makeucp(N,V) {
101: for ( UCV = [], I = 0; I <= N; I++ )
102: UCV = cons(uc(),UCV);
103: for ( L = UCV, I = P = 0; I <= N; I++, L = cdr(L) )
104: P += car(L)*V^I;
105: return [P,UCV];
106: }
107:
108: def ratint(F,V) {
109: L = ratintsep(F,V);
110: Rat = FIRST(L);
111: if ( !SECOND(L) )
112: return L;
113: else {
114: Pf = ratintpf(SECOND(L),V);
115: for ( T = Pf, S = []; T != []; T = cdr(T) )
116: S = cons(ratintlog(car(T),V),S);
117: return [Rat,S];
118: }
119: }
120:
121: def ratintlog(F,V) {
122: Nm = nm(F); Dn = dn(F);
123: C = uc();
124: R = res(V,ptozp(Nm-C*diff(Dn,V)),ptozp(Dn));
125: Rc = FIRST(SECOND(fctr(R)));
126: if ( deg(Rc,C) == 1 ) {
127: VC = -co(Rc,C,0)/co(Rc,C,1);
128: A = gcd(Nm-VC*diff(Dn,V),Dn);
129: return [VC*log(A),0];
130: } else {
131: Root = newalg(Rc);
132: A = gcda(subst(ptozp(Nm-C*diff(Dn,V)),C,Root),subst(ptozp(Dn),C,Root));
133: return [Root*log(A),defpoly(Root)];
134: }
135: }
136:
137: def ratintsep(F,V) {
138: B = dn(F); A = srem(nm(F),B); P = sdiv(nm(F)-R,B);
139: IP = polyint(P,V);
140: G = gcd(B,diff(B,x));
141: if ( type(G) == 1 )
142: return [IP,red(A/B)];
143: H = sdiv(B,G);
144: N = deg(B,V); M = deg(H,V);
145: CL = makeucp(N-M-1,V); DL = makeucp(M-1,V);
146: C = car(CL); CV = car(cdr(CL));
147: D = car(DL); DV = car(cdr(DL));
148: UCV = append(CV,DV);
149: S = solveuc(A-(diff(C,V)*H-C*sdiv(H*diff(G,V),G)+D*G),V,UCV);
150: C = substv(C,S); D = substv(D,S);
151: return [IP+C/G,red(D/H)];
152: }
153:
154: def polyint(P,V) {
155: if ( !P )
156: return 0;
157: if ( type(P) == 1 )
158: return P*V;
159: for ( I = deg(P,V), T = 0; I >= 0; I-- )
160: T += coef(P,I)/(I+1)*V^(I+1);
161: return T;
162: }
163:
164: def ratintpf(P,V) {
165: NmP = nm(P); DnP = dn(P);
166: DnPf = fctr(DnP);
167: L = length(DnPf) - 1;
168: if ( L == 1 )
169: return [P];
170:
171: Lc = FIRST(car(DnPf)); DnPf = cdr(DnPf);
172: NmP = sdiv(NmP,Lc); DnP = sdiv(DnP,Lc);
173: Nm = newvect(L); Dn = newvect(L);
174: for ( I = 0, F = DnPf; I < L; I++, F = cdr(F) )
175: Dn[I] = FIRST(car(F));
176:
177: for ( I = 0, U = -NmP, Vl = []; I < L; I++ ) {
178: CL = makeucp(deg(Dn[I],V)-1,V);
179: Nm[I] = FIRST(CL); Vl = append(Vl,SECOND(CL));
180: U += sdiv(DnP,Dn[I])*Nm[I];
181: }
182:
183: S = solveuc(U,V,Vl);
184: for ( I = 0, F = []; I < L; I++ )
185: if ( T = substv(Nm[I],S) )
186: F = cons(T/Dn[I],F);
187: return F;
188: }
189:
190: def solveuc(P,V,L) {
191: for ( N = deg(P,V), E = [], I = N; I >= 0; I-- )
192: if ( C = coef(P,I) )
193: E = cons(C,E);
194: EVG = eqsimp(E,L);
195: if ( FIRST(EVG) == [] )
196: return THIRD(EVG);
197: else {
198: S = solve(FIRST(EVG),SECOND(EVG));
199: for ( T = S, G = THIRD(EVG); G != []; G = cdr(G) ) {
200: VV = car(G);
201: T = cons([FIRST(VV),substv(SECOND(VV),S)],T);
202: }
203: return T;
204: }
205:
206: }
207:
208: #if 0
209: def append(A,B) {
210: return A == [] ? B : cons(car(A),append(cdr(A),B));
211: }
212: #endif
213:
214: def eqsimp(E,Vs) {
215: for ( Got = []; ; ) {
216: if ( (VV = searchmonic(E,Vs)) == [] )
217: return [E,Vs,Got];
218: V = FIRST(VV); Val = SECOND(VV);
219: Vs = subtract(Vs,V);
220: for ( T = []; E != []; E = cdr(E) )
221: if ( S = subst(car(E),V,Val) )
222: T = cons(S,T);
223: E = T;
224: for ( T = [VV]; Got != []; Got = cdr(Got) ) {
225: VV1 = car(Got);
226: T = cons([FIRST(VV1),subst(SECOND(VV1),V,Val)],T);
227: }
228: Got = T;
229: }
230: }
231:
232: def searchmonic(E,Vs) {
233: for ( ; E != []; E = cdr(E) )
234: for ( P = car(E), T = Vs; T != []; T = cdr(T) ) {
235: V = car(T); C = diff(P,V);
236: if ( C == 1 )
237: return [V,-(P-V)];
238: else if ( C == -1 )
239: return [V,P+V];
240: }
241: return [];
242: }
243:
244: def subtract(S,E) {
245: for ( T = []; S != []; S = cdr(S) )
246: if ( car(S) == E )
247: return append(T,cdr(S));
248: else
249: T = cons(car(S),T);
250: }
251: end$