c++和c#sin函数大值的不同结果
c#
当我使用大数时，我遇到了 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?¹? 真值的一部分。（我使用Sollya计算了这个近似值，中间精度为 165 位。）
但是，通过询问带有签名的 C# 函数public static double Sin (double a)或带有签名的 C++ 函数，您不会得到这个答案double sin(double)。为什么？6.2831853071795856?×?10?? 不是 IEEE 754 binary64 或“double”浮点数，因此您最多只能了解附近浮点数的 sin 是什么。在最近的浮点数，和你通常会通过打字获取6.2831853071795856E+45到一个程序，是6283185307179585571582855233194226059181031424，它不同于6.2831853071795856？×？10？由 28417144766805773940818968576 ? 2.84?×?10²?。
IEEE 754 binary64 浮点数 6283185307179585571582855233194226059181031424 是 6.2831853071795856?×?10?? 的一个很好的近似值 在相对误差中（它与真实值的差异小于 10?¹? 部分），但绝对误差 ~2.84?×?10²? 是远远超出期限2罪（和隔靴搔痒的整数倍）。 所以，你问一个双重功能得到的答案将没有相似的问题你的源代码似乎问：你写sin(6.2831853071795856E+45)，而不是罪（6283185307179585600000000000000000000000000000）充其量你会得到罪（6283185307179585571582855233194226059181031424），大约是0.8248163906169679（再次, 正负 10?¹? 真值的一部分）。
这不是浮点的错，这出错了。 浮点算术可以很好地保持相对误差很小——一个好的数学库可以很容易地使用 binary64 浮点数来计算你提出的问题的好答案。如果您的标尺没有超过 10 的数字，那么您对 sin 的输入误差也可能来自一个小的测量误差？渐变。错误可能来自某种近似错误，例如通过使用截断级数来评估给您 sin 输入的任何函数，无论您使用哪种算术来计算该输入。 问题是您要求在较小的相对误差对应于绝对误差的点处评估周期函数（sin） 远远超出函数的时期。
所以，如果你发现自己试图回答 6.2831853071795856?×?10?? 也就是说，您可能做错了事——天真地使用双浮点数学库例程不会帮助您回答问题。但是使这个问题更加复杂的是，您的 C# 和 C++ 实现都无法返回接近 sin(6283185307179585571582855233194226059181031424) 真实值的任何值：
-
Math.Sin的C# 文档宣传该域可能存在与机器相关的限制。
最有可能的是，您使用的是 Intel CPU，并且您的 C# 实现只是执行Intel x87
fsin指令，根据Intel 手册，该指令仅限于域 [?2?³, 2?³] 中的输入，而您的则超出了 2¹ ?²。正如您所观察到的，此域之外的输入是逐字返回的，即使它们对于正弦函数来说是完全无意义的值。将输入输入到肯定有效的范围的一种快速而肮脏的方法是编写：
Math.Sin(Math.IEEERemainder(6.2831853071795856E+45, 2*Math.PI))这样，您就不会误用 Math.Sin 库例程，因此答案至少应该在 [?1,1] 中作为正弦值。您可以安排在 C/C++ 中使用
sin(fmod(6.2831853071795856E+45, 2*M_PI)). 但是你可能会得到接近 0.35680453559729486 的结果，这也是错误的——见下文关于参数减少的内容。 -
sin但是，您正在使用的 C++ 实现已被破坏；在 C++ 标准中对域没有这样的限制，并且使用广泛可用的高质量软件来计算参数缩减模并在缩减域上计算 sin，没有理由搞砸它（即使这不是一个好要问的问题！）。我不知道仅仅通过观察输出是什么错误，但很可能是在参数减少步骤中：因为 sin( + 2) = sin() 和 sin(?) = ?sin()，如果你想要计算任意实数的 sin() 就足以计算 sin() 其中 = + 2 位于 [?,] 中，对于某个整数 。参数减少是给定的计算任务。
典型的基于 x87 的 sin 实现使用该
fldpi指令以 64 位精度以 binary80 格式加载近似值，然后使用该指令fprem1将该近似值的模减少到 。这种近似不是很好：在内部，英特尔架构近似于0x0.c90fdaa22168c234cp + 2 = 3.1415926535897932384585988507819109827323700301349163055419921875与精度66位，并且fldpi然后将其四舍五入到binary80浮点数0x0.c90fdaa22168c235p + 2 = 3.14159265358979323851280895940618620443274267017841339111328125只有64的位精确。相比之下，典型的数学库，例如古老的 fdlibm，通常使用超过 100 位精度的近似值来减少参数 modulo ，这就是为什么 fdlibm 导数能够非常准确地计算 sin(62831853071795855715828552331942260314284)1
然而，具有明显的x87计算
fldpi/fprem1/fsin提供有关？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) 的近似值。在这种情况下，由于 是一个整数（实际上对于任何数量级至少为 2?² 的 binary64 浮点数，它们都是整数），您将得到恰好 0。
不幸的是，.NET 还没有这些函数（截至 2021-02-04），C 或 C++ 标准数学库也不包含它们，尽管您可以轻松找到sinPi 和 cosPi 的示例代码。
当然，您的问题仍然存在这样一个问题，即使在非常输入时评估 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)
请阅读有关泰勒级数的内容（它们与三角函数有关）
还尝试使用GNUplot来“可视化”sin函数（对于合理的数字）。
数学sin上每个有限实数的the都在 -1 和 1 之间。所以 Math.Sin(6.2831853071795856E+45) = 6.2831853071795856E+45 是错误的。
还可以阅读像这样和那样的会议论文。
- @HimBromBeere Every calculation whose result varies randomly over the entire result range for arguments n and n+epsilon is senseless.
- @HimBromBeere The problem is that the "precision" of a number like `6.2831853071795856E+45` is very very small, much under the 6.28 that is the "precision" needed for `sin` (in 6.28 you make a full "round" of sin). The next number after `6.2831853071795856E+45` in double has a "distance" > 1E30
回答
作为旁注，.NET 的 sin/cos/tan 计算方式发生了变化。它发生在 2016 年 6 月 2 日，此提交是在 .NET Core 上进行的。正如作者在floatdouble.cpp中所写：
AMD64 Windows 上的 Sin、Cos 和 Tan 之前是通过调用 x87 浮点代码（fsin、fcos、fptan）在 vmamd64JitHelpers_Fast.asm 中实现的，因为 CRT 助手太慢了。情况不再如此，所有平台都使用 CRT 调用。
请注意，CRT 是 C 语言运行时。这解释了为什么较新版本的 .NET Core在同一平台上使用时会产生与 C++ 相同的结果。他们使用的 CRT 与其他 C++ 程序使用的相同。
感兴趣的人的“旧”代码的最后一个版本：1和2。
目前尚不清楚 .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操作码的“官方”答案，所以它是英特尔（大约1980年）关于6.2831853071795856E+45的罪过多少的“官方”答案，所以如果你使用英特尔，你必须相信，否则你就是叛徒！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", ExactSpelling = true, SetLastError = true)]
[return: MarshalAs(UnmanagedType.Bool)]
private static extern bool VirtualProtect(IntPtr lpAddress, IntPtr dwSize, uint flAllocationType, out uint lpflOldProtect);
[DllImport("kernel32.dll", ExactSpelling = true, SetLastError = true)]
[return: MarshalAs(UnmanagedType.Bool)]
private static extern bool VirtualFree(IntPtr lpAddress, IntPtr dwSize, 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, 0x44, 0x24, 0x08, // fld QWORD PTR [rsp+0x8]
0xD9, 0xFE, // fsin
0xDD, 0x5C, 0x24, 0x08, // fstp QWORD PTR [rsp+0x8]
0xF2, 0x0F, 0x10, 0x44, 0x24, 0x08, // movsd xmm0,QWORD PTR [rsp+0x8]
0xC3, // ret
} :
new byte[]
{
0xDD, 0x44, 0x24, 0x04, // fld QWORD PTR [esp+0x4]
0xD9, 0xFE, // fsin
0xC2, 0x08, 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, 0, buf, body.Length);
// Change the access of the allocated memory from R/W to Execute
uint oldProtection;
bool result = VirtualProtect(buf, (IntPtr)body.Length, 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 位）。我解决了将程序集操作码放在内存块中并将内存标记为可执行的问题。
Post scriptum：请注意，没有人再使用 x87 指令集，因为它很慢。