1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
//! Sealed scalar-float trait used throughout `semiflow-core`.
//!
//! [`SemiflowFloat`] bundles the exact bounds required by every generic Chernoff
//! implementation. Only `f32` and `f64` implement the trait (explicit `impl`
//! blocks; no blanket impl) so the set of accepted scalars cannot grow without
//! a deliberate ADR. [`Dual<F>`](crate::dual::Dual) is the deliberate third
//! member, authorized by ADR-0133.
//!
//! See [`docs/adr/0025-generic-over-float.md`](../docs/adr/) and
//! `contracts/semiflow-core.math.md` for motivation and bound derivation.
use core::{
fmt::{Debug, Display},
ops::{AddAssign, DivAssign, MulAssign, SubAssign},
};
use num_traits::Float;
// ---------------------------------------------------------------------------
// Trait
// ---------------------------------------------------------------------------
/// Sealed scalar-float trait: the set of types accepted by generic `semiflow-core` types.
///
/// ## Sealed design (ADR-0025, v0.9.0 Block D pilot)
///
/// Only `f32` and `f64` implement `SemiflowFloat` — explicit `impl` blocks,
/// no blanket impl. This prevents accidental instantiation with `Complex<f64>`
/// or other numeric types before the relevant ADR lands.
///
/// ## Bound rationale
///
/// The supertrait list mirrors what every Chernoff kernel actually needs:
/// - [`num_traits::Float`] — `sqrt`, `abs`, `is_finite`, `floor`, `exp`, `ln`, …
/// - `AddAssign + SubAssign + MulAssign + DivAssign` — in-place BLAS operations in [`crate::State`].
/// - `Send + Sync + Copy + 'static` — grid and function types stored in structs,
/// shared across thread boundaries when `parallel` is enabled.
/// - `Debug + Display` — error messages include the offending value.
/// - `PartialOrd` — comparisons in validation helpers (`tau < 0`, `a(x) < 0`).
///
/// ## SIMD note
///
/// `f64` uses AVX2/NEON SIMD paths (Catmull-Rom, K-kernel) when the `simd`
/// feature is enabled. `f32` uses scalar-only paths; a dedicated `f32x8`
/// intrinsic path is deferred to a future ADR.
///
/// ## Example
///
/// ```rust
/// use semiflow::float::SemiflowFloat;
/// fn sum_two<F: SemiflowFloat>(a: F, b: F) -> F { a + b }
///
/// // Both concrete float types work:
/// assert_eq!(sum_two(1.0_f64, 2.0_f64), 3.0_f64);
/// assert_eq!(sum_two(1.0_f32, 2.0_f32), 3.0_f32);
/// ```
#[allow(clippy::module_name_repetitions)]
pub trait SemiflowFloat:
Float
+ AddAssign
+ SubAssign
+ MulAssign
+ DivAssign
+ Send
+ Sync
+ Copy
+ Debug
+ Display
+ PartialOrd
+ 'static
{
}
impl SemiflowFloat for f32 {}
impl SemiflowFloat for f64 {}
// ---------------------------------------------------------------------------
// Small numeric helpers used by generic kernels
// ---------------------------------------------------------------------------
/// Return the additive identity (`0.0`) for `F`.
// future generic helper; companion to `one` and `two` which are actively used
#[allow(dead_code)]
#[inline]
pub(crate) fn zero<F: SemiflowFloat>() -> F {
F::zero()
}
/// Return the multiplicative identity (`1.0`) for `F`.
#[inline]
pub(crate) fn one<F: SemiflowFloat>() -> F {
F::one()
}
/// Return `2.0` as `F`.
#[inline]
pub(crate) fn two<F: SemiflowFloat>() -> F {
let o = one::<F>();
o + o
}
/// Return `0.5` as `F`.
#[inline]
pub(crate) fn half<F: SemiflowFloat>() -> F {
one::<F>() / two::<F>()
}
/// Convert an `f64` literal to `F`.
///
/// Uses `num_traits::cast::ToPrimitive` + `from_f64`. Panics in debug if
/// the conversion fails; in release the unwrap degrades to zero (the
/// `num_traits` contract for out-of-range).
#[inline]
pub(crate) fn from_f64<F: SemiflowFloat>(v: f64) -> F {
F::from(v).unwrap_or_else(F::zero)
}
/// `π^{-D/2}` computed without `powf`, so the value is bit-identical everywhere.
///
/// `f64::powf` lowers to the platform `pow`, which is not correctly rounded and
/// whose glibc IFUNC variants differ by CPU. That matters here because this is a
/// GLOBAL normalisation multiplying every output of the d-D Gauss-Hermite
/// kernels: a 1-ULP difference in the prefactor moves every value, which is what
/// made `G_BINDING_SMOLYAK_PARITY`'s 0-ULP golden pass on five CI runners and
/// fail on the sixth (ADR-0191 AMENDMENT 5).
///
/// Every operation here is IEEE-754 correctly rounded — multiplication, `sqrt`
/// and division are all required to be — so the result is reproducible across
/// platforms and libm versions.
///
/// Reproducibility is the whole claim; accuracy is NOT. Measured against a
/// 60-digit reference on glibc x86-64: at `D = 2, 4, 6` this and `powf` agree
/// bit-for-bit, and at `D = 3, 5` they differ by 1 ULP with `powf` on the
/// correctly rounded side. A future reader tempted to call this "more accurate"
/// should re-measure per `D` per platform first — the point is that this one
/// gives the same bits everywhere, not that it gives better ones.
pub(crate) fn inv_pi_pow_half<F: SemiflowFloat>(d: usize) -> F {
let pi = core::f64::consts::PI;
let mut acc = 1.0_f64;
for _ in 0..(d / 2) {
acc *= pi;
}
if d % 2 == 1 {
acc *= libm::sqrt(pi);
}
from_f64::<F>(1.0 / acc)
}