Rokusa (Pure)

README.md

Some useful code written in ion

aggregations.ion

// @opportunity @perf: reduce data dependency by processing n-elems at a time SIMD style.

func sum_array(first: float*, num: usize, stride_size: usize): float
{
    #assert(_is_aligned(stride_size, alignof(*first)));
    res: float;
    for (i := 0; i<num; i++) {
        res += *_deref_float(first, stride_size, i);
    }
    return res;
}

func mean_array(first: float*, num: usize, stride_size: usize): float
{
    return sum_array(first, num, stride_size) / float(num);
}

func root_mean_square_array(first: float*, num: usize, stride_size: usize): float
{
    #assert(_is_aligned(stride_size, alignof(*first)));
    res: float;
    for (i:=0; i<num; i++) {
        x := *_deref_float(first, stride_size, i);
        res += x * x;
    }
    return libc.sqrt(res / float(num));
}

func minmax_array(first: float*, num: usize, stride_size: usize): Range_Float
{
    #assert(_is_aligned(stride_size, alignof(*first)));
    res := minmax_unit();
    for (i:=0; i<num; i++) {
        x := *_deref_float(first, stride_size, i);
        minmax_inplace(&res, x);
    }
    return res;
}

func minmax_unit(): Range_Float { return aabb1_cover_unit(); }

func minmax_inplace(range: Range_Float*, x: float): Range_Float
{
    *range = aabb1_cover_pt(*range, {x});
    return *range;
}



func _deref_float(first: float*, stride_size: usize, index: usize): float*
{
    dest_byte_ptr := (:uint8*)first;
    return (:float*)(dest_byte_ptr + index*stride_size);
}

func _is_aligned(offset: uintptr, alignment: usize): bool
{
    return 0 == offset & (alignment - 1);
}

audio_noise.ion

// pink noise: aka 1/f noise, has a frequency rolloff of 3dB per octave.

struct PinkPK3
{
    b: float[7];
    output: float;
}

// white noise filtered by an approximated 3dB/octave filter
func pink_noise_pk3_44100hz_advance(pk3: PinkPK3*, white_noise: float)
{
    // Filter by Paul Kellet (pk3 = (Black))
    // [email protected]
    //
    // Used in CSound for instance.
    //
    // Filter to make pink noise from white  (updated March 2000)
    // ------------------------------------
    //
    // This is an approximation to a -10dB/decade filter using a weighted sum
    // of first order filters. It is accurate to within +/-0.05dB above 9.2Hz
    // (44100Hz sampling rate). Unity gain is at Nyquist, but can be adjusted
    // by scaling the numbers at the end of each line.
    //
    // If 'white' consists of uniform random numbers, such as those generated
    // by the rand() function, 'pink' will have an almost gaussian level
    // distribution.
    s := pk3;
    s.b[0] = 0.99886 * s.b[0] + white_noise * 0.0555179;
    s.b[1] = 0.99332 * s.b[1] + white_noise * 0.0750759;
    s.b[2] = 0.96900 * s.b[2] + white_noise * 0.1538520;
    s.b[3] = 0.86650 * s.b[3] + white_noise * 0.3104856;
    s.b[4] = 0.55000 * s.b[4] + white_noise * 0.5329522;
    s.b[5] = -0.7616 * s.b[5] - white_noise * 0.0168980;
    s.output = s.b[0] + s.b[1] + s.b[2] + s.b[3] + s.b[4] + s.b[5] + s.b[6] + white_noise * 0.5362;
    // gain, to keep output within [-1,1] (tried for about 530M samples)
    s.output *= 0.094229373;
    s.b[6] = white_noise * 0.115926;
}

bits.ion

func mask_uint32(offset: int): uint32
{
    return offset >= 32? ~uint32(0) : ((1<<offset) - 1);
}

func mask_uint64(offset: int): uint64
{
    return  offset >= 64? ~uint64(0) : ((uint64(1)<<offset) - 1);
}

func mask_uint16(offset: int): uint16
{
    return offset >= 16? ~uint16(0) : ((uint16(1)<<offset) - 1);
}

func bits_uint32(x: uint32, start: int, num_bits: int): uint32
{
    #assert(start >= 0 && num_bits >= 0);
    #assert(start + num_bits <= 32);
    return (x>>start) & mask_uint32(num_bits);
}

func bits_uint64(x: uint64, start: int, num_bits: int): uint64
{
    #assert(start >= 0 && num_bits >= 0);
    #assert(start + num_bits <= 64);
    return (x>>start) & mask_uint64(num_bits);
}

func bits_uint16(x: uint16, start: int, num_bits: int): uint16
{
    #assert(start >= 0 && num_bits >= 0);
    #assert(start + num_bits <= 16);
    return (x>>start) & mask_uint16(num_bits);
}

func rotl_uint32(x: uint32, offset: int): uint32
{
    #assert(offset >= 0 && offset < 32);
    return (x << offset) | (x >> (32 - offset));
}

