//! Table-based logarithm and exponential for f64 `pow` and `root`, from Arm Optimized Routines' `pow`. The logarithm reduces x by a centre c whose reciprocal has few bits, so that z/c - 1 needs no division, and looks up ln c. The exponential reduces by multiples of ln 2/256 and looks up 2^(j/256) as a float and the tail of its rounding, which leaves a short polynomial. Both steps carry a double word. Each item reads its own table entries, so LLVM vectorizes the loops that call them only where the target has gathers, as AVX2 does. On NEON and in wasm they run one item at a time. `log2_single` and `exp2_single` are the logarithm and exponential of Arm's scalar `powf`, for `Pow`'s and `Root`'s scalar forms at f32 and f16. See `THIRD-PARTY.md`.
use crate::{dd::Double, exp::clamp, float::sealed::Core};
use std::f64::consts::LOG2_E;
/// `ix - off`, k as an f64 and the bits of z, for the float with bits `ix` = z × 2^k and z from the float with bits `off` to twice it. The table index is in the bits of `ix - off` below its top 12, which hold k. They're read with a bias of 2048, so that no shift needs a sign, and k converts through `Core::from_low`.
#[inline]
fn split(ix: u64, off: u64) -> (u64, f64, u64) {
let tmp = ix.wrapping_sub(off);
let k = f64::from_low(tmp.wrapping_add(2048 << 52) >> 52, 2048.0);
(tmp, k, ix.wrapping_sub(tmp & (0xfff << 52)))
}
/// The bits of 0x1.69555p-1. z lies from it to twice it, which keeps ln z small near x = 1.
const OFF: u64 = 0x3fe6955500000000;
/// 1/c for each of the 128 parts of z's range. Each has at most 8 significant bits, so z/c - 1 is exact.
const INVC: [f64; 128] = [
1.4140625, 1.40625, 1.3984375, 1.390625, 1.3828125, 1.375, 1.3671875, 1.359375, 1.3515625, 1.34375, 1.3359375, 1.3359375, 1.328125, 1.3203125, 1.3125,
1.3046875, 1.296875, 1.2890625, 1.2890625, 1.28125, 1.2734375, 1.265625, 1.2578125, 1.25, 1.25, 1.2421875, 1.234375, 1.2265625, 1.2265625, 1.21875,
1.2109375, 1.203125, 1.203125, 1.1953125, 1.1875, 1.1875, 1.1796875, 1.171875, 1.171875, 1.1640625, 1.15625, 1.1484375, 1.1484375, 1.140625, 1.140625,
1.1328125, 1.125, 1.125, 1.1171875, 1.109375, 1.109375, 1.1015625, 1.1015625, 1.09375, 1.0859375, 1.0859375, 1.078125, 1.078125, 1.0703125, 1.0625, 1.0625,
1.0546875, 1.0546875, 1.046875, 1.046875, 1.0390625, 1.0390625, 1.03125, 1.03125, 1.0234375, 1.015625, 1.015625, 1.0078125, 1.0078125, 1.0, 1.0, 0.9921875,
0.984375, 0.9765625, 0.96875, 0.9609375, 0.953125, 0.9453125, 0.94140625, 0.93359375, 0.92578125, 0.91796875, 0.9140625, 0.90625, 0.8984375, 0.89453125,
0.88671875, 0.8828125, 0.875, 0.87109375, 0.86328125, 0.859375, 0.8515625, 0.84765625, 0.83984375, 0.8359375, 0.83203125, 0.82421875, 0.8203125, 0.8125,
0.80859375, 0.8046875, 0.80078125, 0.79296875, 0.7890625, 0.78515625, 0.78125, 0.7734375, 0.76953125, 0.765625, 0.76171875, 0.7578125, 0.75390625,
0.74609375, 0.7421875, 0.73828125, 0.734375, 0.73046875, 0.7265625, 0.72265625, 0.71875, 0.71484375, 0.7109375,
];
/// ln c rounded to a multiple of 2^-43, which keeps k × `LN2.0` + ln c exact.
const LOGC: [f64; 128] = [
-0.3464667673462145, -0.34092658697056777, -0.3353555419211034, -0.3297532863724655, -0.32411946865420305, -0.31845373111855224, -0.31275571000389846,
-0.3070250352949415, -0.3012613305781997, -0.2954642128938758, -0.28963329258306203, -0.28963329258306203, -0.28376817313062475, -0.27786845100342816,
-0.2719337154836694, -0.2659635484970977, -0.25995752443691345, -0.25391520998095984, -0.25391520998095984, -0.2478361639045943, -0.2417199368871934,
-0.23556607131274632, -0.22937410106487732, -0.22314355131425145, -0.22314355131425145, -0.21687393830063684, -0.21056476910735, -0.2042155414286526,
-0.2042155414286526, -0.19782574332987224, -0.19139485299967873, -0.18492233849406148, -0.18492233849406148, -0.1784076574728033, -0.17185025692663203,
-0.17185025692663203, -0.1652495728952772, -0.15860503017665906, -0.15860503017665906, -0.15191604202584585, -0.1451820098444614, -0.13840232285906495,
-0.13840232285906495, -0.131576357788731, -0.131576357788731, -0.12470347850091912, -0.11778303565643, -0.11778303565643, -0.11081436634026431,
-0.10379679368168127, -0.10379679368168127, -0.09672962645856842, -0.09672962645856842, -0.089612158689647, -0.08244366921110213, -0.08244366921110213,
-0.07522342123763792, -0.07522342123763792, -0.06795066190852594, -0.06062462181648698, -0.06062462181648698, -0.053244514518837605, -0.053244514518837605,
-0.045809536031242715, -0.045809536031242715, -0.038318864302141264, -0.038318864302141264, -0.03077165866670839, -0.03077165866670839,
-0.023167059281490765, -0.015504186535963527, -0.015504186535963527, -0.0077821404420319595, -0.0077821404420319595, 0.0, 0.0, 0.007843177461040796,
0.01574835696817445, 0.023716526617363343, 0.03174869831457272, 0.039845908547249564, 0.048009219186383234, 0.056239718322899535, 0.06038051098892083,
0.06871389254808946, 0.07711730334438016, 0.08559193033545398, 0.08985632912185793, 0.09844007281321865, 0.10709813555638448, 0.11145544092528326,
0.12022742699821265, 0.12464244520731427, 0.13353139262449076, 0.13800567301939282, 0.14701474296180095, 0.15154989812720032, 0.16068238169043525,
0.16528009093906348, 0.17453941635187675, 0.17920142945774842, 0.18388527877016259, 0.1933193110035063, 0.19806991376208316, 0.20763936477828793,
0.21245865121420593, 0.21730127569003344, 0.22216746534115828, 0.2319714654378231, 0.2369097470783572, 0.24187253642048745, 0.2468600779315011,
0.2569104137850218, 0.26197371574153294, 0.2670627852490952, 0.27217788591576664, 0.27731928541618345, 0.28248725557466514, 0.29290401643288533,
0.29815337231912054, 0.3034304294199046, 0.30873548164959175, 0.3140688276249648, 0.31943077076641657, 0.3248216194012912, 0.33024168687052224,
0.3356912916381134, 0.3411707574027787,
];
/// The rest of ln c.
const LOGC_TAIL: [f64; 128] = [
5.929407345889625e-15, -2.544157440035963e-14, -3.443525940775045e-14, -2.500123826022799e-15, -8.929337133850617e-15, 1.7625431312172662e-14,
1.5688303180062087e-15, 2.9655274673691784e-14, 3.7923164802093147e-14, 3.993416384387844e-14, 1.9352855826489123e-14, 1.9352855826489123e-14,
-1.9852665484979036e-14, -2.814323765595281e-14, 2.7643769993528702e-14, -4.025092402293806e-14, -1.2621729398885316e-14, -3.600176732637335e-15,
-3.600176732637335e-15, 1.3029797173308663e-14, 4.8230289429940886e-14, -2.0592242769647135e-14, 3.149265065191484e-14, 4.169796584527195e-14,
4.169796584527195e-14, 2.2477465222466186e-14, 3.6507188831790577e-16, -3.827767260205414e-14, -3.827767260205414e-14, -4.7641388950792196e-14,
4.9278276214647115e-14, 4.9485167661250996e-14, 4.9485167661250996e-14, -1.5003333854266542e-14, -2.7194441649495324e-14, -2.7194441649495324e-14,
-2.99659267292569e-14, 2.0472357800461955e-14, 2.0472357800461955e-14, 3.879296723063646e-15, -3.6506824353335045e-14, -5.4183331379008994e-14,
-5.4183331379008994e-14, 1.1729485484531301e-14, 1.1729485484531301e-14, -3.811763084710266e-14, 4.654729747598445e-14, 4.654729747598445e-14,
-2.5799991283069902e-14, 3.7700471749674615e-14, 3.7700471749674615e-14, 1.7306161136093256e-14, 1.7306161136093256e-14, -4.012913552726574e-14,
2.7541708360737882e-14, 2.7541708360737882e-14, 5.0396178134370583e-14, 5.0396178134370583e-14, 1.8195060030168815e-14, 5.213620639136504e-14,
5.213620639136504e-14, 2.532168943117445e-14, 2.532168943117445e-14, -5.148849572685811e-14, -5.148849572685811e-14, 4.6652946995830086e-15,
4.6652946995830086e-15, -4.529814257790929e-14, -4.529814257790929e-14, -4.361324067851568e-14, -1.7274567499706107e-15, -1.7274567499706107e-15,
-2.298941004620351e-14, -2.298941004620351e-14, 0.0, 0.0, -1.4902732911301337e-14, -3.527980389655325e-14, -4.730054772033249e-14, 7.580310369375161e-15,
-4.9893776716773285e-14, -2.262629393030674e-14, -2.345674491018699e-14, -1.3352588834854848e-14, -3.765296820388875e-14, 5.1128335719851986e-14,
-5.046674438470119e-14, 3.1218748807418837e-15, 3.3871241029241416e-14, -1.7376727386423858e-14, 3.957125899799804e-14, -5.2849453521890294e-14,
-3.767012502308738e-14, 3.1859736349078334e-14, 5.0900642926060466e-14, 8.710783796122478e-15, 6.157896229122976e-16, 3.821577743916796e-14,
3.9440046718453496e-14, 2.2924522154618074e-14, -3.742530094732263e-14, -2.5223102140407338e-14, -1.0320443688698849e-14, 1.0634128304268335e-14,
-4.3425422595242564e-14, -1.2527395755711364e-14, -5.204008743405884e-14, -3.979844515951702e-15, -4.7955860343296286e-14, 5.015686013791602e-16,
-7.252318953240293e-16, 2.4688324156011588e-14, 5.465121253624792e-15, 4.102651071698446e-14, -4.996736502345936e-14, 4.903580708156347e-14,
5.089628039500759e-14, 1.1782016386565151e-14, 4.727452940514406e-14, -4.4204083338755686e-14, 1.548345993498083e-14, 2.1522127491642888e-14,
1.1054030169005386e-14, -5.534326352070679e-14, -5.351646604259541e-14, 5.4612144489920215e-14, 2.8136969901227338e-14, -1.156568624616423e-14,
];
/// ln 2, split into a part whose product with any exponent is exact, and the rest.
const LN2: (f64, f64) = (0.6931471805598903, 5.497923018708371e-14);
/// -1/2, then the coefficients of ln(1 + r) - r + r²/2 in the grouping that `ln` evaluates, for |r| up to 1/128.
const A: [f64; 7] = [-0.5, -0.6666666666666679, 0.5000000000000007, 0.7999999995323976, -0.6666666663487739, -1.142909628459501, 1.0000415263675542];
/// ln x as a double word with about 15 bits beyond an f64's, for a positive x, subnormals included. 0 gives -∞, and ∞ and NaN give themselves. It's inlined always: LLVM would otherwise leave it out of line, which would stop the loops that call it vectorizing.
#[inline(always)]
pub(crate) fn ln<const FUSED: bool>(x: f64) -> Double<f64> {
let ix = if x < f64::MIN_POSITIVE { (x * 4503599627370496.0).to_bits().wrapping_sub(52 << 52) } else { x.to_bits() };
let y = ln_bits::<FUSED>(ix);
let y = if x.is_nan() | (x == f64::INFINITY) { Double::new(x, 0.0) } else { y };
if x == 0.0 { Double::new(f64::NEG_INFINITY, 0.0) } else { y }
}
/// `ln` for a positive normal x, with none of its special cases. Other arguments give a meaningless result.
#[inline(always)]
pub(crate) fn ln_normal<const FUSED: bool>(x: f64) -> Double<f64> { ln_bits::<FUSED>(x.to_bits()) }
/// ln of the float with bits `ix`, read as a positive normal float. `ln` passes a subnormal's bits scaled into that range.
#[inline(always)]
fn ln_bits<const FUSED: bool>(ix: u64) -> Double<f64> {
let (tmp, k, iz) = split(ix, OFF);
let i = (tmp >> 45) as usize & 127;
let (z, invc) = (f64::from_bits(iz), INVC[i]);
// Without FMA, z splits into a high part with 21 bits and the rest, so that rhi and rhi² are exact.
let (r, rhi, rlo) = if FUSED { let r = z.mla::<true>(invc, -1.0); (r, r, 0.0) } else {
let zhi = f64::from_bits(iz.wrapping_add(1 << 31) & (!0 << 32));
let (rhi, rlo) = (zhi * invc - 1.0, (z - zhi) * invc);
(rhi + rlo, rhi, rlo)
};
let t1 = k.mla::<FUSED>(LN2.0, LOGC[i]);
let t2 = t1 + r;
let lo1 = k.mla::<FUSED>(LN2.1, LOGC_TAIL[i]);
let lo2 = t1 - t2 + r;
let (ar, ar2) = (A[0] * r, r * (A[0] * r));
let (hi, lo3, lo4) = if FUSED { let hi = t2 + ar2; (hi, ar.mla::<true>(r, -ar2), t2 - hi + ar2) } else {
let arhi = A[0] * rhi;
let hi = t2 + rhi * arhi;
(hi, rlo * (ar + arhi), t2 - hi + rhi * arhi)
};
let p = r.mla::<FUSED>(A[6], A[5]).mla::<FUSED>(ar2, r.mla::<FUSED>(A[4], A[3])).mla::<FUSED>(ar2, r.mla::<FUSED>(A[2], A[1])) * (r * ar2);
let lo = lo1 + lo2 + lo3 + lo4 + p;
Double::new(hi + lo, hi - (hi + lo) + lo)
}
/// 2^(j/256) for j from 0 to 255, each with j × 2^44 subtracted from its bits, so that adding k × 2^44 for the whole k = 256e + j gives the bits of 2^(k/256) while the exponent stays in range.
const SBITS: [u64; 256] = [
0x3ff0000000000000, 0x3feffb1afa5abcbf, 0x3feff63da9fb3335, 0x3feff168143b0281, 0x3fefec9a3e778061, 0x3fefe7d42e11bbcc, 0x3fefe315e86e7f85,
0x3fefde5f72f654b1, 0x3fefd9b0d3158574, 0x3fefd50a0e3c1f89, 0x3fefd06b29ddf6de, 0x3fefcbd42b72a836, 0x3fefc74518759bc8, 0x3fefc2bdf66607e0,
0x3fefbe3ecac6f383, 0x3fefb9c79b1f3919, 0x3fefb5586cf9890f, 0x3fefb0f145e46c85, 0x3fefac922b7247f7, 0x3fefa83b23395dec, 0x3fefa3ec32d3d1a2,
0x3fef9fa55fdfa9c5, 0x3fef9b66affed31b, 0x3fef973028d7233e, 0x3fef9301d0125b51, 0x3fef8edbab5e2ab6, 0x3fef8abdc06c31cc, 0x3fef86a814f204ab,
0x3fef829aaea92de0, 0x3fef7e95934f312e, 0x3fef7a98c8a58e51, 0x3fef76a45471c3c2, 0x3fef72b83c7d517b, 0x3fef6ed48695bbc0, 0x3fef6af9388c8dea,
0x3fef672658375d2f, 0x3fef635beb6fcb75, 0x3fef5f99f8138a1c, 0x3fef5be084045cd4, 0x3fef582f95281c6b, 0x3fef54873168b9aa, 0x3fef50e75eb44027,
0x3fef4d5022fcd91d, 0x3fef49c18438ce4d, 0x3fef463b88628cd6, 0x3fef42be3578a819, 0x3fef3f49917ddc96, 0x3fef3bdda27912d1, 0x3fef387a6e756238,
0x3fef351ffb82140a, 0x3fef31ce4fb2a63f, 0x3fef2e85711ece75, 0x3fef2b4565e27cdd, 0x3fef280e341ddf29, 0x3fef24dfe1f56381, 0x3fef21ba7591bb70,
0x3fef1e9df51fdee1, 0x3fef1b8a66d10f13, 0x3fef187fd0dad990, 0x3fef157e39771b2f, 0x3fef1285a6e4030b, 0x3fef0f961f641589, 0x3fef0cafa93e2f56,
0x3fef09d24abd886b, 0x3fef06fe0a31b715, 0x3fef0432edeeb2fd, 0x3fef0170fc4cd831, 0x3feefeb83ba8ea32, 0x3feefc08b26416ff, 0x3feef96266e3fa2d,
0x3feef6c55f929ff1, 0x3feef431a2de883b, 0x3feef1a7373aa9cb, 0x3feeef26231e754a, 0x3feeecae6d05d866, 0x3feeea401b7140ef, 0x3feee7db34e59ff7,
0x3feee57fbfec6cf4, 0x3feee32dc313a8e5, 0x3feee0e544ede173, 0x3feedea64c123422, 0x3feedc70df1c5175, 0x3feeda4504ac801c, 0x3feed822c367a024,
0x3feed60a21f72e2a, 0x3feed3fb2709468a, 0x3feed1f5d950a897, 0x3feecffa3f84b9d4, 0x3feece086061892d, 0x3feecc2042a7d232, 0x3feeca41ed1d0057,
0x3feec86d668b3237, 0x3feec6a2b5c13cd0, 0x3feec4e1e192aed2, 0x3feec32af0d7d3de, 0x3feec17dea6db7d7, 0x3feebfdad5362a27, 0x3feebe41b817c114,
0x3feebcb299fddd0d, 0x3feebb2d81d8abff, 0x3feeb9b2769d2ca7, 0x3feeb8417f4531ee, 0x3feeb6daa2cf6642, 0x3feeb57de83f4eef, 0x3feeb42b569d4f82,
0x3feeb2e2f4f6ad27, 0x3feeb1a4ca5d920f, 0x3feeb070dde910d2, 0x3feeaf4736b527da, 0x3feeae27dbe2c4cf, 0x3feead12d497c7fd, 0x3feeac0827ff07cc,
0x3feeab07dd485429, 0x3feeaa11fba87a03, 0x3feea9268a5946b7, 0x3feea84590998b93, 0x3feea76f15ad2148, 0x3feea6a320dceb71, 0x3feea5e1b976dc09,
0x3feea52ae6cdf6f4, 0x3feea47eb03a5585, 0x3feea3dd1d1929fd, 0x3feea34634ccc320, 0x3feea2b9febc8fb7, 0x3feea23882552225, 0x3feea1c1c70833f6,
0x3feea155d44ca973, 0x3feea0f4b19e9538, 0x3feea09e667f3bcd, 0x3feea052fa75173e, 0x3feea012750bdabf, 0x3fee9fdcddd47645, 0x3fee9fb23c651a2f,
0x3fee9f9298593ae5, 0x3fee9f7df9519484, 0x3fee9f7466f42e87, 0x3fee9f75e8ec5f74, 0x3fee9f8286ead08a, 0x3fee9f9a48a58174, 0x3fee9fbd35d7cbfd,
0x3fee9feb564267c9, 0x3feea024b1ab6e09, 0x3feea0694fde5d3f, 0x3feea0b938ac1cf6, 0x3feea11473eb0187, 0x3feea17b0976cfdb, 0x3feea1ed0130c132,
0x3feea26a62ff86f0, 0x3feea2f336cf4e62, 0x3feea3878491c491, 0x3feea427543e1a12, 0x3feea4d2add106d9, 0x3feea589994cce13, 0x3feea64c1eb941f7,
0x3feea71a4623c7ad, 0x3feea7f4179f5b21, 0x3feea8d99b4492ed, 0x3feea9cad931a436, 0x3feeaac7d98a6699, 0x3feeabd0a478580f, 0x3feeace5422aa0db,
0x3feeae05bad61778, 0x3feeaf3216b5448c, 0x3feeb06a5e0866d9, 0x3feeb1ae99157736, 0x3feeb2fed0282c8a, 0x3feeb45b0b91ffc6, 0x3feeb5c353aa2fe2,
0x3feeb737b0cdc5e5, 0x3feeb8b82b5f98e5, 0x3feeba44cbc8520f, 0x3feebbdd9a7670b3, 0x3feebd829fde4e50, 0x3feebf33e47a22a2, 0x3feec0f170ca07ba,
0x3feec2bb4d53fe0d, 0x3feec49182a3f090, 0x3feec674194bb8d5, 0x3feec86319e32323, 0x3feeca5e8d07f29e, 0x3feecc667b5de565, 0x3feece7aed8eb8bb,
0x3feed09bec4a2d33, 0x3feed2c980460ad8, 0x3feed503b23e255d, 0x3feed74a8af46052, 0x3feed99e1330b358, 0x3feedbfe53c12e59, 0x3feede6b5579fdbf,
0x3feee0e521356eba, 0x3feee36bbfd3f37a, 0x3feee5ff3a3c2774, 0x3feee89f995ad3ad, 0x3feeeb4ce622f2ff, 0x3feeee07298db666, 0x3feef0ce6c9a8952,
0x3feef3a2b84f15fb, 0x3feef68415b749b1, 0x3feef9728de5593a, 0x3feefc6e29f1c52a, 0x3feeff76f2fb5e47, 0x3fef028cf22749e4, 0x3fef05b030a1064a,
0x3fef08e0b79a6f1f, 0x3fef0c1e904bc1d2, 0x3fef0f69c3f3a207, 0x3fef12c25bd71e09, 0x3fef16286141b33d, 0x3fef199bdd85529c, 0x3fef1d1cd9fa652c,
0x3fef20ab5fffd07a, 0x3fef244778fafb22, 0x3fef27f12e57d14b, 0x3fef2ba88988c933, 0x3fef2f6d9406e7b5, 0x3fef33405751c4db, 0x3fef3720dcef9069,
0x3fef3b0f2e6d1675, 0x3fef3f0b555dc3fa, 0x3fef43155b5bab74, 0x3fef472d4a07897c, 0x3fef4b532b08c968, 0x3fef4f87080d89f2, 0x3fef53c8eacaa1d6,
0x3fef5818dcfba487, 0x3fef5c76e862e6d3, 0x3fef60e316c98398, 0x3fef655d71ff6075, 0x3fef69e603db3285, 0x3fef6e7cd63a8315, 0x3fef7321f301b460,
0x3fef77d5641c0658, 0x3fef7c97337b9b5f, 0x3fef81676b197d17, 0x3fef864614f5a129, 0x3fef8b333b16ee12, 0x3fef902ee78b3ff6, 0x3fef953924676d76,
0x3fef9a51fbc74c83, 0x3fef9f7977cdb740, 0x3fefa4afa2a490da, 0x3fefa9f4867cca6e, 0x3fefaf482d8e67f1, 0x3fefb4aaa2188510, 0x3fefba1bee615a27,
0x3fefbf9c1cb6412a, 0x3fefc52b376bba97, 0x3fefcac948dd7274, 0x3fefd0765b6e4540, 0x3fefd632798844f8, 0x3fefdbfdad9cbe14, 0x3fefe1d802243c89,
0x3fefe7c1819e90d8, 0x3fefedba3692d514, 0x3feff3c22b8f71f1, 0x3feff9d96b2a23d9,
];
/// 2^(j/256) divided by its entry in `SBITS`, less 1, as the bits of an f64, from Arm's `exp_data.c`.
const TAIL: [u64; 256] = [
0x0000000000000000, 0xbc84e82fc61851ac, 0x3c9b3b4f1a88bf6e, 0xbc82985dd8521d32, 0xbc7160139cd8dc5d, 0x3c651e617061bfbd, 0xbc905e7a108766d1,
0x3c845fad437fa426, 0x3c8cd2523567f613, 0xbc954529642b232f, 0xbc8bce8023f98efa, 0x3c8293708ef5c32e, 0x3c60f74e61e6c861, 0xbc95b9280905b2a4,
0x3c90a3e45b33d399, 0x3c84f31f32c4b7e7, 0x3c979aa65d837b6d, 0x3c9407fb30d06420, 0x3c8eb51a92fdeffc, 0xbc9a5d04b3b9911b, 0x3c3ebe3d702f9cd1,
0xbc937a01f0739546, 0xbc6a033489906e0b, 0x3c8b8268b04ef0a5, 0xbc9556522a2fbd0e, 0xbc9ac46e44a2ebcc, 0xbc5080ef8c4eea55, 0xbc65704e90c9f860,
0xbc91c923b9d5f416, 0xbc897cea57e46280, 0x3c80d3e3e95c55af, 0x3c56f01429e2b9d2, 0xbc801b15eaa59348, 0x3c6e653b2459034b, 0xbc8f1ff055de323d,
0x3c92cc7ea345b7dc, 0x3c8b898c3f1353bf, 0x3c957bfb2876ea9e, 0xbc96d99c7611eb26, 0x3c8cdc1873af2155, 0x3c9aecf73e3a2f60, 0xbc9493684653a131,
0xbc8fe782cb86389d, 0xbc98e2899077520a, 0x3c8a6f4144a6c38d, 0x3c9120fcd4f59273, 0x3c807a05b0e4047d, 0x3c89b788c188c9b8, 0x3c968efde3a8a894,
0x3c877afbca90ef84, 0x3c875e18f274487d, 0x3c91512f082876ee, 0x3c80472b981fe7f2, 0x3c9a02f0c7d75ec6, 0xbc96b87b3f71085e, 0xbc803297e78260bf,
0x3c82f7e16d09ab31, 0xbc95b77e5ccd9fbf, 0xbc3d219b1a6fbffa, 0xbc91e75c40b4251e, 0x3c8b3782720c0ab4, 0x3c98a911f1f7785a, 0x3c6e149289cecb8f,
0xbc61e7c998db7dbb, 0x3c834d754db0abb6, 0x3c85425c11faadf4, 0x3c864201e2ac744c, 0xbc979517a03e2847, 0x3c8fdd395dd3f84a, 0xbc800e2a46da4bee,
0xbc86a3803b8e5b04, 0xbc87430803972b34, 0xbc924aedcc4b5068, 0xbc954de30ae02d94, 0xbc9907f81b512d8e, 0xbc94f2487e1c03ec, 0xbc71d1e83e9436d2,
0x3c914a5432fcb2f4, 0xbc991919b3ce1b15, 0x3c79c3bba5562a2f, 0x3c859f48a72a4c6d, 0xbc85a71612e21658, 0xbc9312607a28698a, 0x3c86421f6f1d24d6,
0xbc58a78f4817895b, 0xbc9348a6815fce65, 0xbc7c2c9b67499a1b, 0x3c835c43984d9871, 0x3c4363ed60c2ac11, 0xbc632afc8d9473a0, 0x3c9666093b0664ef,
0xbc95fc5e44de020e, 0x3c6ecce1daa10379, 0xbc7ea0148327c42f, 0x3c93ff8e3f0f1230, 0xbc7a843ad1a88022, 0x3c7690cebb7aafb0, 0x3c892ca3bf144e63,
0x3c931dbdeb54e077, 0xbc902c99b04aa8b0, 0xbc8f94340071a38e, 0x3c73e34f67e67118, 0xbc87deccdc93a349, 0xbc75a3b1197ba0f0, 0xbc78dec6bd0f385f,
0x3c81bd2888075068, 0xbc861246ec7b5cf6, 0xbc896be8ae89ef8f, 0x3c93350518fdd78e, 0xbc88e6ac90348602, 0x3c7b98b72f8a9b05, 0xbc91af7f1365c3ac,
0x3c9063e1e21c5409, 0xbc943a3540d1898a, 0x3c34c7855019c6ea, 0xbc951f58ddaa8090, 0x3c9432e62b64c035, 0xbc82e1648e50a17c, 0xbc8ce44a6199769f,
0x3c95f30eda98a575, 0xbc8c33c53bef4da8, 0x3c917ecda8a72159, 0xbc845378892be9ae, 0xbc9345f3cee1ae6e, 0xbc93cedd78565858, 0xbc85c33fdf910406,
0x3c5710aa807e1964, 0x3c81079ab5789604, 0xbc93b3efbf5e2228, 0x3c727df161cd7778, 0xbc6a12ad8734b982, 0x3c93f9924a05b767, 0xbc6367efb86da9ee,
0xbc87557939a8b5ef, 0xbc80dc3d54e08851, 0x3c51ed2f56fa9d1a, 0xbc781f647e5a3ecf, 0xbc88e67a9006c909, 0xbc86ee4ac08b7db0, 0x3c86597566977ac8,
0xbc8619321e55e68a, 0x3c92c0b7028a5c3a, 0x3c909ccb5e09d4d3, 0x3c8a30faf49cc78c, 0xbc7b32dcb94da51d, 0xbc92dad3519d7b5b, 0x3c94ecfd5467c06b,
0x3c87d51410fd15c2, 0x3c65ebe1abd66c55, 0xbc760a3629969871, 0xbc88a1c52fb3cf42, 0x3c8b18c6e3fdef5d, 0xbc9369b6f13b3734, 0x3c90ec1ddcb1390a,
0xbc805e843a19ff1e, 0xbc522cea4f3afa1e, 0xbc94d450d872576e, 0x3c7c88549b958471, 0x3c90ad675b0e8a00, 0x3c931143962f7877, 0x3c8db72fc1f0eab4,
0x3c93e9e96f112479, 0xbc65b6609cc5e7ff, 0xbc8dac42a4a38df0, 0x3c7bf68359f35f44, 0x3c8b99dd98b1ed84, 0xbc93091fa71e3d83, 0xbc7885ad50cbb750,
0xbc5da9b88b6c1e29, 0xbc82d5e85f3e0301, 0xbc6c23f97c90b959, 0xbc51669428996971, 0xbc92434322f4f9aa, 0x3c71f2b2c1c4c014, 0xbc85ca6cd7668e4b,
0xbc9294f304f166b6, 0x3c71affc2b91ce27, 0xbc8a1e58414c07d3, 0x3c6dd235e10a73bb, 0xbc79740b58a20091, 0xbc87c50422622263, 0x3c9165830a2b96c2,
0x3c8b1c86e3e231d5, 0xbc903d5cbe27874b, 0xbc91bbd1d3bcbb15, 0x3c5986178980fce0, 0x3c90cc319cee31d2, 0xbc89472975b1f2a5, 0x3c8469846e735ab3,
0x3c7d8157a34b7e7f, 0xbc82dfcd978e9db4, 0x3c8c8a4e231ebb7d, 0x3c8c1a7792cb3387, 0xbc888c8d11a142e5, 0xbc907b8f4ad1d9fa, 0x3c889c2ea41433c7,
0xbc55c3d956dcaeba, 0xbc7274aedac8ff80, 0xbc90a40e3da6f640, 0x3c85c620ce76df06, 0xbc68d6f438ad9334, 0xbc8fda52e1b51e41, 0xbc91eee26b588a35,
0xbc32141a7b3e2cd8, 0x3c74ffd70a5fddcd, 0xbc302899507554e5, 0xbc91bdfbfa9298ac, 0xbc80dda2d4c0010c, 0x3c736eae30af0cb3, 0xbc8a007daadf8d68,
0x3c8ee3325c9ffd94, 0x3c836909391181d3, 0x3c84e08fd10959ac, 0xbc811cd7dbdf9547, 0x3c63cdaf384e1a67, 0xbc7ac28b7bef6621, 0x3c676b2c6c921968,
0xbc7030587207b9e1, 0xbc808a1883ccb5d2, 0xbc8cc734592af7fc, 0xbc8fad5d3ffffa6f, 0x3c87752a44f587e8, 0xbc900dae3875a949, 0x3c85b66fefeef52e,
0x3c74a385a63d07a7, 0x3c5159d9d908a96e, 0xbc82919e2040220f, 0x3c8c254d16117a68, 0x3c8e5a50d5c192ac, 0xbc8d8c329fbd0e03, 0x3c843a59ac016b4b,
0xbc8ea6e6fbd5f2a6, 0xbc82d52107b43e1f, 0xbc63e8e3eab2cbb4, 0xbc892ab93b470dc9, 0xbc8b7966cd0d2cd9, 0x3c74b604603a88d3, 0xbc776caa4c2ff1cf,
0x3c83c5ec519d7271, 0xbc81d5fc525d9940, 0xbc8ff7128fd391f0, 0x3c855cd8aaea3d21, 0xbc8dae98e223747d, 0x3c8269947c2bed4a, 0x3c8ec3bc41aa2008,
0xbc83b6137e9afe9e, 0x3c842b94c3a9eb32, 0xbc69fa74878ba7c7, 0x3c8a64a931d185ee, 0x3c901f3a75ee0efe, 0xbc8e37bae43be3ed, 0xbc516a9ce6ed84fa,
0x3c77893b4d91cd9d, 0xbc699c7db2effc76, 0x3c5305c14160cc89, 0x3c64b458677f9840,
];
/// ln 2/256, split into a part with 32 significant bits, whose product with any k here is exact, and the rest.
const LN2_256: (f64, f64) = (0.0027076061742263846, -1.6409824502660487e-13);
/// The coefficients of (e^r - 1 - r)/r² in the grouping that `exp` evaluates, for |r| up to ln 2/512.
const C: [f64; 3] = [0.49999999999999756, 0.1666666799409282, 0.04166667054490923];
/// e^(hi + lo), for a double word whose lo is far below hi. lo joins the reduced argument, except where `clamp` changes hi and the result is 0 or ∞ whatever lo is. The result is 2^e × 2^(j/256) × e^r, with 2^(j/256) carried as `SBITS` and its `TAIL`. 2^e multiplies in two factors, so a result below the normal floats rounds once, in the second.
#[inline(always)]
pub(crate) fn exp<const FUSED: bool>(hi: f64, lo: f64) -> f64 {
let x = clamp(hi, 1100.0);
let (ki, r) = reduce::<FUSED>(x);
let j = ki & 255;
let t = f64::from_bits(SBITS[j as usize] + (j << 44));
let y = t.mla::<FUSED>(series::<FUSED>(j, if x == hi { r + lo } else { r }), t);
// e + 2048, and its halves as exponent fields.
let e = (ki >> 8) - ((1 << 43) - 2048);
let h = e >> 1;
y * f64::from_bits((h - 1) << 52) * f64::from_bits((e - h - 1) << 52)
}
/// The largest |hi| that `exp_normal` takes. Below it, e^hi is a normal float.
pub(crate) const NORMAL: f64 = 512.0;
/// `exp` for |hi| below `NORMAL`, with no clamp. 2^e joins the exponent bits of 2^(j/256), so no multiply scales the result. Other arguments give a meaningless result.
#[inline(always)]
pub(crate) fn exp_normal<const FUSED: bool>(hi: f64, lo: f64) -> f64 {
let (ki, r) = reduce::<FUSED>(hi);
let j = ki & 255;
// `ki << 44` puts j and e in their bit positions and shifts out the 2^51.
let t = f64::from_bits(SBITS[j as usize].wrapping_add(ki << 44));
t.mla::<FUSED>(series::<FUSED>(j, r + lo), t)
}
/// The bits of `SHIFT` + k, whose low 52 hold 2^51 + k, and r, for x = k ln 2/256 + r.
#[inline(always)]
fn reduce<const FUSED: bool>(x: f64) -> (u64, f64) {
let s = x.mla::<FUSED>(256.0 * LOG2_E, f64::SHIFT);
let k = s - f64::SHIFT;
(s.to_bits() & ((1 << 52) - 1), k.mla::<FUSED>(-LN2_256.1, k.mla::<FUSED>(-LN2_256.0, x)))
}
/// e^r - 1 plus the relative rounding error of 2^(j/256) as `SBITS` carries it.
#[inline(always)]
fn series<const FUSED: bool>(j: u64, r: f64) -> f64 {
let r2 = r * r;
f64::from_bits(TAIL[j as usize]) + r2.mla::<FUSED>(C[0], r) + r * r2 * r.mla::<FUSED>(C[2], C[1])
}
/// The bits of 0x1.66p-1 as an f32. `log2_single` reduces x by a centre c of one of 16 subintervals from it to twice it.
const OFF_SINGLE: u32 = 0x3f330000;
/// 1/c for each of `log2_single`'s subintervals.
const INVC_SINGLE: [f64; 16] = [
1.398907162146528, 1.3403141896637998, 1.286432210124115, 1.2367150214269895, 1.1906977166711752, 1.1479821020556429, 1.1082251448272158,
1.0711297413057381, 1.036437278977283, 1.0, 0.9492859795739057, 0.8951049428609004, 0.8476821620351103, 0.8050314851692001, 0.7664671008843108,
0.731428603316328,
];
/// log2 c for each of `log2_single`'s subintervals.
const LOGC_SINGLE: [f64; 16] = [
-0.48430022186289673, -0.42257122959194704, -0.3633754347673556, -0.30651309567405577, -0.25180720160537634, -0.19910014943794563, -0.14825100623281615,
-0.09913323807318392, -0.051632812977629436, 0.0, 0.07508531937943004, 0.15987125980713107, 0.2384046664317681, 0.31288288605863257, 0.38370422656453185,
0.451211048935815,
];
/// The coefficients of log2(1 + r), highest degree first, for r = z/c - 1.
const A_SINGLE: [f64; 5] = [0.288457581109214, -0.36092606229713164, 0.480898481472577, -0.7213474675006291, 1.4426950408774342];
/// log2 x in f64, for the bits `ix` of a positive normal f32, within about 2^-33 relative. It's Arm Optimized Routines' scalar `powf` logarithm. Each item reads its own table entries, so only scalar code calls it. It reads its polynomial through `loaded`, which makes LLVM load each coefficient with one instruction. Otherwise LLVM builds each one on every call from four moves and a transfer between registers, which costs about a sixth of the time of `Pow`'s scalar form at f32.
#[inline(always)]
pub(crate) fn log2_single<const FUSED: bool>(ix: u32) -> f64 {
let tmp = ix.wrapping_sub(OFF_SINGLE);
let i = (tmp >> 19) as usize % 16;
let top = tmp & 0xff800000;
let z = f64::from(f32::from_bits(ix.wrapping_sub(top)));
let r = z.mla::<FUSED>(INVC_SINGLE[i], -1.0);
let y0 = LOGC_SINGLE[i] + f64::from(top as i32 >> 23);
let r2 = r * r;
let a = loaded(&A_SINGLE);
let y = r.mla::<FUSED>(a[0], a[1]);
let p = r.mla::<FUSED>(a[2], a[3]);
y.mla::<FUSED>(r2 * r2, p.mla::<FUSED>(r2, r.mla::<FUSED>(a[4], y0)))
}
/// The bits of 2^(i/32), less i << 47, for each i below 32. Adding k << 47 for k = 32 m + i gives the bits of 2^(k/32).
const EXP2_SINGLE: [u64; 32] = [
0x3ff0000000000000, 0x3fefd9b0d3158574, 0x3fefb5586cf9890f, 0x3fef9301d0125b51, 0x3fef72b83c7d517b, 0x3fef54873168b9aa, 0x3fef387a6e756238,
0x3fef1e9df51fdee1, 0x3fef06fe0a31b715, 0x3feef1a7373aa9cb, 0x3feedea64c123422, 0x3feece086061892d, 0x3feebfdad5362a27, 0x3feeb42b569d4f82,
0x3feeab07dd485429, 0x3feea47eb03a5585, 0x3feea09e667f3bcd, 0x3fee9f75e8ec5f74, 0x3feea11473eb0187, 0x3feea589994cce13, 0x3feeace5422aa0db,
0x3feeb737b0cdc5e5, 0x3feec49182a3f090, 0x3feed503b23e255d, 0x3feee89f995ad3ad, 0x3feeff76f2fb5e47, 0x3fef199bdd85529c, 0x3fef3720dcef9069,
0x3fef5818dcfba487, 0x3fef7c97337b9b5f, 0x3fefa4afa2a490da, 0x3fefd0765b6e4540,
];
/// The coefficients of 2^r - 1, highest degree first, for |r| up to 1/64.
const C_SINGLE: [f64; 3] = [0.05550361559341535, 0.2402284522445722, 0.6931471806916203];
/// 2^x in f64, for |x| below 126, where an f32 result is normal, within about 2^-34 relative. It's Arm Optimized Routines' scalar `powf` exponential. Each item reads its own table entry, so only scalar code calls it. It reads its polynomial as `log2_single` does.
#[inline(always)]
pub(crate) fn exp2_single<const FUSED: bool>(x: f64) -> f64 {
const SHIFT: f64 = 1.5 * 4503599627370496.0 / 32.0;
let s = x + SHIFT;
let r = x - (s - SHIFT);
let ki = s.to_bits();
let t = f64::from_bits(EXP2_SINGLE[ki as usize % 32].wrapping_add(ki << 47));
let c = loaded(&C_SINGLE);
let y = r.mla::<FUSED>(c[0], c[1]).mla::<FUSED>(r * r, r.mla::<FUSED>(c[2], 1.0));
y * t
}
/// The constants `c`, read from memory as loads that LLVM can't fold into the code.
#[inline(always)]
fn loaded<const N: usize>(c: &'static [f64; N]) -> [f64; N] {
// SAFETY: `c` refers to a valid, aligned array that lives for the whole program.
unsafe { std::ptr::read_volatile(c) }
}