softposit/p16e1/math/
sin_pi.rs

1use super::P16E1;
2
3impl P16E1 {
4    pub const fn sin_pi(self) -> Self {
5        let ui_a = self.to_bits();
6
7        let mut f = ui_a as u64;
8
9        let mut sign = f & 0x8000;
10        if sign != 0 {
11            f = 0x10000 - f; // 2's complement if negative
12        }
13        if f > 31743 {
14            // input value is an integer?
15            if f == 0x8000 {
16                return Self::NAR; // sinpi(NaR) is NaR
17            } else {
18                return Self::ZERO; // sinpi of an integer is zero
19            }
20        }
21        if f == 0 {
22            // sinpi(0) = 0
23            return Self::ZERO;
24        }
25        let mut s: i32;
26        if (f & 0x4000) != 0 {
27            // decode regime
28            s = 16;
29            while (f & 0x2000) != 0 {
30                f <<= 1;
31                s += 2;
32            }
33        } else {
34            s = 14;
35            while (f & 0x2000) == 0 {
36                f <<= 1;
37                s -= 2;
38            }
39        }
40        if (f & 0x1000) != 0 {
41            s += 1; // decode exponent
42        }
43        f = (f & 0x0FFF) | 0x1000; // get 12-bit fraction and restore hidden bit
44        f = if s < 0 { f >> -s } else { f << s };
45        f &= 0x_1FFF_FFFF; // fixed-point with 28-bit fraction
46        let mut s = f >> 27; // the quadrant is the multiple of 1/2
47        f &= 0x_07FF_FFFF; // input value modulo 1/2
48        if (s & 2) != 0 {
49            sign ^= 0x8000; // quadrants 2 and 3 flip the sign
50        }
51        if f == 0 {
52            return Self::from_bits(if (s & 1) != 0 {
53                (sign as u16) | 0x4000
54            } else {
55                0
56            });
57        }
58        if (s & 1) != 0 {
59            f = 0x_0800_0000 - f;
60        }
61        f = poly(f);
62        s = 1; // convert 28-bit fixed-point to a posit
63        while (f & 0x_0800_0000) == 0 {
64            f <<= 1;
65            s += 1;
66        }
67        let bit = s & 1;
68        s = (s >> 1) + 14 + bit;
69        if bit == 0 {
70            f &= 0x_07FF_FFFF; // encode exponent bit
71        }
72        f |= 0x_1000_0000; // encode regime termination bit
73        let bit = 1_u64 << (s - 1);
74        if ((f & bit) != 0) && (((f & (bit - 1)) != 0) || ((f & (bit << 1)) != 0)) {
75            // round to nearest, tie to even
76            f += bit;
77        }
78        f >>= s;
79        Self::from_bits((if sign != 0 { 0x10000 - f } else { f }) as u16)
80    }
81}
82
83#[inline]
84const fn poly(f: u64) -> u64 {
85    if f < 0x_000A_5801 {
86        return (f * 102_943) >> 15; // linear approximation suffices
87    }
88    let fs = f >> 11;
89    let fsq = (fs * fs) >> 8;
90    let mut s = (fsq * 650) >> 25;
91    s = (fsq * (9_813 - s)) >> 23;
92    s = (fsq * (334_253 - s)) >> 23;
93    s = (fsq * (5_418_741 - s)) >> 22;
94    (fs * (52_707_180 - s)) >> 13
95}
96
97#[test]
98fn test_sin_pi() {
99    use rand::Rng;
100    let mut rng = rand::thread_rng();
101    for _ in 0..crate::NTESTS16 {
102        let p_a: P16E1 = rng.gen();
103        let f_a = f64::from(p_a);
104        let p = p_a.sin_pi();
105        let f = (f_a * core::f64::consts::PI).sin();
106        let expected = P16E1::from(f);
107        if p.is_zero() {
108            continue;
109        }
110        assert_eq!(p, expected);
111    }
112}
softposit/p16e1/math/sin_pi.rs

softposit/p16e1/math/
sin_pi.rs