func rotl_uint64(x: uint64, offset: int): uint64
{
    #assert(offset >= 0 && offset < 64);
    return (x << offset) | (x >> (64 - offset));
}

func rotr_uint32(x: uint32, offset: int): uint32
{
    #assert(offset >= 0 && offset < 32);
    return (x >> offset) | (x << (32 - offset));
}

func rotr_uint64(x: uint64, offset: int): uint64
{
    #assert(offset >= 0 && offset < 64);
    return (x >> offset) | (x << (64 - offset));
}

func lowest_bit_only_uint32(x: uint32): uint32
{
    y : uint32 = -int32(x);
    return x & y;
}

func lowest_bit_only_uint64(x: uint64): uint64
{
    y : uint64 = -int64(x);
    return x & y;
}

func mask_below_uint32(x: uint32, num_bits: int): uint32
{
    return x & ~mask_uint32(num_bits);
}

func mask_below_uint64(x: uint64, num_bits: int): uint64
{
    return x & ~mask_uint64(num_bits);
}

// bitcount: also known as:
// popcnt for x86 chips
// popcount for arm chips

func bitcount_uint32(x: uint32): int
{
    ca: int;
    cb: int;
    a := x;
    b := x>>1;
    n := 32;
    while (n && (a || b)) {
        ca += a & 0b1;
        cb += b & 0b1;
        a >>= 2;
        b >>= 2;
        n -= 2;
    }
    return ca + cb;
}

func bitcount_uint64(x: uint64): int
{
    ca: int;
    cb: int;
    a := x;
    b := x>>1;
    n := 64;
    while (n && (a || b)) {
        ca += a & 0b1;
        cb += b & 0b1;
        a >>= 2;
        b >>= 2;
        n -= 2;
    }
    return ca + cb;
}

covers.ion

// Covers / axis aligned bounding boxes.
// (Inclusive bounded ranges)
//
// Operations:
// - `cover` creates a box from two points, adding points to boxes or adding boxes together.
//   The unit of the cover operation is a kind of invalid box, a kind of zero box.
// - `dilate`/`erode` resize a box
// 

// Range
struct AxisAlignedBoundingBox1D
{
    min, max: float;
}
typedef aabb1 = AxisAlignedBoundingBox1D;
typedef Range_Float = aabb1;

func aabb1_add_inplace(a: aabb1*, b: float): aabb1 { *a = aabb1_add(*a, b); return *a; }
func aabb1_cover_pt_inplace(a: aabb1*, b: pt1): aabb1 { *a = aabb1_cover_pt(*a, b); return *a; }

func aabb1_pt_size(a: pt1, b: float): aabb1
{
    return {
        min = a.x,
        max = a.x + b,
    };
}

func pt1_cover(a: pt1, b: pt1): aabb1
{
    return a.x<b.x? { min=a.x, max=b.x } : { min=b.x, max=a.x };
}

func aabb1_add(a: aabb1, b: float): aabb1
{
    return { min = a.min + b, max = a.max + b };
}

func aabb1_cover_unit(): aabb1 { return { FLOAT_MAX, FLOAT_MIN }; }

func aabb1_cover_pt(a: aabb1, b: pt1): aabb1
{
    return {
        min = min_float(a.min, b.x),
        max = max_float(a.max, b.x),
    };
}

struct AxisAlignedBoundingBox2D
{
    union {
        struct { 
            min_x, min_y: float;
            max_x, max_y: float;
        }
        struct {
            min_xy: float[2];
            max_xy: float[2];
        }
        ltrb: float[4];
    }
}
#static_assert(sizeof(AxisAlignedBoundingBox2D) == 2*2*sizeof(float))

typedef aabb2 = AxisAlignedBoundingBox2D;
func aabb2_add_inplace(a: aabb2*, b: vec2): aabb2 { *a = aabb2_add(*a, b); return *a; }
func aabb2_cover_pt_inplace(a: aabb2*, p: pt2): aabb2 { *a = aabb2_cover_pt(*a, p); return *a; }
func aabb2_cover_inplace(a: aabb2*, b: aabb2): aabb2 { *a = aabb2_cover(*a, b); return *a; }
func aabb2_sub_minkowski_inplace(a: aabb2*, b: aabb2): aabb2 { *a = aabb2_sub_minkowski(*a, b); return *a; }
func aabb2_dilate_inplace(a: aabb2*, b: vec2): aabb2 { *a = aabb2_dilate(*a, b); return *a; }
func aabb2_erode_inplace(a: aabb2*, b: vec2): aabb2 { *a = aabb2_erode(*a, b); return *a; }
func aabb2_cover_unit(): aabb2 { return { min_x=FLOAT_MAX, min_y=FLOAT_MAX, max_x=FLOAT_MIN, max_y=FLOAT_MIN };}

func pt2_box_cover(a: pt2, b: pt2): aabb2
{
    mmx := pt1_cover({a.vec.x}, {b.vec.x});
    mmy := pt1_cover({a.vec.y}, {b.vec.y});
    return {
        min_x = mmx.min,
        min_y = mmy.min,
        max_x = mmx.max,
        max_y = mmy.max,
    };
}

