3

I have some code that should find approximate value of sin(15°) by Taylor series definition and compare it to the built-in function sin.

I have different results.

  • X - radian value
  • R - current value of series sum
  • Eps - accuracy value
  • S - sign
  • F - factorial value of N
  • XN - value of X to the power of N
  • N - power, 1,3,5 ...
% base
sin_taylor(_,R,Eps,S,FN,XN,_) :-
    abs(S*XN/FN) < Eps,
    R = 0, !.

% step
sin_taylor(X, R, Eps, S, FN, XN, N) :-
    S1 = S*(-1),
    N1 = N+1,
    FN1 = FN*(2*N+2)*(2*N+3),
    XN1 = XN*X*X,
    sin_taylor(X,R1,Eps,S1,FN1,XN1,N1),
    R is R1+S*XN/FN.

% auxiliary predicate to supply parameters to the main one
sin(X,R,Eps) :-
    sin_taylor(X,R1,Eps,-1,1,X*X,1),
    R is 1+R1.

Results in the console:

?- X is 15*(pi / 180), sin(X,R,0.0001).
   X = 0.2617993877991494,
   R = 0.931695959721973.

?- X is 15*(pi / 180), R0 is sin(X).
   X = 0.2617993877991494,
   R0 = 0.25881904510252074.
Will Ness
  • 70,110
  • 9
  • 98
  • 181
Michael
  • 339
  • 2
  • 11
  • @RobertBaron, I guess it's evident: what am I doing wrong? What's the mistake? – Michael Aug 23 '19 at 11:56
  • Sorry, @RobertBaron, I thought I did my best to explain every variable I used. I think that I made a mistake in recursion step part of code, but I can't figure out where exactly. If my explanation of variables is not so clear, please, tell me what should be explained thoroughly so you can help me. Thanks – Michael Aug 23 '19 at 13:29
  • 4
    There is a problem with expressions like `XN1=XN*X*X` as these are not computation on the right hand and assignments on the left hand but unifications. The variable XN1 and the tree `*(XN,*(X,X))` are unified. I don't see how this works ... `is` should be used. – David Tonhofer Aug 23 '19 at 22:19
  • 3
    @DavidTonhofer: it works (somewhat) because expressions are passed down the recursion, and the final is/2 compute them all together – CapelliC Aug 24 '19 at 08:51

2 Answers2

1

Cannot recognize the correct nth-term in your code... where is the factorial computation hidden? Remember that early time optimizations are roots of every evils in SW engineering :)

A straightforward translation of the solution listed in Wikipedia lead me to this code (well, should cite Maclaurin and Taylor in predicate names...)

:- module(sin_taylor,
          [sin/3
          ,rad_deg/2
          ,fact/2
          ]).

% help predicate to give parameters to the main one
sin(X,R,Eps):-
    syn_taylor(X,Eps,0,R).

syn_taylor(X,Eps,N,R) :-
    S is -1**N,
    T is 2*N+1,
    fact(T,D),
    E is X**T,
    Q is S*E/D,
    (   Q < Eps
    ->  R = Q
    ;   M is N+1,
        syn_taylor(X,Eps,M,R1),
        R is Q+R1
    ).

rad_deg(R,D) :-
    var(R) -> R is D*(pi / 180).
    % tbd compute D from R

fact(N,F) :- N>0 -> N1 is N-1, fact(N1,F1), F is N*F1 ; F=1.

that gives reasonable results:

?- rad_deg(X,15),sin(X,R,0.0001).
X = 0.2617993877991494,
R = 0.2588088132736575.

