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/rand.spice b/std/math/rand.spice index 6618fc856..d6e7728a1 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 RAND_STATE = 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 = RAND_STATE; + x ^= x >> 12ul; + x ^= x << 25ul; + x ^= x >> 27ul; + RAND_STATE = 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) { + RAND_STATE = 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/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..02438780e --- /dev/null +++ b/test/test-files/std/math/rand-extra/source.spice @@ -0,0 +1,61 @@ +import "std/math/rand"; + +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 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!"); +}