正弦函数解析（sin.rs）

news2025/1/4 10:09:26

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);
}