sin.rs文件提供了sin函数的实现,它计算并返回一个浮点数(f64类型)的正弦值。这个函数首先处理了一些特殊情况,如极小的值、无穷大和NaN(非数字),然后使用rem_pio2函数将输入参数x归约到[-π/2, π/2]的范围内,并根据归约的结果和n的值(n是x除以π/2的商的整数部分,对4取模的结果),选择调用k_sin或k_cos函数来计算最终的正弦值。
一、代码中关键部分的解释
1. 特殊值处理
- 如果x的绝对值非常小(小于2^-26),函数会直接返回x,因为在这个范围内,sin(x)的值非常接近x。
- 如果x是无穷大或NaN,函数会返回NaN。
2. 参数归约
使用rem_pio2函数将x归约到[-π/2, π/2]范围内,并返回三个值:n(x除以π/2的商的整数部分),y0(x的归约值的主要部分),和y1(x的归约值的次要部分,用于提高精度)。
3. 选择正确的函数来计算正弦值
根据n的值(对4取模),选择调用k_sin或k_cos函数,并可能取反。这是因为正弦和余弦函数在[-π/2, π/2]范围内是对称的,并且可以通过旋转坐标系来相互转换。
4. 测试
提供了一个测试函数test_near_pi,它测试了当x接近π时的正弦值计算是否正确。注意,由于浮点数的表示误差,测试值可能不是精确的π,而是一个接近π的浮点数。
5. force_eval!宏:
这个宏在代码中用于确保某些表达式被评估,即使编译器可能认为它们对最终结果没有影响。如浮点数的精度限制、下溢等。这在处理极小值或强制触发浮点异常时特别有用。但是它是不安全的,会因为特定硬件产生错误,不建议使用。代码如下:
macro_rules! force_eval {
($e:expr) => {
unsafe { ::core::ptr::read_volatile(&$e) }
};
}
6. 位操作和浮点表示:
代码使用了位操作来提取浮点数的指数和尾数部分,以及通过特定的位模式来构造特定的浮点数(如x1p120,它表示2^120)。
请注意,k_sin和k_cos函数以及rem_pio2函数的具体实现在这段代码中没有给出。这些函数是计算正弦和余弦值的核心,rem_pio2函数负责将输入参数归约到基本区间内。
此外,#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]属性是一个条件编译属性,它在启用了测试和assert_no_panic特性时,确保sin函数不会引发恐慌(panic)。这对于编写健壮的库函数非常重要,因为它们不应该在正常情况下崩溃。
二、源代码
use super::{k_cos, k_sin, rem_pio2};
// sin(x)
// 返回x的正弦函数。
//
// 核心函数:
// k_sin ... [-pi/4,pi/4]上的正弦函数
// k_cos ... [-pi/4,pi/4]上的余弦函数
// rem_pio2 ... argument reduction routine
//
// 方法。
//设S、C分别表示[-PI/4、+PI/4]上的sin、cos。在[-pi/4,+pi/4]中,将参数x减小到y1+y2=x-k*pi/2,并设n=k mod 4。
//我们有
//
// n sin(x) cos(x)
// --------------------------------
// 0 S C
// 1 C -S
// 2 -S -C
// 3 -C S
// --------------------------------
//
//特殊情况:
//设trig为“sin”或“cos”中的任何一个。
//trig(+-INF)是NaN,带有信号;
//trig(NaN)是NaN;
//准确度:
//TRIG(x)返回几乎四舍五入的TRIG(x)
/// x的正弦(f64)。(sine)
///
/// x单位是弧度。
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
pub fn sin(x: f64) -> f64 {
let x1p120 = f64::from_bits(0x4770000000000000); // x1p120 == 2^120
/* High word of x. */
let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;//去除符号位,即提取指数11位和尾数前20位
/* |x| ~< pi/4 */
if ix <= 0x3fe921fb {//等同于不超过后32位都是1的值(即0.7853984832763671),pi/4的值为0.785398163397448309615660845819875721_f64,会有略大于pi/4的数放进来。
if ix < 0x3e500000 {// 2^-26,即1.4901161193847656e-8,考虑后32位可能都是1,即1.4901175404702368e-8
/* |x| < 2**-26 */
/* 如果x!=0,则提高不精确性0,如果低于正常值,则下溢*/
if ix < 0x00100000 {//即非约数或0,
force_eval!(x / x1p120);//制造一个下溢
} else {
force_eval!(x + x1p120);//小数加大数,结果小数被舍去,是否有报错信息就不清楚了
}
return x;//以上程序说明,|x|足够小时,即小于1.4901175404702368e-8时返回 x
}
return k_sin(x, 0.0, 0);//-45°至45°排除小值后调用该函数。
}
/* sin(无穷大 或 非数时返回非数 */
if ix >= 0x7ff00000 {
return x - x;//返回NaN值
}
/* 将输入参数归约到一个周期内的值来简化计算 */
let (n, y0, y1) = rem_pio2(x);//调用归约函数
match n & 3 {
0 => k_sin(y0, y1, 1),
1 => k_cos(y0, y1),
2 => -k_sin(y0, y1, 1),
_ => -k_cos(y0, y1),
}
}
#[test]
fn test_near_pi() {
let x = f64::from_bits(0x400921fb000FD5DD); // 3.141592026217707
let sx = f64::from_bits(0x3ea50d15ced1a4a2); // 6.273720864039205e-7
let result = sin(x);
#[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
let result = force_eval!(result);
assert_eq!(result, sx);
}