func aabb2_intersects(a: aabb2, p: pt2): bool
{
    if (p.vec.x < a.min_x || p.vec.x > a.max_x) { return false; }
    if (p.vec.y < a.min_y || p.vec.y > a.max_y) { return false; }
    return true;
}

func aabb2_min_pt(a: aabb2*): pt2 const*
{
    pt2_ptr : pt2*;
    f := &a.min_xy[0]; 
    memcpy(&pt2_ptr, &f, sizeof(pt2_ptr));
    return pt2_ptr;
}

func aabb2_max_pt(a: aabb2*): pt2 const*
{
    pt2_ptr : pt2*;
    f := &a.max_xy[0]; 
    memcpy(&pt2_ptr, &f, sizeof(pt2_ptr));
    return pt2_ptr;
}

func aabb2_center_pt(a: aabb2): pt2
{
    return {{
        x = 0.5*(a.max_xy[0] + a.min_xy[0]),
        y = 0.5*(a.max_xy[1] + a.min_xy[1]),
    }};
}

func aabb2_add(a: aabb2, b: vec2): aabb2
{
    return {
        min_x = a.min_x + b.x,
        min_y = a.min_y + b.y,
        max_x = a.max_x + b.x,
        max_y = a.max_y + b.y,
    };
}

// box dilation
func aabb2_sub_minkowski(a: aabb2, b: aabb2): aabb2
{
    return {
        min_x = a.min_x - b.max_x,
        min_y = a.min_y - b.max_y,
        max_x = a.max_x - b.min_x,
        max_y = a.max_y - b.min_y,
    };
}

func aabb2_dilate(a: aabb2, b: vec2): aabb2
{
    #assert(b.x >= 0 && b.y >= 0);
    return {
        min_x = a.min_x - b.x,
        min_y = a.min_y - b.y,
        max_x = a.max_x + b.x,
        max_y = a.max_y + b.y,
    };
}

func aabb2_erode(a: aabb2, b: vec2): aabb2
{
    #assert(b.x >= 0 && b.y >= 0);
    hx := min_float(0.5*(a.max_x - a.min_x), b.x);
    hy := min_float(0.5*(a.max_y - a.min_y), b.y);
    return {
        min_x = a.min_x + hx,
        min_y = a.min_y + hy,
        max_x = a.max_x - hx,
        max_y = a.max_y - hy,
    };
}

func aabb2_cover_origin_size(o: pt2, size: vec2): aabb2
{
    return pt2_box_cover(o, pt2_add(o, pt2{size}));
}

func aabb2_cover(a: aabb2, b: aabb2): aabb2
{
    return {
        min_x = min_float(a.min_x, b.min_x),
        min_y = min_float(a.min_y, b.min_y),
        max_x = max_float(a.max_x, b.max_x),
        max_y = max_float(a.max_y, b.max_y),
    };
}

func aabb2_cover_pt(a: aabb2, p: pt2): aabb2
{
    return {
        min_x = min_float(a.min_x, p.vec.x),
        min_y = min_float(a.min_y, p.vec.y),
        max_x = max_float(a.max_x, p.vec.x),
        max_y = max_float(a.max_y, p.vec.y),
    };
}

domains.ion

import libc {
    fabs,
    round,
    tanh
}

func clamp(x: int, f: int, l: int): int { return min(max(x, f), l); }
func clamp_float(x: float, f: float, l: float): float { return min(max(x, f), l); }

func clamp_int8(x: int8, f: int8, l: int8): int8 { return min(max(x, f), l); }
func clamp_int16(x: int16, f: int16, l: int16): int16 { return min(max(x, f), l); }
func clamp_int32(x: int32, f: int32, l: int32): int32 { return min(max(x, f), l); }
func clamp_int64(x: int64, f: int64, l: int64): int64 { return min(max(x, f), l); }

func quantize_float(v: float, denominator: float): float
{
    return round(v*denominator)/denominator;
}

// (-inf,inf) to [-1,1]
func map_unit_from_all(x: float): float
{
  return x / (1.0 + fabs(x));
}

// (-inf,inf) to [-1,1]
func map_unit_from_all_tanh(x: float): float
{
  return tanh(x);
}

// [0,inf) to [0,1]
func map_positive_unit_from_positives(x: float): float
{
  return x / (1 + x);
}

// [-INT_MAX,INT_MAX] to [-1,1] (INT_MIN will return a value outside of the range)
func map_unit_from_int32(x: int32): float
{
  // compute in doubles, which can safely represent all 32bit integers
  return float(double(x) / double(INT_MAX));
}

floats.ion

// ieee754 floats
// @url: https://en.wikipedia.org/wiki/IEEE_754

import libc { memcpy } // memcpy is used here as a bit-casting intrinsic

#static_assert(sizeof(float) == sizeof(uint32))

