API reference — numerics

Contents

API reference — numerics#

Device-legal scalar math, the working-type indirection expression trees compute through, and the two modules built on top of it: banded emulated double precision (for FP64-less GPUs) and compensated/atomic accumulation.

Namespace aether::math#

The device-legal scalar dispatch facade — the SAME call spelling works over a native float/double and, where certified, a banded operand (see Banded/emulated real arithmetic). fma and the banded-constrained sincos overload are part of this facade too, but their declarations carry a trailing C++20 requires-clause that this toolchain’s Doxygen (1.9.1) cannot parse into a clean signature — a documented upstream limitation (docs/doxygen_allowlist.txt), not a gap in their own docstrings.

The elementwise set, all aether::math::NAME on float/double, reachable from aether/aether.h. Device compiles call the CUDA math library function of the same C name; host compiles call std::. “vector” marks the functions that AETHER_HOST_VECTOR_MATH routes (for double, in a GCC x86-64 host translation unit) to the faithfully rounded packet math so scalar loops vectorise; the rest stay scalar on the host.

Group

Functions

Host vector route

exponential, logarithm

exp, exp2, exp10, expm1, log, log2, log10, log1p, pow

vector

roots

sqrt, rsqrt, cbrt, hypot

vector

trigonometric

sin, cos, sincos, tan, asin, acos, atan, atan2

vector

hyperbolic

sinh, cosh, tanh, asinh, acosh, atanh

vector

rounding

floor, ceil, trunc, round (half away from zero, as in C), rint (half to even, numpy’s round/rint)

vector

sign, selection

abs, copysign, fmax, fmin, min, fdim

vector / inline

other

sign, clip, fmod, remainder, fma, isnan, isinf, isfinite

scalar

special

erf, erfc

scalar

Where the semantics differ from C or numpy: remainder follows numpy (result has the sign of the divisor), not C’s IEEE remainder; fmod is C’s (sign of the dividend), as is numpy’s fmod. round is C’s; numpy’s round is rint. sign(NaN) is 0 (numpy returns NaN). clip propagates NaN from any argument, as numpy’s does.

Device against host: the rounding, sign/selection, fmod, remainder, clip, fma and classification functions give the same bits; the transcendental ones agree within the CUDA Math API’s documented error plus glibc’s (a few ULP).

Not provided yet: lgamma, tgamma, digamma, Bessel functions, and integer-specific operations.

template<class T>
inline T aether::math::sign(T x)#

Absolute value.

…on a banded operand: three XORs against the leading limb’s own sign bit, exact and total, returning the working carrier.

|a| with the sign of b. copysign is one of the five certified banded entry points, and a facade entry that exists for the emulated scalar but not for double would be a mode-divergence trap of exactly the kind this facade exists to prevent.

…on a banded operand: ONE XOR of two sign bits, then the mask.

Floor.

Maximum of two values (floating-point fmax semantics — a NaN operand is ignored when the other is not NaN).

…on a banded operand. The comparison is the certified sub’s leading limb, never a lexicographic limb compare (which is WRONG — see banded::detail::fmax’s own counterexample).

Minimum of two values (floating-point fmin semantics).

…on a banded operand.

Minimum of two values — alias entry point alongside

fmax/fmin; same dispatch as fmin.

See also

fmax.

…on a banded operand — the same certified fmin, under the alias spelling, so a consumer migrating by spelling swap keeps both names.

a raised to the power b.

…on a banded operand: pow = exp(b*log(a)) over BandExpLog.h’s ported exp/log, certified over banded::detail::bandPowAdmits

.

Square root.

See also

banded::detail::BandedFacade::pow.

Sign of x: +1 for x > 0, -1 for x < 0, 0 otherwise (x == +0, x == -0, or NaN — a NaN compares false to both > and < zero, so it falls through to 0 the same way IEEE zero does; documented rather than special-cased). No device/host split needed — pure comparison against T{0}, identical on both routes (unlike fma/isfinite above, which route to different device/host intrinsics).

template<class T>
inline bool aether::math::isfinite(T x)#

Whether x is finite (not Inf, not NaN). Supersedes aether/expr/Reduce.h’s local detail::isfinite_ shim with the real dispatch header.

