c# - c++ 和 c# sin 函数大值的不同结果

Question

当我使用大数时，我遇到了 C# 中 Math.Sin 函数的这种奇怪行为；例如:
C#:.Net 4.7.2:Math.Sin(6.2831853071795856E+45) = 6.2831853071795856E+45
C++:sin(6.2831853071795856E+45) = -0.089650623841643268
任何想法如何获得与 C++ 相同的结果？
C# 示例:

double value = 6.2831853071795856E+45;    
Console.WriteLine(Math.Sin(value));

C++ 示例:

double x = 6.2831853071795856E+45;
double result;
    
result = sin(x);
cout << "sin(x) = " << result << endl;

Answer 1

这两个答案都非常错误——但你提出的问题很可能会让你陷入困境。
sin(6.2831853071795856 × 10⁴⁵)的真值约为0.09683996046341126；此近似值与真实值的差异小于真实值的 10⁻¹⁶ 部分。 (我使用具有 165 位中间精度的 Sollya 计算了这个近似值。)
但是，通过询问带有签名 public static double Sin (double a) 的 C# 函数，您不会得到这个答案。 ，或带有签名 double sin(double) 的 C++ 函数.为什么？ 6.2831853071795856 × 10⁴⁵ 不是 IEEE 754 binary64 或“double”浮点数，因此您最多只能了解附近浮点数的 sin 是什么。最接近的浮点数，通常输入 6.2831853071795856E+45 会得到什么进入一个程序，是6283185307179585571582855233194226059181031424，它与6.2831853071795856 × 10⁴⁵1457612⁵1457612⁵1457612⁵1457612⁵1457612⁵1487612⁵1457612⁵1487612⁵1457612⁵1457612⁵1457612⁵1487612⁵1487612⁵9
该IEEE 754 binary64浮点数6283185307179585571582855233194226059181031424是一个很好的近似6.2831853071795856×10⁴⁵在相对误差(它相差小于10⁻¹⁷份真值的)，但绝对误差〜2.84×10²⁸远远超出了周期2π罪(远不及𝜋的整数倍)。 因此，您通过询问双重函数得到的答案与您的源代码似乎提出的问题没有相似之处:您在哪里编写 sin(6.2831853071795856E+45) ，而不是罪(6283185307179585600000000000000000000000000000)充其量你会得到罪(6283185307179585571582855233194226059181031424)，大约是0.8248163906169679(再次，加上或减去10个-16厘米部分的真正值(value))。
这不是浮点的错，这出错了。 浮点算术可以很好地保持相对误差很小——一个好的数学库可以很容易地使用 binary64 浮点数来计算你提出的问题的好答案。如果您的标尺没有超过 10⁴⁵ 的等级，那么您对 ​​sin 的输入误差也可能来自一个小的测量误差。错误可能来自某种近似错误，例如通过使用截断级数来评估给您 sin 输入的任何函数，无论您使用哪种算术来计算该输入。 问题是您要求在一个小的相对误差对应于远远超出函数周期的绝对误差的点处评估周期函数 (sin)。

因此，如果您发现自己试图回答 6.2831853071795856 × 10⁴⁵ 的罪过是什么的问题，那么您可能做错了什么——天真地使用双浮点数学库例程并不能帮助您回答问题。但是使这个问题更加复杂的是，您的 C# 和 C++ 实现都无法返回接近 sin(6283185307179585571582855233194226059181031424) 真实值的任何值:

C# documentation for Math.Sin宣传该域可能存在依赖于机器的限制。
最有可能的是，您使用的是 Intel CPU，并且您的 C# 实现只是执行 Intel x87 fsin指令，根据 Intel manual仅限于域 [−2⁶³, 2⁶³] 中的输入，而您的输入超出了 2¹⁵²。正如您所观察到的，该域之外的输入是逐字返回的，即使它们对于正弦函数来说是完全无意义的值。
将输入输入到肯定有效的范围的一种快速而肮脏的方法是编写:

