% From lim@scorpio Thu Jun 17 14:09:36 1993
% adapted for CHRs by thom fruehwirth 930617
% fourier very slow, maybe loops with beale/1, opt1/1 loops, opt2/1 works

rmin(E):- M=:=E. % thom: no optimisation
rmax(E):- (-M)=:=E. % thom: no optimisation

%
beale([X1,X2,X3,X4,X5,X6,X7]) :-
    X1 + 1/4 * X4 - 8 * X5 - X6 + 9 * X7 =:= 0,
    X2 + 1/2 * X4 - 12 * X5 - 1/2 * X6 + 3 * X7 =:= 0,
    X3 + X6 =:= 1,
    X1 >= 0,
    X2 >= 0,
    X3 >= 0,
    X4 >= 0,
    X5 >= 0,
    X6 >= 0,
    X7 >= 0,
    rmin( - 3/4 * X4 + 20 * X5 - 1/2 * X6 + 6* X7).

% topb(L,n)
topb(L,N) :-
    topb1([],L,0,N).

topb1(Li,Li,I,I) :- !.
topb1(Li,Lo,I,N) :-
    insertb(P,Li,Lt),
    J is I+1,
    putcons(Lt,1,J),
    topb1(Lt,Lo,J,N).

insertb(P,[],[P]).
insertb(P,[A|B],[P,A|B]).
insertb(P,[A|B],[A|C]) :-
    insertb(P,B,C).

putcons(_,M,N) :-
    M > N,
    !.
putcons([P|R],M,N) :-
    M0 is M - 1,
%   bwriteln(P > M0/N),
%   bwriteln(P < M/N),
    P > M0/N,
    P < M/N,
    M1 is M + 1,
    putcons(R,M1,N).

bwriteln(X) :-
    writeln(X).
bwriteln(X) :-
    writeln(delete(X)).

% may loop
fib(0,1).
fib(1,1).
fib(N,Z) :-
    Z =:= X1 + X2,
    N1 =:= N-1,
    N2 =:= N-2,
    fib(N1,X1),
    fib(N2,X2).


laplace([_, _]) :- !.
laplace([H1, H2, H3|T]):-
        laplace_vec(H1, H2, H3),
        laplace([H2, H3|T]).

laplace_vec([_, _], [_, _], [_, _]) :- !.
laplace_vec([_TL, T, TR|T1], [ML, M, MR|T2], [_BL, B, BR|T3]):-
        B + T + ML + MR - 4 * M =:= 0,
        laplace_vec([T, TR|T1], [M, MR|T2], [B, BR|T3]).

%
laplace5(M) :-
    M = [
         [0,0,0,0,0],
         [100,R,S,T,100],
         [100,U,V,W,100],
         [100,X,Y,Z,100],
         [100,100,100,100,100]
         ],
    laplace(M).

% [chipc]: laplace7(X).
%
% X = [[0, 0, 0, 0, 0, 0, 0], [100, (5260/99), (63625/1716), (42545/1287), (63625/1716), (5260/99), 100], [100, (388405/5148), (2050/33), (149485/2574), (2050/33), (388405/5148), 100], [100, (1125/13), (2025/26), (75), (2025/26), (1125/13), 100], [100, (477845/5148), (2900/33), (221765/2574), (2900/33), (477845/5148), 100], [100, (9590/99), (162425/1716), (120805/1287), (162425/1716), (9590/99), 100], [100, 100, 100, 100, 100, 100, 100]]

%
laplace7(M) :-
    M = [
         [0,0,0,0,0,0,0],
         [100,R11,R12,R13,R14,R15,100],
         [100,R21,R22,R23,R24,R25,100],
         [100,R31,R32,R33,R34,R35,100],
         [100,R41,R42,R43,R44,R45,100],
         [100,R51,R52,R53,R54,R55,100],
         [100,100,100,100,100,100,100]
        ],
    laplace(M).

% [chipc]: chipOpt(X,Y,Z).
% 
% X = (8/5)
% Y = (6/5)
% Z = (14/5)

%
chipOpt(X1,X2,X3) :-
    X1 + 2 * X2 =< 4,
    3 * X1 + X2 =< 6,
    X3 =:= X1 + X2,
    rmax(X3).

% chipfact(n,1,N)
chipfact(X,Y,M) :-
    X =:= 0,
    !,
    Y =:= M.
