我探索了如何使用几种不同的方法实现余弦。其中一个实现几乎比math.h快3倍，如果您可以使用4位小数精度。
您有没有想过，用您最喜欢的编程语言编写的数学库是如何实现三角函数的，比如余弦？这是一个很常见的函数，你可以在任何一个数学库中找到，所以它必须相当简单，对吧？哦，不。它肯定不是。
这一切都始于我的朋友和同事斯蒂芬·马兹博士，当时他正在研究一个操作系统的内核，我建议他在屏幕上画余弦函数。我经常使用余弦作为图形应用程序的余弦。问题：他的内核不能使用C标准库(再见，math.h！)。它的目标是RISC-V体系结构(与Intel fcos指令完全不同！)。
您必须认识到，我既不是数学家，也不是系统编程大师。事实上，在我大学学习数学和计算机科学的10年里，不知何故，我从来没有上过三角学课程。因此，我将带您通过我在C语言中研究和实现余弦的障碍。我的目标是：
我们将探索几种计算余弦的方法和一些优化。因为我不是C忍者，所以我会避免任何精心设计的诡计或微优化(但如果你知道一些，请告诉我！)。事实上，对于性能稍差的代码，我经常选择可读性更好的代码。我们将在进行过程中进行基准测试，以了解我们行动的影响。所有代码和基准都可以在GitHub存储库中找到。
在我所有的编程时间里，只有一种情况下我使用余弦：游戏，游戏，游戏。
这是我在几乎所有游戏中都使用的基本代码，用于向给定方向移动对象。无论是从自上而下的角度看玩家，还是从屏幕上飞来飞去的投射物。
当我开始这个超越的兔子洞时，我发现了近似余弦的泰勒级数方法。
这似乎是数学库中计算余弦的事实方法，至少在较高级别上是这样。给出更多的项，近似值就会更精确。
我脑海中浮现的一个想法是使用查找表。它是一个预先计算的值数组，可用于在给定一些输入的情况下查找最接近的余弦值。在计算能力更加有限的年代，这种情况并不少见。我在GitHub上找不到任何使用表作为trig函数的著名项目，但我确信它们确实存在。
CORDIC是另一个在我的搜索中不断出现的词。它是一种迭代方法，仅使用加法、减法、位移位和小查找表即可计算正弦和余弦等函数。它通常是在硬件中实现的，可以追溯到20世纪50年代末，或者在经常在低端CPU或微控制器上运行的软件中实现，比如计算器中的那些。这种方法在过去相当流行，并且被Intel8087数学协处理器、TI-85计算器和HP-48G计算器使用。然而，我找不到任何关于今天是否经常使用它的参考资料。有关更多详细信息，请参阅维基百科的文章或描述该方法的原始论文，或者查看用C编写的实现。我不会将我的方法与之进行比较，但我有点好奇它是如何站得住脚的。原文中的数字：
然后是Intel CPU fcos指令。该指令在给定FPU寄存器中存储的-2^63到2^63范围内以弧度表示的浮点值的情况下计算余弦，并将其替换为结果。我不确定现代英特尔CPU是否仍然使用CORDIC方法，如果是，它是在硬件中实现还是在软件中实现。在分解了几个使用余弦的程序之后，我找不到一个真正使用fcos指令的程序。虽然速度很快，但其他人已经记录了指令的不准确性，尤其是当参数接近pi/2的倍数时。有关更多详细信息，请参阅此非官方文档或英特尔官方参考。
现在，任何理智的人都会使用C'；的数学。出于日常目的，你可能应该这样做。但请记住，我们不能使用标准库中的任何内容。而且，那也不会有任何乐趣。我确实比较了数学在我电脑上的余弦函数和Wolfram Alpha的精确度。我发现Math.h的精确度高达15位数，这远远超过了我所需要的。为了了解标准库如何计算余弦，我查看了几个C标准库实现的源代码。Glibc、Newlibc、MUSL等。尽管它们看起来像是在使用泰勒级数，但这些对我来说有点太难理解了。它们彼此都非常不同，通常依赖于几个密集的函数，到处都是幻数，有预计算值的表，并且使用大量的位欺骗。有人花了很多时间把这些弄快。
这是我在MUSL中尝试浏览相关代码时的屏幕截图。从cos()到__cos()到__rem_pio2()。
随着我逐步实现不同的余弦计算方法，我将从两个角度对它们进行比较：运行时和准确性。对于运行时，每个函数使用一定范围的输入值执行1亿次，并使用time.h的时钟函数对其进行计时。为了准确起见，它取我的函数的结果与Math.h；的结果在一定输入范围内的差值，然后返回最坏的情况。例如，精确值0.0002意味着在最坏的情况下，对于大范围输入的一个输入，我的实现与math.h相差0.0002。
//以秒为单位测量多次执行的时间。数值越小越好。Double Runtime(DOUBLE(*FUNC)(DOUBLE)){CLOCK_T START=CLOCK()；FOR(int i=0；i<；100000000；i++){易失性DOUBLE C=FUNC(I/10000.0)；(Void)c；}RETURN(CLOCK()-START)/(DOUBLE)CLOCKS_Per_Sec；}//找出与math.h相比精度最差的情况。数值越小越好。双精度(Double(*func)(Double)){Double w=-1；Double Start=0；Double Stop=const_2pi；Double Step=0.0000001；For(Double i=Start；I<；Stop；i+=Step){Double c=absd(func(I)-cos(I))；if(c&>；w){w=c；}}返回w；}
该基准测试是用clang 11.0.3编译的，运行在13英寸的2018 MacBook Pro上，配备2.7 GHz i7 CPU和16 GB RAM。
您可以在GitHub资源库中找到所有基准代码。感谢Marz博士重写了它，使其具有易于使用的界面。
当我开始使用输入值(如0.1%和0.235)测试它时，这似乎出人意料地准确。当我把它画在math.h旁边时，我的热情就消退了。
洋红色的线条是math.h，绿色是我的函数。它看起来相当精确，介于-pi和+pi之间，但是随后就爆炸了。
DOUBLE COS_TALYLE_TEXAL_6TERMS_NAIVE(DOUBLE X){Return 1-((x*x)/(2))+((x*x)/(24))-((x*x)/(720))+((x*x)/(40320))-((x*x)/(。3628800))+((x*x)/(479001600))；}。
在这一点上，我突然想到了3个可能的改进：缩小输入范围，减少冗余计算的数量，以及不断添加更多的术语。
我尝试的第一个优化是缩小范围。输入值越大，该方法的精度就越低。因为余弦每2pi重复一次，所以我们只想做x=x%(2*pi)；。然而，在C语言中，模运算符不能处理浮点数，所以我们自己做了一个。
#定义modd(x，y)((X)-(Int)((X)/(Y))*(Y))DOUBLE COS_泰勒_TEXAL_6TERMS_2pi(DOUBLE X){x=modd(x，const_2pi)；返回1-((x*x)/(2))+((x*x)/(24))-((x*x)/(720))+((x*x)/(40320))-((x*x)/(3628800))+((。X*x)/(479001600))；}
这对于pi以上的值更好，但是直到2pi它仍然非常不准确。我们可以做得更好，因为余弦值等于π的每个倍数，除了符号反转。要做到这一点，我们可以以2pi取模，如果值大于pi，则减去pi，计算泰勒级数，然后应用正确的符号。所以我们实际上只计算了从0到π的余弦。
DOUBLE COS_TELLER_TEXAL_6TERMS_pi(DOUBLE X){x=modd(x，const_2pi)；char sign=1；if(x>；const_Pi){x-=const_Pi；sign=-1；}回车符号*(1-((x*x)/(2))+((x*x)/(24))-((x*x)/(720))+((x*x)/(40320))-((x*x)/(。3628800))+((x*x)/(479001600))；}。
下一次优化涉及到删除一些冗余计算。您会注意到代码中到处都是x*x。我所做的只是减少了一些具有双xx=x*x；的乘法。
DOUBLE COS_泰勒_TEXAL_6TERMS(DOUBLE X){x=modd(x，const_2pi)；char sign=1；if(x>；const_Pi){x-=const_Pi；sign=-1；}DOUBLE xx=x*x；回车符号*(1-((Xx)/(2))+((xx*xx)/(24))-((xx*xx*xx)/(720))+((xx*xx)/(40320))-((xx*xx)/(3628800))+((xx*xx)/(479001600)；}。
这是一场巨大的表演胜利！我也试了双倍，但没有看到太大的不同，所以我继续前进。
我仍然不确定要用多少个术语。所以我尝试了10个术语，看看它是如何提高准确性的：
DOUBLE COS_泰勒_TEXAL_10TERMS(DOUBLE X){x=modd(x，const_2pi)；char sign=1；if(x>；const_Pi){x-=const_Pi；sign=-1；}DOUBLE xx=x*x；回车符号*(1-((Xx)/(2))+((xx*xx)/(24))-((xx*xx*xx)/(720))+((xx*xx)/(40320))-((xx*xx)/(3628800))+((xx*xx)/(479001600))-((xx*。Xx*xx)/(87178291200))+((xx*xx)/(20922789888000))-((xx*xx)/(6402373705728000))+((xx*xx)/(2432902008176640000))；}。
此时，在查看图表时，10个术语线与math.h重叠。进步了！当将最坏情况的精确度与math.h进行比较时，这4项是个笑话。6个术语的偏差是0.0001，这比我需要的更准确，而10个术语的偏差只有0.00000000007。喔！
然而，更多的术语是以高昂的运行时成本为代价的。从基准来看，天真的4个术语只用了0.4秒，6个术语用了0.94秒，10个术语用了1.46秒。同时，math.h只需要大约1.04秒，精度更高。
在向马兹博士展示了我的进步后，他做了一些代数魔术，并把他的改良版本送给了我。他的方法通过存储消除了大量冗余的计算，并且还有一个额外的好处，那就是允许您指定所需的词条数量。对于特定的应用程序，这可能很方便，这样您就可以将不同程度的精度/速度作为参数。
Double cos_Taylor_running_yterm(Double x，int y){int div=(Int)(x/const_Pi)；x=x-(div*const_Pi)；char sign=1；if(div%2！=0)sign=-1；Double Result=1.0；Double Intern=1.0；Double Num=x*x；for(int i=1；i<；=y；i++){Double Comp=2.0*i；Double den=comp。IF(i%2==0)RESULT+=INTER；ELSE RESULT-=INTER；}RETURN SIGN*RESULT；}。
为了进行基准测试，我没有使用带有第二个参数的这个版本。相反，我复制了该函数并对循环进行了硬编码，使其适用于常量值(如6、10和16)。
增加更多的期限肯定会带来递减的回报。在16个条件下，最坏情况的准确率降低了0.0000000000000009，但是运行时基准测试需要2.57秒。它一点也不慢，但和数学比起来…。它是。
我想尝试的另一个选项是查找表。其想法是预先计算一组值，并将它们硬编码到一个数组中。在计算机出现之前，查找表就已经存在很久了，所以这并不是一种新技术。在这种情况下，我希望放弃一点内存将带来巨大的运行时优势，同时仍然足够准确。
为了生成查找表，Marz博士编写了一个Python脚本，该脚本生成一个C头文件，其中包含一个数组，其中每个元素都是使用math.h计算得出的余弦值。非常聪明！
从数学导入cos，pidef main(f，精度，名称)：f.write("；Double%s[]={\n"；%name)j=0 p=0.0而True：f.write("；{：.20f}，"；.format(cos(P)j+=1 p+=精度，如果p>；2*pi：中断f.write("；1.0}；\n"；)f。Const int%s_size=%d；\n"；%(名称，j+1)if__name__=='；__main__'；：main(open("；cotable_1.h"；，"；w"；)，1.0，"；cotable_1"；)main(open("；cotable_0_1.h"；，"；w"；)main("；cotable_0_1.h"；，"；w"；)main(打开("；cotable_0_1.h"；，"；w"；)。)，0.001，"；COSTRATE_0_01"；)Main(OPEN("；COSTABLE_0_001.h"；，"；W"；)Main(OPEN("；COSTABLE_0_001.h"；，"；w"；)，0.10，"；COSTRATE_0_001"；)Main(OPEN("；COSTABLE_0_001.h"；，"；W"；)Main(OPEN。COSTABLE_0_0001.h"；，"；w"；)，0.0001，"；COSTABLE_0_0001"；)
我们想用不同的精度测试我们的桌子。我们生成了包含8个值、64个值、630值、6285值和62833个值的表。成本以增加可执行文件的形式出现。1.0和0.1表并不明显，但是其他表分别增加了大约5KB、50KB和500KB的可执行文件大小。
DOUBLE ABSD(DOUBLE A){*(UNSIGNED LONG*)&；a)&；=~(1UL<；<；63)；return a；}DOUBLE COS_TABLE_0_01(DOUBLE X){x=ABSD(X)；x=modd(x，const_2pi)；return Cotable_0_01[(Int)(x*100+0.5)]；}。
这些桌子似乎在精确度上达到了很好的平衡。最小表的最坏情况精度为0.49，因此不可用。但是，表大小每增加一位，精度就会增加1位：0.049、0.0049、0.00049和0.000049。每个表的运行时测试大约为0.38。真快啊！(实际上，我们将运行时降到了大约0.33，但是代码很难看。)。
当然，查找表很棒，但是对于表中没有的值，我们可以做得更好。引言，线性插值。对于取两个值之间的加权平均值来说，这只是一个听起来很酷的术语。现在，当输入值不在表中时，我们将根据哪个表项更接近来计算近似值。代码：
#定义LERP(w，v1，v2)((1-(W))*(V1)+(W)*(V2))DOUBLE COS_TABLE_0_01_LERP(DOUBLE X){x=absd(X)；x=modd(x，const_2pi)；Double i=x*100.0；int index=(Int)i；返回LERP(i-index，/*权重*/CoStability_0_01[索引]，/*下限*/CoStability_0_01[索引+1]/*上限*/)；}。
以下是我们的函数与math.h(越低越好！)在最坏情况下的精确度的比较：
下面是用于计算100,000,000个值的运行时函数的比较(越低越好！)：
那么我推荐使用什么呢？Math.h，如果可能的话。这些函数都不是特别慢，而且大多数都足够精确。但是下次我制作一个严重依赖于触发器函数的游戏时，我会使用0_001表。