const FLOAT_MAX = 3.40282346e+38;
const FLOAT_MIN = -3.40282346e+38;
const FLOAT_EXPONENT_MIN = 1;
const FLOAT_EXPONENT_MAX = 254;
const FLOAT_MANTISSA_MIN = 0;
const FLOAT_MANTISSA_MAX = 0x7f_ffff;

const DOUBLE_MAX = +1.3482698511467367e+308d;
const DOUBLE_MIN = -1.3482698511467367e+308d;
const DOUBLE_EXPONENT_MIN = 1;
const DOUBLE_EXPONENT_MAX = 2046;
const DOUBLE_MANTISSA_MIN = 0;
const DOUBLE_MANTISSA_MAX = 0x7_ffff_ffff_ffff;

// Some transcendental numbers:
const FLOAT_TWOPI = 6.28318530717958647693;
const FLOAT_PI = 3.14159265358979323846;
const FLOAT_HALFPI = 1.57079632679489661923;
const FLOAT_SQRT2 = 1.41421356237;

struct UnpackedFloat32
{
    mantissa: uint32;
    exponent: uint32;
    is_negative: bool;
}

struct UnpackedFloat64
{
    mantissa: uint64;
    exponent: uint64;
    is_negative: bool;
}

enum FloatGroup
{
    FloatGroup_Zero,
    FloatGroup_Infinite,
    FloatGroup_NaNs,
    FloatGroup_Denormals,
    FloatGroup_Normals,
}

func classify_float32(v: float): FloatGroup
{
    uv := unpack_float32(v);
    if (uv.exponent == 0) {
        return uv.mantissa == 0? FloatGroup_Zero:FloatGroup_Denormals;
    } else if (uv.exponent == 0xff) {
        return uv.mantissa == 0? FloatGroup_Infinite:FloatGroup_NaNs;
    }
    return FloatGroup_Normals;
}

func classify_float64(v: double): FloatGroup
{
    uv := unpack_float64(v);
    if (uv.exponent == 0) {
        return uv.mantissa == 0? FloatGroup_Zero:FloatGroup_Denormals;
    } else if (uv.exponent == 2047) {
        return uv.mantissa == 0? FloatGroup_Infinite:FloatGroup_NaNs;
    }
    return FloatGroup_Normals;
}

func is_denormal_float32(v: float): bool
{
    return classify_float32(v) == FloatGroup_Denormals;
}

func unpack_float32(v: float): UnpackedFloat32
{
    b32: uint32;
    memcpy(&b32, &v, sizeof(b32));
    return {
        mantissa = bits_uint32(b32, 0, 23),
        exponent = bits_uint32(b32, 23, 8),
        is_negative = bits_uint32(b32, 31, 1),
    };
}

func pack_float32(uv: UnpackedFloat32): float
{
    #assert(uv.exponent <= mask_uint32(8));
    #assert(uv.mantissa <= mask_uint32(23));
    b32: uint32 = uv.mantissa | uv.exponent<<23 | uint32(uv.is_negative)<<31;
    v: float;
    memcpy(&v, &b32, sizeof(v));
    return v;
}

func unpack_float64(v: double): UnpackedFloat64
{
    b64: uint64;
    memcpy(&b64, &v, sizeof(b64));
    return {
        mantissa = bits_uint64(b64, 0, 52),
        exponent = bits_uint64(b64, 52, 11),
        is_negative = bits_uint64(b64, 63, 1),
    };
}

func pack_float64(uv: UnpackedFloat64): double
{
    #assert(uv.exponent <= mask_uint64(11));
    #assert(uv.mantissa <= mask_uint64(52));
    b64: uint64 = uv.mantissa | uv.exponent<<52 | uint64(uv.is_negative)<<63;
    v: double;
    memcpy(&v, &b64, sizeof(v));
    return v;
}

integrations.ion

// integration methods

// verlet integration
//
 // Integrate in many timesteps the acceleration and velocity of a
 // particle, in order to find its new position, velocity and
 // acceleration.
 //
 // More to the point, it computes x(t+timestep), v(t)
 // from:
 //  - x(t)
 //   - x(t-timestep) -- you must keep two steps of position memory
 //   - a(t)
 //
 // x(t+timestep) is returned, and out_v will contain v(t) if non null.
 //
 // @param x the dimension to integrate. x[0] corresponds to the
 // current value. It is update with the new current value.
 // @param x_1[0] to the previous value. It is updated to reflect
 // the new previous value.
 // @param a is the acceleration to apply
 // @param out_v is optional
 func verlet_integrate_inplace(timestep: float, a: float, x: float*, x_1: float*, out_v: float*): float
 {
    sq := timestep * timestep;
//	new_x := sq*(a + 2.0*x[0]/sq - x_1[0]/sq);
    new_x := sq*a + (2.0*x[0] - x_1[0]);

	if (out_v) {
		*out_v = (new_x - x_1[0]) / (2.0*timestep);
	}

	x_1[0] = x[0];
	x[0] = new_x;


	return new_x;
 }

linear_feedback_shift_registers.ion

