问题描述
当我使用大数时，我在 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;
解决方法
这两个答案都非常错误——但您提出的问题可能会给您带来麻烦。
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⁴⁵通过28417144766805773940818968576≈2.84×10²⁸
IEEE 754 binary64 浮点数 62831853071795855715828552331942226059181031424 是 6.28318530717958561 的一个很好的近似值，但误差小于 rel ⁴1⁴1⁴1⁴⁴1⁴⁴绝对误差 ~2.84 × 10²⁸远远超出sin 的周期 2?（远不及 ? 的整数倍）。 所以你通过询问双函数得到的答案将不存在相似的问题你的源代码似乎问：你写sin(6.2831853071795856E+45)，而不是罪（6283185307179585600000000000000000000000000000）充其量你会得到罪（6283185307179585571582855233194226059181031424），大约是0.8248163906169679（再次，正负10 -16厘米的部位真实值）。
这不是浮点数的错。 浮点数算法可以很好地保持相对误差很小——一个好的数学库可以很容易地使用 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?? 位于 [−?,?] 中，对于某个整数 ?。参数减少是给定?的计算?的任务。
sin 的典型基于 x87 的实现使用
fldpi指令以 64 位精度加载二进制 80 格式的 ? 的近似值，然后使用fprem1将该近似值的模减少到 ?。这种近似不是很好：在内部，英特尔架构近似于π通过0x0.c90fdaa22168c234cp + 2 = 3.1415926535897932384585988507819109827323700301349163055419921875与精度66位，和fldpi，则其四舍五入为binary80浮点数0x0.c90fdaa22168c235p + 2 = 3.14159265358979323851280895940618620443274267017841339111328125，只有64位精度。相比之下，典型的数学库，例如古老的 fdlibm，通常使用精度超过 100 位的近似值来减少参数取模 ?，这就是为什么 fdlibm 导数能够非常准确地计算 sin(62831855715828552361942423619428 /p>
然而，使用
fldpi/fprem1/fsin进行明显的 x87 计算给出了大约 -0.8053589558881794，而将 x87 单元设置为 binary64 算术（53 位精度）时的结果相同binary80 算术（64 位精度），或者首先使用 binary64 逼近 ?，给出大约 0.35680453559729486。很明显，您的 C++ 数学库正在做其他事情来为一个糟糕的问题提供错误的答案！
从你输入的数字来看，我猜你可能一直想看看当你试图评估 2? 的大倍数的罪时会发生什么：6.2831853071795856 × 10⁴⁵与 2?10×⁵有一个小的相对误差。当然，这样的 sin 总是为零，但也许你更普遍地想要计算函数 sin(2??)，其中 ? 是一个浮点数。
浮点运算标准 IEEE 754-2019 推荐（但不强制）操作 sinPi、cosPi、tanPi 等，其中 sinPi(?) = sin(?⋅?)。如果您的数学库支持它们，您可以使用 sinPi(2*t) 来获得 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 后的余数很容易做到。
,大数的sin，如6.2831853071795856E+45；没有意义。请记住，Pi 大约为 3.14，您的浮动数字可能是 IEEE 754（有关更多信息，请参阅 http://floating-point-gui.de/）
计算出的数字可能完全错误。
考虑在您的 C 或 C++ 代码中使用静态分析工具，例如 Fluctuat。或者像 CADNA
这样的动态工具根据经验，sin、cos、tan 等三角函数不应与大数（例如，超过几千乘以 PI ℼ = 3.14 的绝对值）。
否则，请使用像 GMPlib 这样的 bignum 库。它会使计算速度降低数千倍，但它可以为您提供对 sin (pi * 1.e+45)
请阅读有关 Taylor series 的内容（它们与三角函数有关）
还尝试使用 GNUplot 来“可视化”sin 函数（对于合理 数字）。
数学上每个有限实数的 sin 都在 -1 和 1 之间。所以 Math.Sin(6.2831853071795856E+45) = 6.2831853071795856E+45 是错误的。
还阅读会议论文，如 this one 和 that one。
,作为旁注，.NET 的 sin/cos/tan 计算方式发生了变化。它发生在 2016 年 6 月 2 日，this commit 在 .NET Core 上。正如作者在floatdouble.cpp中所写：
AMD64 Windows 上的 Sin、Cos 和 Tan 之前在 vm\amd64\JitHelpers_Fast.asm 中实现 通过调用 x87 浮点代码（fsin、fcos、fptan），因为 CRT 助手太慢了。这 不再是这种情况，所有平台都使用 CRT 调用。
请注意，CRT 是 C 语言运行时。这解释了为什么较新版本的 .NET Core 在在同一平台上使用时会产生与 C++ 相同的结果。他们使用的 CRT 与其他 C++ 程序使用的相同。
尚不清楚 .NET Framework 是否/何时继承了这些更改。从一些测试看来，
.NET Core >= 2.0（未测试之前的版本）：Math.Sin(6.2831853071795856E+45) == 0.824816390616968
.NET Framework 4.8（32 位和 64 位）：Math.Sin(6.2831853071795856E+45) == 6.28318530717959E+45
所以它没有继承它们。
有趣的附录
做了一些检查，Math.Sin(6.2831853071795856E+45) == 6.28318530717959E+45是x87程序集的fsin操作码的“官方”答案，所以它是Intel（大约1980年）的“官方”答案6.2831853071795856E+45的罪有多大，所以如果你用的是Intel，你必须相信Math.Sin(6.2831853071795856E+45) == 6.28318530717959E+45，否则你就是叛徒！ ???
如果您想再次检查：
public static class TrigAsm
{
[DllImport("kernel32.dll",ExactSpelling = true,SetLastError = true)]
private static extern IntPtr VirtualAlloc(IntPtr lpAddress,IntPtr dwSize,uint flAllocationType,uint flProtect);
[DllImport("kernel32.dll",SetLastError = true)]
[return: MarshalAs(UnmanagedType.Bool)]
private static extern bool VirtualProtect(IntPtr lpAddress,out uint lpflOldProtect);
[DllImport("kernel32.dll",SetLastError = true)]
[return: MarshalAs(UnmanagedType.Bool)]
private static extern bool VirtualFree(IntPtr lpAddress,uint dwFreeType);
private const uint PAGE_READWRITE = 0x04;
private const uint PAGE_EXECUTE = 0x10;
private const uint MEM_COMMIT = 0x1000;
private const uint MEM_RELEASE = 0x8000;
[SuppressUnmanagedCodeSecurity]
[UnmanagedFunctionPointer(CallingConvention.StdCall)]
public delegate double Double2Double(double d);
public static readonly Double2Double Sin;
static TrigAsm()
{
// Opcoes generated with https://defuse.ca/online-x86-assembler.htm
byte[] body = Environment.Is64BitProcess ?
new byte[]
{
0xF2,0x0F,0x11,0x44,0x24,0x08,// movsd QWORD PTR [rsp+0x8],xmm0
0xDD,// fld QWORD PTR [rsp+0x8]
0xD9,0xFE,// fsin
0xDD,0x5C,// fstp QWORD PTR [rsp+0x8]
0xF2,0x10,// movsd xmm0,QWORD PTR [rsp+0x8]
0xC3,// ret
} :
new byte[]
{
0xDD,0x04,// fld QWORD PTR [esp+0x4]
0xD9,// fsin
0xC2,0x00,// ret 0x8
};
IntPtr buf = IntPtr.Zero;
try
{
// We VirtualAlloc body.Length bytes,with R/W access
// Note that from what I've read,MEM_RESERVE is useless
// if the first parameter is IntPtr.Zero
buf = VirtualAlloc(IntPtr.Zero,(IntPtr)body.Length,MEM_COMMIT,PAGE_READWRITE);
if (buf == IntPtr.Zero)
{
throw new Win32Exception();
}
// Copy our instructions in the buf
Marshal.Copy(body,buf,body.Length);
// Change the access of the allocated memory from R/W to Execute
uint oldProtection;
bool result = VirtualProtect(buf,PAGE_EXECUTE,out oldProtection);
if (!result)
{
throw new Win32Exception();
}
// Create a delegate to the "function"
Sin = (Double2Double)Marshal.GetDelegateForFunctionPointer(buf,typeof(Double2Double));
buf = IntPtr.Zero;
}
finally
{
// There was an error!
if (buf != IntPtr.Zero)
{
// Free the allocated memory
bool result = VirtualFree(buf,IntPtr.Zero,MEM_RELEASE);
if (!result)
{
throw new Win32Exception();
}
}
}
}
}
在 .NET 中你不能有内联汇编（Visual Studio 的内联汇编只有 32 位）。我解决了将汇编操作码放在内存块中并将内存标记为可执行的问题。
发布脚本：请注意，没有人再使用 x87 指令集，因为它很慢。
,f(x) = sin(x) 在正负 1 之间振荡，周期为 2 * pi，即 sin(x) = sin(N * 2 * pi + x)。
暂时忘记编程语言并忘记 sin() 并考虑科学记数法 1.23E+45 中的值是什么意思。它可以表示从 1225000...(41 个零) 到 1234999...(42 个 9) 的任何值，即 1E43 个可能的整数（如果允许分数，则可以表示无限数量的值）。
我在这里通过使用十进制术语并将有效数字精确等同于严格定义来快速和宽松地玩。在此处查找链接以了解更多信息：https://en.wikipedia.org/wiki/Precision。我也希望我的指数算术是正确的，但无论如何它的含义应该很清楚。
回到 sin()、C# 和 C++：由于纯数学函数输出重复每 2 * pi (~6.28) 和双精度浮点数（IEEE 754 64 位双精度 https://en.wikipedia.org/wiki/IEEE_754）在 C# 和 C++ 中都有只有 15-17 个有效十进制数字，sin() 函数的输出可能没有意义，因为输入值“缺失”（不精确）至少 30 个有效/相关十进制数字。。 >