Annotation of OpenXM_contrib2/asir2000/lib/ratint, Revision 1.4
1.3 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.4 ! noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.4 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/lib/ratint,v 1.3 2000/08/21 08:31:43 noro Exp $
1.3 noro 49: */
1.1 noro 50: /*
51: * rational function integration (Trager's algorithm)
52: * load("gr"); load("sp");
53: * toplevel : ratint(F,V)
54: * returns a complicated list s.t.
55: * [<rational part>, [[root*log(poly),defpoly],...].
56: */
57:
58: #define FIRST(a) car(a)
59: #define SECOND(a) car(cdr(a))
60: #define THIRD(a) car(cdr(cdr(a)))
61:
62:
63: def substv(P,Sl)
64: {
65: for ( A = P; Sl != []; Sl = cdr(Sl) )
66: A = subst(A,FIRST(car(Sl)),SECOND(car(Sl)));
67: return A;
68: }
69:
70: def co(X,V,D)
71: {
72: for ( I = 0; I < D; I++ )
73: X = diff(X,V);
74: return sdiv(subst(X,V,0),fac(D));
75: }
76:
77: def solve(El,Vl)
78: /*
79: * El : list of linear forms
80: * Vl : list of variable
81: */
82: {
83: N = length(El); M = length(Vl);
84: Mat = newmat(N,M+1);
85: W = newvect(M+1); Index = newvect(N); Vs = newvect(M);
86: for ( I = 0, Tl = Vl; I < M; Tl = cdr(Tl), I++ )
87: Vs[I] = car(Tl);
88: for ( I = 0, Tl = El; I < N; Tl = cdr(Tl), I++ ) {
89: ltov(car(Tl),Vl,W);
90: for ( J = 0; J <= M; J++ )
91: Mat[I][J] = W[J];
92: }
93: Tl = solvemain(Mat,Index,N,M); L = car(Tl); D = car(cdr(Tl));
94: if ( L < 0 )
95: return [];
96: for ( I = L - 1, S = []; I >= 0; I-- ) {
97: for ( J = Index[I]+1, A = 0; J < M; J++ ) {
98: A += Mat[I][J]*Vs[J];
99: }
100: S = cons([Vs[Index[I]],-red((A+Mat[I][M])/D)],S);
101: }
102: return S;
103: }
104:
105: def solvemain(Mat,Index,N,M)
106: /*
107: * Mat : matrix of size Nx(M+1)
108: * Index : vector of length N
109: */
110: {
111: for ( J = 0, L = 0, D = 1; J < M; J++ ) {
112: for ( I = L; I < N && !Mat[I][J]; I++ );
113: if ( I == N )
114: continue;
115: Index[L] = J;
116: for ( K = 0; K <= M; K++ ) {
117: T = Mat[I][K]; Mat[I][K] = Mat[L][K]; Mat[L][K] = T;
118: }
119: for ( I = L + 1, V = Mat[L][J]; I < N; I++ )
120: for ( K = J, U = Mat[I][J]; K <= M; K++ )
121: Mat[I][K] = sdiv(Mat[I][K]*V-Mat[L][K]*U,D);
122: D = V; L++;
123: }
124: for ( I = L; I < N; I++ )
125: for ( J = 0; J <= M; J++ )
126: if ( Mat[I][J] )
127: return -1;
128: for ( I = L - 2, W = newvect(M+1); I >= 0; I-- ) {
129: for ( J = 0; J <= M; J++ )
130: W[J] = 0;
131: for ( G = I + 1; G < L; G++ )
132: for ( H = Index[G], U = Mat[I][H]; H <= M; H++ )
133: W[H] += Mat[G][H]*U;
134: for ( J = Index[I], U = Mat[I][J]; J <= M; J++ )
135: Mat[I][J] = sdiv(Mat[I][J]*D-W[J],U);
136: }
137: return [L,D];
138: }
139:
140: def ltov(P,VL,W)
141: {
142: for ( I = 0, L = VL; L != []; L = cdr(L), I++ ) {
143: W[I] = co(P,car(L),1); P -= W[I]*car(L);
144: }
145: W[I] = P;
146: }
147:
148: def makeucp(N,V) {
149: for ( UCV = [], I = 0; I <= N; I++ )
150: UCV = cons(uc(),UCV);
151: for ( L = UCV, I = P = 0; I <= N; I++, L = cdr(L) )
152: P += car(L)*V^I;
153: return [P,UCV];
154: }
155:
156: def ratint(F,V) {
157: L = ratintsep(F,V);
158: Rat = FIRST(L);
159: if ( !SECOND(L) )
160: return L;
161: else {
162: Pf = ratintpf(SECOND(L),V);
163: for ( T = Pf, S = []; T != []; T = cdr(T) )
164: S = cons(ratintlog(car(T),V),S);
165: return [Rat,S];
166: }
167: }
168:
169: def ratintlog(F,V) {
170: Nm = nm(F); Dn = dn(F);
171: C = uc();
172: R = res(V,ptozp(Nm-C*diff(Dn,V)),ptozp(Dn));
173: Rc = FIRST(SECOND(fctr(R)));
174: if ( deg(Rc,C) == 1 ) {
175: VC = -co(Rc,C,0)/co(Rc,C,1);
176: A = gcd(Nm-VC*diff(Dn,V),Dn);
177: return [VC*log(A),0];
178: } else {
179: Root = newalg(Rc);
180: A = gcda(subst(ptozp(Nm-C*diff(Dn,V)),C,Root),subst(ptozp(Dn),C,Root));
181: return [Root*log(A),defpoly(Root)];
182: }
183: }
184:
185: def ratintsep(F,V) {
186: B = dn(F); A = srem(nm(F),B); P = sdiv(nm(F)-R,B);
187: IP = polyint(P,V);
188: G = gcd(B,diff(B,x));
189: if ( type(G) == 1 )
190: return [IP,red(A/B)];
191: H = sdiv(B,G);
192: N = deg(B,V); M = deg(H,V);
193: CL = makeucp(N-M-1,V); DL = makeucp(M-1,V);
194: C = car(CL); CV = car(cdr(CL));
195: D = car(DL); DV = car(cdr(DL));
196: UCV = append(CV,DV);
197: S = solveuc(A-(diff(C,V)*H-C*sdiv(H*diff(G,V),G)+D*G),V,UCV);
198: C = substv(C,S); D = substv(D,S);
199: return [IP+C/G,red(D/H)];
200: }
201:
202: def polyint(P,V) {
203: if ( !P )
204: return 0;
205: if ( type(P) == 1 )
206: return P*V;
207: for ( I = deg(P,V), T = 0; I >= 0; I-- )
208: T += coef(P,I)/(I+1)*V^(I+1);
209: return T;
210: }
211:
212: def ratintpf(P,V) {
213: NmP = nm(P); DnP = dn(P);
214: DnPf = fctr(DnP);
215: L = length(DnPf) - 1;
216: if ( L == 1 )
217: return [P];
218:
219: Lc = FIRST(car(DnPf)); DnPf = cdr(DnPf);
220: NmP = sdiv(NmP,Lc); DnP = sdiv(DnP,Lc);
221: Nm = newvect(L); Dn = newvect(L);
222: for ( I = 0, F = DnPf; I < L; I++, F = cdr(F) )
223: Dn[I] = FIRST(car(F));
224:
225: for ( I = 0, U = -NmP, Vl = []; I < L; I++ ) {
226: CL = makeucp(deg(Dn[I],V)-1,V);
227: Nm[I] = FIRST(CL); Vl = append(Vl,SECOND(CL));
228: U += sdiv(DnP,Dn[I])*Nm[I];
229: }
230:
231: S = solveuc(U,V,Vl);
232: for ( I = 0, F = []; I < L; I++ )
233: if ( T = substv(Nm[I],S) )
234: F = cons(T/Dn[I],F);
235: return F;
236: }
237:
238: def solveuc(P,V,L) {
239: for ( N = deg(P,V), E = [], I = N; I >= 0; I-- )
240: if ( C = coef(P,I) )
241: E = cons(C,E);
242: EVG = eqsimp(E,L);
243: if ( FIRST(EVG) == [] )
244: return THIRD(EVG);
245: else {
246: S = solve(FIRST(EVG),SECOND(EVG));
247: for ( T = S, G = THIRD(EVG); G != []; G = cdr(G) ) {
248: VV = car(G);
249: T = cons([FIRST(VV),substv(SECOND(VV),S)],T);
250: }
251: return T;
252: }
253:
254: }
255:
256: #if 0
257: def append(A,B) {
258: return A == [] ? B : cons(car(A),append(cdr(A),B));
259: }
260: #endif
261:
262: def eqsimp(E,Vs) {
263: for ( Got = []; ; ) {
264: if ( (VV = searchmonic(E,Vs)) == [] )
265: return [E,Vs,Got];
266: V = FIRST(VV); Val = SECOND(VV);
267: Vs = subtract(Vs,V);
268: for ( T = []; E != []; E = cdr(E) )
269: if ( S = subst(car(E),V,Val) )
270: T = cons(S,T);
271: E = T;
272: for ( T = [VV]; Got != []; Got = cdr(Got) ) {
273: VV1 = car(Got);
274: T = cons([FIRST(VV1),subst(SECOND(VV1),V,Val)],T);
275: }
276: Got = T;
277: }
278: }
279:
280: def searchmonic(E,Vs) {
281: for ( ; E != []; E = cdr(E) )
282: for ( P = car(E), T = Vs; T != []; T = cdr(T) ) {
283: V = car(T); C = diff(P,V);
284: if ( C == 1 )
285: return [V,-(P-V)];
286: else if ( C == -1 )
287: return [V,P+V];
288: }
289: return [];
290: }
291:
292: def subtract(S,E) {
293: for ( T = []; S != []; S = cdr(S) )
294: if ( car(S) == E )
295: return append(T,cdr(S));
296: else
297: T = cons(car(S),T);
298: }
299: 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