?- rad_deg(X,15),R0 is sin(X).
X = 0.2617993877991494,
R0 = 0.25881904510252074.
CapelliC
  • 59,646
  • 5
  • 47
  • 90
  • avoiding the direct calculation of factorials in Taylor series terms is not an optimization, though. :) it is a correctness issue. the floating point error goes haywire for all arguments except those very close to zero. an error term of 1e-13 vs 1e3 (and even 1e-3) is a qualitative, not a quantitative change. – Will Ness Aug 24 '19 at 09:40
  • @WillNess: are you referring to the fact that the factorial is used in denominator ? I don't get it entirely... I was referring to the lack of explicit factorial in OP' code, that could be the main culprit of his error - but I'm not sure. I haven't debugged neither mine or his code.... – CapelliC Aug 24 '19 at 09:57
  • the further from 0, the more terms will be needed to get to a preset precision (tolerance; error) of the calculation; the higher the term the bigger the factorial so it could cause the 16 significant digits runoff, and overall error of the floating point calculation will become intolerable. it has been my experience that we can avoid this by calculating T(n+1) from T(n) instead (and not directly) and get to the areas much further from 0 while inside the allowed error. if we want to calculate the result to its 10th digit, say, we can't afford the factorial to get bigger than 10^6 in magnitude. – Will Ness Aug 24 '19 at 12:22
  • so, yes. :) the direct formula has us calculating two values, E and D, potentially very large in magnitude, then dividing the two. Instead, finding the ratio of T(n+1)/T(n) will cause calculation of values much smaller in magnitude. (cf. https://en.wikipedia.org/wiki/Precision_(computer_science); https://en.wikipedia.org/wiki/Round-off_error; https://en.wikipedia.org/wiki/Numerical_analysis; https://en.wikipedia.org/wiki/Approximation; etc.) – Will Ness Aug 24 '19 at 12:27
1

You get better results, less numerical errors, if you start with the smallest summand and then use a Horner schema. But the problem is, how to figure out which index k should give you the smallest summand ak?

For sin/1 the summand is an = (-1)n * x(2*n+1) / (2n+1)!.

For x = u * 10(-v), where u is the mantissa of x and v is the negated exponent of x, we see that:

    |an| < |x(2*n+1)| = 10((log10(|u|) - v) * (2*n+1))

So if we want p digits, we could try as smallest summand the index

    k = ceiling(p / (2 * (log10(|u|) - v)))

This works best if we already have |x| =< 1/2 as in the question. The smallest summand index k is computed here:

mp_sin(X, P, Y) :-
   K is integer(ceiling(requested(P)/(2*dec_log10(X)))),
   init_sin(K, X, P, U),
   mp_sin(U, X, P, Y).

The predicate init_sin/4 and mp_sin/4 then compute sine backwards from the smallest summand to the biggest summand a0, using Horner schema to get the sum as well at the same time:

mp_sin((0, S), _, _, S) :- !.
mp_sin(U, X, P, Y) :-
   next_sin(U, X, P, V),
   mp_sin(V, X, P, Y).

init_sin(K, X, P, (K, S)) :-
   (K mod 2 =:= 0 -> V = X; V = -X),
   mp_math(V/(2*K+1), P, S).

next_sin((L, T), X, P, (K, S)) :-
   K is L-1,
   (K mod 2 =:= 0 -> V = X; V = -X),
   mp_math((T*X*X/(2*K+2)+V)/(2*K+1), P, S).

The internal predicate mp_math/3 is here used for BigDecimal arithmetic, so that we can also compute sin/1 with a precision of p = 100, which is not possible with usual floats. Here are some results:

Jekejeke Prolog 4, Runtime Library 1.4.1 (20 August 2019)

?- use_module(library(decimal/multi)).
% 7 consults and 0 unloads in 222 ms.
Yes

?- X is mp(15*(pi/180),5), R is mp(sin(X),5).
X = 0d0.26181,
R = 0d0.25883

?- X is mp(15*(pi/180),102), R is mp(sin(X),102).
X = 0d0.2617993877991494365385536152
732919070164307832812588184145787160
256513671905174165523362354451764223
32,
R = 0d0.2588190451025207623488988376
240483283490689013199305138140032073
150569747488019969223679746942496655
21

Since all computations are done in the given precision, when you compute backwards with Horner schema, you need very few bits. Less numerical errors means you can drive the computation with less bits. But you need a little extra precision nevertheless.

To get a p = 4 result I was using p = 5, and to get a p = 100 result I was using p = 102. Did not yet have time to find a heuristic for this extra precision, and make it transparent for the end-user. Also often we could use a smaller k, so this is still work in progress.