home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / lifeos2.zip / LIFE-1.02 / EXAMPLES / GAUSS.LF < prev    next >
Text File  |  1996-06-04  |  16KB  |  471 lines
  1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  2. %                               %
  3. %         GAUSS V. 1.1          %   (pre-alpha version, kind of)
  4. %                               %
  5. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  6.  
  7. % Author: Christophe Bonnet
  8. % (c) Copyright 1992 Digital Equipment Corporation
  9. % All Rights Reserved
  10.  
  11.  
  12. %%%%%%%%%%%%%%%%%%%%%%%%%
  13. % How to use the solver ?
  14. %
  15. % First,needless to say, load it. You'll then have to start the solver.
  16. % In the most current case, at the moment when you start the solver, you
  17. % want to give it immediately some equations to solve. Here is the
  18. % way to do it :
  19. % starter([equation1,equation2...],Channel)?
  20. %
  21. %   The "Channel" must be an uninstantiated variable. It will play the
  22. % role of a constraint queue, the mandatory passage for sending new
  23. % equations to the solver.
  24. %   The equations are of any of the form described below, in the
  25. % normalization module. But don't bother to look through all the code,
  26. % here is the authorized formulae :
  27. %      => All formulae contructed with LIFE variables, numbers, and the
  28. %         "+", "*", "-" and "~" symbols. The three formers have their
  29. %         natural syntax (i.e. binary, infix, left-associative and the
  30. %         such) and meaning. Note that they can't be used as unary
  31. %         operators, although it's a very common usage. The "~" sign is
  32. %         a representation of the unary minus.
  33. %      => Two such expressions, each on one side of an equal sign.
  34. %
  35. % Note : The former are supposed to be equal to zero. The latter
  36. % behave in the obvious way.
  37. %
  38. % At this point, the reader/user may be puzzled by the fact that a
  39. % linear equation solver accepts any type of polynomial equations as
  40. % input. To be honest, the program behave in the following way :
  41. %  - Even if your equation doesn't look linear, it will try to expand,
  42. %  reduce and simplify it, in case all high degree terms could be
  43. %  eliminated.
  44. %  - If, after all its try, the normalizer can't turn your work into sum
  45. %  of couples (I mean product) number/variable, it will fail, with a
  46. %  polite note saying that your entry wasn't linear.
  47. %
  48. % I all has worked fine, the solver will treat your equations, and then
  49. % you'll be back to the usual LIFE environement, with strange values
  50. % assigned to your previously free variables. Note that (and, I know,
  51. % that's a pity) these values will be the only output of the solver.
  52. %
  53. %  - The "Channel" variable : You don't have to worry about is value. It
  54. %  should be something like q(@~). The tlda, as any experimented life
  55. %  user will know, means that a function has residuate on this variable.
  56. %  Don't worry : It is our solver itself, quietly waiting for you to
  57. %  give him more equations.
  58. %  - On the opposite, the values of the variables that has occured in
  59. %  your equation is of great interest :  There you'll find the value, in
  60. %  the mathematical sense, of your variable. The psi-term that is
  61. %  associated with each variable looks like the following :
  62. %
  63. %       X : vr(3,val=>[5*Y,-2*Z,765],wlist=>[]) 
  64. %
  65. % The sort "vr" for "variable" is just here on programming purpose, as
  66. % is the identifier (the first field, with value 3 in our exemple).
  67. % The attribute wlist won't give you very interesting hints either,
  68. % unless you worry about the conception of th program. The more
  69. % significant field is, for the user, the "val" field. If the variable
  70. % has a definite numeric value (deduced from the equation system), it
  71. % will be stored here, in a one element list ([42], for exemple). Else,
  72. % it will provide you with a linear function of other variables, which
  73. % is, in our former example, f(Y,Z)= 5Y - 2Z + 765. This function will
  74. % automatically be replaced by a more precise one, or eventually a
  75. % numeric value, when you'll give further information to the solver.
  76.  
  77.  
  78. %%% feeding the best
  79. %
  80. % The "Channel" variable will allow you to give additional equations,
  81. % that will be solved in the context defined by previous equations (i.e,
  82. % the partial or definitive values your variables have now). You could
  83. % do that "by hand", of course, by directly modifying the content of the
  84. % channel variable, but we won't encourage such a practice. We've made
  85. % for you two cute functions which theoretically do this job properly.
  86. % Their syntax is :
  87. %
  88. %         feeder(equation,Channel)?
  89. %         fast_feeder([eq1,eq2,...],Channel)?
  90. % They work exactly the same way, except that feeder wwill handle a
  91. % single equation while fast_feeder take a list of equation. Use
  92. % whichever you want.
  93. % By the way, you're perfectly free to put new variables in the
  94. % equations. The counter "varcount" will assure them to have different Id.
  95. %
  96. % What else to say ? Well, I think I've told the main points, on the
  97. % user point of view I mean.  So... enjoy !
  98.  
  99.  
  100. %            Christophe.
  101. %
  102. % PS : I hope you won't mind my spelling. Don't worry, I won't run for
  103. % Vice-Presidency of the USA :-)
  104. % By the way, It's quit as bad in my own language. But with a syntax much
  105. % less fantaisist...
  106. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  107.  
  108.  
  109. %%%%%%%%%%%%%%%%%
  110. % changes from 1.0 to 1.1 :
  111. %
  112. % - use of genint for numbering the variables.
  113. % - variables now have a "name" field, which is not used now, but will be..
  114.  
  115. module("gauss") ?
  116. public( starter,feeder,fast_feeder,vr,~) ?
  117.  
  118.  
  119. collector_starter(Q:q(@)) :- @=constraint_collector(Q).
  120. collector_starter(Q) :- failure_handler(Q).
  121.  
  122. constraint_collector(Q:q([E|S]))-> Cont | 
  123.         solve_constraint(E),
  124.     Q <- q(S),
  125.     Cont = constraint_collector(Q).
  126.  
  127. failure_handler(@) :- nl, write("=} I'm sorry, but I'm dying !"),nl,fail.
  128.  
  129. solve_constraint(gauss(L)) :- solve_eq(L).
  130.  
  131. non_strict(feeder,fast_feeder,starter)?
  132.  
  133. feeder(Eq,Q) :- 
  134.     Cstr=gauss(normalize_and_test(Eq)),
  135.     Q=q([Cstr|@]).
  136.  
  137. fast_feeder([],@).
  138. fast_feeder([E|Tl],Q) :- feeder(E,Q), fast_feeder(Tl,Q).
  139.  
  140. starter(EqList,Queue) :- 
  141.     Queue :== @,
  142.     collector_starter(Queue),
  143.     fast_feeder(EqList,Queue).
  144.  
  145.  
  146.  
  147. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  148. % NORMALIZING ARBITRARY POLYNOMIAL FORMS
  149.  
  150. %   Before to give the equations to the solver, we want to normalize them,
  151. % that is, in our case, translating them to an ordered sum of
  152. % coefficient*variable products. The somewhat obuscated form of the code below
  153. % is due to the fact that we want to give the user the freedom of writing ugly
  154. % things. By "ugly things", I mean any arithmetic expression that we can write
  155. % using *,- and + (and ~, who is just a convenient symbol for the unary minus
  156. % --- overloading of the minus sign  bring us a pretty heap of trouble...).
  157. % It's clear that such expressions won't always be linear, but our politic is
  158. % to accept any polynomial expression as an input, and throw them away if the
  159. % normalization module haven't found that they really were nonlinear.
  160.  
  161. % examples of correct expressions (due to Peter Van Roy), who appeared to
  162. % be linear after a trip through our normalizer :
  163.  
  164. % (X+1)*(X-1)-X*X
  165.  
  166. % X*Y+45*(Y+W+3*(W+Q+23))-Y*X
  167.  
  168. % (X*X+Y*Y)-(X+Y)*(Y+X)+X*Y*(1+1)
  169.  
  170. % (A+B+C+D+E+F+G)*(1+2+3+4+3+Z)-Z*(G+F+E+D+C+B+A)
  171.  
  172.  
  173. % note that expression without equal sign are treated as if there were an "= 0"
  174. % just behind them. 
  175.  
  176.  
  177. %%%%%%%%%%%%%%%%%%%%%%%%%
  178. %  The normalize function
  179.  
  180. op(500,fx,~)?
  181.  
  182. non_strict(normalize,norm_first_pass,norm_dist,kill_zeroes,non_zero_term)?
  183.  
  184. normalize(E)-> Expr | 
  185.         norm_first_pass(E,F),
  186.     G= norm_dist(F),!,
  187.     kill_zeroes(G,Expr).
  188.  
  189. kill_zeroes([],[0]).
  190. kill_zeroes([A:real],[A]).
  191. kill_zeroes([A*V|Tl1],[A*V|Tl2]) :- 
  192.     non_zero_term(A),!,
  193.     kill_zeroes(Tl1,Tl2).
  194. kill_zeroes([@|Tl1],Tl2) :-  kill_zeroes(Tl1,Tl2).
  195.  
  196.  
  197. non_zero_term(0) :- !,fail.
  198. non_zero_term(real) :- !.
  199. non_zero_term(A*vr) :- non_zero_term(A).
  200.  
  201.  
  202. %  normalize_and_test    check if th expression (I mean, the equation)
  203. %  is of degre 1, i.e that, after simplification, there is no variables
  204. %  product left in it.
  205.  
  206. normalize_and_test(E)-> Expr | 
  207.         norm_first_pass(E,F),
  208.     G= norm_dist(F),
  209.     kill_zeroes(G,Expr),!,
  210.     deg_check(Expr).
  211.  
  212. non_strict(deg_check) ?
  213. deg_check([real]) :- !.
  214. deg_check([real*vr|Tail]) :- !,deg_check(Tail).
  215. deg_check(@) :- write("=} Your equation is definitely non linear. We can't do
  216.  anything for you..."),!,fail.
  217.              
  218. %%%%%%%%%%%%
  219. % This first part of the normalizer will create the adequate data structures
  220. % for any free variable she will meet, and suppress the minus (unary and
  221. % binary) signs. It will also multiply a lot of variables by one, thus ensuring
  222. % that we have no variable in our expression that aren't the left factor of a
  223. % product. 
  224.  
  225. %dynamic(varcount)?
  226. %varcount -> 0.
  227. %nvar -> vr(X,val=>none,wlist=>[])|X=varcount+1,setq(varcount,X).
  228. nvar -> vr(Id:genint,name=>strcon("_",int2str(Id)),val=>none,wlist=>[]).
  229.  
  230. named_nvar(N:string) -> vr(Id:genint,name=>N,val=>none,wlist=>[]).
  231.  
  232. norm_first_pass(X,1*X) :- 
  233.     X :== @,!,
  234.     X=nvar.
  235. norm_first_pass(V:vr,1*V) :- !.
  236. norm_first_pass(A:real,A) :- !.
  237. norm_first_pass( ~X , Z*X2) :-
  238.     !,Z= -1,
  239.     norm_first_pass(X,X2).
  240. norm_first_pass(X=Y,X2+Z*Y2) :-
  241.     !,
  242.     Z= -1,
  243.     norm_first_pass(X,X2),
  244.     norm_first_pass(Y,Y2).
  245. norm_first_pass(X+Y,X2+Y2) :-
  246.     !,norm_first_pass(X,X2),
  247.     norm_first_pass(Y,Y2).
  248. norm_first_pass(X*Y,X2*Y2) :-
  249.     !,norm_first_pass(X,X2),
  250.     norm_first_pass(Y,Y2).
  251. norm_first_pass(X-Y,X2+Z*Y2) :-
  252.     !,Z= -1,
  253.     norm_first_pass(X,X2),
  254.     norm_first_pass(Y,Y2).
  255.  
  256.  
  257. %%%%%%%%%%%%%%%%%%%%%%%%
  258. % this part (which has been entirely rewritten since the original version of
  259. % this  program -- they are now mainly written in a functionnal way) deal with
  260. % embedded sums and products. The order of the term being vital for the solver,
  261. % our "philosophy" (may Socrates pardon me this offense...)  has been to
  262. % create  and modify the data structure (i.e. mainly lists that represent
  263. % equations)  such that they always stay sorted.
  264.  
  265.  
  266. norm_dist(A:real) -> [A].
  267. norm_dist(T:(real*vr)) -> [T].
  268. norm_dist(T1*T2) -> merge_lists(mult_list_list(norm_dist(T1),norm_dist(T2))).
  269. norm_dist(T1+T2) -> special_merge(norm_dist(T1),norm_dist(T2)).
  270.  
  271. non_strict(varless)?
  272. varless(vr(X),vr(Y)) -> X<Y.
  273. varless(@*vr(X),B) -> varless(vr(X),B).
  274. varless(A,@*vr(Y)) -> varless(A,vr(Y)).
  275. varless(@,real) -> true.
  276. varless(@,@) -> false.
  277.  
  278. non_strict(termless_pred) ?
  279. termless_pred(real,real,false).
  280. termless_pred(@,real,true).
  281. termless_pred(real,@,false).
  282. termless_pred(T1*V,T2*V,B) :- termless_pred(T1,T2,B).
  283. termless_pred(@*vr(X),@*vr(Y),B) :- B = X < Y.
  284.  
  285. termless(T1,T2) -> B | termless_pred(T1,T2,B).
  286.  
  287. non_strict(mult_term_var) ?
  288. mult_term_var(V:vr,A:real) -> `(A*V).
  289. mult_term_var(V1:vr,A*V2) -> cond(varless(V1,V2),`(A*V2*V1),
  290.                   `(T*V2)) | 
  291.         T=mult_term_var(V1,A).
  292.  
  293. non_strict(mult_term_real) ?
  294. mult_term_real(A:real,B:real) -> A*B.
  295. mult_term_real(A,T*V) -> `(T2*V) | T2=mult_term_real(A,T).
  296.  
  297. non_strict(mult_list) ?
  298. mult_list(A:real,L) -> mult_list_real(A,L).
  299. mult_list(V:vr,L) -> mult_list_var(V,L).
  300. mult_list(A*V:vr,L) -> T |
  301.         F = mult_list_var(V,L),
  302.         T = mult_list(A,F).
  303.  
  304. non_strict(mult_list_real) ?
  305. mult_list_real(@,[]) -> [].
  306. mult_list_real(A,[B:real|Tl]) -> [A*B | mult_list_real(A,Tl)].
  307. mult_list_real(A,[V:vr|Tl]) -> [`(A*V) | mult_list_real(A,Tl)].
  308. mult_list_real(A,[T*V|Tl]) -> [`(X*V) | mult_list_real(A,Tl)] | 
  309.         X = mult_term_real(A,T).
  310.  
  311. mult_list_var(@,[]) -> [].
  312. mult_list_var(V,L:[A|Tl]) -> cond(varless(A,V),
  313.                   [mult_term_var(V,A)| mult_list_var(V,Tl)],
  314.                   [`(A*V)| mult_list_var(V,Tl)]).
  315.  
  316. mult_list_list([],@)-> [].
  317. mult_list_list([X|Tl], Liste) -> [mult_list(X,Liste)|
  318.                   mult_list_list(Tl,Liste)].
  319.  
  320. special_merge([],L) -> L.
  321. special_merge(L,[]) -> L.
  322. special_merge(L1:[A|Tl1],L2:[B|Tl2]) ->
  323.     cond(termless(A,B),
  324.          [A|special_merge(Tl1,L2)],
  325.          (cond(termless(B,A),
  326.           [B|special_merge(L1,Tl2)],
  327.           [add_terms(A,B)|special_merge(Tl1,Tl2)]))).
  328.     %| write(termless(A,B)," | ",termless(B,A),"  ",A," + ",B),nl)).
  329.  
  330. non_strict(add_terms)?
  331. add_terms(A:real,B:real) -> A+B.
  332. add_terms(A*V:vr,B*V) -> `(T*V) | T=add_terms(A,B).
  333.  
  334. merge_lists([]) -> []. 
  335. merge_lists([X|Tl]) ->  special_merge(X,merge_lists(Tl)).
  336.  
  337. non_strict(solve_eq,mult_eq,div_eq,do_subst,wake_up)?
  338.  
  339. solve_eq([0]) :- !.                        %succes
  340. solve_eq([real]) :- !,fail.                    %failure
  341. solve_eq([A*X:vr(val=>none)|L]) :- !,inst_var(X,L,A).            %"pivot"
  342. solve_eq([A*vr(val=>V)|L]) :- ML=mult_eq(V,A),            %substit
  343.                   slist_merge(ML,L,NL), 
  344.                   solve_eq(NL).
  345.  
  346. div_eq([A:real],B) -> [-A/B].
  347. div_eq([A*X|Tl],B) -> [`(Z*X)|div_eq(Tl,B)] | Z= -A/B.
  348.  
  349. mult_eq([A:real],B) -> [A*B].
  350. mult_eq([A*X|Tl],B) -> [`(Z*X)|mult_eq(Tl,B)] | Z=A*B.
  351.  
  352. inst_var(X:vr(val=>V,wlist=>W),Eq,A) :- 
  353.     do_subst(Eq,NewEq),
  354.     V <- div_eq(NewEq,A),
  355.     I= -1,
  356.     wake_up(W,[I*X|V]),
  357.     update_vars(X,V),
  358.     W <- [].
  359.  
  360. do_subst(X:[real],X) :- !.
  361. do_subst([T:(A*vr(val=>none))|L],[T|NL]) :- 
  362.     !,do_subst(L,NL).
  363. do_subst([T:(A*vr(val=>V))|Tl],NV) :- 
  364.     ML=mult_eq(V,A),
  365.     slist_merge(ML,Tl,NL),
  366.     do_subst(NL,NV).
  367.  
  368. wake_up([],@) :- !.
  369. wake_up([gauss_wup(lhs_var=>X,coef=>A)|WL],Val) :-
  370. %    nl,write(X),nl,write(WL),nl,
  371.     recompute(X,A,Val),
  372.     wake_up(WL,Val).
  373.  
  374. recompute(X:vr(Id,name=>N,val=>Oldval),Coef,Sval) :-
  375.     ML=mult_eq(Sval,Coef),
  376.     merge_and_clean(Oldval,ML,Newval,Upd,X),
  377.     update_vars(X,Upd),
  378.     X <- vr(Id,name=>N,val=>Newval,wlist=>[]).
  379.  
  380. non_strict(update_vars)?
  381.  
  382. update_vars(@,[]) :- !.
  383. update_vars(@,[real]) :- !.
  384. update_vars(X,[A*Y:vr(Id,name=>N,wlist=>WL)|Tl]) :-
  385.     Y <- vr(Id,name=>N,
  386.         val => none,
  387.         wlist => [gauss_wup(lhs_var=>X, coef=>root_sort(A)) | WL]),
  388.     update_vars(X,Tl).
  389.  
  390. non_strict(slist_merge,merge_and_clean)?
  391.  
  392. slist_merge([],L,L) :- !.
  393. slist_merge(L,[],L) :- !.
  394. slist_merge([A:real],[B:real],[C]) :-  !,C=A+B.
  395. slist_merge([A|L1],L2:[B|@],[A|L3]) :-
  396.     varless(A,B), !,
  397.     slist_merge(L1,L2,L3).
  398. slist_merge(L1:[A|@],[B|L2],[B|L3]) :-
  399.     varless(B,A),!,
  400.     slist_merge(L1,L2,L3).
  401. slist_merge([(A:real)*(V:vr(Id))|L1],[(B:real)*vr(Id)|L2],Z) :-
  402.     C= A+B,
  403.     Z=cond(null(C),L3,[`(C*V)|L3]),
  404.     slist_merge(L1,L2,L3).
  405.  
  406. merge_and_clean([],L,L,L,@) :- !.
  407. merge_and_clean(L,[],L,[],@) :- !.
  408.  
  409. merge_and_clean([A:real],[B:real],[C],[],@) :-  !,C=A+B.
  410.  
  411. merge_and_clean([A|L1],L2:[B|@],[A|L3],Upd,Lhs) :-
  412.     varless(A,B), !,
  413.     merge_and_clean(L1,L2,L3,Upd,Lhs).
  414.  
  415. merge_and_clean(L1:[A|@],[B|L2],[B|L3],[B|Upd],Lhs) :-
  416.     varless(B,A),!,
  417.     merge_and_clean(L1,L2,L3,Upd,Lhs).
  418.  
  419. merge_and_clean([(A:real)*(V:vr(Id))|L1],
  420.         [(B:real)*vr(Id,wlist=>WL)|L2],L3,Upd,Lhs) :-
  421.     A+B=0,!,
  422.     clean_wlist(WL,Lhs),
  423.     merge_and_clean(L1,L2,L3,Upd,Lhs).
  424.  
  425. merge_and_clean([(A:real)*(V:vr(Id))|L1],
  426.         [(B:real)*vr(Id,wlist=>WL)|L2],[T|L3],[T|Upd],Lhs) :-
  427.     C= A+B,
  428.     T=`(C*V),
  429.     clean_wlist(WL,Lhs),
  430.     merge_and_clean(L1,L2,L3,Upd,Lhs).
  431.  
  432. null(0)->true.
  433. null -> false.
  434.  
  435. % les deux premiers arguments de s_merge sont les slists a combiner,
  436. % le troisieme est le resultat, le 4eme la liste des sous_listes du resultat
  437. % dont le premier element vient de la premiere liste. En d'autres termes, c'est
  438. % la liste des cellules de listes que l'on a creees.
  439.  
  440. s_merge([A|L1],L2:[B|@],L3:[A|L2],[L3|PtL]) :-
  441.     varless(A,B), !,
  442.     s_merge2(L1,L3,PtL).
  443. s_merge(L1,L2,L2,PtL) :- s_merge2(L1,[@|L2],PtL).
  444.  
  445. s_merge2([],@,[]) :- !.    
  446. s_merge2(L,[@|X:[]],L) :- !,X <- L.
  447. s_merge2([A:real],[@,B:real],[]) :- !,B <- A+B.
  448. s_merge2([A|L1],LP:[U|L2:[B|@]],[L3|PtL]) :-
  449.     varless(A,B), !,
  450.     LP <- [U|L3:[A|L2]],
  451.     s_merge2(L1,L3,PtL).
  452.  
  453. non_strict(s_merge2)?
  454.  
  455. s_merge2(L1:[A|@],[@|L2:[B|@]],PtL) :-
  456.     varless(B,A), !,
  457.     s_merge2(L1,L2,PtL).
  458. s_merge2([(A:real)*vr(Id,wlist=>WL)|L1],LP:[@|L2:[(B:real)*vr(Id)|@]],PtL) :-
  459.     A+B = 0, !,
  460.     clean_wlist(WL,L2),
  461.     del(L2),
  462.     s_merge2(L1,LP,PtL).
  463. s_merge2([(A:real)*vr(Id)|L1],[@|L2:[(B:real)*vr(Id)|@]],PtL) :-
  464.     B <- A+B,
  465.     s_merge2(L1,L2,PtL).
  466.  
  467. clean_wlist([],@).
  468. clean_wlist(WL:[gauss_wup(lhs_var=>V)|Tail],V) :- !,WL <- Tail,
  469.                           clean_wlist(WL,V).
  470. clean_wlist([@|WL],V) :- clean_wlist(WL,V).
  471.