【问题标题】:Yet another fast trigonometry另一个快速三角函数
【发布时间】:2018-12-14 19:51:33
【问题描述】:
float sinx(float x)
{
    static const float a[] = {-.1666666664,.0083333315,-.0001984090,.0000027526,-.0000000239};
    float xsq = x*x;
    float temp = x*(1 + a[0]*xsq + a[1]*xsq*xsq + a[2]* xsq*xsq*xsq+a[3]*xsq*xsq*xsq*xsq+ a[4]*xsq*xsq*xsq*xsq*xsq);
    return temp;
}

这些常数是如何计算的?如何使用这种方法计算costan? 我可以扩展它以获得更高的精度吗?我想我需要添加更多常量?


上述“快速”正弦的误差与等次泰勒多项式的关系图。

【问题讨论】:

  • 请注意,当增加度数(常数数量)时,|x|>1 的泰勒近似实际上变得更糟。这可以通过使用 sin(x) 的周期为 2π 来解决
  • 好的,那么我对这些常量没问题。我刚刚绘制了错误图,精度很好。现在我想弄清楚tancos 是否可以用同样的方式计算。在我必须在 8 位慢速 MCU 上计算的公式中,我有 tancos。除法很慢,所以最好避免它。
  • @Pablo 您不需要除法,您可以通过预先计算项将其转换为乘法,就像在上面的代码中,除以 6 被乘以 0.16666 代替。 .
  • 如果你知道除数,这是真的。我说的是tan=sin/cos,你不知道除数。除非我错过了你的观点。
  • @Pablo tan 的展开式是 x+x^3/3+(2 x^5)/15+...,您不需要将其计算为 sin/cos .

标签: c floating-point precision trigonometry numerical-methods


【解决方案1】:

在写这篇文章的时候几乎所有的答案都提到了函数 sin 的泰勒展开,但是如果函数的作者是认真的,他不会使用泰勒系数。泰勒系数往往会产生一个多项式近似,它在接近零时比必要的要好,并且在远离零时越来越差。目标通常是在 -π/2…π/2 等范围内获得一致良好的近似值。对于多项式近似,这可以通过应用the Remez algorithm 来获得。一个接地气的解释是thispost。

通过该方法获得的多项式系数接近泰勒系数,因为两个多项式都试图逼近相同的函数,但多项式可能是 对于相同数量的操作更精确,或者在相同(统一)的近似质量下涉及更少的操作。

我无法仅通过查看它们来判断您问题中的系数是完全泰勒系数还是由 Remez 算法获得的略有不同的系数,但即使不是,它也可能是应该使用的。

最后,写(1 + a[0]*xsq + a[1]*xsq*xsq + a[2]* xsq*xsq*xsq+a[3]*xsq*xsq*xsq*xsq+ a[4]*xsq*xsq*xsq*xsq*xsq) 的人需要阅读更好的多项式评估方案,例如Horner's

1 + xsq*(a[0]+ xsq*(a[1] + xsq*(a[2] + xsq*(a[3] + xsq*a[4])))) 使用 N 次乘法而不是 N2/2。

【讨论】:

  • 它在 ENS Lyon 实施和改进(稍微微调系数,以便浮点错误倾向于改善而不是恶化结果):perso.ens-lyon.fr/florent.de.dinechin/recherche/publis/… 他们想要发布它,但他们担心它太难用了。 :)
  • 你能用霍纳规则对原始函数进行一些修正吗?
  • @Pablo 我已经用xsq 中 4 次多项式的完整霍纳评估更新了我的问题。这应该与您问题中的代码几乎相同,更有效。
  • @Pablo sinf() 的一个很好的单精度实现在这里。附加系数无济于事,因为精度受到结果的单精度格式的限制。 github.com/JuliaLang/openlibm/blob/master/src/k_sinf.c
  • 看起来很有趣,我将在我的 MCU 上进行基准测试。实际上我只需要 tan 和 cos,因为我使用 Rhumb 公式计算地理距离。我尝试使用 LUT,然后使用合适的精度表太大而无法放入闪存。原罪需要太多的循环。我一直在寻找合适的解决方案。最糟糕的是我的浮点数是 4 个字节。我可能需要将我找到的任何解决方案转换为定点。彻头彻尾的噩梦。
【解决方案2】:

他们是-1/61/120-1/5040 ..等等。

