thermite_special/
bernoulli.rs

1//! The Bernoulli sequence as vectors: `$B_0, B_1, B_2, B_3, \ldots$`, splatted.
2//!
3//! This is the public face of the Bernoulli tables. The raw per-format tables live in
4//! [`tables::bernoulli`](crate::tables::bernoulli) and hold only the even-index numbers
5//! from `$B_2$` up, because those are the ones that need a table; [`BernoulliSequence`]
6//! puts the head and the zeros back so the whole sequence can be iterated.
7//!
8//! ```
9//! use thermite::prelude::*;
10//! use thermite_special::bernoulli::BernoulliMath;
11//!
12//! type V = Vector<f64>;
13//!
14//! let b: Vec<f64> = V::bernoulli_numbers(-0.5)   // B_1 = -1/2, the caller's choice
15//!     .take(7)
16//!     .map(|v| v.extract::<0>())
17//!     .collect();
18//!
19//! assert_eq!(b, [1.0, -0.5, 1.0 / 6.0, 0.0, -1.0 / 30.0, 0.0, 1.0 / 42.0]);
20//! ```
21//!
22//! # `$B_1$` is yours to pick
23//!
24//! It is the one Bernoulli number the two conventions disagree on (`$-1/2$` from the
25//! generating function `$x/(e^x - 1)$`, `$+1/2$` from `$x/(1 - e^{-x})$`), so the table
26//! does not carry it and this iterator takes it as an argument instead. Pass whichever
27//! your formula assumes. Nothing else in the sequence changes with the choice.
28//!
29//! # Where it stops
30//!
31//! After the last `$B_{2n}$` representable in the element type: `$B_{258}$` for `f64`
32//! (259 items) and `$B_{64}$` for `f32` (65 items). `$|B_{2n}|$` grows factorially, so
33//! there is nothing beyond it to yield. See the
34//! [table docs](crate::tables::bernoulli) for the boundary in full.
35//!
36//! The iterator is [`ExactSizeIterator`], so `len()` gives that count up front.
37
38use core::iter::FusedIterator;
39
40use crate::RealPrimalMath;
41
42pub use crate::tables::bernoulli::{BernoulliNumbers, bernoulli_b2n};
43
44/// The Bernoulli sequence `$B_0, B_1, B_2, \ldots$` as splatted vectors, including the
45/// zero-valued odd terms.
46///
47/// Built by [`BernoulliMath::bernoulli_numbers`] or [`BernoulliSequence::new`]. See the
48/// [module docs](self) for the `$B_1$` convention and where the sequence ends.
49#[derive(Debug, Clone, Copy)]
50pub struct BernoulliSequence<V: RealPrimalMath<Element: BernoulliNumbers>> {
51    /// The subscript of the next number to yield, so `idx` IS `n` in `$B_n$`, not an
52    /// index into the underlying table, which skips `$B_0$`, `$B_1$` and the odd zeros.
53    idx: usize,
54    /// The caller's `$B_1$`. Held for the whole iteration rather than consumed at step
55    /// two, which is what lets the state be a bare counter.
56    b1: V::Element,
57}
58
59impl<V: RealPrimalMath<Element: BernoulliNumbers>> BernoulliSequence<V> {
60    /// Starts at `$B_0$`, with `b1` as the value of `$B_1$`.
61    #[inline]
62    #[must_use]
63    pub const fn new(b1: V::Element) -> Self {
64        Self { idx: 0, b1 }
65    }
66}
67
68impl<V: RealPrimalMath<Element: BernoulliNumbers>> Iterator for BernoulliSequence<V> {
69    type Item = V;
70
71    #[inline]
72    fn next(&mut self) -> Option<V> {
73        let table = <V::Element as BernoulliNumbers>::B2N;
74
75        let n = self.idx;
76        if n > 2 * table.len() {
77            // Deliberately does NOT advance, so the iterator is fused and `idx` cannot
78            // run away past the end.
79            return None;
80        }
81        self.idx = n + 1;
82
83        Some(match n {
84            0 => V::ONE,
85            1 => V::splat(self.b1),
86            // Every odd Bernoulli number past B_1 is zero. The table skips them, which is
87            // why `idx` is the subscript and the table index is derived, never the reverse.
88            _ if n % 2 == 1 => V::ZERO,
89            _ => V::splat(table[n / 2 - 1]),
90        })
91    }
92
93    #[inline]
94    fn size_hint(&self) -> (usize, Option<usize>) {
95        let end = 2 * <V::Element as BernoulliNumbers>::B2N.len() + 1;
96        let remaining = end.saturating_sub(self.idx);
97        (remaining, Some(remaining))
98    }
99}
100
101impl<V: RealPrimalMath<Element: BernoulliNumbers>> ExactSizeIterator for BernoulliSequence<V> {}
102
103impl<V: RealPrimalMath<Element: BernoulliNumbers>> FusedIterator for BernoulliSequence<V> {}
104
105/// Bernoulli numbers for real primal vectors.
106///
107/// Blanket-implemented for every [`RealPrimalMath`] vector whose element carries the
108/// tables, so bringing the trait into scope is all that is needed.
109pub trait BernoulliMath: RealPrimalMath<Element: BernoulliNumbers> + Sized {
110    /// The sequence `$B_0, B_1, B_2, \ldots$` as splatted vectors, zeros included.
111    ///
112    /// `b1` is the value of `$B_1$`, which the tables deliberately do not choose for you.
113    /// See the [module docs](self).
114    #[inline]
115    #[must_use]
116    fn bernoulli_numbers(b1: Self::Element) -> BernoulliSequence<Self> {
117        BernoulliSequence::new(b1)
118    }
119}
120
121impl<V: RealPrimalMath<Element: BernoulliNumbers>> BernoulliMath for V {}
thermite_special/bernoulli.rs

thermite_special/
bernoulli.rs
