1#![allow(dead_code, clippy::cast_possible_truncation, clippy::cast_lossless)]
8
9use crate::algorithms::DoubleWord;
10use core::num::Wrapping;
11
12pub use self::{reciprocal_2_mg10 as reciprocal_2, reciprocal_mg10 as reciprocal};
13
14#[doc = crate::algorithms::unstable_warning!()]
16#[inline(always)]
19#[must_use]
20pub fn reciprocal_ref(d: u64) -> u64 {
21 debug_assert!(d >= (1 << 63));
22 let r = u128::MAX / u128::from(d);
23 debug_assert!(r >= (1 << 64));
24 debug_assert!(r < (1 << 65));
25 r as u64
26}
27
28#[doc = crate::algorithms::unstable_warning!()]
30#[inline(always)]
52#[must_use]
53pub fn reciprocal_mg10(d: u64) -> u64 {
54 const ZERO: Wrapping<u64> = Wrapping(0);
55 const ONE: Wrapping<u64> = Wrapping(1);
56
57 static TABLE: [u16; 256] = [
59 2045, 2037, 2029, 2021, 2013, 2005, 1998, 1990, 1983, 1975, 1968, 1960, 1953, 1946, 1938,
60 1931, 1924, 1917, 1910, 1903, 1896, 1889, 1883, 1876, 1869, 1863, 1856, 1849, 1843, 1836,
61 1830, 1824, 1817, 1811, 1805, 1799, 1792, 1786, 1780, 1774, 1768, 1762, 1756, 1750, 1745,
62 1739, 1733, 1727, 1722, 1716, 1710, 1705, 1699, 1694, 1688, 1683, 1677, 1672, 1667, 1661,
63 1656, 1651, 1646, 1641, 1636, 1630, 1625, 1620, 1615, 1610, 1605, 1600, 1596, 1591, 1586,
64 1581, 1576, 1572, 1567, 1562, 1558, 1553, 1548, 1544, 1539, 1535, 1530, 1526, 1521, 1517,
65 1513, 1508, 1504, 1500, 1495, 1491, 1487, 1483, 1478, 1474, 1470, 1466, 1462, 1458, 1454,
66 1450, 1446, 1442, 1438, 1434, 1430, 1426, 1422, 1418, 1414, 1411, 1407, 1403, 1399, 1396,
67 1392, 1388, 1384, 1381, 1377, 1374, 1370, 1366, 1363, 1359, 1356, 1352, 1349, 1345, 1342,
68 1338, 1335, 1332, 1328, 1325, 1322, 1318, 1315, 1312, 1308, 1305, 1302, 1299, 1295, 1292,
69 1289, 1286, 1283, 1280, 1276, 1273, 1270, 1267, 1264, 1261, 1258, 1255, 1252, 1249, 1246,
70 1243, 1240, 1237, 1234, 1231, 1228, 1226, 1223, 1220, 1217, 1214, 1211, 1209, 1206, 1203,
71 1200, 1197, 1195, 1192, 1189, 1187, 1184, 1181, 1179, 1176, 1173, 1171, 1168, 1165, 1163,
72 1160, 1158, 1155, 1153, 1150, 1148, 1145, 1143, 1140, 1138, 1135, 1133, 1130, 1128, 1125,
73 1123, 1121, 1118, 1116, 1113, 1111, 1109, 1106, 1104, 1102, 1099, 1097, 1095, 1092, 1090,
74 1088, 1086, 1083, 1081, 1079, 1077, 1074, 1072, 1070, 1068, 1066, 1064, 1061, 1059, 1057,
75 1055, 1053, 1051, 1049, 1047, 1044, 1042, 1040, 1038, 1036, 1034, 1032, 1030, 1028, 1026,
76 1024,
77 ];
78
79 debug_assert!(d >= (1 << 63));
80 let d = Wrapping(d);
81
82 let d0 = d & ONE;
83 let d9 = d >> 55;
84 let d40 = ONE + (d >> 24);
85 let d63 = (d + ONE) >> 1;
86 let v0 = Wrapping(*unsafe { TABLE.get_unchecked((d9.0 - 256) as usize) } as u64);
88 let v1 = (v0 << 11) - ((v0 * v0 * d40) >> 40) - ONE;
89 let v2 = (v1 << 13) + ((v1 * ((ONE << 60) - v1 * d40)) >> 47);
90 let e = ((v2 >> 1) & (ZERO - d0)) - v2 * d63;
91 let v3 = (mul_hi(v2, e) >> 1) + (v2 << 31);
92 let v4 = v3 - muladd_hi(v3, d, d) - d;
93
94 v4.0
95}
96
97#[doc = crate::algorithms::unstable_warning!()]
99#[inline(always)]
106#[must_use]
107pub fn reciprocal_2_mg10(d: u128) -> u64 {
108 debug_assert!(d >= (1 << 127));
109 let (d0, d1) = d.split();
110
111 let mut v = reciprocal(d1);
112 let (mut p, overflow) = d1.wrapping_mul(v).overflowing_add(d0);
113 if overflow {
114 v = v.wrapping_sub(1);
115 if p >= d1 {
116 v = v.wrapping_sub(1);
117 p = p.wrapping_sub(d1);
118 }
119 p = p.wrapping_sub(d1);
120 }
121 let (t0, t1) = u128::mul(v, d0).split();
122
123 let (p, overflow) = p.overflowing_add(t1);
124 if overflow {
125 v = v.wrapping_sub(1);
126 if u128::join(p, t0) >= d {
127 v = v.wrapping_sub(1);
128 }
129 }
130 v
131}
132
133#[inline(always)]
134#[must_use]
135fn mul_hi(a: Wrapping<u64>, b: Wrapping<u64>) -> Wrapping<u64> {
136 Wrapping(u128::mul(a.0, b.0).high())
137}
138
139#[inline(always)]
140#[must_use]
141fn muladd_hi(a: Wrapping<u64>, b: Wrapping<u64>, c: Wrapping<u64>) -> Wrapping<u64> {
142 Wrapping(u128::muladd(a.0, b.0, c.0).high())
143}
144
145#[cfg(test)]
146mod tests {
147 use super::*;
148 use proptest::proptest;
149
150 #[test]
151 fn test_reciprocal() {
152 proptest!(|(n: u64)| {
153 let n = n | (1 << 63);
154 let expected = reciprocal_ref(n);
155 let actual = reciprocal_mg10(n);
156 assert_eq!(expected, actual);
157 });
158 }
159
160 #[test]
161 fn test_reciprocal_2() {
162 assert_eq!(reciprocal_2_mg10(1 << 127), u64::MAX);
163 assert_eq!(reciprocal_2_mg10(u128::MAX), 0);
164 assert_eq!(
165 reciprocal_2_mg10(0xd555_5555_5555_5555_5555_5555_5555_5555),
166 0x3333_3333_3333_3333
167 );
168 assert_eq!(
169 reciprocal_2_mg10(0xd0e7_57b0_2171_5fbe_cba4_ad0e_825a_e500),
170 0x39b6_c5af_970f_86b3
171 );
172 assert_eq!(
173 reciprocal_2_mg10(0xae5d_6551_8a51_3208_a850_5491_9637_eb17),
174 0x77db_09d1_5c3b_970b
175 );
176 }
177}