[数值计算] 三角函数的泰勒级数批处理简单实现:正弦,余弦

neorobin

cmd 不支持三角函数, 如果仅要求4位左右有效数字, 可以用三角函数的泰勒级数前几项来近似计算

正弦,余弦 定义如下:





cmd 同样不支持小数, 而 正弦, 余弦 的值域在 [0,1], 我们把 sin(x) 的自变量和结果值都放大来实现.

令 Zn = 10000^n, 例如: Z1=10^4, Z2=10^8, Z3=10^12, Z5=10^20,...

把 sin(x) 级数展开式中 的自变量和结果值都放大 Z1=10000 倍,
10000*sin(x) = x*Z1 - (x*Z1)^3/(3!*Z2) + (x*Z1)^5/(5!*Z4) - (x*Z1)^7/(7!*Z6) + ...


在 32 位系统中, 5!*Z4 已经超出整数范围 [-2^31, 2^31-1], 而且 x 的值仅达 0.13 时就可以让
(x*Z1)^3 出现上溢, 对展开式第 3 项变形, 得

(x*Z1)^5/(5!*Z4) = X/a*X/b*X/c*X/d*X/(5!*Z4/(a*b*c*d));    (其中, X = x*Z1)

这样以乘除交替进行, 从而在一定范围内避免或减少计算溢出

为避免在 |x| = PI/2 时产生上溢, 以及 |x| 接近 0 时减少下溢, 列出不等式组

PP^2/a < MAX 且 PP^3/a/b < MAX 且 PP^4/a/b/c < MAX 且 PP^5/a/b/c/d <MAX 且 deg/a > 1 且 deg^2/a/b > 1 且 0 < c < MAX 且 0 < d < MAX 且
5!*z4/a/b/c/d < MAX          (其中, PP=PI/2*Z1=15708, MAX=2^31-1, deg=PI/180*Z1)

解上不等式组, 在解集范围内为 5 个系数  a, b, c, d, 5!*z4/a/b/c/d  选取了一组值为

a = 1; b = 1875; c = 15625; d = 16000; 5!*z4/a/b/c/d=2560000

其他各项依此类得, 仅取展开式前 4 项, 可以得到如下

SIN=(X-X*X/1875*X/320000+X*X/1875*X/15625*X/16000*X/2560000-X*X/1875*X/15360*X/15625*X/15625*X/16000*X/44800000)

对余弦函数, 取展开式前 5 项, 同样方法可以得到

COS=(10000-X*X/20000+X*X/1875*X/15625*X/819200-X*X/1875*X/15360*X/15625*X/15625*X/10240000+X*X/1875*X/15360*X/15625*X/15625*X/16000*X/15625*X/229376000)

即使如此, 当角度值大于一定值时, 仍将发生上溢而产生错误结果.
根据 正弦, 余弦函数都具有最小周期 2PI, 容易把任意角转化为 [0,2PI) 内的角来计算, 更进一步, [0,2PI) 内的角还可以转化为 [-PI/2,PI/2]内的角来计算:


角 t 在[0,2pi)内,

0 <= t < pi/2          t -> t         in [-pi/2, 0)
不作变换

pi/2 <= t < 3pi/2   t -> pi-t      in (-pi/2, pi/2]
sin(t) = sin(pi-t),    cos(t) = -cos(pi-t)

3pi/2 <= t < 2pi    t -> -2pi+t  in [-pi/2, 0)
sin(t) = sin(-2pi+t),  cos(t) = cos(-2pi+t)


这样 对任意角, 通过转化, 角度模值控制在 [-pi/2, pi/2] 内, 就可以避免出现上溢出错.

下面是测试代码

@echo off 

set "SIN=(A-A*A/1875*A/320000+A*A/1875*A/15625*A/16000*A/2560000-A*A/1875*A/15360*A/15625*A/15625*A/16000*A/44800000)"
set "COS=(10000-A*A/20000+A*A/1875*A/15625*A/819200-A*A/1875*A/15360*A/15625*A/15625*A/10240000+A*A/1875*A/15360*A/15625*A/15625*A/16000*A/15625*A/229376000)"

REM 以 Pi/8 为步长在 [0, 2Pi) 区间类的 15 个测试角度
set "List=0,3927,7854,11781,15708,19635,23562,27489,31416,35343,39270,43197,47124,51051,54978,58905"
set /a "p=31416, p2=62832, pn2=-62832, p#2=15708, p3#2=47124, p3#2_=p3#2-1"
set R_Sin=
set R_Cos=

for %%a in (%List%) do (
        set /a "t=%%a"
        set /a "s1=(t-p#2^t-p3#2)>>31, s3=p3#2_-t>>31, t=(-t&s1)+(t&~s1)+(p&s1)+(pn2&s3), r1=%SIN:A=t%"
        set /a "t=%COS:A=t%, r2=(-t&s1)+(t&~s1)"        
        setlocal EnableDelayedExpansion
        set R_Sin=!R_Sin!!r1!,
        set R_Cos=!R_Cos!!r2!,
)
set R_
pause
复制代码

输出结果

R_Cos=10000,9238,7071,3826,0,-3826,-7071,-9238,-10000,-9238,-7071,-3826,0,3826,7071,9238,
R_Sin=0,3827,7071,9239,9999,9239,7071,3827,0,-3827,-7071,-9239,-9999,-9239,-7071,-3827,
复制代码

CrLf

你和 plp 都是高斯

523066680

这趋势，真的要搞一个批处理的2.5D引擎吗……

luckboy45

虽然不明白有什么用，但我还是不明觉厉，给你顶起来！~