← Back to branch summary

~burner/xsb/debianized-xsb

« back to all changes in this revision

Viewing changes to examples/chr/gauss.chr

  • Committer: Michael R. Head
  • Date: 2006-09-06 22:11:55 UTC
  • Revision ID: burner@n23-20060906221155-7e398d23438a7ee4
Add the files from the 3.0.1 release package

Show diffs side-by-side

added added

removed removed

Lines of Context:
examples/chr/gauss.chr
 
1
%% Ported to hProlog by Tom Schrijvers, tom.schrijvers@cs.kuleuven.ac.be
 
2
%%
 
3
%% $Id: gauss.chr,v 1.2 2003/10/03 14:50:37 tschrijvers Exp $
 
4
 
 
5
%% Example usage:
 
6
%%
 
7
%%      ?- [gauss].
 
8
%%      ?- X + Y #= 1, X - Y #= 1, label.
 
9
%%      
 
10
%%      Y = -0.00000000000000000e+00
 
11
%%      X = 1.00000000000000000e+00
 
12
%%      Yes
 
13
%%      ?- X + Y #= 1, X - Y #= 0, label.
 
14
%%      
 
15
%%      Y = X = 5.00000000000000000e-01
 
16
%%      Yes
 
17
%%
 
18
 
 
19
%% solving linear polynomial equations by variable elimination a la Gauss
 
20
%% thom fruehwirth 910610,911213,920124,930602,931223, 980311
 
21
%% 961107 christian holzbaur for SICStus CHR 
 
22
%% complete for equalities, leaves equalities implicit, slow
 
23
%% may loop if variables of the equations are unified
 
24
 
 
25
:- chr_module(gauss).
 
