Recursive program to compute Taylor approx of cosi

Question

I'm still pretty new to Prolog, and i'm not sure why this code isn't working. I believe it is most likely a problem with the base case or in the last 3 lines of the recursive case. Everything else works just fine.

This program determines cosine calculated with series approximation,

to do so it needs to calculate the factorial of 2K, also -1 ^ K, and then uses these 2 calculations in the final equation (this is done in % Recursive Case).

% Factorial from class
fact(0, 1).
fact(N, F) :- 
    N > 0,
    N1 is N-1,
    fact(N1, F1),
    F is F1 * N.

% Calculate -1 ^ K
signCnt(0,1).
signCnt(K,S) :- 
    K > 0,
    K1 is K - 1,
    signCnt(K1,S1),
    S is S1 * -1.

% Base case
cosN(N,_,_,0).

% Recursive case
cosN(K,N,X,Y) :- K < N,
    signCnt(K,S),
    K2 is 2 * K,
    fact(K2,F),
    Yk is (S * X**K2)/F,
    K1 is K + 1,
    cosN(K1,N,X,Y1),
    Y is Y1 + Yk.

cosN(N,X,Y) :- 
    N>0,
    cosN(0,N,X,Y).

Inputs should be in the form

?- cosN(25,pi,Y).

with an expected output of

Y = -1.0 ;
false.

however, it doesn't go through the recursion properly and the output ends up looking like this:

where 5 and pi could be anything so long as pi remains in form pi (i.e. pi/2, pi/3), also there shouldn't be any additional lines added, as we were given a line number restriction. Lines should be edited/replaced. Anything to point me in the right direction would be greatly appreciated too.

(Thank you to Guy Coder for help formatting)

---

Edit by Guy Coder

Some test cases using SWI-Prolog

:- begin_tests(cosine_approximation).

factorial_test_case_generator(0,1).
factorial_test_case_generator(1,1).
factorial_test_case_generator(2,2).
factorial_test_case_generator(3,6).
factorial_test_case_generator(4,24).
factorial_test_case_generator(5,120).
factorial_test_case_generator(6,720).
factorial_test_case_generator(7,5040).
factorial_test_case_generator(8,40320).
factorial_test_case_generator(20,2432902008176640000).

test('factorial',[nondet,forall(factorial_test_case_generator(N,Factorial))]) :-
    fact(N,Factorial).

signCnt_test_case_generator(0,1).
signCnt_test_case_generator(1,-1).
signCnt_test_case_generator(2,1).
signCnt_test_case_generator(3,-1).
signCnt_test_case_generator(4,1).
signCnt_test_case_generator(5,-1).

test('signCnt',[nondet,forall(signCnt_test_case_generator(N,Sign))]) :-
    signCnt(N,Sign).

:- end_tests(cosine_approximation).

Example run:

?- make.
% c:/users/eric/documents/projects/prolog/so_question_161 compiled 0.00 sec, 5 clauses
% PL-Unit: cosine_approximation .......... done
% All 10 tests passed
true.

Answer 1

Just an addendum

I had to think about whether the program indeed sums small floats together into successively larger floats and not small floats onto larger floats (which might render the result more imprecise than it needs to), but it does.

Although it's inelegant to recompute the factorial fully at each element of the Taylor series and to not use -1 * (k mod 2) to obtain (-1)^k directly, instead going through recursion.

Here is the call diagram for orientation:

Addendum 2: Code for more efficient computation

So I availed myself to some time to perform the exercise of writing a cos approximation that just recurses on itself and carries all the ancillary information for computing the terms and and the sum.

% ===
% Entry point!
% Evaluate the Taylor series for cos(z) at "z" (not too far from 0, probably
% less than 1). The terms (sum elements) for index values 0..K are computed
5 and added. (K >= 0)
% ===

taylor_cos(Res,Z,Kmax,Verbose) :- 
   Zf is Z*1.0, % make a float
   float(Zf),
   integer(Kmax),Kmax >= 0,
   Zsq is Zf*Zf,
   at_element_k(Res,0,Kmax,Zsq,_,_,Verbose).

% The value computed is always the first one

even(K) :- integer(K), (K mod 2) =:= 0. % eval left & compare numerically
odd(K)  :- integer(K), (K mod 2) =:= 1. % eval left & compare numerically

% Compute (-1)^k, k an integer >= 0.
% Computed value is on first place in predicate argument list.

