From a751dbde16d7fe3c9139945e0d424d32782b0993 Mon Sep 17 00:00:00 2001 From: Auberer Date: Mon, 15 Jun 2026 17:48:07 +0200 Subject: [PATCH 1/3] Extend std/math with classification, special functions and RNG Add the missing elementary and numeric helpers to the math module: - Float classification: isNan, isInf, isFinite, signbit, copysign, inf(), nan() - Double-returning rounding: truncD, floorD, ceilD, roundD; plus fmod, remainder, modf, frexp, ldexp - Exp/log extras: exp2, expm1, log1p - Inverse hyperbolic: asinh, acosh, atanh - Special functions: tgamma, lgamma (Lanczos), erf, erfc (Abramowitz-Stegun) - Utility: signum, lerp, fma, isqrt, array overloads of max/min - Constants: TAU, PHI, INV_PI, INV_SQRT_2 Replace the C rand() dependency with a deterministic, seedable xorshift64* generator (seed, randULong, randDouble, randDouble(min,max), randBool, randGaussian) and reimplement randInt on top of it, making randomness portable and reproducible. Cover all additions with assert-based reference tests under test/test-files/std/math. Co-Authored-By: Claude Opus 4.8 --- std/math/const.spice | 8 +- std/math/fct.spice | 495 ++++++++++++++++++ std/math/rand.spice | 95 +++- test/test-files/std/math/const-extra/cout.out | 1 + .../std/math/const-extra/source.spice | 9 + .../std/math/fct-classification/cout.out | 1 + .../std/math/fct-classification/source.spice | 29 + .../std/math/fct-exp-log-extra/cout.out | 1 + .../std/math/fct-exp-log-extra/source.spice | 23 + .../std/math/fct-inverse-hyperbolic/cout.out | 1 + .../math/fct-inverse-hyperbolic/source.spice | 25 + .../std/math/fct-rounding-double/cout.out | 1 + .../std/math/fct-rounding-double/source.spice | 62 +++ test/test-files/std/math/fct-special/cout.out | 1 + .../std/math/fct-special/source.spice | 30 ++ test/test-files/std/math/fct-utility/cout.out | 1 + .../std/math/fct-utility/source.spice | 42 ++ test/test-files/std/math/rand-extra/cout.out | 1 + .../std/math/rand-extra/source.spice | 59 +++ 19 files changed, 881 insertions(+), 4 deletions(-) create mode 100644 test/test-files/std/math/const-extra/cout.out create mode 100644 test/test-files/std/math/const-extra/source.spice create mode 100644 test/test-files/std/math/fct-classification/cout.out create mode 100644 test/test-files/std/math/fct-classification/source.spice create mode 100644 test/test-files/std/math/fct-exp-log-extra/cout.out create mode 100644 test/test-files/std/math/fct-exp-log-extra/source.spice create mode 100644 test/test-files/std/math/fct-inverse-hyperbolic/cout.out create mode 100644 test/test-files/std/math/fct-inverse-hyperbolic/source.spice create mode 100644 test/test-files/std/math/fct-rounding-double/cout.out create mode 100644 test/test-files/std/math/fct-rounding-double/source.spice create mode 100644 test/test-files/std/math/fct-special/cout.out create mode 100644 test/test-files/std/math/fct-special/source.spice create mode 100644 test/test-files/std/math/fct-utility/cout.out create mode 100644 test/test-files/std/math/fct-utility/source.spice create mode 100644 test/test-files/std/math/rand-extra/cout.out create mode 100644 test/test-files/std/math/rand-extra/source.spice diff --git a/std/math/const.spice b/std/math/const.spice index bb9dd432e..9226b3758 100644 --- a/std/math/const.spice +++ b/std/math/const.spice @@ -13,4 +13,10 @@ public const double PI_4 = 0.78539816339744830962; public const double SQRT_2 = 1.41421356237309504880; public const double SQRT_E = 1.64872127070012814685; public const double SQRT_PI = 1.7724538509055160273; -public const double SQRT_PHI = 1.27201964951406896425; \ No newline at end of file +public const double SQRT_PHI = 1.27201964951406896425; + +public const double TAU = 6.28318530717958647692; // 2 * PI +public const double PHI = 1.61803398874989484820; // Golden ratio + +public const double INV_PI = 0.31830988618379067154; // 1 / PI +public const double INV_SQRT_2 = 0.70710678118654752440; // 1 / SQRT_2 \ No newline at end of file diff --git a/std/math/fct.spice b/std/math/fct.spice index c7fc0bad3..9ed655d22 100644 --- a/std/math/fct.spice +++ b/std/math/fct.spice @@ -558,3 +558,498 @@ public inline f almostEqual(double input1, double input2) { public inline f almostEqual(double input1, double input2, double epsilon) { return abs(input1 - input2) <= epsilon; } + +/** + * Return the sign of the input as -1, 0 or 1 (of the same type as the input). + * + * @param x Input number + * @return -1 if x is negative, 1 if x is positive, 0 otherwise + */ +public inline f signum(Numeric x) { + if x > cast(0) { return cast(1); } + if x < cast(0) { return cast(-1); } + return cast(0); +} + +/** + * Linearly interpolate between a and b by the factor t. + * Note: t is not clamped, so values outside of [0, 1] extrapolate. + * + * @param a Start value (t = 0) + * @param b End value (t = 1) + * @param t Interpolation factor + * @return Interpolated value + */ +public inline f lerp(double a, double b, double t) { + return a + (b - a) * t; +} + +/** + * Multiply-add convenience helper, computed as a * b + c. + * Note: This is not a hardware fused multiply-add and does not guarantee single rounding. + * + * @param a First factor + * @param b Second factor + * @param c Addend + * @return a * b + c + */ +public inline f fma(double a, double b, double c) { + return a * b + c; +} + +/** + * Calculate the maximum value of a non-empty array. + * + * @param values Array of values + * @param size Number of elements in the array + * @return Maximum value of the array + */ +public f max(Numeric[] values, unsigned long size) { + Numeric res = values[0]; + for unsigned long i = 1ul; i < size; i++ { + if values[i] > res { res = values[i]; } + } + return res; +} + +/** + * Calculate the minimum value of a non-empty array. + * + * @param values Array of values + * @param size Number of elements in the array + * @return Minimum value of the array + */ +public f min(Numeric[] values, unsigned long size) { + Numeric res = values[0]; + for unsigned long i = 1ul; i < size; i++ { + if values[i] < res { res = values[i]; } + } + return res; +} + +/** + * Calculate the integer square root (largest whole number r with r * r <= x). + * Note: For negative inputs 0 is returned. + * + * @param x Input number + * @return Integer square root of input + */ +public f isqrt(WholeNumber x) { + if x <= cast(0) { return cast(0); } + WholeNumber lo = cast(1); + WholeNumber hi = x; + WholeNumber res = cast(1); + while lo <= hi { + WholeNumber mid = lo + (hi - lo) / cast(2); + // mid * mid <= x, rewritten as mid <= x / mid to avoid overflow + if mid <= x / mid { + res = mid; + lo = mid + cast(1); + } else { + hi = mid - cast(1); + } + } + return res; +} + +// ------------------------------------------------- Double rounding -------------------------------------------------- + +/** + * Truncate towards zero, returning a double. + * Unlike trunc, this also works for magnitudes outside of the int range. + * + * @param x Input number + * @return Truncated value as a double + */ +public f truncD(double x) { + // Values with a magnitude of 2^63 or more are already integral + if abs(x) >= 9223372036854775808.0 { return x; } + return cast(cast(x)); +} + +/** + * Round down to the nearest whole number, returning a double. + * + * @param x Input number + * @return Largest whole number <= x as a double + */ +public f floorD(double x) { + double t = truncD(x); + return (x < 0.0 && t != x) ? t - 1.0 : t; +} + +/** + * Round up to the nearest whole number, returning a double. + * + * @param x Input number + * @return Smallest whole number >= x as a double + */ +public f ceilD(double x) { + double t = truncD(x); + return (x > 0.0 && t != x) ? t + 1.0 : t; +} + +/** + * Round to the nearest whole number (halves away from zero), returning a double. + * + * @param x Input number + * @return Rounded value as a double + */ +public f roundD(double x) { + return x >= 0.0 ? floorD(x + 0.5) : ceilD(x - 0.5); +} + +/** + * Split the input into its integral and fractional parts. + * The integral part (truncated towards zero) is stored in intPart. + * + * @param x Input number + * @param intPart Output reference receiving the integral part + * @return Fractional part of the input (same sign as the input) + */ +public f modf(double x, double& intPart) { + intPart = truncD(x); + return x - intPart; +} + +/** + * Calculate the floating point remainder of x / y (with the sign of x). + * Note: For y == 0 NaN is returned. + * + * @param x Dividend + * @param y Divisor + * @return Remainder of x / y + */ +public f fmod(double x, double y) { + if y == 0.0 { return nan(); } + return x - truncD(x / y) * y; +} + +/** + * Calculate the IEEE remainder of x / y, where the quotient is rounded to the nearest integer. + * Note: For y == 0 NaN is returned. + * + * @param x Dividend + * @param y Divisor + * @return IEEE remainder of x / y + */ +public f remainder(double x, double y) { + if y == 0.0 { return nan(); } + return x - roundD(x / y) * y; +} + +/** + * Decompose the input into a normalized mantissa in [0.5, 1) and a power of two, + * such that x == mantissa * 2^exponent. The exponent is stored in exponent. + * + * @param x Input number + * @param exponent Output reference receiving the exponent + * @return Normalized mantissa + */ +public f frexp(double x, int& exponent) { + exponent = 0; + if x == 0.0 || isNan(x) || isInf(x) { return x; } + double m = abs(x); + while m >= 1.0 { + m /= 2.0; + exponent++; + } + while m < 0.5 { + m *= 2.0; + exponent--; + } + return x < 0.0 ? -m : m; +} + +/** + * Multiply the input by an integer power of two, computing x * 2^exponent. + * + * @param x Input number + * @param exponent Power of two + * @return x * 2^exponent + */ +public f ldexp(double x, int exponent) { + double res = x; + int e = exponent; + while e > 0 { + res *= 2.0; + e--; + } + while e < 0 { + res /= 2.0; + e++; + } + return res; +} + +// ----------------------------------------------- Float classification ----------------------------------------------- + +/** + * Construct a positive infinity value. + * + * @return Positive infinity + */ +public f inf() { + result = 0.0; + unsigned long bits = 0x7ff0000000000000ul; + unsafe { + sCopyUnsafe(cast(&bits), cast(&result), sizeof(bits)); + } +} + +/** + * Construct a quiet not-a-number (NaN) value. + * + * @return NaN + */ +public f nan() { + result = 0.0; + unsigned long bits = 0x7ff8000000000000ul; + unsafe { + sCopyUnsafe(cast(&bits), cast(&result), sizeof(bits)); + } +} + +/** + * Check whether the input is a not-a-number (NaN) value. + * + * @param x Input number + * @return true if x is NaN + */ +public inline f isNan(double x) { + // NaN is the only value that does not compare equal to itself + return x != x; +} + +/** + * Check whether the input is positive or negative infinity. + * + * @param x Input number + * @return true if x is infinite + */ +public inline f isInf(double x) { + // Infinity is the only non-zero value that is unchanged by doubling + return x != 0.0 && x * 2.0 == x; +} + +/** + * Check whether the input is a finite number (neither NaN nor infinite). + * + * @param x Input number + * @return true if x is finite + */ +public inline f isFinite(double x) { + return !isNan(x) && !isInf(x); +} + +/** + * Check whether the sign bit of the input is set (true for negative values and -0.0). + * + * @param x Input number + * @return true if the sign bit is set + */ +public f signbit(double x) { + unsigned long bits = 0ul; + unsafe { + sCopyUnsafe(cast(&x), cast(&bits), sizeof(x)); + } + return (bits >> 63ul) != 0ul; +} + +/** + * Copy the sign of one value onto the magnitude of another. + * + * @param magnitude Value providing the magnitude + * @param sign Value providing the sign + * @return A value with the magnitude of magnitude and the sign of sign + */ +public inline f copysign(double magnitude, double sign) { + return signbit(sign) ? -abs(magnitude) : abs(magnitude); +} + +// ------------------------------------------------- Exp / log extras ------------------------------------------------- + +/** + * Calculate 2 raised to the power of the input. + * + * @param x Input number + * @return 2 raised to the power of input + */ +public inline f exp2(double x) { + return exp(x * LN2); +} + +/** + * Calculate e^x - 1, accurately even for inputs close to zero. + * + * @param x Input number + * @return e raised to the power of input, minus one + */ +public f expm1(double x) { + if abs(x) >= 0.5 { return exp(x) - 1.0; } + // Series: e^x - 1 = x + x^2/2! + x^3/3! ... + double res = 0.0; + double term = x; + double n = 1.0; + while res + term != res { + res += term; + n += 1.0; + term *= x / n; + } + return res; +} + +/** + * Calculate the natural logarithm of (1 + x), accurately even for inputs close to zero. + * Note: For inputs <= -1 the result 0.0 is returned, as the result is not a real number. + * + * @param x Input number + * @return Natural logarithm of (1 + x) + */ +public f log1p(double x) { + if x <= -1.0 { return 0.0; } + if abs(x) >= 0.5 { return ln(1.0 + x); } + // Series: ln(1 + x) = x - x^2/2 + x^3/3 - ... + double res = 0.0; + double term = x; + double n = 1.0; + double sign = 1.0; + double inc = term / n; + while res + inc != res { + res += inc; + term *= x; + n += 1.0; + sign = -sign; + inc = sign * term / n; + } + return res; +} + +// --------------------------------------------- Inverse hyperbolic functions --------------------------------------------- + +/** + * Calculate the inverse hyperbolic sine of the input. + * + * @param x Input number + * @return Inverse hyperbolic sine of input + */ +public inline f asinh(double x) { + return x < 0.0 ? -ln(-x + sqrt(x * x + 1.0)) : ln(x + sqrt(x * x + 1.0)); +} + +/** + * Calculate the inverse hyperbolic cosine of the input. + * Note: For inputs < 1 the result 0.0 is returned, as the result is not a real number. + * + * @param x Input number (>= 1) + * @return Inverse hyperbolic cosine of input + */ +public inline f acosh(double x) { + if x < 1.0 { return 0.0; } + return ln(x + sqrt(x * x - 1.0)); +} + +/** + * Calculate the inverse hyperbolic tangent of the input. + * Note: For inputs outside of (-1, 1) the result 0.0 is returned, as the result is not a real number. + * + * @param x Input number in range (-1, 1) + * @return Inverse hyperbolic tangent of input + */ +public inline f atanh(double x) { + if x <= -1.0 || x >= 1.0 { return 0.0; } + return 0.5 * ln((1.0 + x) / (1.0 - x)); +} + +// ------------------------------------------------- Special functions ------------------------------------------------- + +/** + * Calculate the gamma function of the input, using the Lanczos approximation. + * + * @param x Input number + * @return Gamma function of input + */ +public f tgamma(double x) { + // Reflection formula for the left half plane: Gamma(x) = PI / (sin(PI*x) * Gamma(1-x)) + if x < 0.5 { + // sin expects degrees, so PI*x radians correspond to 180*x degrees + return PI / (sin(180.0 * x) * tgamma(1.0 - x)); + } + double[9] c = [ + 0.99999999999980993, + 676.5203681218851, + -1259.1392167224028, + 771.32342877765313, + -176.61502916214059, + 12.507343278686905, + -0.13857109526572012, + 9.9843695780195716e-6, + 1.5056327351493116e-7 + ]; + double xx = x - 1.0; + double a = c[0]; + double t = xx + 7.5; // g (= 7) + 0.5 + for int i = 1; i < 9; i++ { + a += c[i] / (xx + cast(i)); + } + // sqrt(2*PI) = SQRT_2 * SQRT_PI + return SQRT_2 * SQRT_PI * pow(t, xx + 0.5) * exp(-t) * a; +} + +/** + * Calculate the natural logarithm of the absolute value of the gamma function, + * using the Lanczos approximation in logarithmic space (avoids overflow for large inputs). + * + * @param x Input number + * @return Natural logarithm of the absolute value of the gamma function of input + */ +public f lgamma(double x) { + // Reflection: ln|Gamma(x)| = ln(PI) - ln|sin(PI*x)| - ln|Gamma(1-x)| + if x < 0.5 { + return ln(PI) - ln(abs(sin(180.0 * x))) - lgamma(1.0 - x); + } + double[9] c = [ + 0.99999999999980993, + 676.5203681218851, + -1259.1392167224028, + 771.32342877765313, + -176.61502916214059, + 12.507343278686905, + -0.13857109526572012, + 9.9843695780195716e-6, + 1.5056327351493116e-7 + ]; + double xx = x - 1.0; + double a = c[0]; + double t = xx + 7.5; + for int i = 1; i < 9; i++ { + a += c[i] / (xx + cast(i)); + } + // 0.5 * ln(2*PI) = ln(SQRT_2 * SQRT_PI) + return ln(SQRT_2 * SQRT_PI) + (xx + 0.5) * ln(t) - t + ln(a); +} + +/** + * Calculate the error function of the input, using the Abramowitz and Stegun approximation (7.1.26). + * Note: The maximum absolute error of the approximation is about 1.5e-7. + * + * @param x Input number + * @return Error function of input + */ +public f erf(double x) { + bool negative = x < 0.0; + double t = 1.0 / (1.0 + 0.3275911 * abs(x)); + double poly = ((((1.061405429 * t - 1.453152027) * t + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t; + double y = 1.0 - poly * exp(-x * x); + return negative ? -y : y; +} + +/** + * Calculate the complementary error function of the input (1 - erf(x)). + * + * @param x Input number + * @return Complementary error function of input + */ +public inline f erfc(double x) { + return 1.0 - erf(x); +} diff --git a/std/math/rand.spice b/std/math/rand.spice index 6618fc856..2fb30a68c 100644 --- a/std/math/rand.spice +++ b/std/math/rand.spice @@ -1,4 +1,42 @@ -ext f rand(); +// Imports +import "std/math/fct"; + +// Internal generator state (xorshift64*). Seeded with a fixed non-zero default so that +// an unseeded program still produces a well-defined, reproducible sequence. +unsigned long randState = 0x2545f4914f6cdd1dul; + +/** + * Advance the internal xorshift64* generator and return the next raw 64-bit value. + * + * @return Next raw pseudo-random value + */ +f nextRaw() { + unsigned long x = randState; + x ^= x >> 12ul; + x ^= x << 25ul; + x ^= x >> 27ul; + randState = x; + return x * 0x2545f4914f6cdd1dul; +} + +/** + * Seed the pseudo-random number generator. + * Note: A seed of 0 is replaced by a fixed non-zero constant, as the generator cannot escape the all-zero state. + * + * @param seedValue Seed value + */ +public p seed(unsigned long seedValue) { + randState = seedValue == 0ul ? 0x9e3779b97f4a7c15ul : seedValue; +} + +/** + * Generate a raw pseudo-random unsigned long, uniformly distributed over the whole 64-bit range. + * + * @return Random unsigned long + */ +public f randULong() { + return nextRaw(); +} /** * Generates a random integer between min and max @@ -9,5 +47,56 @@ ext f rand(); */ public f randInt(int min, int max) { if min > max { return max; } - return rand() % (max - min + 1) + min; -} \ No newline at end of file + unsigned long range = cast(max - min) + 1ul; + return min + cast(nextRaw() % range); +} + +/** + * Generate a random double, uniformly distributed in the half-open interval [0, 1). + * + * @return Random double in [0, 1) + */ +public f randDouble() { + // Use the upper 53 bits (the mantissa width of a double) and scale by 2^-53 + return cast(nextRaw() >> 11ul) * (1.0 / 9007199254740992.0); +} + +/** + * Generate a random double, uniformly distributed in the half-open interval [min, max). + * Note: If min > max, the function always returns max + * + * @param min Lower bound (inclusive) + * @param max Upper bound (exclusive) + * @return Random double in [min, max) + */ +public f randDouble(double min, double max) { + if min > max { return max; } + return min + randDouble() * (max - min); +} + +/** + * Generate a random boolean value with equal probability for true and false. + * + * @return Random boolean + */ +public f randBool() { + return (nextRaw() & 1ul) == 1ul; +} + +/** + * Generate a normally distributed random double with the given mean and standard deviation, + * using the Box-Muller transform. + * + * @param mean Mean of the distribution + * @param stdDev Standard deviation of the distribution + * @return Normally distributed random double + */ +public f randGaussian(double mean, double stdDev) { + double u1 = randDouble(); + double u2 = randDouble(); + // Guard against ln(0), as randDouble() can return exactly 0.0 + if u1 < 1.0e-300 { u1 = 1.0e-300; } + double magnitude = stdDev * sqrt(-2.0 * ln(u1)); + // cos expects degrees, so 2*PI*u2 radians correspond to 360*u2 degrees + return mean + magnitude * cos(360.0 * u2); +} diff --git a/test/test-files/std/math/const-extra/cout.out b/test/test-files/std/math/const-extra/cout.out new file mode 100644 index 000000000..cebe8db0e --- /dev/null +++ b/test/test-files/std/math/const-extra/cout.out @@ -0,0 +1 @@ +All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/const-extra/source.spice b/test/test-files/std/math/const-extra/source.spice new file mode 100644 index 000000000..67dc73210 --- /dev/null +++ b/test/test-files/std/math/const-extra/source.spice @@ -0,0 +1,9 @@ +import "std/math/const"; + +f main() { + assert TAU > 6.283185 && TAU < 6.283186; + assert PHI > 1.618033 && PHI < 1.618034; + assert INV_PI > 0.318309 && INV_PI < 0.318310; + assert INV_SQRT_2 > 0.707106 && INV_SQRT_2 < 0.707107; + printf("All assertions passed!"); +} diff --git a/test/test-files/std/math/fct-classification/cout.out b/test/test-files/std/math/fct-classification/cout.out new file mode 100644 index 000000000..cebe8db0e --- /dev/null +++ b/test/test-files/std/math/fct-classification/cout.out @@ -0,0 +1 @@ +All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-classification/source.spice b/test/test-files/std/math/fct-classification/source.spice new file mode 100644 index 000000000..f0cf2071a --- /dev/null +++ b/test/test-files/std/math/fct-classification/source.spice @@ -0,0 +1,29 @@ +import "std/math/fct"; + +f main() { + // inf / nan construction and detection + assert isInf(inf()); + assert isInf(-inf()); + assert !isFinite(inf()); + assert isNan(nan()); + assert !isNan(inf()); + assert !isInf(nan()); + assert isFinite(0.0); + assert isFinite(-123.456); + assert !isFinite(nan()); + assert !isNan(1.0); + assert !isInf(1.0e300); // large but finite + + // signbit + assert signbit(-1.0); + assert !signbit(1.0); + assert signbit(-inf()); + assert !signbit(inf()); + + // copysign + assert copysign(3.0, -2.0) == -3.0; + assert copysign(-3.0, 2.0) == 3.0; + assert copysign(3.0, 5.0) == 3.0; + + printf("All assertions passed!"); +} diff --git a/test/test-files/std/math/fct-exp-log-extra/cout.out b/test/test-files/std/math/fct-exp-log-extra/cout.out new file mode 100644 index 000000000..cebe8db0e --- /dev/null +++ b/test/test-files/std/math/fct-exp-log-extra/cout.out @@ -0,0 +1 @@ +All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-exp-log-extra/source.spice b/test/test-files/std/math/fct-exp-log-extra/source.spice new file mode 100644 index 000000000..88cf236ec --- /dev/null +++ b/test/test-files/std/math/fct-exp-log-extra/source.spice @@ -0,0 +1,23 @@ +import "std/math/fct"; + +f main() { + // exp2 + assert almostEqual(exp2(0.0), 1.0); + assert almostEqual(exp2(10.0), 1024.0, 1.0e-3); + assert almostEqual(exp2(-1.0), 0.5); + + // expm1 (accurate near zero) + assert almostEqual(expm1(0.0), 0.0); + assert almostEqual(expm1(1.0), 1.718281828); + assert almostEqual(expm1(0.0001), 0.0001000050); + assert almostEqual(expm1(-0.5), -0.393469340); + + // log1p (accurate near zero) + assert almostEqual(log1p(0.0), 0.0); + assert almostEqual(log1p(1.0), 0.693147180); + assert almostEqual(log1p(0.0001), 0.0000999950); + assert log1p(-1.0) == 0.0; // out of domain convention + assert log1p(-2.0) == 0.0; + + printf("All assertions passed!"); +} diff --git a/test/test-files/std/math/fct-inverse-hyperbolic/cout.out b/test/test-files/std/math/fct-inverse-hyperbolic/cout.out new file mode 100644 index 000000000..cebe8db0e --- /dev/null +++ b/test/test-files/std/math/fct-inverse-hyperbolic/cout.out @@ -0,0 +1 @@ +All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-inverse-hyperbolic/source.spice b/test/test-files/std/math/fct-inverse-hyperbolic/source.spice new file mode 100644 index 000000000..94cf05402 --- /dev/null +++ b/test/test-files/std/math/fct-inverse-hyperbolic/source.spice @@ -0,0 +1,25 @@ +import "std/math/fct"; + +f main() { + // asinh + assert almostEqual(asinh(0.0), 0.0); + assert almostEqual(asinh(1.0), 0.881373587); + assert almostEqual(asinh(-1.0), -0.881373587); + assert almostEqual(asinh(2.0), 1.443635475); + + // acosh + assert almostEqual(acosh(1.0), 0.0); + assert almostEqual(acosh(2.0), 1.316957897); + assert almostEqual(acosh(10.0), 2.993222846); + assert acosh(0.5) == 0.0; // out of domain convention + + // atanh + assert almostEqual(atanh(0.0), 0.0); + assert almostEqual(atanh(0.5), 0.549306144); + assert almostEqual(atanh(-0.5), -0.549306144); + assert atanh(1.0) == 0.0; // out of domain convention + assert atanh(-1.0) == 0.0; + assert atanh(2.0) == 0.0; + + printf("All assertions passed!"); +} diff --git a/test/test-files/std/math/fct-rounding-double/cout.out b/test/test-files/std/math/fct-rounding-double/cout.out new file mode 100644 index 000000000..cebe8db0e --- /dev/null +++ b/test/test-files/std/math/fct-rounding-double/cout.out @@ -0,0 +1 @@ +All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-rounding-double/source.spice b/test/test-files/std/math/fct-rounding-double/source.spice new file mode 100644 index 000000000..d8e503934 --- /dev/null +++ b/test/test-files/std/math/fct-rounding-double/source.spice @@ -0,0 +1,62 @@ +import "std/math/fct"; + +f main() { + // truncD + assert truncD(1.9) == 1.0; + assert truncD(-1.9) == -1.0; + assert truncD(5.0) == 5.0; + assert truncD(1.0e20) == 1.0e20; // already integral, outside long range + + // floorD + assert floorD(1.9) == 1.0; + assert floorD(-1.1) == -2.0; + assert floorD(2.0) == 2.0; + assert floorD(-2.0) == -2.0; + + // ceilD + assert ceilD(1.1) == 2.0; + assert ceilD(-1.9) == -1.0; + assert ceilD(2.0) == 2.0; + + // roundD + assert roundD(1.4) == 1.0; + assert roundD(1.5) == 2.0; + assert roundD(-1.5) == -2.0; + assert roundD(-1.4) == -1.0; + + // modf + double ip = 0.0; + double frac = modf(3.75, ip); + assert ip == 3.0; + assert almostEqual(frac, 0.75); + double ip2 = 0.0; + double frac2 = modf(-2.25, ip2); + assert ip2 == -2.0; + assert almostEqual(frac2, -0.25); + + // fmod + assert almostEqual(fmod(5.3, 2.0), 1.3); + assert almostEqual(fmod(-5.3, 2.0), -1.3); + assert isNan(fmod(1.0, 0.0)); + + // remainder + assert almostEqual(remainder(5.3, 2.0), -0.7); + assert almostEqual(remainder(5.0, 3.0), -1.0); + + // frexp + int e = 0; + double m = frexp(8.0, e); + assert e == 4; + assert almostEqual(m, 0.5); + int e2 = 0; + double m2 = frexp(0.0, e2); + assert e2 == 0; + assert m2 == 0.0; + + // ldexp + assert almostEqual(ldexp(0.5, 4), 8.0); + assert almostEqual(ldexp(1.5, 2), 6.0); + assert almostEqual(ldexp(3.0, -1), 1.5); + + printf("All assertions passed!"); +} diff --git a/test/test-files/std/math/fct-special/cout.out b/test/test-files/std/math/fct-special/cout.out new file mode 100644 index 000000000..cebe8db0e --- /dev/null +++ b/test/test-files/std/math/fct-special/cout.out @@ -0,0 +1 @@ +All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-special/source.spice b/test/test-files/std/math/fct-special/source.spice new file mode 100644 index 000000000..ea6a1e6af --- /dev/null +++ b/test/test-files/std/math/fct-special/source.spice @@ -0,0 +1,30 @@ +import "std/math/fct"; + +f main() { + // tgamma + assert almostEqual(tgamma(1.0), 1.0, 1.0e-6); + assert almostEqual(tgamma(2.0), 1.0, 1.0e-6); + assert almostEqual(tgamma(3.0), 2.0, 1.0e-6); + assert almostEqual(tgamma(4.0), 6.0, 1.0e-5); + assert almostEqual(tgamma(5.0), 24.0, 1.0e-4); + assert almostEqual(tgamma(0.5), 1.772453851, 1.0e-6); // sqrt(PI) + assert almostEqual(tgamma(1.5), 0.886226925, 1.0e-6); // sqrt(PI) / 2 + + // lgamma + assert almostEqual(lgamma(1.0), 0.0, 1.0e-6); + assert almostEqual(lgamma(2.0), 0.0, 1.0e-6); + assert almostEqual(lgamma(5.0), 3.178053830, 1.0e-5); // ln(24) + assert almostEqual(lgamma(0.5), 0.572364942, 1.0e-6); // ln(sqrt(PI)) + + // erf + assert almostEqual(erf(0.0), 0.0, 1.0e-6); + assert almostEqual(erf(1.0), 0.842700793, 1.0e-6); + assert almostEqual(erf(-1.0), -0.842700793, 1.0e-6); + assert almostEqual(erf(0.5), 0.520499878, 1.0e-6); + + // erfc + assert almostEqual(erfc(0.0), 1.0, 1.0e-6); + assert almostEqual(erfc(1.0), 0.157299207, 1.0e-6); + + printf("All assertions passed!"); +} diff --git a/test/test-files/std/math/fct-utility/cout.out b/test/test-files/std/math/fct-utility/cout.out new file mode 100644 index 000000000..cebe8db0e --- /dev/null +++ b/test/test-files/std/math/fct-utility/cout.out @@ -0,0 +1 @@ +All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-utility/source.spice b/test/test-files/std/math/fct-utility/source.spice new file mode 100644 index 000000000..421a5e98f --- /dev/null +++ b/test/test-files/std/math/fct-utility/source.spice @@ -0,0 +1,42 @@ +import "std/math/fct"; + +f main() { + // signum + assert signum(5) == 1; + assert signum(-5) == -1; + assert signum(0) == 0; + assert signum(-3.5) == -1.0; + assert signum(3.5) == 1.0; + assert signum(0.0) == 0.0; + assert signum(-7l) == -1l; + + // fma + assert almostEqual(fma(2.0, 3.0, 4.0), 10.0); + assert almostEqual(fma(-2.0, 3.0, 1.0), -5.0); + + // lerp + assert almostEqual(lerp(0.0, 10.0, 0.0), 0.0); + assert almostEqual(lerp(0.0, 10.0, 1.0), 10.0); + assert almostEqual(lerp(0.0, 10.0, 0.5), 5.0); + assert almostEqual(lerp(2.0, 4.0, 0.25), 2.5); + + // isqrt + assert isqrt(0) == 0; + assert isqrt(1) == 1; + assert isqrt(16) == 4; + assert isqrt(17) == 4; + assert isqrt(99) == 9; + assert isqrt(100) == 10; + assert isqrt(-5) == 0; + assert isqrt(1000000000000l) == 1000000l; + + // array max / min + int[6] ints = [ 3, 7, 1, 9, 2, 5 ]; + assert max(ints, len(ints)) == 9; + assert min(ints, len(ints)) == 1; + double[4] doubles = [ 2.5, -1.5, 3.25, 0.0 ]; + assert almostEqual(max(doubles, len(doubles)), 3.25); + assert almostEqual(min(doubles, len(doubles)), -1.5); + + printf("All assertions passed!"); +} diff --git a/test/test-files/std/math/rand-extra/cout.out b/test/test-files/std/math/rand-extra/cout.out new file mode 100644 index 000000000..cebe8db0e --- /dev/null +++ b/test/test-files/std/math/rand-extra/cout.out @@ -0,0 +1 @@ +All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/rand-extra/source.spice b/test/test-files/std/math/rand-extra/source.spice new file mode 100644 index 000000000..dfcad044b --- /dev/null +++ b/test/test-files/std/math/rand-extra/source.spice @@ -0,0 +1,59 @@ +import "std/math/rand"; +import "std/math/fct"; + +f main() { + // Reproducibility: the same seed yields the same sequence + seed(42ul); + double a1 = randDouble(); + double a2 = randDouble(); + seed(42ul); + assert randDouble() == a1; + assert randDouble() == a2; + + // randDouble range [0, 1) + seed(123ul); + for int i = 0; i < 1000; i++ { + double r = randDouble(); + assert r >= 0.0 && r < 1.0; + } + + // randDouble(min, max) range + for int i = 0; i < 1000; i++ { + double r = randDouble(-5.0, 5.0); + assert r >= -5.0 && r < 5.0; + } + assert randDouble(5.0, 1.0) == 1.0; // min > max convention + + // randInt range [min, max] + for int i = 0; i < 1000; i++ { + int n = randInt(1, 6); + assert n >= 1 && n <= 6; + } + assert randInt(10, 1) == 1; // min > max convention + + // randULong reproducibility + seed(7ul); + unsigned long u = randULong(); + seed(7ul); + assert randULong() == u; + + // randBool produces both values across many draws + seed(999ul); + bool sawTrue = false; + bool sawFalse = false; + for int i = 0; i < 100; i++ { + if randBool() { + sawTrue = true; + } else { + sawFalse = true; + } + } + assert sawTrue && sawFalse; + + // randGaussian returns finite values + for int i = 0; i < 1000; i++ { + assert isFinite(randGaussian(0.0, 1.0)); + } + + printf("All assertions passed!"); +} From 2a99d655669f467e89c1c1434cfb76c1c5e6a3fb Mon Sep 17 00:00:00 2001 From: Auberer Date: Mon, 15 Jun 2026 18:07:50 +0200 Subject: [PATCH 2/3] Drop fct.spice additions and their tests; keep constants and RNG Revert all std/math/fct.spice additions and remove their reference tests, along with the const.spice test. Keep the new math constants and the seedable RNG. Rework the randGaussian check to not depend on the removed isFinite helper. Co-Authored-By: Claude Opus 4.8 --- std/math/fct.spice | 495 ------------------ test/test-files/std/math/const-extra/cout.out | 1 - .../std/math/const-extra/source.spice | 9 - .../std/math/fct-classification/cout.out | 1 - .../std/math/fct-classification/source.spice | 29 - .../std/math/fct-exp-log-extra/cout.out | 1 - .../std/math/fct-exp-log-extra/source.spice | 23 - .../std/math/fct-inverse-hyperbolic/cout.out | 1 - .../math/fct-inverse-hyperbolic/source.spice | 25 - .../std/math/fct-rounding-double/cout.out | 1 - .../std/math/fct-rounding-double/source.spice | 62 --- test/test-files/std/math/fct-special/cout.out | 1 - .../std/math/fct-special/source.spice | 30 -- test/test-files/std/math/fct-utility/cout.out | 1 - .../std/math/fct-utility/source.spice | 42 -- .../std/math/rand-extra/source.spice | 12 +- 16 files changed, 7 insertions(+), 727 deletions(-) delete mode 100644 test/test-files/std/math/const-extra/cout.out delete mode 100644 test/test-files/std/math/const-extra/source.spice delete mode 100644 test/test-files/std/math/fct-classification/cout.out delete mode 100644 test/test-files/std/math/fct-classification/source.spice delete mode 100644 test/test-files/std/math/fct-exp-log-extra/cout.out delete mode 100644 test/test-files/std/math/fct-exp-log-extra/source.spice delete mode 100644 test/test-files/std/math/fct-inverse-hyperbolic/cout.out delete mode 100644 test/test-files/std/math/fct-inverse-hyperbolic/source.spice delete mode 100644 test/test-files/std/math/fct-rounding-double/cout.out delete mode 100644 test/test-files/std/math/fct-rounding-double/source.spice delete mode 100644 test/test-files/std/math/fct-special/cout.out delete mode 100644 test/test-files/std/math/fct-special/source.spice delete mode 100644 test/test-files/std/math/fct-utility/cout.out delete mode 100644 test/test-files/std/math/fct-utility/source.spice diff --git a/std/math/fct.spice b/std/math/fct.spice index 9ed655d22..c7fc0bad3 100644 --- a/std/math/fct.spice +++ b/std/math/fct.spice @@ -558,498 +558,3 @@ public inline f almostEqual(double input1, double input2) { public inline f almostEqual(double input1, double input2, double epsilon) { return abs(input1 - input2) <= epsilon; } - -/** - * Return the sign of the input as -1, 0 or 1 (of the same type as the input). - * - * @param x Input number - * @return -1 if x is negative, 1 if x is positive, 0 otherwise - */ -public inline f signum(Numeric x) { - if x > cast(0) { return cast(1); } - if x < cast(0) { return cast(-1); } - return cast(0); -} - -/** - * Linearly interpolate between a and b by the factor t. - * Note: t is not clamped, so values outside of [0, 1] extrapolate. - * - * @param a Start value (t = 0) - * @param b End value (t = 1) - * @param t Interpolation factor - * @return Interpolated value - */ -public inline f lerp(double a, double b, double t) { - return a + (b - a) * t; -} - -/** - * Multiply-add convenience helper, computed as a * b + c. - * Note: This is not a hardware fused multiply-add and does not guarantee single rounding. - * - * @param a First factor - * @param b Second factor - * @param c Addend - * @return a * b + c - */ -public inline f fma(double a, double b, double c) { - return a * b + c; -} - -/** - * Calculate the maximum value of a non-empty array. - * - * @param values Array of values - * @param size Number of elements in the array - * @return Maximum value of the array - */ -public f max(Numeric[] values, unsigned long size) { - Numeric res = values[0]; - for unsigned long i = 1ul; i < size; i++ { - if values[i] > res { res = values[i]; } - } - return res; -} - -/** - * Calculate the minimum value of a non-empty array. - * - * @param values Array of values - * @param size Number of elements in the array - * @return Minimum value of the array - */ -public f min(Numeric[] values, unsigned long size) { - Numeric res = values[0]; - for unsigned long i = 1ul; i < size; i++ { - if values[i] < res { res = values[i]; } - } - return res; -} - -/** - * Calculate the integer square root (largest whole number r with r * r <= x). - * Note: For negative inputs 0 is returned. - * - * @param x Input number - * @return Integer square root of input - */ -public f isqrt(WholeNumber x) { - if x <= cast(0) { return cast(0); } - WholeNumber lo = cast(1); - WholeNumber hi = x; - WholeNumber res = cast(1); - while lo <= hi { - WholeNumber mid = lo + (hi - lo) / cast(2); - // mid * mid <= x, rewritten as mid <= x / mid to avoid overflow - if mid <= x / mid { - res = mid; - lo = mid + cast(1); - } else { - hi = mid - cast(1); - } - } - return res; -} - -// ------------------------------------------------- Double rounding -------------------------------------------------- - -/** - * Truncate towards zero, returning a double. - * Unlike trunc, this also works for magnitudes outside of the int range. - * - * @param x Input number - * @return Truncated value as a double - */ -public f truncD(double x) { - // Values with a magnitude of 2^63 or more are already integral - if abs(x) >= 9223372036854775808.0 { return x; } - return cast(cast(x)); -} - -/** - * Round down to the nearest whole number, returning a double. - * - * @param x Input number - * @return Largest whole number <= x as a double - */ -public f floorD(double x) { - double t = truncD(x); - return (x < 0.0 && t != x) ? t - 1.0 : t; -} - -/** - * Round up to the nearest whole number, returning a double. - * - * @param x Input number - * @return Smallest whole number >= x as a double - */ -public f ceilD(double x) { - double t = truncD(x); - return (x > 0.0 && t != x) ? t + 1.0 : t; -} - -/** - * Round to the nearest whole number (halves away from zero), returning a double. - * - * @param x Input number - * @return Rounded value as a double - */ -public f roundD(double x) { - return x >= 0.0 ? floorD(x + 0.5) : ceilD(x - 0.5); -} - -/** - * Split the input into its integral and fractional parts. - * The integral part (truncated towards zero) is stored in intPart. - * - * @param x Input number - * @param intPart Output reference receiving the integral part - * @return Fractional part of the input (same sign as the input) - */ -public f modf(double x, double& intPart) { - intPart = truncD(x); - return x - intPart; -} - -/** - * Calculate the floating point remainder of x / y (with the sign of x). - * Note: For y == 0 NaN is returned. - * - * @param x Dividend - * @param y Divisor - * @return Remainder of x / y - */ -public f fmod(double x, double y) { - if y == 0.0 { return nan(); } - return x - truncD(x / y) * y; -} - -/** - * Calculate the IEEE remainder of x / y, where the quotient is rounded to the nearest integer. - * Note: For y == 0 NaN is returned. - * - * @param x Dividend - * @param y Divisor - * @return IEEE remainder of x / y - */ -public f remainder(double x, double y) { - if y == 0.0 { return nan(); } - return x - roundD(x / y) * y; -} - -/** - * Decompose the input into a normalized mantissa in [0.5, 1) and a power of two, - * such that x == mantissa * 2^exponent. The exponent is stored in exponent. - * - * @param x Input number - * @param exponent Output reference receiving the exponent - * @return Normalized mantissa - */ -public f frexp(double x, int& exponent) { - exponent = 0; - if x == 0.0 || isNan(x) || isInf(x) { return x; } - double m = abs(x); - while m >= 1.0 { - m /= 2.0; - exponent++; - } - while m < 0.5 { - m *= 2.0; - exponent--; - } - return x < 0.0 ? -m : m; -} - -/** - * Multiply the input by an integer power of two, computing x * 2^exponent. - * - * @param x Input number - * @param exponent Power of two - * @return x * 2^exponent - */ -public f ldexp(double x, int exponent) { - double res = x; - int e = exponent; - while e > 0 { - res *= 2.0; - e--; - } - while e < 0 { - res /= 2.0; - e++; - } - return res; -} - -// ----------------------------------------------- Float classification ----------------------------------------------- - -/** - * Construct a positive infinity value. - * - * @return Positive infinity - */ -public f inf() { - result = 0.0; - unsigned long bits = 0x7ff0000000000000ul; - unsafe { - sCopyUnsafe(cast(&bits), cast(&result), sizeof(bits)); - } -} - -/** - * Construct a quiet not-a-number (NaN) value. - * - * @return NaN - */ -public f nan() { - result = 0.0; - unsigned long bits = 0x7ff8000000000000ul; - unsafe { - sCopyUnsafe(cast(&bits), cast(&result), sizeof(bits)); - } -} - -/** - * Check whether the input is a not-a-number (NaN) value. - * - * @param x Input number - * @return true if x is NaN - */ -public inline f isNan(double x) { - // NaN is the only value that does not compare equal to itself - return x != x; -} - -/** - * Check whether the input is positive or negative infinity. - * - * @param x Input number - * @return true if x is infinite - */ -public inline f isInf(double x) { - // Infinity is the only non-zero value that is unchanged by doubling - return x != 0.0 && x * 2.0 == x; -} - -/** - * Check whether the input is a finite number (neither NaN nor infinite). - * - * @param x Input number - * @return true if x is finite - */ -public inline f isFinite(double x) { - return !isNan(x) && !isInf(x); -} - -/** - * Check whether the sign bit of the input is set (true for negative values and -0.0). - * - * @param x Input number - * @return true if the sign bit is set - */ -public f signbit(double x) { - unsigned long bits = 0ul; - unsafe { - sCopyUnsafe(cast(&x), cast(&bits), sizeof(x)); - } - return (bits >> 63ul) != 0ul; -} - -/** - * Copy the sign of one value onto the magnitude of another. - * - * @param magnitude Value providing the magnitude - * @param sign Value providing the sign - * @return A value with the magnitude of magnitude and the sign of sign - */ -public inline f copysign(double magnitude, double sign) { - return signbit(sign) ? -abs(magnitude) : abs(magnitude); -} - -// ------------------------------------------------- Exp / log extras ------------------------------------------------- - -/** - * Calculate 2 raised to the power of the input. - * - * @param x Input number - * @return 2 raised to the power of input - */ -public inline f exp2(double x) { - return exp(x * LN2); -} - -/** - * Calculate e^x - 1, accurately even for inputs close to zero. - * - * @param x Input number - * @return e raised to the power of input, minus one - */ -public f expm1(double x) { - if abs(x) >= 0.5 { return exp(x) - 1.0; } - // Series: e^x - 1 = x + x^2/2! + x^3/3! ... - double res = 0.0; - double term = x; - double n = 1.0; - while res + term != res { - res += term; - n += 1.0; - term *= x / n; - } - return res; -} - -/** - * Calculate the natural logarithm of (1 + x), accurately even for inputs close to zero. - * Note: For inputs <= -1 the result 0.0 is returned, as the result is not a real number. - * - * @param x Input number - * @return Natural logarithm of (1 + x) - */ -public f log1p(double x) { - if x <= -1.0 { return 0.0; } - if abs(x) >= 0.5 { return ln(1.0 + x); } - // Series: ln(1 + x) = x - x^2/2 + x^3/3 - ... - double res = 0.0; - double term = x; - double n = 1.0; - double sign = 1.0; - double inc = term / n; - while res + inc != res { - res += inc; - term *= x; - n += 1.0; - sign = -sign; - inc = sign * term / n; - } - return res; -} - -// --------------------------------------------- Inverse hyperbolic functions --------------------------------------------- - -/** - * Calculate the inverse hyperbolic sine of the input. - * - * @param x Input number - * @return Inverse hyperbolic sine of input - */ -public inline f asinh(double x) { - return x < 0.0 ? -ln(-x + sqrt(x * x + 1.0)) : ln(x + sqrt(x * x + 1.0)); -} - -/** - * Calculate the inverse hyperbolic cosine of the input. - * Note: For inputs < 1 the result 0.0 is returned, as the result is not a real number. - * - * @param x Input number (>= 1) - * @return Inverse hyperbolic cosine of input - */ -public inline f acosh(double x) { - if x < 1.0 { return 0.0; } - return ln(x + sqrt(x * x - 1.0)); -} - -/** - * Calculate the inverse hyperbolic tangent of the input. - * Note: For inputs outside of (-1, 1) the result 0.0 is returned, as the result is not a real number. - * - * @param x Input number in range (-1, 1) - * @return Inverse hyperbolic tangent of input - */ -public inline f atanh(double x) { - if x <= -1.0 || x >= 1.0 { return 0.0; } - return 0.5 * ln((1.0 + x) / (1.0 - x)); -} - -// ------------------------------------------------- Special functions ------------------------------------------------- - -/** - * Calculate the gamma function of the input, using the Lanczos approximation. - * - * @param x Input number - * @return Gamma function of input - */ -public f tgamma(double x) { - // Reflection formula for the left half plane: Gamma(x) = PI / (sin(PI*x) * Gamma(1-x)) - if x < 0.5 { - // sin expects degrees, so PI*x radians correspond to 180*x degrees - return PI / (sin(180.0 * x) * tgamma(1.0 - x)); - } - double[9] c = [ - 0.99999999999980993, - 676.5203681218851, - -1259.1392167224028, - 771.32342877765313, - -176.61502916214059, - 12.507343278686905, - -0.13857109526572012, - 9.9843695780195716e-6, - 1.5056327351493116e-7 - ]; - double xx = x - 1.0; - double a = c[0]; - double t = xx + 7.5; // g (= 7) + 0.5 - for int i = 1; i < 9; i++ { - a += c[i] / (xx + cast(i)); - } - // sqrt(2*PI) = SQRT_2 * SQRT_PI - return SQRT_2 * SQRT_PI * pow(t, xx + 0.5) * exp(-t) * a; -} - -/** - * Calculate the natural logarithm of the absolute value of the gamma function, - * using the Lanczos approximation in logarithmic space (avoids overflow for large inputs). - * - * @param x Input number - * @return Natural logarithm of the absolute value of the gamma function of input - */ -public f lgamma(double x) { - // Reflection: ln|Gamma(x)| = ln(PI) - ln|sin(PI*x)| - ln|Gamma(1-x)| - if x < 0.5 { - return ln(PI) - ln(abs(sin(180.0 * x))) - lgamma(1.0 - x); - } - double[9] c = [ - 0.99999999999980993, - 676.5203681218851, - -1259.1392167224028, - 771.32342877765313, - -176.61502916214059, - 12.507343278686905, - -0.13857109526572012, - 9.9843695780195716e-6, - 1.5056327351493116e-7 - ]; - double xx = x - 1.0; - double a = c[0]; - double t = xx + 7.5; - for int i = 1; i < 9; i++ { - a += c[i] / (xx + cast(i)); - } - // 0.5 * ln(2*PI) = ln(SQRT_2 * SQRT_PI) - return ln(SQRT_2 * SQRT_PI) + (xx + 0.5) * ln(t) - t + ln(a); -} - -/** - * Calculate the error function of the input, using the Abramowitz and Stegun approximation (7.1.26). - * Note: The maximum absolute error of the approximation is about 1.5e-7. - * - * @param x Input number - * @return Error function of input - */ -public f erf(double x) { - bool negative = x < 0.0; - double t = 1.0 / (1.0 + 0.3275911 * abs(x)); - double poly = ((((1.061405429 * t - 1.453152027) * t + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t; - double y = 1.0 - poly * exp(-x * x); - return negative ? -y : y; -} - -/** - * Calculate the complementary error function of the input (1 - erf(x)). - * - * @param x Input number - * @return Complementary error function of input - */ -public inline f erfc(double x) { - return 1.0 - erf(x); -} diff --git a/test/test-files/std/math/const-extra/cout.out b/test/test-files/std/math/const-extra/cout.out deleted file mode 100644 index cebe8db0e..000000000 --- a/test/test-files/std/math/const-extra/cout.out +++ /dev/null @@ -1 +0,0 @@ -All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/const-extra/source.spice b/test/test-files/std/math/const-extra/source.spice deleted file mode 100644 index 67dc73210..000000000 --- a/test/test-files/std/math/const-extra/source.spice +++ /dev/null @@ -1,9 +0,0 @@ -import "std/math/const"; - -f main() { - assert TAU > 6.283185 && TAU < 6.283186; - assert PHI > 1.618033 && PHI < 1.618034; - assert INV_PI > 0.318309 && INV_PI < 0.318310; - assert INV_SQRT_2 > 0.707106 && INV_SQRT_2 < 0.707107; - printf("All assertions passed!"); -} diff --git a/test/test-files/std/math/fct-classification/cout.out b/test/test-files/std/math/fct-classification/cout.out deleted file mode 100644 index cebe8db0e..000000000 --- a/test/test-files/std/math/fct-classification/cout.out +++ /dev/null @@ -1 +0,0 @@ -All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-classification/source.spice b/test/test-files/std/math/fct-classification/source.spice deleted file mode 100644 index f0cf2071a..000000000 --- a/test/test-files/std/math/fct-classification/source.spice +++ /dev/null @@ -1,29 +0,0 @@ -import "std/math/fct"; - -f main() { - // inf / nan construction and detection - assert isInf(inf()); - assert isInf(-inf()); - assert !isFinite(inf()); - assert isNan(nan()); - assert !isNan(inf()); - assert !isInf(nan()); - assert isFinite(0.0); - assert isFinite(-123.456); - assert !isFinite(nan()); - assert !isNan(1.0); - assert !isInf(1.0e300); // large but finite - - // signbit - assert signbit(-1.0); - assert !signbit(1.0); - assert signbit(-inf()); - assert !signbit(inf()); - - // copysign - assert copysign(3.0, -2.0) == -3.0; - assert copysign(-3.0, 2.0) == 3.0; - assert copysign(3.0, 5.0) == 3.0; - - printf("All assertions passed!"); -} diff --git a/test/test-files/std/math/fct-exp-log-extra/cout.out b/test/test-files/std/math/fct-exp-log-extra/cout.out deleted file mode 100644 index cebe8db0e..000000000 --- a/test/test-files/std/math/fct-exp-log-extra/cout.out +++ /dev/null @@ -1 +0,0 @@ -All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-exp-log-extra/source.spice b/test/test-files/std/math/fct-exp-log-extra/source.spice deleted file mode 100644 index 88cf236ec..000000000 --- a/test/test-files/std/math/fct-exp-log-extra/source.spice +++ /dev/null @@ -1,23 +0,0 @@ -import "std/math/fct"; - -f main() { - // exp2 - assert almostEqual(exp2(0.0), 1.0); - assert almostEqual(exp2(10.0), 1024.0, 1.0e-3); - assert almostEqual(exp2(-1.0), 0.5); - - // expm1 (accurate near zero) - assert almostEqual(expm1(0.0), 0.0); - assert almostEqual(expm1(1.0), 1.718281828); - assert almostEqual(expm1(0.0001), 0.0001000050); - assert almostEqual(expm1(-0.5), -0.393469340); - - // log1p (accurate near zero) - assert almostEqual(log1p(0.0), 0.0); - assert almostEqual(log1p(1.0), 0.693147180); - assert almostEqual(log1p(0.0001), 0.0000999950); - assert log1p(-1.0) == 0.0; // out of domain convention - assert log1p(-2.0) == 0.0; - - printf("All assertions passed!"); -} diff --git a/test/test-files/std/math/fct-inverse-hyperbolic/cout.out b/test/test-files/std/math/fct-inverse-hyperbolic/cout.out deleted file mode 100644 index cebe8db0e..000000000 --- a/test/test-files/std/math/fct-inverse-hyperbolic/cout.out +++ /dev/null @@ -1 +0,0 @@ -All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-inverse-hyperbolic/source.spice b/test/test-files/std/math/fct-inverse-hyperbolic/source.spice deleted file mode 100644 index 94cf05402..000000000 --- a/test/test-files/std/math/fct-inverse-hyperbolic/source.spice +++ /dev/null @@ -1,25 +0,0 @@ -import "std/math/fct"; - -f main() { - // asinh - assert almostEqual(asinh(0.0), 0.0); - assert almostEqual(asinh(1.0), 0.881373587); - assert almostEqual(asinh(-1.0), -0.881373587); - assert almostEqual(asinh(2.0), 1.443635475); - - // acosh - assert almostEqual(acosh(1.0), 0.0); - assert almostEqual(acosh(2.0), 1.316957897); - assert almostEqual(acosh(10.0), 2.993222846); - assert acosh(0.5) == 0.0; // out of domain convention - - // atanh - assert almostEqual(atanh(0.0), 0.0); - assert almostEqual(atanh(0.5), 0.549306144); - assert almostEqual(atanh(-0.5), -0.549306144); - assert atanh(1.0) == 0.0; // out of domain convention - assert atanh(-1.0) == 0.0; - assert atanh(2.0) == 0.0; - - printf("All assertions passed!"); -} diff --git a/test/test-files/std/math/fct-rounding-double/cout.out b/test/test-files/std/math/fct-rounding-double/cout.out deleted file mode 100644 index cebe8db0e..000000000 --- a/test/test-files/std/math/fct-rounding-double/cout.out +++ /dev/null @@ -1 +0,0 @@ -All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-rounding-double/source.spice b/test/test-files/std/math/fct-rounding-double/source.spice deleted file mode 100644 index d8e503934..000000000 --- a/test/test-files/std/math/fct-rounding-double/source.spice +++ /dev/null @@ -1,62 +0,0 @@ -import "std/math/fct"; - -f main() { - // truncD - assert truncD(1.9) == 1.0; - assert truncD(-1.9) == -1.0; - assert truncD(5.0) == 5.0; - assert truncD(1.0e20) == 1.0e20; // already integral, outside long range - - // floorD - assert floorD(1.9) == 1.0; - assert floorD(-1.1) == -2.0; - assert floorD(2.0) == 2.0; - assert floorD(-2.0) == -2.0; - - // ceilD - assert ceilD(1.1) == 2.0; - assert ceilD(-1.9) == -1.0; - assert ceilD(2.0) == 2.0; - - // roundD - assert roundD(1.4) == 1.0; - assert roundD(1.5) == 2.0; - assert roundD(-1.5) == -2.0; - assert roundD(-1.4) == -1.0; - - // modf - double ip = 0.0; - double frac = modf(3.75, ip); - assert ip == 3.0; - assert almostEqual(frac, 0.75); - double ip2 = 0.0; - double frac2 = modf(-2.25, ip2); - assert ip2 == -2.0; - assert almostEqual(frac2, -0.25); - - // fmod - assert almostEqual(fmod(5.3, 2.0), 1.3); - assert almostEqual(fmod(-5.3, 2.0), -1.3); - assert isNan(fmod(1.0, 0.0)); - - // remainder - assert almostEqual(remainder(5.3, 2.0), -0.7); - assert almostEqual(remainder(5.0, 3.0), -1.0); - - // frexp - int e = 0; - double m = frexp(8.0, e); - assert e == 4; - assert almostEqual(m, 0.5); - int e2 = 0; - double m2 = frexp(0.0, e2); - assert e2 == 0; - assert m2 == 0.0; - - // ldexp - assert almostEqual(ldexp(0.5, 4), 8.0); - assert almostEqual(ldexp(1.5, 2), 6.0); - assert almostEqual(ldexp(3.0, -1), 1.5); - - printf("All assertions passed!"); -} diff --git a/test/test-files/std/math/fct-special/cout.out b/test/test-files/std/math/fct-special/cout.out deleted file mode 100644 index cebe8db0e..000000000 --- a/test/test-files/std/math/fct-special/cout.out +++ /dev/null @@ -1 +0,0 @@ -All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-special/source.spice b/test/test-files/std/math/fct-special/source.spice deleted file mode 100644 index ea6a1e6af..000000000 --- a/test/test-files/std/math/fct-special/source.spice +++ /dev/null @@ -1,30 +0,0 @@ -import "std/math/fct"; - -f main() { - // tgamma - assert almostEqual(tgamma(1.0), 1.0, 1.0e-6); - assert almostEqual(tgamma(2.0), 1.0, 1.0e-6); - assert almostEqual(tgamma(3.0), 2.0, 1.0e-6); - assert almostEqual(tgamma(4.0), 6.0, 1.0e-5); - assert almostEqual(tgamma(5.0), 24.0, 1.0e-4); - assert almostEqual(tgamma(0.5), 1.772453851, 1.0e-6); // sqrt(PI) - assert almostEqual(tgamma(1.5), 0.886226925, 1.0e-6); // sqrt(PI) / 2 - - // lgamma - assert almostEqual(lgamma(1.0), 0.0, 1.0e-6); - assert almostEqual(lgamma(2.0), 0.0, 1.0e-6); - assert almostEqual(lgamma(5.0), 3.178053830, 1.0e-5); // ln(24) - assert almostEqual(lgamma(0.5), 0.572364942, 1.0e-6); // ln(sqrt(PI)) - - // erf - assert almostEqual(erf(0.0), 0.0, 1.0e-6); - assert almostEqual(erf(1.0), 0.842700793, 1.0e-6); - assert almostEqual(erf(-1.0), -0.842700793, 1.0e-6); - assert almostEqual(erf(0.5), 0.520499878, 1.0e-6); - - // erfc - assert almostEqual(erfc(0.0), 1.0, 1.0e-6); - assert almostEqual(erfc(1.0), 0.157299207, 1.0e-6); - - printf("All assertions passed!"); -} diff --git a/test/test-files/std/math/fct-utility/cout.out b/test/test-files/std/math/fct-utility/cout.out deleted file mode 100644 index cebe8db0e..000000000 --- a/test/test-files/std/math/fct-utility/cout.out +++ /dev/null @@ -1 +0,0 @@ -All assertions passed! \ No newline at end of file diff --git a/test/test-files/std/math/fct-utility/source.spice b/test/test-files/std/math/fct-utility/source.spice deleted file mode 100644 index 421a5e98f..000000000 --- a/test/test-files/std/math/fct-utility/source.spice +++ /dev/null @@ -1,42 +0,0 @@ -import "std/math/fct"; - -f main() { - // signum - assert signum(5) == 1; - assert signum(-5) == -1; - assert signum(0) == 0; - assert signum(-3.5) == -1.0; - assert signum(3.5) == 1.0; - assert signum(0.0) == 0.0; - assert signum(-7l) == -1l; - - // fma - assert almostEqual(fma(2.0, 3.0, 4.0), 10.0); - assert almostEqual(fma(-2.0, 3.0, 1.0), -5.0); - - // lerp - assert almostEqual(lerp(0.0, 10.0, 0.0), 0.0); - assert almostEqual(lerp(0.0, 10.0, 1.0), 10.0); - assert almostEqual(lerp(0.0, 10.0, 0.5), 5.0); - assert almostEqual(lerp(2.0, 4.0, 0.25), 2.5); - - // isqrt - assert isqrt(0) == 0; - assert isqrt(1) == 1; - assert isqrt(16) == 4; - assert isqrt(17) == 4; - assert isqrt(99) == 9; - assert isqrt(100) == 10; - assert isqrt(-5) == 0; - assert isqrt(1000000000000l) == 1000000l; - - // array max / min - int[6] ints = [ 3, 7, 1, 9, 2, 5 ]; - assert max(ints, len(ints)) == 9; - assert min(ints, len(ints)) == 1; - double[4] doubles = [ 2.5, -1.5, 3.25, 0.0 ]; - assert almostEqual(max(doubles, len(doubles)), 3.25); - assert almostEqual(min(doubles, len(doubles)), -1.5); - - printf("All assertions passed!"); -} diff --git a/test/test-files/std/math/rand-extra/source.spice b/test/test-files/std/math/rand-extra/source.spice index dfcad044b..02438780e 100644 --- a/test/test-files/std/math/rand-extra/source.spice +++ b/test/test-files/std/math/rand-extra/source.spice @@ -1,5 +1,4 @@ import "std/math/rand"; -import "std/math/fct"; f main() { // Reproducibility: the same seed yields the same sequence @@ -50,10 +49,13 @@ f main() { } assert sawTrue && sawFalse; - // randGaussian returns finite values - for int i = 0; i < 1000; i++ { - assert isFinite(randGaussian(0.0, 1.0)); - } + // randGaussian is reproducible for a given seed + seed(55ul); + double g1 = randGaussian(0.0, 1.0); + double g2 = randGaussian(0.0, 1.0); + seed(55ul); + assert randGaussian(0.0, 1.0) == g1; + assert randGaussian(0.0, 1.0) == g2; printf("All assertions passed!"); } From a696df8c1d2cdca9a5cba8565675bb8d746ab74c Mon Sep 17 00:00:00 2001 From: Auberer Date: Mon, 15 Jun 2026 18:28:03 +0200 Subject: [PATCH 3/3] Rename RNG state global to RAND_STATE for valid type identifier Co-Authored-By: Claude Opus 4.8 --- std/math/rand.spice | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/std/math/rand.spice b/std/math/rand.spice index 2fb30a68c..d6e7728a1 100644 --- a/std/math/rand.spice +++ b/std/math/rand.spice @@ -3,7 +3,7 @@ import "std/math/fct"; // Internal generator state (xorshift64*). Seeded with a fixed non-zero default so that // an unseeded program still produces a well-defined, reproducible sequence. -unsigned long randState = 0x2545f4914f6cdd1dul; +unsigned long RAND_STATE = 0x2545f4914f6cdd1dul; /** * Advance the internal xorshift64* generator and return the next raw 64-bit value. @@ -11,11 +11,11 @@ unsigned long randState = 0x2545f4914f6cdd1dul; * @return Next raw pseudo-random value */ f nextRaw() { - unsigned long x = randState; + unsigned long x = RAND_STATE; x ^= x >> 12ul; x ^= x << 25ul; x ^= x >> 27ul; - randState = x; + RAND_STATE = x; return x * 0x2545f4914f6cdd1dul; } @@ -26,7 +26,7 @@ f nextRaw() { * @param seedValue Seed value */ public p seed(unsigned long seedValue) { - randState = seedValue == 0ul ? 0x9e3779b97f4a7c15ul : seedValue; + RAND_STATE = seedValue == 0ul ? 0x9e3779b97f4a7c15ul : seedValue; } /**