Annotation of OpenXM_contrib2/asir2018/lib/ratint, Revision 1.1
1.1 ! 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
! 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
! 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: *
! 48: * $OpenXM$
! 49: */
! 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 / asir2018 / lib