minus_one_tothe_k( 1,K) :- even(K),!. % ! to make this deterministic
minus_one_tothe_k(-1,K) :- odd(K).    % actually no need to test odd(K)

% Compute (2*k)!, k an integer >= 0, if (2*(k-1))! is known.
% Computed value is on first place in predicate argument list.
% The base case is conceptually jarring as the "prior value" can be anything.
% This is not unlike a function becoming evaluatable because of lazy evaluation.

two_times_k_factorial(1  ,0,_)        :- !.
two_times_k_factorial(Res,K,ResPrior) :- K>0, Res is ResPrior*K*(4*K-2).

% Compute (z^(2*k)), k an integer >= 0, if (z^(2*(k-1))) is known.
% z² is passed too so that we do not need to recompute it again and again.
% Computed value is on first place in predicate argument list.

z_tothe_2k(1,   0, _   ,_)        :- !.
z_tothe_2k(Res, K, Zsq ,ResPrior) :- K>0, Res is ResPrior * Zsq.

% Compute the Taylor series by summing the elements(k) with k in [0..Kmax)
% (so Kmax >= 1).
% When calling this initially, the values for TTKFprior and ZTT2Kprior
% are of no importance. 
% The procedures calls itself recursively to compute element(i), element(i+1)
% etc. based on prior intermediate values. The base case is attained when
% K > Kmax. The sum accumulates in SumFromKmaxBackwards when the recursion
% comes back up the stack.

at_element_k(0.0,K,Kmax,_,_,_,Verbose) :-
   K > Kmax,!,
   ((Verbose = verbose) -> 
   format("past the end as K=~d > Kmax=~d, returning back up the stack\n",[K,Kmax]) ; true).

at_element_k(SumFromKmaxBackwards,K,Kmax,Zsq,TTKFprior,ZTT2Kprior,Verbose) :- 
   minus_one_tothe_k(M1TTK,K),                 % M1TTK = (-1)^K
   two_times_k_factorial(TTKF,K,TTKFprior),    % TTKF  = f(K,TTKFprior)
   z_tothe_2k(ZTT2K,K,Zsq,ZTT2Kprior),         % ZTT2K = f(K,z²,ZTT2Kprior)
   ElementK is M1TTK * ZTT2K / TTKF,           % element_k = M1TTK * (ZTT2K / TTKF)
   ((Verbose = verbose) -> format("element(~d) = ~e\n",[K,ElementK]) ; true),
   KP1 is K+1,
   at_element_k(SumFromKmaxBackwardsPrior,KP1,Kmax,Zsq,TTKF,ZTT2K,Verbose),
   SumFromKmaxBackwards is SumFromKmaxBackwardsPrior + ElementK,
   ((Verbose = verbose) -> format("taylor-series-sum(~d ... ~d) = ~e (added ~e to prior value ~e)\n",
                                  [K,Kmax,SumFromKmaxBackwards, ElementK, SumFromKmaxBackwardsPrior]) ; true).

Run this! The Verbose variable is set to verbose to generate more printout during the Taylor series computation. We compute 11 terms of the series (indexes 0...10).

?- taylor_cos(Res,0.01,10,verbose).
element(0) = 1.000000e+00
element(1) = -5.000000e-05
element(2) = 4.166667e-10
element(3) = -1.388889e-15
element(4) = 2.480159e-21
element(5) = -2.755732e-27
element(6) = 2.087676e-33
element(7) = -1.147075e-39
element(8) = 4.779477e-46
element(9) = -1.561921e-52
element(10) = 4.110318e-59
past the end as K=11 > Kmax=10, returning back up the stack
taylor-series-sum(10 ... 10) = 4.110318e-59 (added 4.110318e-59 to prior value 0.000000e+00)
taylor-series-sum(9 ... 10) = -1.561920e-52 (added -1.561921e-52 to prior value 4.110318e-59)
taylor-series-sum(8 ... 10) = 4.779476e-46 (added 4.779477e-46 to prior value -1.561920e-52)
taylor-series-sum(7 ... 10) = -1.147074e-39 (added -1.147075e-39 to prior value 4.779476e-46)
taylor-series-sum(6 ... 10) = 2.087675e-33 (added 2.087676e-33 to prior value -1.147074e-39)
taylor-series-sum(5 ... 10) = -2.755730e-27 (added -2.755732e-27 to prior value 2.087675e-33)
taylor-series-sum(4 ... 10) = 2.480156e-21 (added 2.480159e-21 to prior value -2.755730e-27)
taylor-series-sum(3 ... 10) = -1.388886e-15 (added -1.388889e-15 to prior value 2.480156e-21)
taylor-series-sum(2 ... 10) = 4.166653e-10 (added 4.166667e-10 to prior value -1.388886e-15)
taylor-series-sum(1 ... 10) = -4.999958e-05 (added -5.000000e-05 to prior value 4.166653e-10)
taylor-series-sum(0 ... 10) = 9.999500e-01 (added 1.000000e+00 to prior value -4.999958e-05)
Res = 0.9999500004166653.