// LFSR: Linear Feedback Shift Register

// @ion @defect register cannot be used as a symbol, it must be escaped when
// generated in C. This is probably true of every keyword in C that fails to
// be escaped.

func lfsr16_advance(reg: uint16*): bool
{
    x := *reg;
    bit: bool = 
        bits_uint16(x, 1, 1) ^
        bits_uint16(x, 2, 1) ^
        bits_uint16(x, 3, 1) ^
        bits_uint16(x, 5, 1);
    *reg = (x>>1) | (bit<<15);
    
    return bit;
}

func lfsr16_current(reg: uint16): bool
{
    return bits_uint16(reg, 15, 1);
}

// @url: https://github.com/Stenzel/newshadeofpink/blob/master/pink52.h
func lfsr64_stenzel_advance(reg: uint64*): bool
{
    x := *reg;
    bit: bool = bits_uint64(x, 63, 1);
    x <<= 1;
    x ^= bit & 0xA800_0000_0000_0001ull;    
    *reg = x;
    return bit;
}

models.ion

struct Line2d_Float
{
    yslope: float;
    y: float;
}

// @precondition: input x MUST not all be equal
// @precondition: num>1
func linear_least_square_2d_array(first_x: float*, first_y: float*, num: usize, stride_x: usize, stride_y: usize): Line2d_Float
{
    #assert(0 == (stride_x & (alignof(*first_x) - 1)));
    #assert(0 == (stride_y & (alignof(*first_y) - 1)));
    #assert(num > 1);

    xsum: float;
    sqxsum: float;
    ysum: float;
    xysum: float;
    for (i:=0; i<num; i++) {
        x := *_deref_float(first_x, stride_x, i);
        xsum += x;
        sqxsum += x*x;
        y := *_deref_float(first_y, stride_y, i);
        ysum += y;
        xysum += x*y;
    }
    scale := 1.0/(xsum*xsum - num*sqxsum);
    #assert(scale != 0.0);
    return {
        yslope = scale*(ysum*xsum - num*xysum),
        y = scale*(xsum*xysum - ysum*sqxsum)
    };
}

numerics.ion

// atan approximation

import libc
{
    fabs,
    copysign,
}

// polynomial approximation of arctangent.
// approximation is valid if z in [-1, 1]
// @url: https://www.dsprelated.com/showarticle/1052.php#tabs1-comments
func atan_unit_approx(z: float): float
{
     m := z<0.0? -1.0:1.0;
     z *= m;
     y := 0.141499*z*z*z*z + 
        -0.343315*z*z*z +
        -0.016224*z*z +
        1.003839*z +
        -0.000158;
     return m*y;
}

// polynomial approximation of atan2
// @url: https://www.dsprelated.com/showarticle/1052.php#tabs1-comments
func atan2_approx(y: float, x: float): float
{
     ay := fabs(y);
     ax := fabs(x);
     invert := ay > ax;
     z := invert? ax/ay : ay/ax;            // [0,1]
     th := atan_unit_approx(z);                  // [0,π/4]
     if(invert) { th = FLOAT_HALFPI - th; } // [0,π/2]
     if(x < 0)  { th = FLOAT_PI - th; }     // [0,π]
     th = copysign(th, y);                  // [-π,π]
     return th;
}

// Pseudo Sine function.
//
// Originally from Stefan Stenzel, CEO of Waldorf
// @url: https://namoseley.wordpress.com/2017/01/03/pseudo-sinusoidal-waveform-ii/
// https://twitter.com/5tenzel/status/816225338330140673

// approximation of sin(PI*x) when x in [-1,1]
func waldorf_cycle(x: float): float
{
    abs_x := x < 0.? -x:x;
    return x*4.0*(1 - abs_x);
}

order.ion

func min(a: int, b: int): int { return a<b? a:b; }
func max(a: int, b: int): int { return b<a? a:b; }

func min_float(a: float, b: float): float { return a<b? a:b; }
func max_float(a: float, b: float): float { return b<a? a:b; }

func min_int8(a: int8, b: int8): int8 { return a<b? a:b; }
func max_int8(a: int8, b: int8): int8 { return b<a? a:b; }
func min_int16(a: int16, b: int16): int16 { return a<b? a:b; }
func max_int16(a: int16, b: int16): int16 { return b<a? a:b; }
func min_int32(a: int32, b: int32): int32 { return a<b? a:b; }
func max_int32(a: int32, b: int32): int32 { return b<a? a:b; }
func min_int64(a: int64, b: int64): int64 { return a<b? a:b; }
func max_int64(a: int64, b: int64): int64 { return b<a? a:b; }

func min_double(a: double, b: double): double { return a<b? a:b; }
func max_double(a: double, b: double): double { return b<a? a:b; }

random_numbers.ion

// @def{rng}: random number generator

// @url: http://pcg-random.org
// based on minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org
// Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)
struct PCG32
{
    state: uint64; // any value possible
    inc: uint64; // @precondition{(inc & 1) == 1}: inc must be odd
}

