Module Tolk.Transcendental

Transcendental functions as arithmetic.

A target without instructions for exp2, log2, sin or sqrt computes them from additions, multiplications, comparisons, selections, casts and bit reinterpretations. The functions below build those computations as nodes: polynomial approximations on a reduced argument, with the reduction and the reconstruction done on the bits of the floating-point representation.

exp2, log2 and sin are defined in Dtype.t.Float32 and Dtype.t.Float64. They return IEEE's special values at infinities, NaNs and zeros, and the sine of a zero is that zero, its sign kept. Elsewhere each result is within one unit in the last place of the correctly rounded one: at most one value of its type away from it, on the ordered line of the type's values. The other floats compute them in Dtype.t.Float32 and round once (patterns), which keeps them within one unit too.

Errors. A function given a node of another type than it is defined for raises Invalid_argument.

Bit manipulation

#val exponent_bias : Dtype.t -> int

exponent_bias dt is the bias of the exponent field of the float dt: 2^e-1 - 1 for e exponent bits (Dtype.finfo), and one more for Dtype.t.Fp8e4m3fnuz and Dtype.t.Fp8e5m2fnuz.

Raises Invalid_argument if dt is not in Dtype.floats.

#val shl : Ops.t -> int -> Ops.t

shl x n is x * 2ⁿ: x shifted left by n bits if x is an integer.

Raises Invalid_argument if n is negative.

#val shr : Ops.t -> int -> Ops.t

shr x n is x // 2ⁿ, the division rounding towards negative infinity: x shifted right by n bits if x is an integer.

Raises Invalid_argument if n is negative.

#val rintk : Ops.t -> Ops.t

rintk d is the float d rounded to the nearest integer, halves away from zero, as the signed integer of d's width: Dtype.t.Int16, Dtype.t.Int32 or Dtype.t.Int64.

#val pow2if : Ops.t -> Dtype.t -> Ops.t

pow2if q dt is 2^q as a float, for q an integer in the range of normal exponents of that float: Dtype.t.Float64 for an Dtype.t.Int64 q, Dtype.t.Float32 for an Dtype.t.Int32 q, and dt for an Dtype.t.Int16 q.

#val frexp : Ops.t -> Ops.t * Ops.t

frexp v is (m, e) with |v| = |m| * 2^e and 0.5 <= |m| < 1, for a normal v. m has v's type, and v's sign in Dtype.t.Float16 and Dtype.t.Float32; in Dtype.t.Float64 it is positive. e has the unsigned integer type of v's width, and wraps around for |v| < 0.5.

Argument reductions

A reduction brings an angle d >= 0 near zero, where a polynomial approximates the sine, and says which multiple of a half-turn or of a quarter-turn it removed.

#val payne_hanek_reduction : Ops.t -> Ops.t * Ops.t

payne_hanek_reduction d is (r, q) with d = q * pi/2 + r and |r| <= pi/4, for 1 <= d, infinities excluded. r has d's type and q is a Dtype.t.Int32 whose value modulo 4 is the quadrant of d. It multiplies d's whole mantissa by the bits of 1/(2pi) at its exponent, to 128 bits of the fraction, and rounds the remainder once from them, which keeps r within an ulp of the exact remainder whatever d's magnitude.

#val cody_waite_reduction : Ops.t -> Ops.t * Ops.t

cody_waite_reduction d is (r, q) with d = q * pi/2 + r and |r| about pi/4 at most, for |d| <= 39800. r has d's type and q is a Dtype.t.Int32 whose value modulo 4 is the quadrant of d. It subtracts q * pi/2 in several parts whose products are exact, which is precise only while q is small.

Functions

#val xsin : ?fast:bool -> ?switch_over:float -> Ops.t -> Ops.t

xsin ~fast ~switch_over d is the sine of d, NaN at infinities and NaN. Angles below switch_over (default 30., at most 39800.) in magnitude are reduced with cody_waite_reduction, and larger ones with payne_hanek_reduction, unless fast (default false), which assumes every angle is below switch_over and builds only the first. fast holds when d's bounds (Ops.vmin, Ops.vmax) lie below switch_over in magnitude.

#val xexp2 : Ops.t -> Ops.t

xexp2 d is 2^d: inf where it overflows d's type, 0 where it underflows past the subnormals, and NaN at NaN.

#val xlog2 : Ops.t -> Ops.t

xlog2 d is the base-2 logarithm of d, subnormals included: inf at inf, -inf at both zeros, and NaN below zero and at NaN.

#val xpow : Ops.t -> Ops.t -> Ops.t

xpow base exponent is base to the power exponent, computed as exp2 (exponent * log2 |base|). A negative base gives NaN for an exponent that is not an integer, except for -inf, and the result's sign follows the exponent's parity otherwise. Any base to the power 0 is 1, NaN and infinities included.

#val xpow_int : Ops.t -> Ops.t -> Dtype.t -> Ops.t option

xpow_int base exponent dt is base to the power exponent in the integer type dt, wrapping at its width as repeated products do. A negative exponent gives 0, but 1 for a base of 1, and -1 or 1 by the exponent's parity for -1 in a signed dt. It is None if dt is not in Dtype.ints.

Raises Invalid_argument if exponent is a weak expression whose bounds no 64-bit integer type holds (Ops.fits_no_int): compiled code would compute it wrapped, and a power reads every bit of its exponent, so a wrong exponent is a wrong value. A literal exponent is exact at any size.

Rewriting

#val patterns : force:bool -> Op.Set.t -> (unit, Ops.t) Ops.Pattern_matcher.t

patterns ~force ops rewrites each Op.t.Exp2, Op.t.Log2 and Op.t.Sin whose operation is not in ops, or every one if force, to xexp2, xlog2 and xsin. On the other floats, Dtype.t.Float16, Dtype.t.Bfloat16 and the 8-bit ones, the operation is computed in Dtype.t.Float32 and cast back. An Op.t.Sqrt under the same condition becomes xpow d 0.5.

Arguments equal to those of an earlier call give that call's matcher.