template<class T>
inline bool aether::math::isnan(T x)#

Whether x is NaN.

template<class T>
inline bool aether::math::isinf(T x)#

Whether x is +inf or -inf.

template<class T>
inline T aether::math::clip(T x, T lo, T hi)#

x limited to [lo, hi], with numpy’s clip semantics: minimum(maximum(x, lo), hi) where maximum/minimum propagate NaN (unlike fmax/fmin), so a NaN in any argument gives NaN, and lo > hi gives hi. Pure compares, identical on host and device.

template<class T>
inline T aether::math::remainder(T a, T b)#

…floor on a banded operand: forwards to the ported BandRound.h::floor via BandedFacade.

Least integer >= x.

…ceil on a banded operand: -floor(-x), sign-bit flips only.

Nearest integer, half away from zero (C’s round, not rint’s half-to-even; numpy’s round is rint).

…round on a banded operand: Shewchuk grow-expansion + early-stop floor cascade, ported verbatim.

Nearest integer in the current rounding mode (half to even by default; numpy’s rint/round).

Round toward zero.

…trunc on a banded operand: copysign(floor(|x|), x).

Positive difference: a > b ? a - b : +0.

…fdim on a banded operand: forwards to Band.h’s own fdim.

Floating-point remainder of a/b, sign of the dividend.

…fmod on a banded operand: forwards to BandRound.h::fmod — exact on span <= 72 (bandFmodAdmits), ported verbatim.

Floored remainder of a/b with numpy’s remainder semantics: the result has the sign of the DIVISOR b (a - floor(a/b)*b, computed exactly from fmod), a zero result is copysign(0, b), and remainder(-1, inf) == inf as in numpy. NOT C’s IEEE remainder (round-to-nearest quotient); use fmod for the truncated (sign of the dividend) form. Scalar on the host.

template<class T>
inline T aether::math::sin(T x)#

Sine of x. Device FP64 path routes to detail::cwSin (Cody- Waite below 2^31, full-range Payne-Hanek above, NaN for Inf/NaN, STACK=0 — see detail/BoundedTrig.h); device FP32 routes to the native sinf; host routes to std::sin.

template<class T>
inline T aether::math::cos(T x)#

Cosine of x. Device FP64 path routes to detail::cwCos; device FP32 routes to the native cosf; host routes to std::cos.

template<class T>
inline T aether::math::rsqrt(T x)#

…sqrt on a banded operand: forwards to the certified RsqrtCore.h::sqrt_ via BandedFacade.

Reciprocal square root, 1/sqrt(x). Device routes to the hardware SFU intrinsic (rsqrtf/::rsqrt); host computes 1/sqrt(x) (under AETHER_HOST_VECTOR_MATH, the faithfully rounded packet rsqrt for double) — NOT claimed bit-identical between the routes: the seed itself differs by construction, exactly as the certified Band family’s own rsqrt is not.

template<class T>
inline T aether::math::rsqrtCube(T x)#

…rsqrt on a banded operand: forwards to BandRoot.h’s rsqrtIeee (the rsqrt(+Inf)=+0 fix) via BandedFacade.

x^(-3/2). No CUDA/std:: intrinsic exists; the scalar (non- Band) leg composes this file’s own rsqrt rather than duplicating a seed/refine schedule — unlike the certified Band overload (RsqrtCore.h::rsqrtCube), which fuses to avoid materializing the root, there is no certified-bound reason to fuse a plain float/double composition.

aether::WorkingType#

template<class T>
struct WorkingType#

The working (on-register) carrier for storage scalar T, plus the two conversions between them. PRIMARY = identity; specialize for a scalar whose arithmetic lives on a different type.

A specialization must provide the same three names:

  • using type = W; the working carrier

  • static type toWorking(const T&) storage → working (the region ENTRY)

  • static T fromWorking(const type&) working → storage (the PACK terminal)

Both conversions are AETHER_DEVICEHOST(). A specialization whose fromWorking is not constexpr is fine — every caller in expr/ is a template, so it is simply never constant-evaluated for that instantiation.

Banded emulated real (aether/banded/banded.h)#

Not included from the aether/aether.h umbrella — see Banded/emulated real arithmetic for why and for the three types’ roles.

struct Band#