The 80-column mind of Stackoverflow is getting a bit on my nerves. We have a gazillion pixels of width on screens nowadays, and they are unused and left white because "Muh Visual Design"!! Anyway...

Now add some code to generate Count test floats evenly distributed between From and To. The generator/4 generates successive values on backtracking. The cos_compare/3 compares what our cos-approximating function computes and what the system computes (which comes somewhere deep down in a library).

generator(X,From,To,1) :- 
   From =< To,
   From_f is From*1.0,
   To_f   is To*1.0,
   X      is (From_f + To_f) / 2.0.

generator(X,From,To,Count) :- 
   integer(Count), 
   Count > 1,
   From =< To,
   From_f  is From*1.0,
   To_f    is To*1.0,
   Delta_f is (To_f - From_f)/(Count * 1.0),
   CountM1 is Count-1, 
   between(0,CountM1,I), 
   X is From_f + Delta_f*I.

cos_compare(Z,Kmax,Verbose) :-
   taylor_cos(Res,Z,Kmax,Verbose),
   Cos is cos(Z),
   Delta is abs(Res-Cos),
   format("For z = ~e, k_max = ~d, difference to real cos = ~e\n", [Z, Kmax, Delta]).

Then let's actually compare 100 values between float -4.0 and float +4.0 and, where we compute 11 terms (indexes 0..11) of the Taylor series at each value:

run(Verbose) :- forall(generator(Z,-4.0,+4.0,100), cos_compare(Z,10,Verbose)).

?- run(quiet).  
For z = -4.000000e+00, k_max = 10, difference to real cos = 1.520867e-08
For z = -3.920000e+00, k_max = 10, difference to real cos = 9.762336e-09
For z = -3.840000e+00, k_max = 10, difference to real cos = 6.209067e-09
For z = -3.760000e+00, k_max = 10, difference to real cos = 3.911487e-09
For z = -3.680000e+00, k_max = 10, difference to real cos = 2.439615e-09
......
For z = 3.680000e+00, k_max = 10, difference to real cos = 2.439615e-09
For z = 3.760000e+00, k_max = 10, difference to real cos = 3.911487e-09
For z = 3.840000e+00, k_max = 10, difference to real cos = 6.209067e-09
For z = 3.920000e+00, k_max = 10, difference to real cos = 9.762336e-09
true.

Not looking so bad.

Addendum 3: Using SWI-Prolog "dicts" to communicate between predicates

I have found that when writing Perl functions, it is often advantageous to short-circuit position-based argument passing and pass a single bunch of name-value pairs, aka 'hashes' instead. This adds a lot of flexibility (named parameters, easy to add parameters, easy to debug, easy to pass parameters to subfunctions etc.)

Let's try this here too.

This is restricted to SWI-Prolog because "dicts" are an SWI-Prolog feature. Code like this renders Prolog's indexing mechanism useless, as now every predicate has exactly the same argument, Dict, so should be relatively slow at runtime.

Just the modified predicates are

taylor_cos(Res,Z,Kmax,Verbose) :-
   Zf is Z*1.0, % make a float
   float(Zf),
   integer(Kmax),Kmax >= 0,
   Zsq is Zf*Zf,
   at_element_k(taylor{  sum     : Res  % the result
                        ,k       : 0
                        ,kmax    : Kmax
                        ,zsq     : Zsq
                        ,ttkf_prior  : _
                        ,ztt2k_prior : _
                        ,verbose : Verbose }).


% ---
% Base case, when k > kmax
% ---

% We map the passed "Dict" to a sub-Dict to grab values.
% As this is "unification", not only "pattern matching" the value for
% sum "0.0" is shared into "Dict".

at_element_k(Dict) :-
   taylor{  sum     : 0.0
           ,k       : K
           ,kmax    : Kmax
           ,verbose : Verbose } :< Dict,

   K > Kmax,  % guard
   !,         % commit
   ((Verbose = verbose) ->
      format("past the end as K=~d > Kmax=~d, returning back up the stack\n",[K,Kmax])
      ; true).

