thermite_special/
polylog.rs

1//! The polylogarithm order.
2//!
3//! [`polylog`](crate::SpecialMath::polylog) takes its order as a [`PolylogOrder`]: a
4//! **scalar**, uniform across the packet, tagged by the class of order it carries.
5//! `$\mathrm{Li}_n$` at whole-number `n` is a table lookup: every coefficient of the unity
6//! series is a tabulated `$\zeta$` value, the leading `$\Gamma(1-s)(-\mu)^{s-1}$` term
7//! collapses into `$H_{n-1} - \ln(-\mu)$`, and the far
8//! field is the Bernoulli-polynomial inversion formula. At arbitrary real `s` the same
9//! series needs a sweep of live `$\zeta(s-k)$` evaluations, a real power, and Wood's
10//! m-th-root identity in the far field. Those are different algorithms with costs an
11//! order of magnitude apart.
12//!
13//! # Why a scalar, not a vector
14//!
15//! Every order-dependent quantity (`$\zeta(s-k)/k!$`, `$k^{-s}$`, `$\Gamma(1-s)$`, the
16//! near-integer brackets) is a _per-call_ scalar precompute, splatted once. A per-lane
17//! order would pay that sweep per lane, and the regime choices that depend on `s` (near an
18//! integer, non-positive) would become masks over arms every lane then has to evaluate.
19//! A caller with several orders runs several calls.
20//!
21//! # The payload types
22//!
23//! `PolylogOrder<E, S>` carries the invoking vector's own element types: `E` is its
24//! `Element` (`f64` on an `f64` vector, `Complex<f64>` on a complex one, a `Dual` element on
25//! a dual one) and `S` is its `Signed` lane element (`i64` on an `f64` vector, `i32` on an
26//! `f32` one), so an order is spelled in the arithmetic of the type it is used with and
27//! nothing is ever converted. The trait signature is
28//! `PolylogOrder<Self::Element, <Self::Signed as GenericVector>::Element>`.
29//!
30//! # Downgrading
31//!
32//! [`simplify`](PolylogOrder::simplify) narrows [`Real`](PolylogOrder::Real) to
33//! [`Integer`](PolylogOrder::Integer) when the value is _exactly_ whole. Nothing narrower
34//! is attempted: unlike the Bessel order there is no half-integer shortcut (Wood's
35//! half-integer series is the general one with tabulated constants). A near-integer
36//! order is a **correctness** hazard for the general arm rather than a cost choice, which
37//! is why the general kernel fuses the two cancelling poles algebraically instead of
38//! trusting the caller to have snapped.
39
40use thermite::element::{FloatElement, SignedIntegerElement};
41use thermite::math::scalar::Unwrap;
42use thermite::prelude::*;
43
44/// The order `$s$` of a polylogarithm, tagged with the class of order it carries. See the
45/// [module documentation](self).
46#[derive(Debug, Clone, Copy, PartialEq)]
47pub enum PolylogOrder<E, S> {
48    /// `$s = n$`, a whole number of either sign, in the vector's signed lane type. The
49    /// tabulated arm for `$n \ge 1$`, closed forms at `0` and `1`, the reflected series for
50    /// `$n < 0$`.
51    Integer(S),
52
53    /// Arbitrary real `$s$`, in the vector's element type. The general algorithm.
54    Real(E),
55}
56
57impl<E: FloatElement, S: SignedIntegerElement> PolylogOrder<E, S> {
58    /// Narrow [`Real`](Self::Real) to [`Integer`](Self::Integer) when the value is exactly
59    /// a whole number. Never changes the value of `$s$`.
60    ///
61    /// `V` is any real vector over `E` (the one the order is about to be used with) and
62    /// only supplies the float-to-integer lane conversion.
63    #[inline(always)]
64    pub fn simplify<V>(self) -> Self
65    where
66        V: FloatVector<Element = E>,
67        V::Signed: GenericVector<Element = S>,
68    {
69        match self {
70            Self::Integer(_) => self,
71            Self::Real(s) => {
72                let r = FloatElement::round(s);
73                if r == s {
74                    Self::Integer(V::splat(r).to_signed_integer().extract::<0>())
75                } else {
76                    self
77                }
78            }
79        }
80    }
81}
82
83/// The order is a scalar in both the vector and the scalar spelling of `polylog`, so it
84/// crosses the scalar surface unchanged.
85impl<E, S> Unwrap for PolylogOrder<E, S> {
86    type Unwrapped = Self;
87
88    #[inline(always)]
89    fn wrap(value: Self) -> Self {
90        value
91    }
92
93    #[inline(always)]
94    fn unwrap(self) -> Self {
95        self
96    }
97}
thermite_special/polylog.rs

thermite_special/
polylog.rs