或者更确切地说:-1/3!、1/5!-1/7!1/9!...等

sin x 看泰勒系列in here

它的正下方有cos x:

对于cos x,从上图可以看出,常量为-1/2!, 1/4!, -1/6!, 1/8!...

tan x 略有不同:

所以要针对 cosx 进行调整:

float cosx(float x)
{
    static const float a[] = {-.5, .0416666667,-.0013888889,.0000248016,-.0000002756};
    float xsq = x*x;
    float temp = (1 + a[0]*xsq + a[1]*xsq*xsq + a[2]* xsq*xsq*xsq+a[3]*xsq*xsq*xsq*xsq+ a[4]*xsq*xsq*xsq*xsq*xsq);
    return temp;
}

【讨论】:

  • cosa[0]*xsq还是a[0]*x*x*x等?
  • @Pablo 完全符合给定的 cos,从该部分的 sin 版本中删除 x* 就足够了。 sin 是前面的 1 次方,所以这就是为什么要乘以.. x(x² + x⁴ ..) 与 x³ + x⁵ ...
【解决方案3】:

系数与数学函数手册中给出的系数相同。 Abramowitz 和 Stegan (1964),第 76 页,归功于 Carlson 和 Goldstein,Rational approximations of functions,Los Alamos Scientific Laboratory (1955)。

第一个可以在http://www.jonsson.eu/resources/hmf/pdfwrite_600dpi/hmf_600dpi_page_76.pdf找到。

第二个在http://www.osti.gov/bridge/servlets/purl/4374577-0deJO9/4374577.pdf。第 37 页给出:

关于您的第三个问题,“我可以扩展它以获得更高的精度吗?”,http://lol.zoy.org/wiki/doc/maths/remez 有一个可下载的 Remez 算法的 C++ 实现;它提供(未经我检查)sin 的 6 阶多项式的系数:

 error: 3.9e-14
 9.99999999999624e-1
-1.66666666660981e-1
 8.33333330841468e-3
-1.98412650240363e-4
 2.75568408741356e-6
-2.50266363478673e-8
 1.53659375573646e-10

当然,您需要从 float 更改为 double 才能实现任何改进。这也可以回答您关于costan 的第二个问题。

另外,我在 cmets 中看到最终需要一个定点答案。大约 26 年前,我在 8031 汇编器中实现了 32 位定点版本;我会尝试挖掘它,看看它是否有任何有用的东西。

更新:如果您坚持使用 32 位双精度数,那么我认为您将精度提高“一两位”的唯一方法就是忘记浮点数并使用固定点。令人惊讶的是,谷歌似乎没有发现任何东西。以下代码提供概念验证,在标准 Linux 机器上运行:

#include <stdio.h>
#include <math.h>
#include <stdint.h>

// multiply two 32-bit fixed-point fractions (no rounding)
#define MUL32(a, b) ((uint64_t)(a) * (b) >> 32)

// sin32:  Fixed-point sin calculation for first octant, coefficients from
//         Handbook for Computing Elementary Functions, by Lyusternik et al, p. 89.
// input:  0 to 0xFFFFFFFF, giving fraction of octant 0 to PI/8, relative to 2**32
// output: 0 to 0.7071, relative to 2**32

static uint32_t sin32(uint32_t x) { // x in 1st octant, = radians/PI*8*2**32
   uint32_t y, x2 = MUL32(x, x);    // x2 = x * x
   y = 0x000259EB;                  // a7 = 0.000 035 877 1
   y = 0x00A32D1E  - MUL32(x2, y);  // a5 = 0.002 489 871 8
   y = 0x14ABBA77  - MUL32(x2, y);  // a3 = 0.080 745 367 2
   y = 0xC90FDA73u - MUL32(x2, y);  // a1 = 0.785 398 152 4
   return MUL32(x, y);
}

int main(void) {
   int i;
   for (i = 0; i < 45; i += 2) { // 0 to 44 degrees
      const double two32 = 1LL << 32;
      const double radians = i * M_PI / 180;
      const uint32_t octant = i / 45. * two32; // fraction of 1st octant
      printf("%2d  %+.10f  %+.10f  %+.10f    %+.0f\n", i,  
         sin(radians) - sin32(octant) / two32,
         sin(radians) - sinf(radians),
         sin(radians) - (float)sin(radians),
         sin(radians) * two32 - sin32(octant));
   }
   return 0;
}   