Certified 3-slot signed-limb carrier. No sign field, no scale field, NO MEMORY FORM (that is BandedReal/BandCell8, BandCell8.h).

Public Functions

inline void toBits(std::uint32_t &hi, std::uint32_t &lo, std::uint32_t &tail) const#

Reinterpret this Band’s three limbs as raw IEEE-754 FP32 bit patterns — the inverse of fromBits, equally bare.

Public Members

float hi#

leading limb; carries the overall sign natively

float lo#

second limb, |lo| <= ulp(hi) after each certified op

float tail#

residue limb; what lifts the carrier past the df64 ceiling

Public Static Functions

static inline Band fromBits(std::uint32_t hi, std::uint32_t lo, std::uint32_t tail)#

Reinterpret three raw IEEE-754 FP32 bit patterns as a Band’s limbs — a bare per-limb bit_cast, no rounding, no normalization.

Public Static Attributes

static constexpr int certifiedBits = 53#

The certified effective width of this carrier. Declared EXPLICITLY rather than inherited from any default — a default that happens to agree would give the right answer for the wrong reason and would survive a later change to that default.

Friends

inline friend Band operator+(Band a, Band b)#

Banded add.

See also

detail::add.

inline friend Band operator-(Band a)#

Banded unary negate.

See also

detail::neg.

inline friend Band operator-(Band a, Band b)#

Banded sub.

See also

detail::sub.

inline friend Band operator*(Band a, Band b)#

Banded mul.

See also

detail::mul.

inline friend Band operator/(Band a, Band b)#

Banded div — ONE named primitive with its own body and its own battery, never the composition mul(a, recip(b)).

See also

detail::div.

struct BandCell8#

8-byte storage cell: one 64-bit word whose high half IS an IEEE-754 binary32 and whose low half is the next 32 significand bits.

The single member is deliberately a bare word rather than a bitfield pair: bitfield layout is implementation-defined, and every operation this type supports is a shift and a mask anyway.

Public Functions

inline Band band() const#

→ Band. The stored word’s three limbs, about six instructions.

A thin spelling of detail::bandFromCell8 (body supplied out-of-line at the bottom of this file, once that function is complete). Exact on the value: the word holds 56 significand bits and this redistributes them into the format’s canonical 24 + 23 + 9 split without rounding. The limb boundaries are the format’s, not those of whatever Band produced the word.

Public Members

std::uint64_t w#

the packed word: [ binary32 leading limb : 32 more bits ]

Public Static Functions

static inline BandCell8 fromTexel(int2 t)#

Reassemble a cell from a fetched int2 texel. The entire fetch-side decode — two shifts and an OR, no arithmetic.

The word order here is not the one the field names suggest: x is the low word and y the high one, because an int2 texel over eight little-endian bytes maps x -> bytes 0-3, and a uint64_t’s low half lives there. So x carries the unsigned 32-bit continuation and y the leading binary32 with the sign and exponent. Reading them the other way round does not perturb a value — it reinterprets a significand field as an exponent, which is wrong by an arbitrary power of two and frequently by a sign, and which symmetric test data cannot see. Checked by a value round-trip over asymmetric rows with a deliberately swapped arm shown to fail (BandCell8Cert.PunnedLoadStoreRoundTripsBitwise’s texel sibling).

Friends

inline friend bool operator==(BandCell8 a, BandCell8 b)#

Storage identity — the packed words are equal.

A bit comparison, deliberately not a numeric one. The tier-1 codec is canonical (one word per value), so on cells produced by cell8FromBand the two coincide. They part company on the words the format keeps distinct on purpose: +0 (all-zero) versus -0 (the sign bit alone) are unequal here, which is the right answer for “is this buffer still the fill value” and the wrong one for arithmetic. Compare decoded carriers, or BandedReals, when a numeric answer is wanted.

struct BandedReal#

The 8-byte banded storage scalar: a BandCell8 word with numeric semantics.

See also

the file header for what it deliberately lacks.

Public Functions

inline constexpr std::uint64_t toBits() const#

The raw 64-bit word.

inline constexpr BandCell8 cell() const#

The stored cell, undecoded.

inline Band band() const#

…the named spelling. Also the name a coefficient-table seam detects a storage codec by (requires { e.band(); } rather than naming any codec type), so a type without this member silently leaves the certified decode path.