chipfact(X,Y,M) :-
    X1 =:= X - 1,
    M1 =:= X * Y,
    chipfact(X1,M1,M).

% order of magnitude slower than chipfact/3
fact(0,1).
fact(1,1).
fact(N,R) :- 1 < N, N =< R, M =:= N-1, fact(M,T), R =:= N * T.

%
mg(P,T,I,B,MP):-
    T > 0,
    T =< 1,
    B + MP =:= P * (1 + I/100).
mg(P,T,I,B,MP):-
    T > 1,
    I1 =:= I / 100,
    T1 =:= T -1,
    P2 =:= P * (1 + I1) - MP,
    mg(P2, T1, I, B, MP).

mg1(X,Y,Z) :-
    2 =:= T,
    1 =:= I,
    T > 1,
    I1 =:= I / 100,
    T1 =:= T -1,
    P2 =:= P * (1 + I1) - MP,
    T1 > 0,
    T1 =< 1,
    B + MP =:= P2 * (1 + I/100).

% [chipc]: top0(X,Y,Z).
%
% X = _r94 
% Y = (-9462212541120451001/1000000000000000000) * _r94  + (10462212541120451001/1000000000000000000) * _r90 
% Z = (101/100) * _r94  + (-1) * _r90      More? (;) 
%
% ---------------------------------------------------------
% [chipc]: top1(X,Y,Z).
%
% X = _r94 
% Y = (-101/100) * _r94  + (201/100) * _r90 
% Z = (101/100) * _r94  + (-1) * _r90      More? (;) 
% ---------------------------------------------------------
% [chipc]: top(X,Y,Z).
% 
% X = _r94 
% Y = (-2101900399479668244827490915525641902001/100000000000000000000000000000000000000) * _r94  + (2201900399479668244827490915525641902001/100000000000000000000000000000000000000) * _r90 
% Z = (101/100) * _r94  + (-1) * _r90      More? (;) 
%

%
top0(P, B, MP) :- mg(P,10,1,B,MP).
top1(P, B, MP) :- mg(P,2,1,B,MP).
top(P, B, MP) :- mg(P,20,1,B,MP).
  

% Detection of Implied Equalities -------------------------------------------

%
top4([A,B,C,D]) :-
        A=<B, 
	B=<C, 
	C=<D, 
	A>=D. 

%     B = A,
%     C = A,
%     D = A  
/*
L = [A_m108, B_m128, C_m336, D_m660] 
 
Constraints:
eq0([B_m128 * 1, D_m660 * -1], 0)
eq0([C_m336 * 1, D_m660 * -1], 0)
eq0([A_m108 * -1, D_m660 * 1], 0)
*/ 

%
% X = [(1/3), (0), (13/3)]
%
%
opt1([X1,X2,X3]) :-
    X1 + X2 + 2 * X3 =< 9,
    X1 + X2 - X3     =< 2,
   -X1 + X2 + X3     =< 4,
    X1 >= 0,
    X2 >= 0,
    X3 >= 0,
    rmin(X1 + X2 - 4 * X3).

%
% X = [(0), (0)]
%
%
opt2([X1,X2]) :-
    X1 + 2 * X2 =< 4,
    X2 =< 1,
    X1 >= 0,
    X2 >= 0,
    rmin(X1 + X2).

available_res(10).
available_res(14).
available_res(27).
available_res(60).
available_res(100).
available_cell(10).
available_cell(20).

ohm(V,I,R) :-
%    bwriteln(V =:= I * R),
    V =:= I * R.

sum([],Z) :- 
%    bwriteln(Z =:= 0),
    Z =:= 0.
sum([H|T],N) :-
%    bwriteln(N =:= H + M),
    N =:= H + M,
    sum(T,M).

kirchoff(L) :-
    sum(L,0).

% X = [(200/37), (540/37)]     More? (;) 
% 
% X = [(140/37), (600/37)]     More? (;) 
% 
% X = [(540/127), (2000/127)]     More? (;) 

%
ohm_example([V1,V2]) :-
    29/2 < V2, V2 < 65/4,
    available_res(R1),
    available_res(R2),
    available_cell(V),
    ohm(V1,I1,R1), ohm(V2,I2,R2),
    kirchoff([I1,-I2]), kirchoff([-V,V1,V2]).