系数来自 Handbook for Computing Elementary Functions,Lyusternik 等人,p。 89、here

我选择这个特殊函数的唯一原因是它比你原来的系列少了一个词。

结果是:

 0  +0.0000000000  +0.0000000000  +0.0000000000    +0
 2  +0.0000000007  +0.0000000003  +0.0000000012    +3
 4  +0.0000000010  +0.0000000005  +0.0000000031    +4
 6  +0.0000000012  -0.0000000029  -0.0000000011    +5
 8  +0.0000000014  +0.0000000011  -0.0000000044    +6
10  +0.0000000014  +0.0000000050  -0.0000000009    +6
12  +0.0000000011  -0.0000000057  +0.0000000057    +5
14  +0.0000000006  -0.0000000018  -0.0000000061    +3
16  -0.0000000000  +0.0000000021  -0.0000000026    -0
18  -0.0000000005  -0.0000000083  -0.0000000082    -2
20  -0.0000000009  +0.0000000095  -0.0000000107    -4
22  -0.0000000010  -0.0000000007  +0.0000000139    -4
24  -0.0000000009  -0.0000000106  +0.0000000010    -4
26  -0.0000000005  +0.0000000065  -0.0000000049    -2
28  -0.0000000001  -0.0000000032  -0.0000000110    -0
30  +0.0000000005  -0.0000000126  -0.0000000000    +2
32  +0.0000000010  +0.0000000037  -0.0000000025    +4
34  +0.0000000015  +0.0000000193  +0.0000000076    +7
36  +0.0000000013  -0.0000000141  +0.0000000083    +6
38  +0.0000000007  +0.0000000011  -0.0000000266    +3
40  -0.0000000005  +0.0000000156  -0.0000000256    -2
42  -0.0000000009  -0.0000000152  -0.0000000170    -4
44  -0.0000000005  -0.0000000011  -0.0000000282    -2

因此,我们看到这个定点计算比sinf()(float)sin() 精确大约十倍,并且正确到29 位。在 MUL32() 中使用舍入而不是截断只取得了边际改进。

【讨论】:

  • 既然你似乎有这本书,你能说一下作者选择的系数有什么属性吗?我认为 -.1666666664(而不是 -.1666666666 或 -.1666666667)可能反映了系数的 IEEE 754 浮点性质,但如果它们是在 1955 年选择的,肯定不是,所以这可能是统一的毕竟是多项式逼近。
  • 感谢您提供引人入胜的链接。
  • @PascalCuoq:看论文,它看起来确实像一个统一多项式近似,但 Carlson 和 Goldstein 将他们的方法描述为“本质上是一种经验方法”(第 12 页)。
  • 是的,“它 [...] 希望该方法最终可以建立在坚实的数学基础上”让我笑了。维基百科将 Remez 算法的日期定为 1934 年。有相似之处,但不清楚 Carlson 和 Goldstein 是否意识到这一点。
  • @JosephQuinsey 不幸的是 double 和 float 的大小相同 - 4 字节
【解决方案4】:

该函数正在使用其泰勒展开式计算sin 的值:

这些常数是各种 -1/3!, 1/5!等等(参见例如here 了解泰勒级数的其他函数)。

现在,如果您指定了序列的每个项,则 sin(x) 的泰勒展开对于每个 x 都是精确的,但是 AFAIK 有更快、更精确的方法确定软件中的三角函数。

此外,许多处理器提供直接在处理器中实现的此类功能(例如,在 x86 上为它们提供现成的操作码),因此通常无需为这些东西操心。

【讨论】:

  • 我试过 CORDIC 并得到更差的精度结果。使用没有三角函数和慢除法的 8 位 AVR。
  • 根据我的阅读,通常使用 CORDIC 作为第一步,然后使用一个简短的泰勒级数来润色结果。参见例如here.
  • 是的,我也读过。这是值得尝试的事情……对于 CORDIC+Taylor 来说,周期可能太多了。
【解决方案5】:
cos(x) = sqrt(1 - sin^2(x))
tan(x) = sin(x)/cos(x)

Sin(x) = x -x^3/3! + x^5/5! + (-1)^k*x^(2k+1)/(2k+1)! , k = 1, 2, ...

这是无限函数

【讨论】:

  • 我对常量更感兴趣。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2015-04-20
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多