Annotation of OpenXM_contrib2/asir2018/lib/solve, 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: def kenzan(El,Sl)
! 51: {
! 52: for ( Tl = El; Tl != []; Tl = cdr(Tl) ) {
! 53: if ( substv(car(Tl),Sl) ) {
! 54: print("kenzan : error");
! 55: return 0;
! 56: }
! 57: }
! 58: print("kenzan : ok");
! 59: return 1;
! 60: }
! 61:
! 62: def substv(P,Sl)
! 63: {
! 64: for ( A = P; Sl != []; Sl = cdr(Sl) )
! 65: A = subst(A,car(car(Sl)),car(cdr(car(Sl))));
! 66: return A;
! 67: }
! 68:
! 69: def co(X,V,D)
! 70: {
! 71: for ( I = 0; I < D; I++ )
! 72: X = diff(X,V);
! 73: return sdiv(subst(X,V,0),fac(D));
! 74: }
! 75:
! 76: def solve(El,Vl)
! 77: /*
! 78: * El : list of linear forms
! 79: * Vl : list of variable
! 80: */
! 81: {
! 82: N = length(El); M = length(Vl);
! 83: Mat = newmat(N,M+1);
! 84: W = newvect(M+1); Index = newvect(N); Vs = newvect(M);
! 85: for ( I = 0, Tl = Vl; I < M; Tl = cdr(Tl), I++ )
! 86: Vs[I] = car(Tl);
! 87: for ( I = 0, Tl = El; I < N; Tl = cdr(Tl), I++ ) {
! 88: solve_ltov(car(Tl),Vl,W);
! 89: for ( J = 0; J <= M; J++ )
! 90: Mat[I][J] = W[J];
! 91: }
! 92: Tl = solvemain(Mat,Index,N,M); L = car(Tl); D = car(cdr(Tl));
! 93: if ( L < 0 )
! 94: return [];
! 95: for ( I = L - 1, S = []; I >= 0; I-- ) {
! 96: for ( J = Index[I]+1, A = 0; J < M; J++ ) {
! 97: A += Mat[I][J]*Vs[J];
! 98: }
! 99: S = cons([Vs[Index[I]],-red((A+Mat[I][M])/D)],S);
! 100: }
! 101: if ( kenzan(El,S) )
! 102: return S;
! 103: else
! 104: return [];
! 105: return S;
! 106: }
! 107:
! 108: def solvemain(Mat,Index,N,M)
! 109: /*
! 110: * Mat : matrix of size Nx(M+1)
! 111: * Index : vector of length N
! 112: */
! 113: {
! 114: for ( J = 0, L = 0, D = 1; J < M; J++ ) {
! 115: for ( I = L; I < N && !Mat[I][J]; I++ );
! 116: if ( I == N )
! 117: continue;
! 118: Index[L] = J;
! 119: for ( K = 0; K <= M; K++ ) {
! 120: T = Mat[I][K]; Mat[I][K] = Mat[L][K]; Mat[L][K] = T;
! 121: }
! 122: for ( I = L + 1, V = Mat[L][J]; I < N; I++ )
! 123: for ( K = J, U = Mat[I][J]; K <= M; K++ )
! 124: Mat[I][K] = sdiv(Mat[I][K]*V-Mat[L][K]*U,D);
! 125: D = V; L++;
! 126: }
! 127: for ( I = L; I < N; I++ )
! 128: for ( J = 0; J <= M; J++ )
! 129: if ( Mat[I][J] )
! 130: return -1;
! 131: for ( I = L - 2, W = newvect(M+1); I >= 0; I-- ) {
! 132: for ( J = 0; J <= M; J++ )
! 133: W[J] = 0;
! 134: for ( G = I + 1; G < L; G++ )
! 135: for ( H = Index[G], U = Mat[I][H]; H <= M; H++ )
! 136: W[H] += Mat[G][H]*U;
! 137: for ( J = Index[I], U = Mat[I][J]; J <= M; J++ )
! 138: Mat[I][J] = sdiv(Mat[I][J]*D-W[J],U);
! 139: }
! 140: return [L,D];
! 141: }
! 142:
! 143: def solve_ltov(P,VL,W)
! 144: {
! 145: for ( I = 0, L = VL; L != []; L = cdr(L), I++ ) {
! 146: W[I] = co(P,car(L),1); P -= W[I]*car(L);
! 147: }
! 148: W[I] = P;
! 149: }
! 150: end$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>