% ---
% Default case, when k <= kmax
% ---

% We map the passed "Dict" to a sub-Dict to grab values.
% We use ":<" instead of "=" so that, if the passed Dict has more values
% than expected (which can happen during program extension and fiddling),
% "partial unification" can still proceed, "=" would fail. However, no
% values can be missing!
% This gives us also the funny possibility of completely ignoring Kmax in
% the "input Dict", it doesn't appear anywhere and is still passed down
% through the recursive call. Well, it *does* appear because we print it
% out.

at_element_k(Dict) :-
   taylor{  sum         : SumFromKmaxBackwards  % the output value, to be captured by the caller
           ,k           : K                     % index of the current term/element in the Taylor sum
           ,kmax        : Kmax                  % max index for which a term/element will be computed
           ,zsq         : Zsq                   % z², a constant
           ,ttkf_prior  : TTKFprior             % prior "two times k factorial" i.e. (2*(k-1))!
           ,ztt2k_prior : ZTT2Kprior            % prior "z to the 2*k" i.e. z^(2*(k-1))
           ,verbose     : Verbose } :< Dict,    % emit messages about progress if Verbose = verbose

   minus_one_tothe_k(M1TTK,K),                       % compute (-1)^K
   two_times_k_factorial(TTKF,K,TTKFprior),          % compute (2*k)! based on prior value
   z_tothe_2k(ZTT2K,K,Zsq,ZTT2Kprior),               % compute z^(2*k) based on prior value
   ElementK is M1TTK * ZTT2K / TTKF,                 % compute value for Taylor sum term/element at k

   % (isn't there a better way to print conditionally?)

   ((Verbose = verbose) ->
      format("element(~d) = ~e\n",[K,ElementK])
      ; true),

   % create a NextDict from Dict for recursive call

   KP1 is K+1,
   put_dict( _{ sum        : SumFromKmaxBackwardsPrior
               ,k          : KP1
               ,ttkf_prior : TTKF
               ,ztt2k_prior: ZTT2K }, Dict, NextDict),

   % recursive call 
   % (foundational thought: the procedure is really a **channel-doing-computations between the series of dicts**)

   at_element_k(NextDict),

   % on return, complete summing the Taylor series backwards from highest index to the current index k

   SumFromKmaxBackwards is SumFromKmaxBackwardsPrior + ElementK,

   % (more conditional printing)

   ((Verbose = verbose) ->
      format("taylor-series-sum(~d ... ~d) = ~e (added ~e to prior value ~e)\n",
            [K,Kmax,SumFromKmaxBackwards,ElementK,SumFromKmaxBackwardsPrior])
      ; true).

Is it more readable? I suppose so.

Answer 2

Base case was wrong, should've been cosN(N,N,_,0). as K and N must both be equal to N when the program is finished the recursive process.

Test cases:

:- begin_tests(cosine_approximation).

factorial_test_case_generator(0,1).
factorial_test_case_generator(1,1).
factorial_test_case_generator(2,2).
factorial_test_case_generator(3,6).
factorial_test_case_generator(4,24).
factorial_test_case_generator(5,120).
factorial_test_case_generator(6,720).
factorial_test_case_generator(7,5040).
factorial_test_case_generator(8,40320).
factorial_test_case_generator(20,2432902008176640000).

test('factorial',[nondet,forall(factorial_test_case_generator(N,Factorial))]) :-
    fact(N,Factorial).

signCnt_test_case_generator(0,1).
signCnt_test_case_generator(1,-1).
signCnt_test_case_generator(2,1).
signCnt_test_case_generator(3,-1).
signCnt_test_case_generator(4,1).
signCnt_test_case_generator(5,-1).

test('signCnt',[nondet,forall(signCnt_test_case_generator(N,Sign))]) :-
    signCnt(N,Sign).

cosN_test_case_generator(3,pi/2,0.01996895776487828).
cosN_test_case_generator(5,pi,-0.9760222126236076).
cosN_test_case_generator(25,pi,-1.0).
cosN_test_case_generator(10,pi/2,-3.3306690738754696e-15).

test('cosN',[nondet,forall(cosN_test_case_generator(N,X,Y))]) :-
    cosN(N,X,Y).

:- end_tests(cosine_approximation).

Example run:

?- make.
% /Users/oliverclarke/prolog/lab5-quiz compiled 0.00 sec, 3 clauses
% PL-Unit: cosine_approximation .................... done
% All 20 tests passed
true.