libm 0.2.16

libm in pure Rust
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121

//! Architecture-specific support for x86-32 without SSE2
//!
//! We use an alternative implementation on x86, because the
//! main implementation fails with the x87 FPU used by
//! debian i386, probably due to excess precision issues.
//!
//! See https://github.com/rust-lang/compiler-builtins/pull/976 for discussion on why these
//! functions are implemented in this way.

pub fn ceil(mut x: f64) -> f64 {
    unsafe {
        core::arch::asm!(
            "fld qword ptr [{x}]",
            // Save the FPU control word, using `x` as scratch space.
            "fstcw [{x}]",
            // Set rounding control to 0b10 (+∞).
            "mov word ptr [{x} + 2], 0x0b7f",
            "fldcw [{x} + 2]",
            // Round.
            "frndint",
            // Restore FPU control word.
            "fldcw [{x}]",
            // Save rounded value to memory.
            "fstp qword ptr [{x}]",
            x = in(reg) &mut x,
            // All the x87 FPU stack is used, all registers must be clobbered
            out("st(0)") _, out("st(1)") _,
            out("st(2)") _, out("st(3)") _,
            out("st(4)") _, out("st(5)") _,
            out("st(6)") _, out("st(7)") _,
            options(nostack),
        );
    }
    x
}

pub fn floor(mut x: f64) -> f64 {
    unsafe {
        core::arch::asm!(
            "fld qword ptr [{x}]",
            // Save the FPU control word, using `x` as scratch space.
            "fstcw [{x}]",
            // Set rounding control to 0b01 (-∞).
            "mov word ptr [{x} + 2], 0x077f",
            "fldcw [{x} + 2]",
            // Round.
            "frndint",
            // Restore FPU control word.
            "fldcw [{x}]",
            // Save rounded value to memory.
            "fstp qword ptr [{x}]",
            x = in(reg) &mut x,
            // All the x87 FPU stack is used, all registers must be clobbered
            out("st(0)") _, out("st(1)") _,
            out("st(2)") _, out("st(3)") _,
            out("st(4)") _, out("st(5)") _,
            out("st(6)") _, out("st(7)") _,
            options(nostack),
        );
    }
    x
}
/// Implements the exponential functions with `x87` assembly.
///
/// This relies on the instruction `f2xm1`, which computes `2^x - 1` (for
/// |x| < 1). This transcendental instruction is documented to produce results
/// with error below 1ulp (in the native double-extended precision format). This
/// translates to correctly rounded results for f32, but results in f64 may have
/// 1ulp error, which may depend on the hardware.
macro_rules! x87exp {
    ($float_ty:ident, $word_size:literal, $fn_name:ident,  $load_op:literal) => {
        pub fn $fn_name(mut x: $float_ty) -> $float_ty { unsafe {
            core::arch::asm!(
                // Prepare the register stack as
                // ```
                // st(0) = y = x*log2(base)
                // st(1) = 1.0
                // st(2) = round(y)
                // ```
                concat!($load_op, " ", $word_size, " ptr [{x}]"),
                "fld1",
                "fld st(1)",
                "frndint",
                "fxch st(2)",

                // Compare y with round(y) to determine if y is finite and
                // not an integer. If so, compute `exp2(y - round(y))` into
                // st(1). Otherwise skip ahead with `st(1) = 1.0`
                "fucom st(2)",
                "fstsw ax",
                "test ax, 0x4000",
                "jnz 2f",
                "fsub st(0), st(2)", // st(0) = y - round(y)
                "f2xm1",             // st(0) = 2^st(0) - 1.0
                "fadd st(1), st(0)", // st(1) = 1 + st(0) = exp2(y - round(y))
                "2:",

                // Finally, scale by `exp2(round(y))` and clear the stack.
                "fstp st(0)",
                "fscale",
                concat!("fstp ", $word_size, " ptr [{x}]"),
                "fstp st(0)",
                x = in(reg) &mut x,
                out("ax") _,
                out("st(0)") _, out("st(1)") _,
                out("st(2)") _, out("st(3)") _,
                out("st(4)") _, out("st(5)") _,
                out("st(6)") _, out("st(7)") _,
                options(nostack),
            );
            x
        }}
    };
}

x87exp!(f32, "dword", x87_exp2f, "fld");
x87exp!(f64, "qword", x87_exp2, "fld");
x87exp!(f32, "dword", x87_exp10f, "fldl2t\nfmul");
x87exp!(f64, "qword", x87_exp10, "fldl2t\nfmul");
x87exp!(f32, "dword", x87_expf, "fldl2e\nfmul");
x87exp!(f64, "qword", x87_exp, "fldl2e\nfmul");