Skip to content

Instantly share code, notes, and snippets.

@Marc-B-Reynolds
Marc-B-Reynolds / xoroshiro128p_fail.c
Last active April 6, 2025 12:27
test driver that demos "xoroshiro128+" failing PractRand
// compile with whatever then run PractRand:
// ./test | RNG_test stdin64 -tlmin 256KB -tf 2 -tlmax 512GB -seed 0
//****************************************************************************
// verbatim from: https://prng.di.unimi.it/xoroshiro128plus.c
/* Written in 2016-2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
To the extent possible under law, the author has dedicated all copyright
and related and neighboring rights to this software to the public domain
@Marc-B-Reynolds
Marc-B-Reynolds / hastings_atan.sollya
Last active February 8, 2025 18:36
Using sollya to "repeat" Cecil Hastings' 1955 approximations for atan. Works as an example of using a 'seed' approximation & using sollya to refine.
// Repeating Cecil Hastings atan approximations from 1955 using Sollya.
// Other than "amusement factor" these use Hasting's approximations
// to "seed" the sollya calls to 'fpminimax'.
//
// Current world usage could include using an initial approximation
// from some multiprecision or computer algebra system and then
// Sollya to refine the constants to hardware sized (half,single,double,etc)
//
// Additionally the wide choice of the range would cause Sollya to fail
// to converge if we didn't feed it an initial guess.
@Marc-B-Reynolds
Marc-B-Reynolds / atan_spitball.c
Created January 20, 2025 21:11
3 minimal single parameter atan implementations on [-1,1]. Only accurate to about 1/5 a degree.
// Public Domain under http://unlicense.org, see link for details.
// EXCEPT "CORE-MATH" marked section
// {gcc,clang} -O3 -march=native -Wall -Wextra -Wconversion -Wpedantic -Wno-unused-function -fno-math-errno -ffp-contract=off atan_spitball.c -o atan_spitball -lm
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
@Marc-B-Reynolds
Marc-B-Reynolds / mum_max.md
Created January 3, 2025 10:20
MUM max 32 'k' values for

For bit widths on $[8,18]$ brute force all values of $k$ for the MUM step and discover the 32 values that produce the largest output set (codomain). pop is poulation count, runs is the number of runs of 1s (linear).

bits = 8

k size $2^{-8}$ size pop runs k (binary)
0000009e 00000073 0.449219 5 2 1..1111.
00000076 0000006f 0.433594 5 2 .111.11.
00000016 0000006e 0.429688 3 2 ...1.11.
@Marc-B-Reynolds
Marc-B-Reynolds / prng_pop.c
Created December 6, 2024 14:38
Skims two methods: bisection method for random numbers with given popcount & branchfree random bit permutation with bit-scatter or sheep-and-goats op.
// UNLICENSE or whatever other you want that doesn't hold me
// responsible for you choosing to use any of this.
// a sketch of a few niche things:
// 1) generate a random number with a given population count
// 2) perform a random bit permutation of the input
// 3) not quite random versions of (2)
// 4) use (2) or (3) to produce stateful variants of (1)
#include <stdint.h>
@Marc-B-Reynolds
Marc-B-Reynolds / test.md.html
Last active November 13, 2024 04:06
Single file: markdeep + plotly that mostly works

Header

Some bold text before an inline function: $y = x^2$

NOTE: The gist preview is completely whacked. Click on raw for the source.

@Marc-B-Reynolds
Marc-B-Reynolds / unorm8_to_half.sollya
Last active November 9, 2024 15:08
UNORM8 to half float validation
/* -*- mode: c; -*- */
// multiply by recip: 10 aren't correctly rounded
print("multiply by 1/255");
for i from 0 to 255 do
{
cr = HP(i/255);
r = HP(i*0x1.01p-8);
@Marc-B-Reynolds
Marc-B-Reynolds / f32_quadratic_hq.c
Created July 8, 2023 08:16
single precision quadratic via double promotion. avoids spurious over and underflow.
// whatever return type
typedef struct { float h,l; } f32_pair_t;
f32_pair_t f32_quadratic_hq(float A, float B, float C)
{
double a = (double)A;
double b = (double)B;
double c = (double)C;
@Marc-B-Reynolds
Marc-B-Reynolds / quick.c
Last active January 31, 2023 23:30
50 example 64-bit bijections using CRC32-C opcode (brute force validate)
// if 'f' is the 64-bit input CRC32-C (with 'incoming' 32-bit value of zero) then
// there are 50 bijections (invertiable functions) which are a sum of f
// and a simple xorshift. So the follow two forms:
//
// f(x) ^ (x ^ (x << S))
// f(x) ^ (x ^ (x >> S))
//
// and 23 using a rotate instead of a shift:
//
@Marc-B-Reynolds
Marc-B-Reynolds / cbrt.c
Last active February 23, 2023 18:03
bit-exact portable cube root and reciprocal thereof
// Public Domain under http://unlicense.org, see link for details.
// except:
// * core-math function `cr_cbrtf` (see license below)
// * musl flavored fdlib function `fdlibm_cbrtf` (see license below)
// code and test driver for cube root and it's reciprocal based on:
// "Fast Calculation of Cube and Inverse Cube Roots Using
// a Magic Constant and Its Implementation on Microcontrollers"
// Moroz, Samotyy, Walczyk, Cieslinski, 2021
// (PDF: https://www.mdpi.com/1996-1073/14/4/1058)