26
 
 
27
:- export (#=)/2, label/0.
 
28
 
 
29
:- import
 
30
        ground/1,
 
31
        length/2
 
32
   from basics.
 
33
 
 
34
:- import
 
35
        call/1
 
36
   from standard.
 
37
 
 
38
:- import
 
39
        globalize/1,
 
40
        var_compare/3
 
41
   from ordering.
 
42
 
 
43
?- op(700,xfx,#=).
 
44
 
 
45
%% math-utilities.pl ===========================================================
 
46
%% auxiliary predicates for math*.pl constraint solvers
 
47
%% thom fruehwirth 1991-92, revised 930518,931223,940304
 
48
%% 961030 christian holzbaur, SICStus adaption
 
49
 
 
50
%% SETTINGS --------------------------------------------------------------------
 
51
 
 
52
:- op(100,xfx,equals).
 
53
 
 
54
%% with float_precision single or double
 
55
%% adapt precision for zero/1 test
 
56
%% :- ( current_module(eclipse) -> 
 
57
%%       get_flag(float_precision,G)
 
58
%%    ; 
 
59
%%       G = double
 
60
%%    ),
 
61
%% ( G==single -> setval(precision,1.0e-06),setval(mprecision,-1.0e-06)
 
62
%% ; 
 
63
%%   G==double -> setval(precision,1.0e-12),setval(mprecision,-1.0e-12)
 
64
%% ).
 
65
getval(precision,1.0e-06).
 
66
getval(mprecision,-1.0e-06).
 
67
 
 
68
 
 
69
%% 
 
70
%% 
 
71
%% AUXILIARY PREDICATES -------------------------------------------------------
 
72
 
 
73
%%
 
74
%% sort X*K,slack(_)*K with globalized Xs
 
75
%%
 
76
sort1(A,B):-
 
77
        msort(A,C),
 
78
        ( (C=[X*_|_],nonvar(X),X=slack(_))->A=B;B=C). % slacks unordered why?
 
79
 
 
80
msort( L, S) :-
 
81
        length( L, Len),
 
82
        msort( Len, L, [], S).
 
83
 
 
84
msort( 0, L,     L, []) :- ! .
 
85
msort( 1, [X|L], L, [X]) :- ! .
 
86
msort( N, L0, L2, S) :-
 
87
        P is N>>1,
 
88
        Q is N-P,
 
89
        msort( P, L0, L1, Sp),
 
90
        msort( Q, L1, L2, Sq),
 
91
        gmerge( Sp, Sq, S).
 
92
 
 
93
gmerge( [], B, B) :- ! .
 
94
gmerge( A, [], A) :- ! .
 
95
gmerge( [A|As], [B|Bs], Res) :-
 
96
        cmp( R, A, B),
 
97
        gmerge( R, A, As, B, Bs, Res).
 
98
 
 
99
gmerge( =, A, As, _, Bs, [A|Rest]) :- gmerge( As, Bs,     Rest).
 
100
gmerge( <, A, As, B, Bs, [A|Rest]) :- gmerge( As, [B|Bs], Rest).
 
101
gmerge( >, A, As, B, Bs, [B|Rest]) :- gmerge( [A|As], Bs, Rest).
 
102
 
 
103
cmp( R, X, Y) :- var(X), var(Y), !, var_compare( R, X, Y) .
 
104
cmp( R, X, _) :- var(X), !, R = (<).
 
105
cmp( R, _, Y) :- var(Y), !, R = (>).
 
106
cmp( R, X, Y) :-
 
107
        functor( X, Fx, Ax),
 
108
        functor( Y, Fy, Ay),
 
109
        compare( Rr, Ax/Fx, Ay/Fy),
 
110
        ( Rr = (=),
 
111
          Ax > 0 ->
 
112
          cmp_args( 1,Ax, X, Y, R)
 
113
        ; 
 
114
          R = Rr
 
115
        ).
 
116
 
 
117
cmp_args( N,M, _, _, R) :- N>M, !, R = (=).
 
118
cmp_args( N,M, X, Y, R) :-
 
119
        arg( N, X, Ax),
 
120
        arg( N, Y, Ay),
 
121
        cmp( Rr, Ax, Ay),
 
122
        ( Rr = (=) ->
 
123
          N1 is N+1,
 
124
          cmp_args( N1,M, X, Y, R)
 
125
        ; 
 
126
          R = Rr
 
127
        ).
 
128
 
 
129
 
 
130
rev([],L,L).
 
131
rev([X|L1],L2,L3):- rev(L1,[X|L2],L3).
 
132
 
 
133
gdelete( X, [X|L],  L).
 
134
gdelete( Y, [X|Xs], [X|Xt]) :-
 
135
        gdelete( Y, Xs, Xt).
 
136
 
 
137
zero(X) :- nonvar(X), X = slack(S), !, zero( S).
 
138
zero(C):-    
 
139
    float(C) -> 
 
140
        getval(precision,P),
 
141
        getval(mprecision,MP),
 
142
        MP < C,                                         % cope with imprecision
 
143
        C < P      
 
144
    ;              
 
145
        C=:=0.
 
146
 
 
147
unwrap(X,Y) :-
 
148
        (  nonvar(X), X = slack(Y) ->
 
149
                true
 
150
        ;
 
151
                X = Y
 
152
        ).
 
153
 
 
154
is_div( C1, C2, C3) :-
 
155
        unwrap( C1, C11),
 
156
        unwrap( C2, C21),
 
157
        unwrap( C3, C31),
 
158
        is_divu( C11, C21, C31).
 
159
 
 
160
is_divu(C1,C2,C3):- zero(C1),!,C3=0.
 
161
is_divu(C1,C2,C3):- X is -(C1/C2), % minus here to get sign needed in handlers
 
162
        avoid_float(X,C3).
 
163
 
 
164
avoid_float(X,C3):-
 
165
        float(X) -> Y is round(X),Z is X-Y,(zero(Z)-> C3 is floor(Y);C3=X) ; C3=X.
 
166
 
 
167
 
 
168
%% HANDLING SLACK VARIABLES ----------------------------------------------------
 
169
 
 
170
%% COMPUTING WITH POLYNOMIALS -------------------------------------------------
 
171
 
 
172
%% gets rounded constant C from is_div/3
 
173
mult_const(eq0(C1,P1),C,eq0(0 ,[])):- C=:=0,! .
 
174
mult_const(eq0(C1,P1),C,eq0(C1,P1)):- C=:=1,! .
 
175
mult_const(eq0(C1,P1),C2,eq0(C,P)):-
 
176
        ( zero(C1) -> C=0 ; C is C1*C2),
 
177
        mult_const1(P1,C2,P).
 
178
mult_const1([],C,[]).
 
179
mult_const1([Xi*Ci|Poly],C,PolyR):-
 
180
        ( zero(Ci) -> PolyR=NPoly ; NCi is Ci*C,PolyR=[Xi*NCi|NPoly]),
 
181
        mult_const1(Poly,C,NPoly).
 
182
 
 
183
%% gets input from const_mult/3
 
184
add_eq0(eq0(C1,P1),eq0(C2,P2),eq0(C,P0)):-
 
185
        Ci is C1+C2,
 
186
        ( zero(Ci) -> C=0 ; C=Ci),
 
187
        add_eq1(P1,P2,P0).
 
188
%%       sort(P,P0).
 
189
add_eq1([],Poly,Poly):- ! .
 
190
add_eq1(Poly,[],Poly):- ! .
 
191
add_eq1([Xi1*Ci1|Poly1],Poly21,Poly):-
 
192
        gdelete(Xi2*Ci2,Poly21,Poly2),Xi2==Xi1,
 
193
        !,
 
194
        Ci is Ci1+Ci2,
 
195
        ( zero(Ci) -> Poly=Poly3 ; Poly=[Xi1*Ci|Poly3]),
 
196
        add_eq1(Poly1,Poly2,Poly3).
 
197
add_eq1([Xi1*Ci1|Poly1],Poly2,[Xi1*Ci1|Poly3]):-
 
198
        add_eq1(Poly1,Poly2,Poly3).
 
199
 
 
200
 
 
201
 
 
202
normalize(A,B,P2,C1):-
 
203
        normalize1(A-B,P),
 
204
        P=eq0(C1,P1),rev(P1,[],P1R), globalize(P1R),
 
205
        sort1(P1,P2).                                 
 
206
 
 
207
normalize1(V,P) :- var(V),!,
 
208
P=eq0(0,[V*1]).
 
209
normalize1(C,P) :- ground(C),!,
 
210
call(C1 is C), P=eq0(C1,[]).
 
211
normalize1(slack(V),P) :- !,
 
212
P=eq0(0,[slack(V)*1]).
 
213
normalize1((+E),P) :- !,
 
214
normalize1(E,P).
 
215
normalize1((-E),P) :- !,
 
216
normalize1(E,P1),
 
217
mult_const(P1,(-1),P).
 
218
normalize1(A*B,C) :- ground(A),!,
 
219
normalize1(B,BN),
 
220
mult_const(BN,A,C).
 
221
normalize1(B*A,C) :- ground(A),!,
 
222
normalize1(B,BN),
 
223
mult_const(BN,A,C).
 
224
normalize1(B/A,C) :- ground(A),!,
 
225
normalize1(B,BN),
 
226
call(A1 is 1/A),
 
227
mult_const(BN,A1,C). 
 
228
normalize1(A-B,C) :- !,
 
229
normalize1(A,AN),
 
230
normalize1((-B),BN),
 
231
add_eq0(AN,BN,C).
 
232
normalize1(A+B,C) :- !,
 
233
normalize1(A,AN),
 
234
normalize1(B,BN),
 
235
add_eq0(AN,BN,C).
 
236
%% normalize1(E,C) :-
 
237
%% C=eq0(0,[CX*1]),
 
238
%% eqnonlin(CX,E).                                      % add a nonlinear equation constraint
 
239
 
 
240
 
 
241
%% end of file math-utilities.pl -----------------------------------------------
 
242
 
 
243
:- constraints (equals)/2.
 
244
%% Poly equals Const, where Poly is list of monomials Variable*Coefficient 
 
245
 
 
246
eliminate(X) @ [X*Coeff1|P1] equals C1 \ P equals C2 
 
247
        <=> gdelete(X*Coeff2,P,P2) | is_div(Coeff2,Coeff1,C), mult_const(eq0(C1,P1),C,eq0(C1C,P1C)), add_eq0(eq0(C2,P2),eq0(C1C,P1C),eq0(C3,P3)), P3 equals C3.
 
248
 
 
249
 
 
250
:- constraints (#=)/2.    
 
251
%% curly brackets as wrapper to avoid name clash with built-in =:=
 
252
 
 
253
%% split @ f( (C, Cs)) <=> f( C ), f( Cs ).
 
254
 
 
255
normalize @ '#='(A,B) <=> 
 
256
        normalize(A,B,Poly,Const), 
 
257
        Poly equals Const.      
 
258
 
 
259
:- constraints label/0.
 
260
 
 
261
label \ equals([V * C], N) # ID <=> var(V), number(C), number(N) | V is -N / C pragma passive(ID).
 
262
label \ equals([],N) # ID <=> N =:= 0 pragma passive(ID).
 
263
label <=> true.
 
264
 
 
265
/* end of file gauss.chr ------------------------------------------------*/
 
266