Public Static Functions

static inline constexpr BandedReal fromCell(BandCell8 cell)#

Wrap a codec word. Explicit (a named factory, not a converting constructor): a BandCell8 may carry a per-array bias frame and this type is defined at bias 0, so the conversion is a claim, not a coincidence.

static inline constexpr BandedReal fromBits(std::uint64_t w)#

Construct from the raw 64-bit word — the only constexpr construction route (see the file header).

static inline BandedReal fromBand(Band b)#

…the named spelling of the same thing.

Friends

inline friend bool operator==(BandedReal a, BandedReal b)#

Numeric equality: -0 == +0 is true, NaN == NaN is false. (BandCell8::operator== answers the storage question and gives the opposite answer on both — that is the reason this type is not an alias.)

Namespace aether::accum#

Compensated-sum and atomic accumulation policies for reductions.

namespace accum#

Typedefs

template<class Real>
using DefaultPolicyT = typename DefaultPolicy<Real>::T#

Variables

constexpr offset_t npos = offset_max#

The “no slot” sentinel tryAppend() returns on overflow — the same bit pattern as aether::offset_max, an offset_t value no real capacity reaches.

class AppendCounter#
#include <AppendCounter.h>

Owning device-side append counter over a declared capacity. One uint32_t atomic slot counter (a host-pinned Chunk, plus a device Chunk when built with CUDA) and one DeviceFlag for overflow.

Public Functions

inline std::size_t capacity() const#

Declared capacity — the bound tryAppend() enforces.

inline AppendCounterView deviceView() const#

The view a kernel (or host code, in AETHER_CPP_MODE) takes by value.

inline std::size_t rawCount() const#

Non-destructive host-side read of the RAW counter (may exceed capacity() — every overflowing thread still incremented it before discovering the overflow). Downloads first (a no-op in AETHER_CPP_MODE).

inline bool overflowed() const#

Non-destructive peek at the overflow flag (does not clear it).

inline void checkOverflow(const std::string &context = "AppendCounter::checkOverflow")#

Throw if the overflow flag is set, and clear it. Call this after commit(), which never throws, so a rejected overflow does not discard the samples that did fit.

Throws:

aether::Error – if the flag was set.

inline void reset()#

Reset the raw counter to 0 for a fresh step. Does NOT touch the overflow flag.

template<class T, std::size_t... Es>
inline std::size_t commit(Array<T, Es...> &arr)#

Read the raw count back, clamp it to capacity(), adopt that many newly-written spare samples into arr via arr.commit(...), then reset this counter for the next step. Never throws; pair with overflowed()/checkOverflow() to surface an overflow.

Returns:

The clamped count actually committed.

class AppendCounterView#
#include <AppendCounter.h>

Non-owning, AETHER_DEVICEHOST()-safe view of one AppendCounter — what a kernel takes by value.

Public Functions

inline offset_t tryAppend() const#

Claim the next slot. Device: atomicAdd(counter, 1u). Host (incl. AETHER_CPP_MODE, where this same code path runs as plain host code): a relaxed std::atomic<uint32_t>::fetch_add.

Returns:

A slot < capacity (write it), or aether::accum::npos when the draw would overflow — capacity is UNCHANGED in that case (the counter itself may keep climbing past capacity across many overflowing threads; that raw value is what AppendCounter::rawCount() reports back, clamped by commit()), and the overflow flag is set. NEVER returns a slot >= capacity.

struct CompensatedAtomic#
#include <AccumPlane.h>

double plane: hardware atomicAdd + exact 2Sum residual in a companion lane, folded at delivery. Concurrency-safe, order-independent.

template<class Real>
struct DefaultPolicy#

The policy a plane uses when the caller does not name one. Primary template deliberately UNDEFINED.

template<>
struct DefaultPolicy<double>#
template<class Real, std::size_t VD, class Policy>
class LaneView#
#include <AccumPlane.h>

The per-sample accumulator handle returned by PlaneView::operator[]. Holds a reference to the plane view and the sample index only — no buffered state, so it may be created and discarded freely (acc[i] += workVec3 reads the way the kernels already read).

Public Functions

