1% math-utilities.pl =========================================================== 2% thom fruehwirth 1991-92, revised 930518,931223 3% 940304 matching mistake in normalize1/2 corrected 4% 940304 math_portray/2 extended with clauses for math-fourier.chr 5% for math* solvers 6 7 8% SETTINGS --------------------------------------------------------------------- 9 10% for use in is/2: precision, slack variables, simulated infimum, etc. 11 12% Code works with flag prefer_rationals on or off 13% and with float_precision single or double 14 15% adapt precision for zero/1 test 16:- get_flag(float_precision,G), 17 (G==single -> setval(precision,1e-06),setval(mprecision,-1e-06) 18 ; 19 G==double -> setval(precision,1e-12),setval(mprecision,-1e-12) 20 ). 21 22slack(X,X). % :- X>=0. 23 24inf(3.40282e38). 25minf((-3.40282e38)). 26sup(1e-45). 27msup((-1e-45)). 28 29 30 31% PRETTY PRINT --------------------------------------------------------------- 32 33math_portray(equals(P,C),P1=:=0):- zero(C),!, 34 make_poly(P,P1). 35math_portray(equals(P,C),P1=:=C1):-!, 36 MC is (-C), 37 avoid_float(MC,C1), 38 make_poly(P,P1). 39math_portray(eqnonlin(X,(E)),X=:=E):-!. 40% for math-fourier.chr: 41math_portray(eq(P,C,R),Out):- zero(C),!, 42 Out=..[R,P1,0], 43 make_poly(P,P1). 44math_portray(eq(P,C,R),Out):-!, 45 Out=..[R,P1,C1], 46 MC is (-C), 47 avoid_float(MC,C1), 48 make_poly(P,P1). 49 50:- define_macro((equals)/2,math_portray/2,[write]). 51:- define_macro((eqnonlin)/2,math_portray/2,[write]). 52:- define_macro((eq)/3,math_portray/2,[write]). 53 54%make_poly([],0). 55make_poly([X*C],-CX):- call_kernel(C<0),!, 56 C1 is (-C), 57 avoid_float(C1,C2), 58 make_mono(C2,X,CX). 59make_poly([X*C],CX):-!, 60 avoid_float(C,C1), 61 make_mono(C1,X,CX). 62make_poly([X*C|P],P1-CX):- call_kernel(C<0),!, 63 C1 is (-C), 64 avoid_float(C1,C2), 65 make_mono(C2,X,CX), 66 make_poly(P,P1). 67make_poly([X*C|P],P1+CX):- 68 avoid_float(C,C1), 69 make_mono(C1,X,CX), 70 make_poly(P,P1). 71 72make_mono(C,X,CX):- nonvar(X),X=slack(Y),!,make_mono(C,Y,CX). 73make_mono(C,X,CX1):- nonvar(X),number(X),!,CX is C*X,avoid_float(CX,CX1). 74make_mono(1,X,X):-!. 75make_mono(1_1,X,X):-!. 76make_mono(C,X,C*X). 77 78 79% AUXILIARY PREDICATES ------------------------------------------------------- 80 81call_kernel((A,B)) ?- !,call_kernel(A),call_kernel(B). 82call_kernel(Built_in):- call_explicit(Built_in,sepia_kernel). 83 84sort1(A,B):- 85 sort(A,C), 86 ((C=[X*_|_],nonvar(X),X=slack(_))->A=B;B=C). % slacks unordered why? 87 88% globalize variables into meta variables so they can be sorted 89globalize(Vars):- delay(Vars-Flag,true),Flag=fired. 90 91rev([],L,L). 92rev([X|L1],L2,L3):- rev(L1,[X|L2],L3). 93 94extract(X*C2,P0,P) ?- delete(Y*C2,P0,P),X==Y,!. 95 96zero(C):- 97 real(C) -> 98 getval(precision,P), % otherwise call-kernel does not work 99 getval(mprecision,MP), 100 call_kernel(MP < C), % cope with imprecision 101 call_kernel(C < P) 102 ; 103 call_kernel(C=:=0). 104 105nonzero(C):- not zero(C). 106 107 108is_div(C1,C2,C3):- nonvar(C3),C3=slack(C4),!,is_div(C1,C2,C4). % slack-case 109is_div(C1,C2,C3):- zero(C1),!,C3=0. 110is_div(C1,C2,C3):- X is -(C1/C2), % minus here to get sign needed in handlers 111 avoid_float(X,C3). 112 113is_mul(C1,C2,C3):- zero(C1),!,C3=0. 114is_mul(C1,C2,C3):- zero(C2),!,C3=0. 115is_mul(C1,C2,C3):- X is C1*C2, 116 avoid_float(X,C3). 117 118avoid_float(X,C3):- 119 real(X) -> Y is round(X),Z is X-Y,(zero(Z)-> C3 is fix(Y);C3=X) ; C3=X. 120 121 122simplifyable(X*C,P,P1):- delete(X*C,P,P1),ground(X),!. 123 124ground(X):- not nonground(X). 125 126 127% HANDLING SLACK VARIABLES ---------------------------------------------------- 128 129all_slacks([]). 130all_slacks([slack(_)*C|P]) ?- 131 all_slacks(P). 132 133all_slacks([],_). 134all_slacks([slack(_)*C|P],S) ?- 135 sign(C,S), 136 all_slacks(P,S). 137 138sign(C,0):- zero(C),!. 139sign(C,S):- call_kernel(C>0) -> S=1 ; S=(-1). 140 141all_zeroes([]). 142all_zeroes([slack(0)*C|P]) :- 143 all_zeroes(P). 144 145 146% COMPUTING WITH POLYNOMIALS ------------------------------------------------- 147 148% gets rounded constant C from is_div/3 149mult_const(eq0(C1,P1),C,eq0(0 ,[])):- call_kernel(C=:=0),!. 150mult_const(eq0(C1,P1),C,eq0(C1,P1)):- call_kernel(C=:=1),!. 151mult_const(eq0(C1,P1),C2,eq0(C,P)):- 152 (zero(C1) -> C=0 ; C is C1*C2), 153 mult_const1(P1,C2,P). 154 mult_const1([],C,[]). 155 mult_const1([Xi*Ci|Poly],C,PolyR):- 156 (zero(Ci) -> PolyR=NPoly ; NCi is Ci*C,PolyR=[Xi*NCi|NPoly]), 157 mult_const1(Poly,C,NPoly). 158 159% gets input from const_mult/3 160add_eq0(eq0(C1,P1),eq0(C2,P2),eq0(C,P0)):- 161 Ci is C1+C2, 162 (zero(Ci) -> C=0 ; C=Ci), 163 add_eq1(P1,P2,P0). 164% sort(P,P0). 165 add_eq1([],Poly,Poly):-!. 166 add_eq1(Poly,[],Poly):-!. 167 add_eq1([Xi1*Ci1|Poly1],Poly21,Poly):- 168 delete(Xi2*Ci2,Poly21,Poly2),Xi2==Xi1, 169 !, 170 Ci is Ci1+Ci2, 171 (zero(Ci) -> Poly=Poly3 ; Poly=[Xi1*Ci|Poly3]), 172 add_eq1(Poly1,Poly2,Poly3). 173 add_eq1([Xi1*Ci1|Poly1],Poly2,[Xi1*Ci1|Poly3]):- 174 add_eq1(Poly1,Poly2,Poly3). 175 176 177 178normalize(A,B,P2,C1):- 179 normalize1(A-B,P), 180 P=eq0(C1,P1),rev(P1,[],P1R),globalize(P1R), 181 sort1(P1,P2). 182 183 normalize1(V,P) ?- var(V),!, 184 P=eq0(0,[V*1]). 185 normalize1(C,P) ?- ground(C),!, 186 C1 is C,P=eq0(C1,[]). 187 normalize1(slack(V),P) ?- !, 188 P=eq0(0,[slack(V)*1]). 189 normalize1((+E),P) ?-!, 190 normalize1(E,P). 191 normalize1((-E),P) ?-!, 192 normalize1(E,P1), 193 mult_const(P1,(-1),P). 194 normalize1(A*B,C) ?- ground(A),!, 195 normalize1(B,BN), 196 mult_const(BN,A,C). 197 normalize1(B*A,C) ?- ground(A),!, 198 normalize1(B,BN), 199 mult_const(BN,A,C). 200 normalize1(B/A,C) ?- ground(A),!, 201 normalize1(B,BN), 202 A1 is 1/A, 203 mult_const(BN,A1,C). 204 normalize1(A-B,C) ?- !, 205 normalize1(A,AN), 206 normalize1((-B),BN), 207 add_eq0(AN,BN,C). 208 normalize1(A+B,C) ?- !, 209 normalize1(A,AN), 210 normalize1(B,BN), 211 add_eq0(AN,BN,C). 212 normalize1(E,C) ?- 213 C=eq0(0,[CX*1]), 214 eqnonlin(CX,E). % add a nonlinear equation constraint 215 216 217% end of file math-utilities.pl -----------------------------------------------