% X = [10, 27, 20]     More? (;) 
% 
% X = [14, 60, 20]     More? (;) 
% 
% X = [27, 100, 20]     More? (;) 
%
ohm_example1([R1,R2,V]) :-
    29/2 < V2, V2 < 65/4,
    available_res(R1),
    available_res(R2),
    available_cell(V),
    ohm(V1,I1,R1), ohm(V2,I2,R2),
    kirchoff([I1,-I2]), kirchoff([-V,V1,V2]).

%
% X = [(14), (60), (20)]
%
%
ohm1([A,B,C]) :-
    A =:= 14,
    B =:= 60,
    C =:= 20,
    29/2 < V2, V2 < 65/4,
    V1/A - V2/B =:= 0,
    V1 + V2 =:= C.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
example( [X0,X1,X2,X3,X4]) :-
        87*X0  +52*X1  +27*X2  -54*X3  +56*X4 =<  -93,
        33*X0  -10*X1  +61*X2  -28*X3  -29*X4 =<   63,
       -68*X0   +8*X1  +35*X2  +68*X3  +35*X4 =<  -85,
        90*X0  +60*X1  -76*X2  -53*X3  +24*X4 =<  -68,
       -95*X0  -10*X1  +64*X2  +76*X3  -24*X4 =<   33,
        43*X0  -22*X1  +67*X2  -68*X3  -92*X4 =<  -97,
        39*X0   +7*X1  +62*X2  +54*X3  -26*X4 =<  -27,
        48*X0  -13*X1   +7*X2  -61*X3  -59*X4 =<   -2,
        49*X0  -23*X1  -31*X2  -76*X3  +27*X4 =<    3,
       -50*X0  +58*X1   -1*X2  +57*X3  +20*X4 =<    6,
       -13*X0  -63*X1  +81*X2   -3*X3  +70*X4 =<   64,
        20*X0  +67*X1  -23*X2  -41*X3  -66*X4 =<   52,
       -81*X0  -44*X1  +19*X2  -22*X3  -73*X4 =<  -17,
       -43*X0   -9*X1  +14*X2  +27*X3  +40*X4 =<   39,
        16*X0  +83*X1  +89*X2  +25*X3  +55*X4 =<   36,
        +2*X0  +40*X1  +65*X2  +59*X3  -32*X4 =<   13,
       -65*X0  -11*X1  +10*X2  -13*X3  +91*X4 =<   49,
        93*X0  -73*X1  +91*X2   -1*X3  +23*X4 =<  -87.



top2 :- example( [X0,X1,X2,X3,X4]).

%     X3<=-5/4-35/68*X2-2/17*X1+X0-35/68*X4,
%     X3<=68/53-76/53*X2+60/53*X1+90/53*X0+24/53*X4,
%     X3<=-1/2-31/27*X2-7/54*X1-13/18*X0+13/27*X4,
%     X3<=17/22+19/22*X2-2*X1-81/22*X0-73/22*X4,
%     X3<=33/76-16/19*X2+5/38*X1+5/4*X0+6/19*X4,
%     X3>=87+91*X2-73*X1+93*X0+23*X4,
%     X3>=-3/76-31/76*X2-23/76*X1+49/76*X0+27/76*X4,
%     X3<=13/9-14/27*X2+1/3*X1+43/27*X0-40/27*X4,
%     X3<=2/19+1/57*X2-58/57*X1+50/57*X0-20/57*X4 

top3 :- example( [X0,_,_,_,X4]).

%     X0>=477804/40409+6973307/969816*X4,
%     X0>=7357764/4517605-5006476/13552815*X4,
%     X0>=58416/36205-4659804/12418315*X4,
%     X0>=3139326/1972045-745308/1972045*X4,
%     X0>=67158/43105-16394/43105*X4,
%     X0>=1327097/6210451-2619277/6210451*X4,
%     X0<=-688135/1217232-2174029/811488*X4 

% [chipc]: top5(X).
% 
% X = [(0), (0), (0), (1000), (0), (0), (50), (200/9) + (-1/9) * _rp522 , (1000/3), (5350/9) + (1/9) * _rp522 , (4000), (0), (7450/3), (0), (0), (5000), (250), (600), (1000/9) + (-5/9) * _rp522 , (0)]
% yes.

