只是一个附录
我不得不考虑程序是否确实将小浮点数加在一起成连续较大的浮点数,而不是将小浮点数加到较大的浮点数上(这可能会使结果比它需要的更不精确),但确实如此。
尽管在泰勒级数的每个元素处完全重新计算阶乘并且不使用直接-1 * (k mod 2)
获取(-1)^k
而是通过递归来获得是不优雅的。
这是方向的调用图:
附录 2:更高效计算的代码
因此,我抽出一些时间来编写一个cos
近似值,该近似值仅在自身上递归,并带有用于计算项和总和的所有辅助信息。
% ===
% 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).
运行这个!该Verbose
变量设置为verbose
在泰勒级数计算期间生成更多打印输出。我们计算了该系列的 11 个术语(索引 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.
Stackoverflow的80 列思想让我有点紧张。如今,我们在屏幕上有数以亿计的像素宽度,由于“Muh Visual Design”,它们未被使用并留白!!反正...
现在添加一些代码来生成Count
均匀分布在From
和之间的测试浮点数To
。在generator/4
回溯时生成连续值。cos_compare/3
比较我们的近似cos
函数计算的内容和系统计算的内容(位于库的深处)。
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]).
然后让我们实际比较 float-4.0
和 float之间的 100 个值+4.0
,我们在每个值处计算泰勒级数的 11 个项(索引 0..11):
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.
看起来没那么糟糕。
附录 3:使用 SWI-Prolog“dicts”在谓词之间进行通信
我发现在编写 Perl 函数时,缩短基于位置的参数传递并传递一组名称-值对(也称为“哈希”)通常是有利的。这增加了很多灵活性(命名参数,易于添加参数,易于调试,易于将参数传递给子函数等)
让我们在这里也试试这个。
这仅限于 SWI-Prolog,因为“dicts”是SWI-Prolog 的一项功能。像这样的代码使 Prolog 的索引机制变得无用,因为现在每个谓词都有完全相同的参数Dict
,所以在运行时应该相对较慢。
只是修改后的谓词是
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).
它更具可读性吗?我想是这样。