一、问题描述
数值计算中，许多数学函数(如正弦、余弦、平方根等)可以通过泰勒级数展开的方法求得，然而，对于正弦函数与余弦函数，采用泰勒级数展开的方法收敛速度慢，计算速度慢。那么，有没有快速计算正弦函数与余弦函数并且可以估算误差限的算法呢？
QT中提供了快速正弦函数、快速余弦函数的算法，代码片段如下：
#define QT_SINE_TABLE_SIZE 256
#define M_PI (3.14159265358979323846)
inline qreal qFastSin(qreal x)
{
int si = int(x * (0.5 * QT_SINE_TABLE_SIZE / M_PI)); // Would be more accurate with qRound, but slower.
qreal d = x - si * (2.0 * M_PI / QT_SINE_TABLE_SIZE);
int ci = si + QT_SINE_TABLE_SIZE / 4;
si &= QT_SINE_TABLE_SIZE - 1;
ci &= QT_SINE_TABLE_SIZE - 1;
return qt_sine_table[si] + (qt_sine_table[ci] - 0.5 * qt_sine_table[si] * d) * d;
}
inline qreal qFastCos(qreal x)
{
int ci = int(x * (0.5 * QT_SINE_TABLE_SIZE / M_PI)); // Would be more accurate with qRound, but slower.
qreal d = x - ci * (2.0 * M_PI / QT_SINE_TABLE_SIZE);
int si = ci + QT_SINE_TABLE_SIZE / 4;
si &= QT_SINE_TABLE_SIZE - 1;
ci &= QT_SINE_TABLE_SIZE - 1;
return qt_sine_table[si] - (qt_sine_table[ci] + 0.5 * qt_sine_table[si] * d) * d;
}
二、推导过程
QT中提供的快速正弦、快速余弦的算法，是一种典型的以空间换时间的查表法。这里以快速正弦函数为例进行讲解，对于快速余弦函数很容易类比推导得到。
由于
s
i
n
(
x
)
sin(x)
sin(x)为周期为
2
π
2\pi
2π的函数，对于任意实数
x
x
x均可以变换到区间
[
0
,
2
π
]
[0,2\pi]
[0,2π]，假设区间
[
0
,
2
π
]
[0,2\pi]
[0,2π]的正弦表
s
i
n
e
T
a
b
l
e
sineTable
sineTable长度为
N
N
N。
求
s
i
n
(
x
)
sin(x)
sin(x)时，实数
x
x
x所对应的最接近其值的正弦表索引(不一定在
0
至
N
−
1
0 至 N-1
0至N−1之间):
s
i
=
r
o
u
n
d
(
x
N
2
π
)
(1)
s_i=round(\dfrac{xN}{2\pi})\tag 1
si=round(2πxN)(1)
由于调用
r
o
u
n
d
round
round取整函数增加了函数调用的开销，这里用其他方式实现
r
o
u
n
d
round
round函数的功能。若
x
>
0
x>0
x>0，则
s
i
g
n
=
1
sign=1
sign=1;否则，
s
i
g
n
=
−
1
sign=-1
sign=−1。计算
r
o
u
n
d
(
x
N
2
π
)
round(\dfrac{xN}{2\pi})
round(2πxN)等价于计算
(
i
n
t
)
(
x
N
2
π
+
s
i
g
n
∗
0.5
)
(int) (\dfrac{xN}{2\pi}+sign*0.5)
(int)(2πxN+sign∗0.5)。
计算
x
x
x与所对应的最接近其值的正弦表横坐标
x
s
i
x_{s_i}
xsi的偏差：
d
=
x
−
s
i
2
π
N
(2)
d=x-s_i \dfrac{2\pi}{N} \tag 2
d=x−siN2π(2)
由于
c
o
s
(
x
)
=
s
i
n
(
x
+
π
/
2
)
cos(x)=sin(x+\pi/2)
cos(x)=sin(x+π/2)，因而，求
c
o
s
(
x
)
cos(x)
cos(x)时，实数
x
x
x所对应的最接近其值的正弦表索引(不一定在
0
至
N
−
1
0 至 N-1
0至N−1之间):
c
i
=
s
i
+
N
/
4
(3)
c_i=s_i + N/4\tag 3
ci=si+N/4(3)
由于除法运算较慢，上式可以采用移位运算代替：
c
i
=
s
i
+
(
N
>
>
2
)
c_i=s_i +( N>>2)
ci=si+(N>>2)
将
s
i
,
c
i
s_i,c_i
si,ci变换到
0
至
N
−
1
0至N-1
0至N−1之间：
{
s
i
=
s
i
&
(
N
−
1
)
c
i
=
c
i
&
(
N
−
1
)
(4)
\begin{cases} s_i=s_i \& (N-1) \\ c_i=c_i \& (N-1) \\ \tag 4 \end{cases}
{si=si&(N−1)ci=ci&(N−1)(4)
其中，
&
\&
&表示按位与。
由于
s
i
n
(
x
)
=
s
i
n
(
x
s
i
+
d
)
=
s
i
n
(
x
s
i
)
c
o
s
(
d
)
+
c
o
s
(
x
s
i
)
s
i
n
(
d
)
sin(x)=sin(x_{s_i}+d)=sin(x_{s_i})cos(d)+cos(x_{s_i})sin(d)
sin(x)=sin(xsi+d)=sin(xsi)cos(d)+cos(xsi)sin(d)，
d
d
d很小，
s
i
n
(
d
)
,
c
o
s
(
d
)
sin(d),cos(d)
sin(d),cos(d)用等价无穷小代替：
{
s
i
n
(
d
)
∼
d
c
o
s
(
d
)
∼
1
−
d
2
/
2
(5)
\begin{cases} sin(d) \sim d \\ cos(d)\sim1-d^2/2 \\ \tag 5 \end{cases}
{sin(d)∼dcos(d)∼1−d2/2(5)
得到：
s
i
n
(
x
)
=
s
i
n
e
T
a
b
l
e
[
s
i
]
+
(
s
i
n
e
T
a
b
l
e
[
c
i
]
−
0.5
∗
s
i
n
e
T
a
b
l
e
[
s
i
]
∗
d
)
∗
d
(6)
sin(x)=sineTable[s_i] +(sineTable[c_i] - 0.5*sineTable[s_i]*d) *d \tag 6
sin(x)=sineTable[si]+(sineTable[ci]−0.5∗sineTable[si]∗d)∗d(6)
特别的，当需要同时计算
s
i
n
(
x
)
,
c
o
s
(
x
)
sin(x),cos(x)
sin(x),cos(x)时，为了进一步减小计算量，避免重复计算
d
,
s
s
,
c
i
d,s_s,c_i
d,ss,ci，可以将计算
s
i
n
(
x
)
,
c
o
s
(
x
)
sin(x),cos(x)
sin(x),cos(x)写成一个函数。
s
i
n
(
x
)
sin(x)
sin(x)在
0
0
0处的泰勒展开为：
s
i
n
(
x
)
=
x
−
x
3
/
6
+
x
5
/
120
−
x
7
/
5040
+
.
.
.
sin(x)=x-x^3/6 + x^5/120- x^7/5040+...
sin(x)=x−x3/6+x5/120−x7/5040+...。
c
o
s
(
x
)
cos(x)
cos(x)在
0
0
0处的泰勒展开为：
c
o
s
(
x
)
=
1
−
x
2
/
2
+
x
4
/
24
−
x
6
/
720
+
.
.
.
cos(x)=1- x^2/2 + x^4/24- x^6/720+...
cos(x)=1−x2/2+x4/24−x6/720+...。
采用式(5)的等价无穷小表示时，
s
i
n
(
d
)
sin(d)
sin(d)的误差限约为
ϵ
s
=
d
m
a
x
3
/
6
\epsilon_{s}=d_{max}^3/6
ϵs=dmax3/6,
c
o
s
(
d
)
cos(d)
cos(d)的误差限约为
ϵ
c
=
d
m
a
x
4
/
24
\epsilon_{c}=d_{max}^4/24
ϵc=dmax4/24,
d
m
a
x
=
π
/
(
N
−
1
)
d_{max}=\pi/(N-1)
dmax=π/(N−1)。
计算
s
i
n
(
x
)
sin(x)
sin(x)的误差限：
ϵ
s
i
n
e
≤
∣
s
i
n
(
x
s
i
)
∣
m
a
x
∗
ϵ
c
+
∣
c
o
s
(
x
s
i
)
∣
m
a
x
∗
ϵ
s
=
ϵ
c
+
ϵ
s
(7)
\epsilon_{sine}\le |sin(x_{s_i})|_{max} * \epsilon_{c} + |cos(x_{s_i})|_{max} * \epsilon_{s}=\epsilon_{c}+\epsilon_{s} \tag 7
ϵsine≤∣sin(xsi)∣max∗ϵc+∣cos(xsi)∣max∗ϵs=ϵc+ϵs(7)
同理可计算
c
o
s
(
x
)
cos(x)
cos(x)的误差限：
ϵ
c
o
s
i
n
e
≤
ϵ
c
+
ϵ
s
(8)
\epsilon_{cosine}\le \epsilon_{c}+\epsilon_{s} \tag 8
ϵcosine≤ϵc+ϵs(8)
三、c代码
快速正弦函数、快速余弦函数的c语言实现如下：
#define SINE_TABLE_SIZE 256
static const double sineTable[SINE_TABLE_SIZE] = {
0.0,
0.024541228522912288,
0.049067674327418015,
0.073564563599667426,
0.098017140329560604,
0.1224106751992162,
0.14673047445536175,
0.17096188876030122,
0.19509032201612825,
0.2191012401568698,
0.24298017990326387,
0.26671275747489837,
0.29028467725446233,
0.31368174039889152,
0.33688985339222005,
0.35989503653498811,
0.38268343236508978,
0.40524131400498986,
0.42755509343028208,
0.44961132965460654,
0.47139673682599764,
0.49289819222978404,
0.51410274419322166,
0.53499761988709715,
0.55557023301960218,
0.57580819141784534,
0.59569930449243336,
0.61523159058062682,
0.63439328416364549,
0.65317284295377676,
0.67155895484701833,
0.68954054473706683,
0.70710678118654746,
0.72424708295146689,
0.74095112535495911,
0.75720884650648446,
0.77301045336273699,
0.78834642762660623,
0.80320753148064483,
0.81758481315158371,
0.83146961230254524,
0.84485356524970701,
0.85772861000027212,
0.87008699110871135,
0.88192126434835494,
0.89322430119551532,
0.90398929312344334,
0.91420975570353069,
0.92387953251128674,
0.93299279883473885,
0.94154406518302081,
0.94952818059303667,
0.95694033573220894,
0.96377606579543984,
0.97003125319454397,
0.97570213003852857,
0.98078528040323043,
0.98527764238894122,
0.98917650996478101,
0.99247953459870997,
0.99518472667219682,
0.99729045667869021,
0.99879545620517241,
0.99969881869620425,
1.0,
0.99969881869620425,
0.99879545620517241,
0.99729045667869021,
0.99518472667219693,
0.99247953459870997,
0.98917650996478101,
0.98527764238894122,
0.98078528040323043,
0.97570213003852857,
0.97003125319454397,
0.96377606579543984,
0.95694033573220894,
0.94952818059303667,
0.94154406518302081,
0.93299279883473885,
0.92387953251128674,
0.91420975570353069,
0.90398929312344345,
0.89322430119551521,
0.88192126434835505,
0.87008699110871146,
0.85772861000027212,
0.84485356524970723,
0.83146961230254546,
0.81758481315158371,
0.80320753148064494,
0.78834642762660634,
0.7730104533627371,
0.75720884650648468,
0.74095112535495899,
0.72424708295146689,
0.70710678118654757,
0.68954054473706705,
0.67155895484701855,
0.65317284295377664,
0.63439328416364549,
0.61523159058062693,
0.59569930449243347,
0.57580819141784545,
0.55557023301960218,
0.53499761988709715,
0.51410274419322177,
0.49289819222978415,
0.47139673682599786,
0.44961132965460687,
0.42755509343028203,
0.40524131400498992,
0.38268343236508989,
0.35989503653498833,
0.33688985339222033,
0.31368174039889141,
0.29028467725446239,
0.26671275747489848,
0.24298017990326407,
0.21910124015687005,
0.19509032201612861,
0.17096188876030122,
0.1467304744553618,
0.12241067519921635,
0.098017140329560826,
0.073564563599667732,
0.049067674327417966,
0.024541228522912326,
0.0,
-0.02454122852291208,
-0.049067674327417724,
-0.073564563599667496,
-0.09801714032956059,
-0.1224106751992161,
-0.14673047445536158,
-0.17096188876030097,
-0.19509032201612836,
-0.2191012401568698,
-0.24298017990326382,
-0.26671275747489825,
-0.29028467725446211,
-0.31368174039889118,
-0.33688985339222011,
-0.35989503653498811,
-0.38268343236508967,
-0.40524131400498969,
-0.42755509343028181,
-0.44961132965460665,
-0.47139673682599764,
-0.49289819222978393,
-0.51410274419322155,
-0.53499761988709693,
-0.55557023301960196,
-0.57580819141784534,
-0.59569930449243325,
-0.61523159058062671,
-0.63439328416364527,
-0.65317284295377653,
-0.67155895484701844,
-0.68954054473706683,
-0.70710678118654746,
-0.72424708295146678,
-0.74095112535495888,
-0.75720884650648423,
-0.77301045336273666,
-0.78834642762660589,
-0.80320753148064505,
-0.81758481315158382,
-0.83146961230254524,
-0.84485356524970701,
-0.85772861000027201,
-0.87008699110871135,
-0.88192126434835494,
-0.89322430119551521,
-0.90398929312344312,
-0.91420975570353047,
-0.92387953251128652,
-0.93299279883473896,
-0.94154406518302081,
-0.94952818059303667,
-0.95694033573220882,
-0.96377606579543984,
-0.97003125319454397,
-0.97570213003852846,
-0.98078528040323032,
-0.98527764238894111,
-0.9891765099647809,
-0.99247953459871008,
-0.99518472667219693,
-0.99729045667869021,
-0.99879545620517241,
-0.99969881869620425,
-1.0,
-0.99969881869620425,
-0.99879545620517241,
-0.99729045667869021,
-0.99518472667219693,
-0.99247953459871008,
-0.9891765099647809,
-0.98527764238894122,
-0.98078528040323043,
-0.97570213003852857,
-0.97003125319454397,
-0.96377606579543995,
-0.95694033573220894,
-0.94952818059303679,
-0.94154406518302092,
-0.93299279883473907,
-0.92387953251128663,
-0.91420975570353058,
-0.90398929312344334,
-0.89322430119551532,
-0.88192126434835505,
-0.87008699110871146,
-0.85772861000027223,
-0.84485356524970723,
-0.83146961230254546,
-0.81758481315158404,
-0.80320753148064528,
-0.78834642762660612,
-0.77301045336273688,
-0.75720884650648457,
-0.74095112535495911,
-0.724247082951467,
-0.70710678118654768,
-0.68954054473706716,
-0.67155895484701866,
-0.65317284295377709,
-0.63439328416364593,
-0.61523159058062737,
-0.59569930449243325,
-0.57580819141784523,
-0.55557023301960218,
-0.53499761988709726,
-0.51410274419322188,
-0.49289819222978426,
-0.47139673682599792,
-0.44961132965460698,
-0.42755509343028253,
-0.40524131400499042,
-0.38268343236509039,
-0.359895036534988,
-0.33688985339222,
-0.31368174039889152,
-0.2902846772544625,
-0.26671275747489859,
-0.24298017990326418,
-0.21910124015687016,
-0.19509032201612872,
-0.17096188876030177,
-0.14673047445536239,
-0.12241067519921603,
-0.098017140329560506,
-0.073564563599667412,
-0.049067674327418091,
-0.024541228522912448
};
/**********************************************************************************************
Function: fast_sin
Description: 快速正弦函数
Input: 浮点数x
Output: 无
Input_Output: 无
Return: sin(x)
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
double fast_sin(double x)
{
int sign = x > 0.0 ? 1 : -1;
int si = (int)(x * SINE_TABLE_SIZE / (2.0 * PI) + sign * 0.5); //+0.5后取整等价于round函数
double d = x - si * (2.0 * PI / SINE_TABLE_SIZE);
int ci = si + (SINE_TABLE_SIZE >> 2);
si &= (SINE_TABLE_SIZE - 1);
ci &= (SINE_TABLE_SIZE - 1);
return sineTable[si] + (sineTable[ci] - 0.5 * sineTable[si] * d) * d;
}
/**********************************************************************************************
Function: fast_cos
Description: 快速余弦函数
Input: 浮点数x
Output: 无
Input_Output: 无
Return: cos(x)
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
double fast_cos(double x)
{
int sign = x > 0.0 ? 1 : -1;
int ci = (int)(x * SINE_TABLE_SIZE / (2.0 * PI) + sign * 0.5); //+0.5后取整等价于round函数
double d = x - ci * (2.0 * PI / SINE_TABLE_SIZE);
int si = ci + (SINE_TABLE_SIZE >> 2);
si &= (SINE_TABLE_SIZE - 1);
ci &= (SINE_TABLE_SIZE - 1);
return sineTable[si] - (sineTable[ci] + 0.5 * sineTable[si] * d) * d;
}
/**********************************************************************************************
Function: fast_sin_cos
Description: 快速正弦函数和快速余弦函数同时计算
Input: 浮点数x
Output: 正弦值sinX，余弦值cosX
Input_Output: 无
Return: 无
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
void fast_sin_cos(double x, double* sinX, double* cosX)
{
int sign = x > 0.0 ? 1 : -1;
int si = (int)(x * SINE_TABLE_SIZE / (2.0 * PI) + sign * 0.5); //+0.5后取整等价于round函数
double d = x - si * (2.0 * PI / SINE_TABLE_SIZE);
int ci = si + (SINE_TABLE_SIZE >> 2);
si &= (SINE_TABLE_SIZE - 1);
ci &= (SINE_TABLE_SIZE - 1);
*sinX = sineTable[si] + (sineTable[ci] - 0.5 * sineTable[si] * d) * d;
*cosX = sineTable[ci] - (sineTable[si] + 0.5 * sineTable[ci] * d) * d;
}
#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <time.h>
#include <math.h>
#define TEST_COUNT 10000000
int main()
{
int i, N[3];
long startTime, finishTime;
double deltaX, maxErrorSine, maxErrorCosine;
double dmax, epsilon, x, s1, c1, s2, c2;
startTime = clock();
for (i = 0; i < TEST_COUNT; i++)
{
s1 = sin(i);
c1 = cos(i);
}
finishTime = clock();
printf("优化前计算时间: %f s\n", (finishTime - startTime) / 1000.0);
startTime = clock();
for (i = 0; i < TEST_COUNT; i++)
{
fast_sin_cos(i, &s2, &c2);
}
finishTime = clock();
printf("优化后计算时间: %f s\n", (finishTime - startTime) / 1000.0);
maxErrorSine = 0.0;
maxErrorCosine = 0.0;
deltaX = 2.0 * PI / (TEST_COUNT - 1);
for (i = 0; i < TEST_COUNT; i++)
{
x = i * deltaX;
s1 = sin(x);
c1 = cos(x);
fast_sin_cos(x, &s2, &c2);
if (fabs(s1 - s2) > maxErrorSine)
{
maxErrorSine = fabs(s1 - s2);
}
if (fabs(c1 - c2) > maxErrorCosine)
{
maxErrorCosine = fabs(c1 - c2);
}
}
printf("快速正弦函数最大计算误差: %e\n", maxErrorSine);
printf("快速余弦函数最大计算误差: %e\n", maxErrorCosine);
N[0] = 256;
N[1] = 512;
N[2] = 1024;
for (i = 0; i < 3; i++)
{
dmax = PI / (N[i] - 1);
epsilon = dmax * dmax * dmax / 6.0 + dmax * dmax * dmax * dmax / 24.0;
printf("当N=%d时，快速正弦函数与快速余弦函数理论误差限: %e\n", N[i], epsilon);
}
getchar();
return 0;
}
3673
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
