3

我有一些代码应该通过泰勒级数定义找到近似值sin(15°)并将其与内置函数进行比较sin。

我有不同的结果。

  • X- 弧度值
  • R- 系列和的当前值
  • Eps- 准确度值
  • S- 符号
  • F- 阶乘值N
  • XN-X的幂的值N
  • N- 功率，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.

控制台中的结果：

?- 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.
4

2 回答 2

2

如果您从最小的加数开始，然后使用霍纳模式，您将获得更好的结果，更少的数字错误。但问题是，如何确定哪个索引k应该给你最小的 summand ？ak

因为sin/1总和是。an = (-1)n * x(2*n+1) / (2n+1)!

对于，其中是 的尾数并且是 的负指数，我们看到：x = u * 10(-v)uxvx

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

因此，如果我们想要p数字，我们可以尝试将索引作为最小和

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

如果我们已经有|x| =< 1/2问题，这效果最好。k此处计算最小的总和指数：

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).

谓词init_sin/4，mp_sin/4然后从最小 summand 到最大 summand 反向计算正弦，同时使用 Horner 模式得到总和：a0

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).

内部谓词mp_math/3在这里用于 BigDecimal 算术，因此我们也可以sin/1使用 的精度进行计算p = 100，这对于通常的浮点数是不可能的。以下是一些结果：

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

由于所有计算都以给定的精度完成，因此当您使用 Horner 模式进行反向计算时，您只需要很少的位。更少的数字错误意味着您可以用更少的位来驱动计算。但是你仍然需要一点额外的精度。

得到p = 4我正在使用的结果p = 5，并得到p = 100我正在使用的结果p = 102。还没有时间为这种额外的精度找到启发式方法，并使其对最终用户透明。我们通常也可以使用更小的k，所以这仍在进行中。

于 2019-08-27T22:44:27.280 回答
1

无法识别代码中正确的第 n 项...阶乘计算隐藏在哪里？请记住，早期优化是软件工程中所有弊端的根源:)

维基百科中列出的解决方案的直接翻译引导我找到这段代码（好吧，应该在谓词名称中引用 Maclaurin和Taylor ......）

:- 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.

给出合理的结果：

?- 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.
于 2019-08-24T08:28:29.047 回答