template<class E>
inline LaneView &operator+=(const Expression<E, typename E::element_type> &term)#

Accumulate a VD-component term into this sample’s lanes.

template<class E>
inline LaneView &accumulate(const Expression<E, typename E::element_type> &term)#

Named spelling of operator+=.

inline Item<Real, VD> delivered() const#

This sample’s delivered value (see PlaneView::delivered).

inline void zero()#

Zero this sample’s lanes, unconditionally.

template<class Real, std::size_t VD = 3, class Policy = DefaultPolicyT<Real>>
class Plane#
#include <AccumPlane.h>

Owning per-sample accumulation plane: VD lanes of double, plus the companion lane when the policy has one. Host-side type; hand deviceView() (CUDA) or hostView() to the code that accumulates.

Public Functions

inline explicit Plane(std::size_t n)#

Allocate a plane for n samples. Storage is UNINITIALIZED (matches Array’s own “raw allocate” convention) — call zeroHost()/zeroDeviceAsync() before first use.

inline std::size_t samples() const#

Number of samples.

inline ViewT hostView()#

Host-side view.

inline ViewT deviceView()#

Device-side view &#8212; what a kernel takes, by value. AETHER_CPP_MODE: the same single chunk as hostView() (mirrors Array::deviceView()).

inline void upload()#

Upload the host lanes to the device lanes. No-op in AETHER_CPP_MODE.

inline void download()#

Download the device lanes into the host lanes. No-op in AETHER_CPP_MODE.

inline void zeroHost()#

The host mirror of zeroDeviceAsync(): one unconditional pass over the whole host plane.

inline ArrayT &valueArray()#

The underlying value-lane Array (allocation, upload/download).

inline ArrayT &compArray()#

The underlying companion-lane Array; 0 samples unless hasCompanionLane.

Public Static Functions

static inline bool zeroIsAllZeroBytes()#

Is an all-zero byte pattern the policy’s zero element? (Always true here — double’s +0.0 is all-zero-bytes — checked at run time since the answer is about object representation.)

template<class Real, std::size_t VD, class Policy = DefaultPolicyT<Real>>
class PlaneView#
#include <AccumPlane.h>

Non-owning device-bindable view of an accumulation plane (public alias aether::AccumPlaneView), passed by value — the only way a kernel touches the plane.

Template Parameters:
  • Real – the scalar the plane delivers (only double is supported).

  • VD – lanes per sample (3 for a vec3 acceleration).

  • Policy – CompensatedAtomic or Serialized.

Public Functions

inline offset_t size() const#

Number of samples the plane covers.

inline LaneT operator[](const SampleIndex &i)#

plane[i] += term &#8212; the accumulation surface.

template<class E>
inline void accumulate(const SampleIndex &i, const Expression<E, typename E::element_type> &term)#

Accumulate a VD-component term into sample i’s lanes. The lanes are visited in ascending component order, so a host run and a device run accumulate in the same order and produce the same bits.

inline void zero(const SampleIndex &i)#

Zero sample i’s lanes &#8212; UNCONDITIONALLY, including the companion lane. No guard, no owner, no “first producer” flag.

inline ItemT delivered(const SampleIndex &i) const#

Sample i’s accumulated value, in Real. One fold for the whole step; returns a register-resident Item.

template<class OutE>
inline void deliver(OutE &out, const SampleIndex &i) const#

Terminal, written straight into a destination expression: plane.deliver(out, i) is out[i] = plane.delivered(i).

inline const StoreViewT &valueView() const#

The value-lane VIEW — diagnostics/tests, not part of the accumulation surface.

inline const StoreViewT &compView() const#

The companion-lane view; a 0-sample view unless hasCompanionLane.

Public Static Attributes

static constexpr std::size_t VecDims = VD#

Lanes per sample.

static constexpr bool concurrentSafe = TraitsT::concurrentSafe#

See also

PolicyTraits::concurrentSafe.

static constexpr bool hasCompanionLane = TraitsT::hasCompanionLane#

Whether a companion residual lane exists.

template<class Policy>
struct PolicyTraits#

Compile-time description of a policy. Primary template deliberately left UNDEFINED — an unknown policy tag is a hard compile error at the plane, never a silent default.