Math.Sin(Math.IEEERemainder(6.2831853071795856E+45, 2*Math.PI))

这样，您就不会误用 Math.Sin 库例程，因此答案至少应该在 [−1,1] 中，因为正弦应该如此。您可以使用 sin(fmod(6.2831853071795856E+45, 2*M_PI)) 安排在 C/C++ 中获得相同或接近的结果。 .但是您可能会得到接近 0.35680453559729486 的结果，这也是错误的——参见下文关于参数减少的内容。

sin的C++实现但是，您正在使用的只是损坏了；在 C++ 标准中对域没有这样的限制，并且使用广泛可用的高质量软件来计算参数减少模好问题!)。
我不知道仅仅通过观察输出是什么错误，但很可能是在参数减少步骤中:因为 sin(𝑥 + 2𝜋) = sin(𝑥) 和 sin(−𝑥) = −sin(𝑥) ，如果你想计算任意实数的 sin(𝑥) 就足以计算 sin(𝑦) 其中 𝑦 = 𝑥 + 2𝜋𝑘 位于 [−𝜋,𝜋] 中，对于某个整数 𝑘。参数减少是给定 𝑥 计算 𝑦 的任务。
典型的基于 x87 的 sin 实现使用 fldpi指令以 64 位精度加载二进制 80 格式的近似值 𝜋，然后使用 fprem1减少模那个近似于𝜋。这种近似不是很好:在内部，Intel 架构近似 𝜋 0x0.c90fdaa22168c234cp+2 = 3.14159265358979323845859109827323700303965 14,16 16 16 16 5 16 5 16 16 16 5 16 16 5 9 7 位精度为 4,595然后将其四舍五入为 binary80 浮点数 0x0.c90fdaa22168c235p+2 = 3.14159265358979323851280895940618620443274267017841338912513 位，精度为 16 位。
相比之下，典型的数学库，例如古老的 fdlibm，通常使用超过 100 位精度的近似值来减少参数取模 𝜋，这就是为什么 fdlibm 导数能够非常准确地计算 sin(62831853071795855715828552331942826035
然而，明显的 x87 计算与 fldpi/fldpi/fprem1给出大约 -0.8053589558881794，与 x87 单元设置为 binary64 算术(53 位精度)而不是 binary80 算术(64 位精度)相同，或者首先使用 binary64 逼近 𝜋，给出大约 0.356804953556972因此，显然您的 C++ 数学库正在做其他事情来为一个糟糕的问题提供错误的答案!

从你输入的数字来看，我猜你可能一直想看看当你尝试评估 2𝜋 的大倍数时会发生什么:6.2831853071795856 × 10⁴⁵与 2𝜋 ×‴10 有一个小的相对误差。当然，这样的 sin 总是为零，但也许你更普遍地想要计算函数 sin(2𝜋𝑡)，其中 𝑡 是一个浮点数。
浮点运算标准 IEEE 754-2019 推荐(但不强制)操作 sinPi、cosPi、tanPi 等，其中 sinPi(𝑡) = sin(𝜋⋅𝑡)。如果您的数学库支持它们，您可以使用 fsin得到 sin(2𝜋𝑡) 的近似值。在这种情况下，由于 𝑡 是一个整数(实际上对于任何二进制 64 位浮点数，其大小至少为 2⁵²，它们都是整数)，你会得到恰好 0。
不幸的是， .NET does not yet have these functions (截至 2021-02-04)，C 或 C++ 标准数学库也不包含它们，尽管您可以轻松找到 example code for sinPi and cosPi floating around .
当然，您的问题仍然存在这样一个问题，即在非常输入的情况下评估 sinPi 函数也很麻烦，因为即使是很小的相对误差(比如 10⁻¹⁷)仍然意味着绝对误差远远超出函数的周期 2。但问题是不会被以超越数为模的糟糕参数缩减所放大:计算浮点数除以 2 后的余数很容易做到。