// check if rng's preconditions are valid
func pcg32_is_valid(x: PCG32): bool
{
    return x.inc & 1 == 1;
}

// advance rng and return a result
func pcg32_advance(rng: PCG32*): uint32
{
    state := rng.state;
    xorshifted: uint32 = uint32(((state>>18u) ^ state) >> 27u);
    rot: int32 = state>>59u;
    res: uint32 = (xorshifted >> rot) | (xorshifted << (-rot & 31));
    rng.state = state * 6364136223846793005ULL + (rng.inc|1);
    return res;
}

// Generating floating pointer numbers.
//
// You may use @symbol{map_unit_from_int32}, however 
// Taylor Campbell documented that the typical division-based
// method of generating random floats gives a non-uniform
// distribution in:
//
// @url: https://mumble.net/~campbell/2014/04/28/uniform-random-float
// @abtract{
// Uniform random floats: How to generate a double-precision
// floating-point number in [0,1] uniformy at random given a 
// uniform random source of bits.
// }

vecs.ion

import libc { 
    memcpy, // used for bitcasting
    sqrt,
}

// Vecs

typedef vec1 = float;

struct Vec2D
{
    union {
        struct { x, y: float; }
        xy: float[2];
    }
}

typedef vec2 = Vec2D;
func vec2_add_inplace(a: vec2*, b: vec2): vec2 { *a = vec2_add(*a, b); return *a; }
func vec2_sub_inplace(a: vec2*, b: vec2): vec2 { *a = vec2_sub(*a, b); return *a; }
func vec2_mul_scalar_inplace(a: vec2*, s: float): vec2 { *a = vec2_mul_scalar(*a, s); return *a; }
func vec2_div_scalar_inplace(a: vec2*, s: float): vec2 { *a = vec2_div_scalar(*a, s); return *a; }
func vec2_mul_entries_inplace(a: vec2*, b: vec2): vec2 { *a = vec2_mul_entries(*a, b); return *a; }
func vec2_unit_or_zero_inplace(a: vec2*): vec2 { *a = vec2_unit_or_zero(*a); return *a; }

func vec2_add(a: vec2, b: vec2): vec2
{
    return { x = a.x + b.x, y = a.y + b.y };
}

func vec2_sub(a: vec2, b: vec2): vec2
{
    return { x = a.x - b.x, y = a.y - b.y };
}

func vec2_mul_scalar(a: vec2, s: float): vec2
{
    return { x = s*a.x, y = s*a.y };
}

func vec2_div_scalar(a: vec2, s: float): vec2
{
    rs := 1.0/s;
    return { x = rs*a.x, y = rs*a.y };
}

func vec2_mul_entries(a: vec2, b: vec2): vec2
{
    return { x = a.x*b.x, y = a.y*b.y };
}

func vec2_dot(a: vec2, b: vec2): float
{
    return a.x*b.x + a.y*b.y;
}

func vec2_unit_or_zero(v: vec2): vec2
{
    sqnorm := vec2_dot(v, v);
    if (sqnorm == 0.0) { return {}; }
    return vec2_mul_scalar(v, sqrt(sqnorm));
}

struct Vec3D
{
    union {
        struct { x, y, z: float; }
        xyz: float[3];
    }
}

typedef vec3 = Vec3D;

func vec3_from_vec2(v: vec2): vec3 { return { x = v.x, y = v.y }; }
func vec3_add_inplace(a: vec3*, b: vec3): vec3 { *a = vec3_add(*a, b); return *a; }
func vec3_sub_inplace(a: vec3*, b: vec3): vec3 { *a = vec3_sub(*a, b); return *a; }
func vec3_mul_scalar_inplace(a: vec3*, s: float): vec3 { *a = vec3_mul_scalar(*a, s); return *a; }
func vec3_div_scalar_inplace(a: vec3*, s: float): vec3 { *a = vec3_div_scalar(*a, s); return *a; }
func vec3_mul_entries_inplace(a: vec3*, b: vec3): vec3 { *a = vec3_mul_entries(*a, b); return *a; }
func vec3_cross_inplace(a: vec3*, b: vec3): vec3 { *a = vec3_cross(*a, b); return *a; }
func vec3_unit_or_zero_inplace(a: vec3*): vec3 { *a = vec3_unit_or_zero(*a); return *a; }

func vec3_add(a: vec3, b: vec3): vec3
{
    return { x = a.x + b.x, y = a.y + b.y, z = a.z + b.z };
}

func vec3_sub(a: vec3, b: vec3): vec3
{
    return { x = a.x - b.x, y = a.y - b.y, z = a.z - b.z };
}

func vec3_mul_scalar(a: vec3, s: float): vec3
{
    return { x = s*a.x, y = s*a.y, z = s*a.z };
}

func vec3_div_scalar(a: vec3, s: float): vec3
{
    rs := 1.0/s;
    return { x = rs*a.x, y = rs*a.y, z = rs*a.z };
}