template<>
struct PolicyTraits<CompensatedAtomic>#
template<>
struct PolicyTraits<Serialized>#
struct Serialized#
#include <AccumPlane.h>

double plane: direct +=, no atomics, no companion lane. The caller owes serialization of the producers writing one sample.

namespace atomic#

Functions

template<class ViewT>
inline double minAbsReal(ViewT &h, const SampleIndex &i, double val)#

Atomic minimum on the magnitude of a signed real value.

Compares |val| against |h[i]|; if smaller, stores copysign(|val|, h[i]) so the sign already at the slot is preserved.

Template Parameters:

ViewT – a scalar (extents<dyn>) double view.

Returns:

The value previously stored at the slot.

inline double minAbsReal(double *h, double val)#

Raw-pointer overload for callers that already hold a raw address.

template<class ViewT>
inline void orBool(ViewT &h, const SampleIndex &i, bool val)#

Monotone OR-store of true into a one-byte bool slot.

Because the underlying operation is an OR with the constant true, the result is independent of the order concurrent writers execute in — a plain monotone store of 1 suffices on both CUDA and host paths. A false call is a no-op.

Template Parameters:

ViewT – a scalar bool view.

inline void orBool(bool *h, bool val)#

Raw-pointer overload for callers that already hold a raw address.

template<class ViewT>
inline int maxInt(ViewT &h, const SampleIndex &i, int val)#

Atomic maximum on a signed 32-bit integer slot.

Template Parameters:

ViewT – a scalar int view.

Returns:

The value previously stored at the slot.

inline int maxInt(int *h, int val)#

Raw-pointer overload for callers that already hold a raw address.

template<class ViewT>
inline void compensatedSum(ViewT &value, ViewT &comp, const SampleIndex &i, double term)#

Compensated accumulation of a double term into a value slot and its companion residual (comp) slot. Delivered quantity is value[i] + comp[i], independent of accumulation order.

View-only: no raw-pointer overload. aether/accum/AccumPlane.h is its only consumer.

Template Parameters:

ViewT – scalar double views (value lane and companion lane).

namespace detail#

Functions

inline double minAbsReal_(double *address, double val)#

Atomic minimum on the magnitude of a signed real value at address, preserving the sign already stored there.

inline void orBool_(bool *address, bool val)#

Monotone OR-store of true into a one-byte bool slot at address.

Invariant: the only operation is an idempotent set-to-1, performed as a byte-granular store (ST.U8), so it never touches a neighbouring byte and stays inside a one-byte allocation. The result is visible at the kernel boundary. A future need to clear a slot or read it back within a kernel must switch the slot type to 32 bits and use a word atomic.

inline double twoSumErr_(double a, double b, double s)#

Knuth two-sum error term: given s == fl(a + b), returns the residual (a + b) - s exactly, for any ordering of |a|, |b| and any signs (no fast2Sum precondition). Six FP64 flops, branch-free.

inline void compensatedSum_(double *value, double *comp, double term)#

Compensated accumulation of term into a (value, comp) slot pair. Delivered quantity is value[i] + comp[i], independent of accumulation order. Requires atomicAdd(double*, double) (sm_60+; aether’s lowest configured arch is sm_61).

inline int maxInt_(int *address, int val)#

Atomic maximum on a signed 32-bit integer at address.

namespace detail#

Functions

template<class T, std::size_t VD>
inline View<T, extents<dyn>, layout_stride> laneView_(const View<T, extents<VD, dyn>, layout_stride> &full, std::size_t d)#

Slice lane d (d < VD) out of the owning VD-lane view full (View<T, extents<VD,dyn>, layout_stride>, Array’s own view shape) as a scalar View<T, extents<dyn>, layout_stride>: one pointer offset by lane d’s own stride, with the sample-mode stride carried straight through.

template<class Real, class Policy>
struct PolicyOps#

The storage element and the three operations (zero, accumulate, deliver) of one (scalar, policy) pair. Primary template deliberately UNDEFINED.

template<>
struct PolicyOps<double, CompensatedAtomic>#
#include <AccumPlane.h>

double + compensated atomic. Two lanes; delivery folds them.

template<>
struct PolicyOps<double, Serialized>#
#include <AccumPlane.h>

double + direct serialized +=. One lane; delivery is a read.