top5([X11,X12,X13,X14,X15,X21,X22,X23,X24,X25,Y21,Y22,Y23,Y24,Y25,Z21,Z22,Z23,Z24,Z25]) :-   
              X11 + X12 + X13 + X14 + X15 =:= 1000,
              X21 + X22 + X23 + X24 + X25 =:= 1000,
            
              4*X11 + 5*X21 - Y21 - Z21 =< 0,
             -4*X12 - 5*X22 + Y22 + Z22 =:= 0,
             -4*X13 - 5*X23 + Y24 - Y25 + Z24 - Z25 =:= 0,
             -4*X14 - 5*X24 + Y21 - Y22 - Y23 + Y25
                            + Z21 - Z22 - Z23 + Z25 =:= 0,
             -4*X15 - 5*X25 + Y23 - Y24 + Z23 - Z24 =:= 0,
            
              7*X11 + 9*X21 >= 0,
              7*X12 + 9*X22 =< 3000,
              7*X13 + 9*X23 =< 200,
              7*X14 + 9*X24 =< 10000,
              7*X15 + 9*X25 =< 7000,
            
              Z21 =< 5000,
              Z22 =< 250,
              Z23 =< 600,
              Z24 =< 7000,
              Z25 =< 4000,
            
              X11 >= 0,
              X12 >= 0,
              X13 >= 0,
              X14 >= 0,
              X15 >= 0,
              X21 >= 0, X22 >= 0, X23 >= 0, X24 >= 0, X25 >= 0,
            
              Y21 >= 0, Y22 >= 0, Y23 >= 0, Y24 >= 0, Y25 >= 0,
            
              Z21 >= 0, Z22 >= 0, Z23 >= 0, Z24 >= 0, Z25 >= 0,
            
              M =:= 99999,
              - Min =:= 99999 * X11 + 99999 * X21 + 4 * Y21 + 7 * Y22 + 
			3 * Y23 +  8*Y24 + 5*Y25,
              rmax(Min).

%
top5a(List) :-
	List = [X11,X12,X13,X14,X15,X21,X22,X23,X24,X25,Y21,Y22,Y23,Y24,Y25,Z21,Z22,Z23,Z24,Z25],
              X11 + X12 + X13 + X14 + X15 =:= 1000,
              X21 + X22 + X23 + X24 + X25 =:= 1000,
            
              4*X11 + 5*X21 - Y21 - Z21 =< 0,
             -4*X12 - 5*X22 + Y22 + Z22 =:= 0,
             -4*X13 - 5*X23 + Y24 - Y25 + Z24 - Z25 =:= 0,
             -4*X14 - 5*X24 + Y21 - Y22 - Y23 + Y25
                            + Z21 - Z22 - Z23 + Z25 =:= 0,
             -4*X15 - 5*X25 + Y23 - Y24 + Z23 - Z24 =:= 0,
            
              7*X11 + 9*X21 >= 0,
              7*X12 + 9*X22 =< 3000,
              7*X13 + 9*X23 =< 200,
              7*X14 + 9*X24 =< 10000,
              7*X15 + 9*X25 =< 7000,
            
              Z21 =< 5000,
              Z22 =< 250,
              Z23 =< 600,
              Z24 =< 7000,
              Z25 =< 4000,
            
              X11 >= 0,
              X12 >= 0,
              X13 >= 0,
              X14 >= 0,
              X15 >= 0,
              X21 >= 0, X22 >= 0, X23 >= 0, X24 >= 0, X25 >= 0,
            
              Y21 >= 0, Y22 >= 0, Y23 >= 0, Y24 >= 0, Y25 >= 0,
            
              Z21 >= 0, Z22 >= 0, Z23 >= 0, Z24 >= 0, Z25 >= 0,
            
              M =:= 99999,
	      Min = 23450,
              - Min =:= 99999 * X11 + 99999 * X21 + 4 * Y21 + 7 * Y22 + 
			3 * Y23 +  8*Y24 + 5*Y25.

%       M = 99999,
%       Min = 23450,
%       X11 = 0,
%       X12 = 0,
%       X13 = 0,
%       X14 = 1000,
%       X15 = 0,
%       X21 = 0,
%       X22 = 50,
%       X23 = 1850/3-X25,
%       X24 = 1000/3,
%       Y21 = 4000,
%       Y22 = 0,
%       Y23 = 7450/3,
%       Y24 = 0,
%       Y25 = 0,
%       Z21 = 5000,
%       Z22 = 250,
%       Z23 = 600,
%       Z24 = 9250/3-5*X25,
%       Z25 = 0,
%       X25>=5350/9,
%       X25<=1850/3 




% ===========================================================================