func vec3_mul_entries(a: vec3, b: vec3): vec3
{
    return { x = a.x*b.x, y = a.y*b.y, z = a.z*b.z };
}

func vec3_dot(a: vec3, b: vec3): float
{
    return a.x*b.x + a.y*b.y + a.z*b.z;
}

func vec3_cross(a: vec3, b: vec3): vec3
{
    return { 
        x = a.y*b.z - a.z*b.y,
        y = -(a.x*b.z - a.z*b.x),
        z = a.x*b.y - a.y*b.y
    };
}

func vec3_unit_or_zero(v: vec3): vec3
{
    sqnorm := vec3_dot(v, v);
    if (sqnorm == 0.0) { return {}; }
    return vec3_mul_scalar(v, sqrt(sqnorm));
}

struct Vec4D
{
    union {
        struct { x, y, z, w: float; }
        xyzw: float[4];
    }
}

typedef vec4 = Vec4D;

func vec4_from_vec2(v: vec2): vec4 { return { x = v.x, y = v.y}; }
func vec4_from_vec3(v: vec3): vec4 { return { x = v.x, y = v.y, z = v.z}; }
func vec4_add_inplace(a: vec4*, b: vec4): vec4 { *a = vec4_add(*a, b); return *a; }
func vec4_sub_inplace(a: vec4*, b: vec4): vec4 { *a = vec4_sub(*a, b); return *a; }
func vec4_mul_scalar_inplace(a: vec4*, s: float): vec4 { *a = vec4_mul_scalar(*a, s); return *a; }
func vec4_div_scalar_inplace(a: vec4*, s: float): vec4 { *a = vec4_div_scalar(*a, s); return *a; }
func vec4_mul_entries_inplace(a: vec4*, b: vec4): vec4 { *a = vec4_mul_entries(*a, b); return *a; }
func vec4_unit_or_zero_inplace(a: vec4*): vec4 { *a = vec4_unit_or_zero(*a); return *a; }

func vec4_add(a: vec4, b: vec4): vec4
{
    return { x = a.x + b.x, y = a.y + b.y, z = a.z + b.z, w = a.w + b.w };
}

func vec4_sub(a: vec4, b: vec4): vec4
{
    return { x = a.x - b.x, y = a.y - b.y, z = a.z - b.z, w = a.w - b.w };
}

func vec4_mul_scalar(a: vec4, s: float): vec4
{
    return { x = s*a.x, y = s*a.y, z = s*a.z, w = s*a.w };
}

func vec4_div_scalar(a: vec4, s: float): vec4
{
    rs := 1.0/s;
    return { x = rs*a.x, y = rs*a.y, z = rs*a.z, w = rs*a.w };
}

func vec4_mul_entries(a: vec4, b: vec4): vec4
{
    return { x = a.x*b.x, y = a.y*b.y, z = a.z*b.z, w = a.w*b.w };
}

func vec4_dot(a: vec4, b: vec4): float
{
    return a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
}

func vec4_unit_or_zero(v: vec4): vec4
{
    sqnorm := vec4_dot(v, v);
    if (sqnorm == 0.0) { return {}; }
    return vec4_mul_scalar(v, sqrt(sqnorm));
}

struct Matrix3D
{
    n: float[3][3];
}
#static_assert(sizeof(Matrix3D) == 3*sizeof(vec3))
#static_assert(alignof(Matrix3D) == alignof(vec3))

typedef mat3 = Matrix3D;

func mat3_vec_at(m: mat3*, index: int): vec3*
{
    #assert(index < 3);
    vec3_ptr: vec3*;
    f := &m.n[index][0];
    memcpy(&vec3_ptr, &f, sizeof(vec3_ptr)); // @bitcast
    return vec3_ptr;
}

func vec3_mul_mat(m: mat3, v: vec3): vec3
{
    return {
        x = m.n[0][0]*v.x + m.n[1][0] + m.n[2][0],
        y = m.n[0][1]*v.x + m.n[1][1] + m.n[2][1],
        z = m.n[0][2]*v.x + m.n[1][2] + m.n[2][2],
    };
}

// Example use: normal transformation with a rotation matrix
func bivec3_mul_orthogonal_mat(av: bivec3, m: mat3): bivec3
{
    // Orthogonal(M) => Inverse(M) == Transposed(M)
    // av' = av * Transposed(M)
    // <=>
    // Transposed(av') = M*Transposed(av)
    return bivec3_from_vec3(vec3_mul_mat(m, vec3_from_bivec3(av)));
}

struct Matrix4D
{
    n: float[4][4];
}
#static_assert(sizeof(Matrix4D) == 4*sizeof(vec4))
#static_assert(alignof(Matrix4D) == alignof(vec4))

typedef mat4 = Matrix4D;

func mat4_vec_at(m: mat4*, index: int): vec4*
{
    #assert(index < 4);
    vec4_ptr: vec4*;
    f := &m.n[index][0];
    memcpy(&vec4_ptr, &f, sizeof(vec4_ptr)); // @bitcast
    return vec4_ptr;
}

func vec4_mul_mat(m: mat4, v: vec4): vec4
{
    return {
        x = m.n[0][0]*v.x + m.n[1][0] + m.n[2][0],
        y = m.n[0][1]*v.x + m.n[1][1] + m.n[2][1],
        z = m.n[0][2]*v.x + m.n[1][2] + m.n[2][2],
        z = m.n[0][3]*v.x + m.n[1][3] + m.n[2][3],
    };
}

// Points
struct Pt1D
{
    x: float;
}
typedef pt1 = Pt1D;

func pt1_add_inplace(a: pt1*, b: pt1): pt1 { *a = pt1_add(*a, b); return *a; }
func pt1_mul_scalar_inplace(a: pt1*, s: float): pt1 { *a = pt1_mul_scalar(*a, s); return *a; }
func pt1_clamp_inplace(a: pt1*, cover: aabb1): pt1 { *a = pt1_clamp(*a, cover); return *a; }

func pt1_add(a: pt1, b: pt1): pt1 { return { a.x + b.x }; }
func pt1_sub(a: pt1, b: pt1): vec1 { return a.x - b.x; }
func pt1_mul_scalar(a: pt1, s: float): pt1 { return { a.x * s }; }
func pt1_clamp(a: pt1, cover: aabb1): pt1
{
    if (a.x < cover.min) { return { cover.min }; }
    if (a.x > cover.max) { return { cover.max }; }
    return a;
}

struct Pt2D
{
    vec: vec2;
}
typedef pt2 = Pt2D;

func pt2_add_inplace(a: pt2*, b: pt2): pt2 { *a = pt2_add(*a, b); return *a; }
func pt2_mul_scalar_inplace(a: pt2*, s: float): pt2 { *a = pt2_mul_scalar(*a, s); return *a; }

func pt2_add(a: pt2, b: pt2): pt2 { return { vec = vec2_add(a.vec, b.vec) }; }
func pt2_sub(a: pt2, b: pt2): vec2 { return vec2_sub(a.vec, b.vec); }
func pt2_mul_scalar(a: pt2, s: float): pt2 { return { vec = vec2_mul_scalar(a.vec, s) }; }

struct Pt3D
{
    vec: vec3;
}
typedef pt3 = Pt3D;

func pt3_add_inplace(a: pt3*, b: pt3): pt3 { *a = pt3_add(*a, b); return *a; }
func pt3_mul_scalar_inplace(a: pt3*, s: float): pt3 { *a = pt3_mul_scalar(*a, s); return *a; }

func pt3_add(a: pt3, b: pt3): pt3 { return { vec = vec3_add(a.vec, b.vec) }; }
func pt3_sub(a: pt3, b: pt3): vec3 { return vec3_sub(a.vec, b.vec); }
func pt3_mul_scalar(a: pt3, s: float): pt3 { return { vec = vec3_mul_scalar(a.vec, s) }; }

struct Pt4D
{
    vec: vec4;
}
typedef pt4 = Pt4D;
func pt4_from_pt2(p: pt2): pt4 { return { vec = { x = p.vec.x, y = p.vec.y, w = 1.0 } }; }
func pt4_from_pt3(p: pt3): pt4 { return { vec = { x = p.vec.x, y = p.vec.y, z = p.vec.z, w = 1.0 }  }; }

func pt4_add_inplace(a: pt4*, b: pt4): pt4 { *a = pt4_add(*a, b); return *a; }
func pt4_mul_scalar_inplace(a: pt4*, s: float): pt4 { *a = pt4_mul_scalar(*a, s); return *a; }

func pt4_add(a: pt4, b: pt4): pt4 { return { vec = vec4_add(a.vec, b.vec) }; }
func pt4_sub(a: pt4, b: pt4): vec4 { return vec4_sub(a.vec, b.vec); }
func pt4_mul_scalar(a: pt4, s: float): pt4 { return { vec = vec4_mul_scalar(a.vec, s) }; }


// Antivectors


func vec3_wedge(a: vec3, b: vec3): bivec3
{
    return bivec3_from_vec3(vec3_cross(a, b));
}

func vec3_wedge_bivec3(a: vec3, b: bivec3): float
{
    return vec3_dot(a, vec3_from_bivec3(b));
}

// normals
struct Bivec3D
{
    union {
        struct { x, y, z: float; }
        xyz: float[3];
    }
}
typedef bivec3 = Bivec3D;
func bivec3_from_vec3(v: vec3): bivec3 { return { x = v.x, y = v.y, z = v.z }; }
func vec3_from_bivec3(v: bivec3): vec3 { return { x = v.x, y = v.y, z = v.z }; }

/* @todo @incomplete

// two points make a line
func pt3_wedge(a: pt3, b: pt3): bivec4
{
    return {}; // @todo @incomplete
}

// three points make a plane
func pt3_wedge3(a: pt3, b: pt3, c: pt3): trivec4
{
    return {}; // @todo @incomplete
}
*/

struct Bivec4D; // lines

struct Trivec4D; // planes