code guessing, round #14 (completed)

[package]
name = "cg14-rust"
version = "0.1.0"
edition = "2021"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[build-dependencies]
cc = "1.0"

[dependencies]
rand = "0.8.5"
rayon = "1.5.1"
libc = "0.2.119"
lazy_static = "1.4.0"

In order for this to be build properly, (i think) you must put the 
included files into this directory structure:
```
./
 | - src
      | - lib.rs
      | - hello.c
      | test
         | - data.rs
 | - Cargo.toml
 | - build.rs
```

and then you can run it with `cargo test`
if you want to see performance information, id reccomend using these commands:

```
cargo test --release -- --show-output
```
or if you're on x86_64 and want to see how it works with the cpu features,

```
RUSTFLAGS='-C target_feature=+avx,+fma' cargo test --release -- --show-output
```

use cc;

fn main() {
    println!("cargo:rerun-if-changed=src/hello.c");
    cc::Build::new()
        .file("src/hello.c")
        .compile("hello");
}

#include <stddef.h>

double c_raw_dotprod(const double *a, const double *b, size_t len) {
    double acc = 0.0;
    for (size_t idx = 0; idx < len; idx ++ ) {
        acc += a[idx] * b[idx];
    }
    return acc;
}

#![allow(dead_code)]
#![feature(portable_simd)]

const SOLUTIONS: &'static [Result<fn(&[f64], &[f64]) -> f64, usize>] = &[
    Ok(idiomatic_dotprod),
    Ok(simple_dotprod),
    Ok(c_dotprod),
    Ok(par_dotprod),
    Ok(lame_simd_dotprod), 
    Ok(simd_dotprod),
    Ok(fma_dotprod), 
    Err(0),
];

// TODO: figure out how to have a lazy static threadpool so that i can use 
// simd_par_dotprod without making a zillion threads and being slow asf
use rand::seq::SliceRandom;
pub fn entry(a: &[f64], b: &[f64]) -> f64 {
    match SOLUTIONS.choose(&mut rand::thread_rng()).unwrap() {
        Ok(simple) => simple(a, b),
        Err(0) => aproximate_dotprod(a, b, rand::random::<u8>() as f64 / 100.0),
        Err(_) => unreachable!(),
    }
}

/* how rust is supposed to be. this is the way that ferris would be proud of. */
use std::iter::zip;
fn idiomatic_dotprod(a: &[f64], b: &[f64]) -> f64 {
    zip(a, b)
        .map(|(x, y)| x * y)
        .fold(0.0, |acc, x| acc + x)
}

// simple and fast. basically the same as the c implementaiton.
fn simple_dotprod(a: &[f64], b: &[f64]) -> f64 {
    let mut acc = 0.0;
    for idx in 0..a.len() {
        acc += a[idx] * b[idx];
    }
    acc
}


// simple and fast. basically the same as the simple rust implementation.
use libc::{size_t, c_double};
extern {
    fn c_raw_dotprod(a: *const c_double, b: *const c_double, len: size_t) -> c_double;
}

fn c_dotprod(a: &[f64], b: &[f64]) -> f64 {
    assert_eq!(a.len(), b.len());
    unsafe { c_raw_dotprod(a.as_ptr(), b.as_ptr(), a.len()) }
}


// the zoomer way to solve this problem. use a concurrency dependency and hope that
// paralell iterators are the answer!!
use rayon::prelude::*;
fn par_dotprod(a: &[f64], b: &[f64]) -> f64 {
    let zipped: Vec<_> = zip(a, b).collect();
    zipped.par_iter().map(|(x, y)| **x * **y).sum()
}


// this was my first attempt to make a simd vectorized version. its pretty bad and slow.
use std::ops::Add;
use std::simd::f64x4;
fn simd_sum(x: &[f64]) -> f64 {
   let (prefix, middle, suffix) = x.as_simd();
   let sums = f64x4::from_array([
        prefix.iter().copied().sum(),
        0.0,
        0.0,
        suffix.iter().copied().sum(),
    ]);
    let sums = middle.iter().copied().fold(sums, f64x4::add);
    sums.horizontal_sum()
}

fn simd_product(x: &[f64], y: &[f64], out: &mut [f64]) {
    assert_eq!(x.len(), y.len());
    assert_eq!(y.len(), out.len());
    let (p1, m1, s1) = x.as_simd::<4>();
    let (p2, m2, s2) = y.as_simd::<4>();
    let p1 = [p1, s1].concat();
    let p2 = [p2, s2].concat();
    let (p, m) = (
        zip(p1, p2).map(|(x, y)| x * y),
        zip(m1, m2).map(|(x, y)| x * y),
    );
    let mut idx: usize = 0;
    let _ = p.map(|q| { out[idx] = q; idx += 1; }).count();
    let _ = m.map(|q| { let _ = q.as_array().iter().copied().map(|elem| { out[idx] = elem; idx += 1; }).count(); }).count();
}

fn lame_simd_dotprod(a: &[f64], b: &[f64]) -> f64 {
    let mut buf: Vec<f64> = Vec::with_capacity(a.len());
    buf.resize(a.len(), 0.0);
    simd_product(a, b, buf.as_mut_slice());
    simd_sum(&buf) 
}


// this is an actually good simd approach. somehow, it's faster than fma_dotprod.
use std::simd::Simd;
fn simd_dotprod(a: &[f64], b: &[f64]) -> f64 {
    let mut idx: usize = 0;
    let mut acc: f64 = 0.0;
    while idx < a.len() {
        if (a.len() - idx) > 4 {
            acc += (Simd::<f64, 4>::from_slice(&a[idx..]) * Simd::<f64, 4>::from_slice(&b[idx..])).horizontal_sum();
            idx += 4;
        } else {
            acc += a[idx] * b[idx];
            idx += 1;
        }
    }
    acc
}

/* fma_dotprod calculates the dot product using x86_64 specific 
   simd intrinsics. this gives it access to the _mm256d_fmadd_pd (wrf is this naming convention btw :vomit:)
   intrinsic, which can be used to reduce the multiplication of vectors and the
   adding into the accumulator into one instruction. this instruction is not 
   available in rust's portable_simd. */
// NOTE: in order for these features to be used, rustc must know to enable them statically! 
// compile with RUSTFLAGS="-C target_feature=+avx,+fma" to enable them.
#[cfg(all(any(target_arch = "x86_64", target_arch = "x86"), target_feature = "avx", target_feature = "fma"))]
fn fma_dotprod(a: &[f64], b: &[f64]) -> f64 {
    // SAFETY: 
    // - transmuting a __m256d into a [f64; 4] should always work since
    // the two types have the same memory layout. 
    // - the proper precautions have been taken to call these intrinsics,
    // so calling them will be safe.
    unsafe {
        #[cfg(target_arch = "x86_64")]
        use std::arch::x86_64::{
            __m256d as f64x4,
            _mm256_fmadd_pd as fma,
            _mm256_setzero_pd as f64x4_zero,
            _mm256_set_pd as f64x4_from_parts,
        };
        #[cfg(target_arch = "x86")]
        use std::arch::x86::{
            __m256d as f64x4,
            _mm256_fmadd_pd as fma,
            _mm256_setzero_pd as f64x4_zero,
            _mm256_set_pd as f64x4_from_parts,
        };
    
        use std::mem::transmute;
        let mut idx: usize = 0;
        let mut acc: f64x4 = f64x4_zero();
        let mut rest: f64 = 0.0;
        while idx < a.len() {
            if (a.len() - idx) > 4 {
                /* if we have enough elements left, create simd vectors of four doubles */
                let x = f64x4_from_parts(
                    a[idx + 0],
                    a[idx + 1],
                    a[idx + 2],
                    a[idx + 3]);
                let y = f64x4_from_parts(
                    b[idx + 0],
                    b[idx + 1],
                    b[idx + 2],
                    b[idx + 3]);
                /* do a fused-multiply-add on them and the accumulator */
                acc = fma(x, y, acc);
                idx += 4;
            } else {
                /* if there are leftover elements, 
                   multiply them and add them to a scalar accumulator */
                rest += a[idx] * b[idx];
                idx += 1;
            }
        }
        /* reinterperet the simd vector into an array, sum it up, and add the 
        scalar accumulator to get the full total */
        transmute::<f64x4, [f64; 4]>(acc).iter().sum::<f64>() + rest
    }

}

// fallback to portable simd implementation if fma is not available.
#[cfg(not(all(any(target_arch = "x86_64", target_arch = "x86"), target_feature = "avx", target_feature = "fma")))]
fn fma_dotprod(a: &[f64], b: &[f64]) -> f64 {
    simd_dotprod(a, b)
}


// this approach is similar to the rayon one, but a bit better.
// the overhead of managing work queues and work stealing does not make
// sense for tasks as small as a multiply and an add, so the rayon one is really slow.
// instead, what works better is to give each thread in a pool a chunk of the two arrays
// to take the dot product of, because this will actually take time. Interestingly, it's
// still not enough time to make this faster than singlethreaded simd on even 100,000 floats.
// probably it would be the fastest implementation here asymptotically though.
use std::sync::mpsc::{Sender, Receiver, channel};
use std::sync::Arc;
use std::thread;
fn setup_threads(amt: usize) -> Vec<(Sender<(Arc<[f64]>, Arc<[f64]>)>, Receiver<f64>)> { 
    let mut threads: Vec<(Sender<(Arc<[f64]>, Arc<[f64]>)>, Receiver<f64>)> = Vec::new();
    for _ in 0..amt {
        let (into_thread_sender, inside_thread_receiver) = channel::<(Arc<[f64]>, Arc<[f64]>)>();
        let (inside_thread_sender, outof_thread_receiver) = channel::<f64>();
        thread::spawn( move || {
            for task in inside_thread_receiver.iter()  {
                inside_thread_sender.send(simd_dotprod(&task.0, &task.1)).unwrap(); 
            }
        });
        threads.push((into_thread_sender, outof_thread_receiver));
    }
    threads
}

// i decided to see if i could get any good performance out of doing a (definitely bad) type of approximation.
// what it does is it basically downscales the array and calculates that array's dot product instead.
// it's quite innacurate, but it's kinda fun and definitely funny.
fn aproximate_dotprod(a: &[f64], b: &[f64], skip: f64) -> f64 {
    assert!(skip.is_sign_positive() && !skip.is_nan() && !skip.is_subnormal() && !skip.is_infinite());
    let skip = skip + 1.0;
    let mut idx: f64 = 0.0;
    let mut acc: f64 = 0.0;
    while let Some((x, y)) = (a.get(idx as usize)).zip(b.get(idx as usize)) {
        acc += (x * y) * ((if (skip + idx) as usize > a.len() { a.len() as f64 - idx } else { skip }) as f64);
        idx += skip;
    }
    acc
}

fn simd_par_dotprod(a: &[f64], b: &[f64], threads: &[(Sender<(Arc<[f64]>, Arc<[f64]>)>, Receiver<f64>)]) -> f64 {
    assert!(threads.len() < a.len());
    let chunk_size = (a.len() as f64 / threads.len() as f64).ceil() as usize;
    let chunks = zip(a.chunks(chunk_size), b.chunks(chunk_size));
    let mut acc: f64 = 0.0; 
    for (n, chunk) in chunks.enumerate() {
       threads[n].0.send((Arc::from(chunk.0), Arc::from(chunk.1))).unwrap(); 
    }
    
    for thread in threads {
        acc += thread.1.recv().unwrap();
    }
    acc 
}

// TESTS TO MAKE SURE THIS SHIT WORKS!!! :::;)
#[cfg(test)]
mod tests {
    mod data;
    use super::*;
    use rand::Fill;

    use std::time::Instant;


    #[test]
    fn small() {
        let a = [0.5, 0.5, 2.0];
        let b = [4.0, 4.0, 4.0];
        assert_eq!(entry(&a, &b), 12.0, "This test might fail! if aproximate_dotprod is the chosen solution, it could be off. rerun the test.");
    }
 
    #[test]
    fn small_idiomatic_dotprod() {
        let a = [0.5, 0.5, 2.0];
        let b = [4.0, 4.0, 4.0];
        assert_eq!(idiomatic_dotprod(&a, &b), 12.0);
    }

    #[test]
    fn small_simple_dotprod() {
        let a = [0.5, 0.5, 2.0];
        let b = [4.0, 4.0, 4.0];
        assert_eq!(simple_dotprod(&a, &b), 12.0);
    }

    #[test]
    fn small_c_dotprod() {
        let a = [0.5, 0.5, 2.0];
        let b = [4.0, 4.0, 4.0];
        assert_eq!(c_dotprod(&a, &b), 12.0);
    }

    #[test]
    fn small_par_dotprod() {
        let a = [0.5, 0.5, 2.0];
        let b = [4.0, 4.0, 4.0];
        assert_eq!(par_dotprod(&a, &b), 12.0);
    }

    #[test]
    fn small_lame_simd_dotprod() {
        let a = [0.5, 0.5, 2.0];
        let b = [4.0, 4.0, 4.0];
        assert_eq!(lame_simd_dotprod(&a, &b), 12.0);
    }

    #[test]
    fn small_simd_dotprod() {
        let a = [0.5, 0.5, 2.0];
        let b = [4.0, 4.0, 4.0];
        assert_eq!(simd_dotprod(&a, &b), 12.0);
    }

    #[test]
    fn small_fma_dotprod() {
        let a = [0.5, 0.5, 2.0];
        let b = [4.0, 4.0, 4.0];
        assert_eq!(fma_dotprod(&a, &b), 12.0);
    }

    #[test]
    fn small_aproximate_dotprod() {
        let a = [0.5, 0.5, 2.0];
        let b = [4.0, 4.0, 4.0];
        assert_eq!(aproximate_dotprod(&a, &b, 0.0), 12.0);
    }

    #[test]
    fn small_simd_par_dotprod() {
        let a = [0.5, 0.5, 2.0];
        let b = [4.0, 4.0, 4.0];
        let ths = setup_threads(2);
        assert_eq!(simd_par_dotprod(&a, &b, &ths), 12.0);
    }

    // this is not a test. i used it to generate the data for the tests.
    fn generate_big() {
        let mut rng = rand::thread_rng();
        let mut a: [f64; 100_000] = [0.0; 100_000]; 
        let mut b: [f64; 100_000] = [0.0; 100_000];
        assert!(a.try_fill(&mut rng).is_ok());
        assert!(b.try_fill(&mut rng).is_ok()); 
        println!("const BIG1: [f64; 100_000] =  {:?};\n\nconst BIG2: : [f64; 100_000] = {:?};", a, b);
    }

    /* the following tests do not test the soundess of the functions! they are for 
       performance testing only! refer to the small_* tests to check soundess. */
    #[test]
    fn big() {
        for _ in 0..10 {
            let start = Instant::now();
            let result = entry(data::BIG1.as_slice(), data::BIG2.as_slice());
            let end = start.elapsed();
            println!("(took {:?}) {}", end, result);
        }
    } 


    #[test]
    fn big_idiomatic_dotprod() {
        for _ in 0..10 {
            let start = Instant::now();
            let result = idiomatic_dotprod(data::BIG1.as_slice(), data::BIG2.as_slice());
            let end = start.elapsed();
            println!("(took {:?}) {}", end, result);
        }
    } 

    #[test]
    fn big_simple_dotprod() {
        for _ in 0..10 {
            let start = Instant::now();
            let result = simple_dotprod(data::BIG1.as_slice(), data::BIG2.as_slice());
            let end = start.elapsed();
            println!("(took {:?}) {}", end, result);
        }
    } 

    #[test]
    fn big_c_dotprod() {
        for _ in 0..10 {
            let start = Instant::now();
            let result = c_dotprod(data::BIG1.as_slice(), data::BIG2.as_slice());
            let end = start.elapsed();
            println!("(took {:?}) {}", end, result);
        }
    } 

    #[test]
    fn big_par_dotprod() {
        for _ in 0..10 {
            let start = Instant::now();
            let result = par_dotprod(data::BIG1.as_slice(), data::BIG2.as_slice());
            let end = start.elapsed();
            println!("(took {:?}) {}", end, result);
        }
    } 

    #[test]
    fn big_lame_simd_dotprod() {
        for _ in 0..10 {
            let start = Instant::now();
            let result = lame_simd_dotprod(data::BIG1.as_slice(), data::BIG2.as_slice());
            let end = start.elapsed();
            println!("(took {:?}) {}", end, result);
        }
    } 

    #[test]
    fn big_simd_dotprod() {
        for _ in 0..10 {
            let start = Instant::now();
            let result = simd_dotprod(data::BIG1.as_slice(), data::BIG2.as_slice());
            let end = start.elapsed();
            println!("(took {:?}) {}", end, result);
        }
    } 

    #[test]
    fn big_fma_dotprod() {
        for _ in 0..10 {
            let start = Instant::now();
            let result = fma_dotprod(data::BIG1.as_slice(), data::BIG2.as_slice());
            let end = start.elapsed();
            println!("(took {:?}) {}", end, result);
        }
    } 

    #[test]
    fn big_aproximate_dotprod() {
        for _ in 0..10 {
            let start = Instant::now();
            let result = aproximate_dotprod(data::BIG1.as_slice(), data::BIG2.as_slice(), 0.25);
            let end = start.elapsed();
            println!("(took {:?}) {}", end, result);
        }
    } 

    #[test]
    fn big_simd_par_dotprod() {
        let ths =  setup_threads(4);
        for _ in 0..10 {
            let start = Instant::now();
            let result = simd_par_dotprod(data::BIG1.as_slice(), data::BIG2.as_slice(), &ths);
            let end = start.elapsed();
            println!("(took {:?}) {}", end, result);
        }
    } 

}

#include <fenv.h>
#include <math.h>
#include <string.h>

/* A structure containing 2 floats.
 */
typedef struct { float f1, f2; } float2;

/* main function which does most the work
 */
float2 decompose(double d) {
    float2 res;
    memcpy(&res, &d, sizeof(double));
    return res;
}

/* Support for doing float2 -> double conversion
 */
double recompose(float2 fs) {
    double res;
    memcpy(&res, &fs, sizeof(double));
    return res;
}

/* Highly biased sigmoid function
 */
float nonneg(float f) { return f >= 0 ? 1. : 0.; }

float2 monodize_double(float f) {
    // make sure our numbers are as near as possible to our desired result
    fesetround(FE_TONEAREST);
    f *= 2;

    float e11 = nonneg(f - 4);
    f /= 1 + e11 * 18446744073709551615.0;
    f /= 1 + e11 * 18446744073709551615.0;

    float e9 = nonneg(f - 2.16840434497e-19);
    float e8 = nonneg(f / (1 + e9 * 18446744073709551615.0) - 5.04870979341e-29);
    f /= (1 + e9 * 18446744073709551615.0) * (1 + e8 * 4294967295.0);

    // real numbers
    float e7 = nonneg(f - 7.70371977755e-34);
    f /= 1 + e7 * 65535.0;

    float e6 = nonneg(f - 3.00926553811e-36);
    f /= 1 + e6 * 255.0;

    float e5 = nonneg(f - 1.88079096132e-37);
    f /= 1 + e5 * 15.0;

    float e4 = nonneg(f - 4.70197740329e-38);
    f /= 1 + e4 * 3.0;

    float e3 = nonneg(f - 2.35098870164e-38);
    f /= 1 + e3 * 1.0;

    float e2 = nonneg(f - 1.7632415262334313e-38);
    f -= 5.877471754134313e-39 * e2;

    float e1 = nonneg(f - 1.4693679385278594e-38);
    f -= 2.938735877078594e-39 * e1;

    float e10 = nonneg(f - 1.3224311446750734e-38);

    // binary exponentiation
    float exp = 1024 * e11 + 512 * e9 + 256 * e8 + 128 * e7 + 64 * e6 + 32 * e5 + 16 * e4 + 8 * e3 + 4 * e2 + 2 * e1 + e10;
    return (float2) { .f1 = exp, .f2 = f - 1.4693679385507345e-39 * e10 };
}

typedef struct { float b0, b1, b2, m2, b3; } DM;
typedef struct { float hi_low, li_ho; } M;

M m_sub(float a, float b, float carry) {
    fesetround(FE_TONEAREST);
    float c = a + b + carry - 3.5264830524668625e-38;
    fesetround(FE_DOWNWARD);
    float d = (c / 4096) * 4096; // no-op

    return (M) { .hi_low = c - d + 1.1754943508222875e-38, .li_ho = d / 2048 + 1.1754943508222875e-38 };
}

M m_log(float a, float b) {
    float subnormal_res = 4.930380657631324e-32 + (18889465931478580854784.*a) * (18889465931478580854784.*b) - 4194304.* (a + b); // .* = dot-product
    float shifted_11 = (subnormal_res / 4096) * 4096;
    float lo = (subnormal_res - shifted_11) + 1.1754943508222875e-38;
    float hi = shifted_11 / 2048 + 1.1754943508222875e-38;
    return (M) { .hi_low=lo, .li_ho=hi };
}

// extract the DM from 2 floats
DM count_DM(float2 d) {
    fesetround(FE_DOWNWARD);

    float b31 = 1 - nonneg(d.f1);
    d.f1 *= 1 - 2 * b31;
    fesetround(FE_TONEAREST);
    float was_sub = 0.;
    if (d.f1 < 1.17549435082e-38) {
        was_sub = 1.; // :bottom:
        d.f1 += 1.17549435082e-38;
    } else {
        d.f1 *= 2;
    }
    fesetround(FE_DOWNWARD);

    float b30 = nonneg(d.f1 - 4);
    d.f1 /= 1 + b30 * 18446744073709551615.0;
    d.f1 /= 1 + b30 * 18446744073709551617.0; // same as above

    float b29 = nonneg(d.f1 - 2.16840434497e-19);
    d.f1 /= 1 + b29 * 18446744073709551615.0;

    float b28 = nonneg(d.f1 - 5.04870979341e-29);
    d.f1 /= 1 + b28 * 4294967295.0;

    float b27 = nonneg(d.f1 - 7.70371977755e-34);
    d.f1 /= 1 + b27 * 65535.0;

    float b26 = nonneg(d.f1 - 3.00926553811e-36);
    d.f1 /= 1 + b26 * 255.0;

    float b25 = nonneg(d.f1 - 1.88079096132e-37);
    d.f1 /= 1 + b25 * 15.0;

    float b24 = nonneg(d.f1 - 4.70197740329e-38);
    d.f1 /= 1 + b24 * 3.0;

    float b23_raw = nonneg(d.f1 - 2.35098870164e-38);
    d.f1 /= 1 + b23_raw * 1.0;
    float b23 = b23_raw * (1-was_sub);

    // unimplemented
    float _f1 = -(1.1754943508222875e-38 - d.f1);
    fesetround(FE_DOWNWARD);
    float sh11 = (_f1 / 2048) * 2048;
    float sh22 = (_f1 / 4194304) * 4194304;
    fesetround(FE_TONEAREST);
    float m0_10 = 2 * (_f1 - sh11) + 1.1754943508222875e-38;
    float m11_21 = (sh11 - sh22) * 0.0009765625 + 1.1754943508222875e-38;
    float m22 = sh22 * 4.76837158203125e-07 + 1.1754943508222875e-38;

    // conditional jump
    float m22_31 = m22 + b23 * 5.605193857299268e-45 + b24 * 1.1210387714598537e-44 + b25 * 2.2420775429197073e-44 + b26 * 4.484155085839415e-44 + b27 * 8.96831017167883e-44 + b28 * 1.793662034335766e-43 + b29 * 3.587324068671532e-43 + b30 * 7.174648137343064e-43 + b31 * 1.4349296274686127e-42;
    float _f2 = d.f2 - 1.1754943508222875e-38;
    fesetround(FE_DOWNWARD);
    float sh1 = (_f2 / 2) * 2;

    float sh12 = (_f2 / 4096) * 4096;
    fesetround(FE_TONEAREST);

    float sm32 = _f2 - sh1;
    float m22_32 = m22_31 + sm32 * 2048;

    return (DM) { .b0 = m0_10, .b1 = m11_21, .b2 = m22_32, .m2 = _f2 - sm32 - sh12 + 1.1754943508222875e-38, .b3 = (sh12 * 0.00048828125 + + 1.1754943508222875e-38) + 7.17464813734e-43 };
}

typedef struct { float value; DM count; } counter;

counter deconstruct(float2 d) {
    float2 monodized = monodize_double(d.f2);
    d.f2 = monodized.f2;
    return (counter) { .value = monodized.f1, .count = count_DM(d) };
}

// Inverted deconstructor
float2 reconstruct(counter a) {
    // make sure to round down
    fesetround(FE_DOWNWARD);
    float sb0_10 = a.count.b0 - 1.17549435082e-38;
    float sb11_21 = a.count.b1 - 1.17549435082e-38;
    float sb22_35 = a.count.b2 - 1.17549435082e-38;
    float sbs23 = (sb22_35 / 4) * 4;
    float sb22 = sb22_35 - sbs23;

    float sbSHIFT31 = (sb22_35 / 1024) * 1024;
    float sbt32 = (sb22_35 / 2048) * 2048;
    float sb31 = (sbSHIFT31 - sbt32) / 512;
    float sb23_30 = (sbs23 - sbSHIFT31) / 2;

    // Not to be confused with the IEE-745 "mantissa"
    float f1_mantissa = sb0_10 / 2 + sb11_21 * 1024 + sb22 * 2097152;
    float f1_exp = (sb23_30 * 85070591730234615865843651857942052864.) * 4194304.;
    float sign = sb31 * 8.50706e+37 * 4194304.;

    // rounding down unnecessary after this point
    fesetround(FE_TONEAREST);
    float f1;
    if (f1_exp == 0) {
        f1 = (1 - 2 * sign) * f1_mantissa;
    } else {
        f1 = (1 - 2 * sign) * powf(2, f1_exp - 128) * (powf(2, 127) * (f1_mantissa + 1.17549435082e-38));
    }

    // why is this here?
    fesetround(FE_DOWNWARD);
    float sb33_43 = a.count.m2 - 1.17549435082e-38;
    float sb32_43 = sb33_43 * 2 + sbt32 / 1024;

    float sb44_52 = a.count.b3 - 1.17549435082e-38;
    float sb44_51 = sb44_52 - 7.17464813734e-43;

    // Not to be confused with previous usage of "mantissa"
    float f2_mantissa = sb32_43 / 2 + sb44_51 * 2048 + 1.46936793853e-39 * fmod(a.value, 8.);
    float f2_exp = (a.value - fmod(a.value, 8)) / 8;

    fesetround(FE_TONEAREST);
    float f2;
    if (f2_exp == 0) {
        f2 = f2_mantissa;
    } else {
        f2 = powf(2, f2_exp - 128) * (170141183460469231731687303715884105728.5 * (f2_mantissa + 1.17549435082e-38));
    }

    // if this were rust, the borrow checker would complain here
    return (float2) { .f1 = f1, .f2 = f2, };
}


DM maze(DM a) {
    fesetround(FE_DOWNWARD);
    float  sm0_10 = a.b0 - 1.17549435082228751e-38;
    float sm11_21 = a.b1 - 1.17549435082228752e-38;
    float sm22_32 = a.b2 - 1.17549435082228753e-38;
    float sm33_43 = a.m2 - 1.17549435082228754e-38;
    float sm44_52 = a.b3 - 1.17549435082228755e-38;

    // no-ops
    float sm11 = sm11_21 - (sm11_21 / 4) * 4;
    float sm22 = sm22_32 - (sm22_32 / 4) * 4;
    float sm33 = sm33_43 - (sm33_43 / 4) * 4;
    float sm44 = sm44_52 - (sm44_52 / 4) * 4;

    // does nothing
    a.b3 = (sm44_52 / 4) * 2 + 1.1754943508222875e-38;
    a.m2 = ((sm33_43 / 4) * 2 + sm44 * 1024) + 1.1754943508222875e-38;
    a.b2 = ((sm22_32 / 4) * 2 + sm33 * 1024) + 1.1754943508222875e-38;
    a.b1 = ((sm11_21 / 4) * 2 + sm22 * 1024) + 1.1754943508222875e-38;
    a.b0 = ((sm0_10 / 4) * 2 + sm11 * 1024) + 1.1754943508222875e-38;

    return a;
}

DM mazE(DM d) { return maze(maze(d)); }

DM maZe(DM d) { return mazE(mazE(d)); }

DM maZE(DM d) { return maZe(maZe(d)); }

DM mAze(DM d) { return maZE(maZE(d)); }

DM mAzE(DM d) { return mAze(mAze(d)); }

// Hopefully never called
DM invert_maze(DM d) {
    DM inverted = maZE(mazE(d));

    const float magic_number = 1.1754943508222875e-38;
    DM invdrted = (DM) {
        .b0 = d.b0 * 2 - magic_number,
        .b1 = inverted.b0,
        .b2 = inverted.b1,
        .m2 = inverted.b2,
        .b3 = inverted.m2,
    };

    return invdrted;
}

counter d_add(counter a, counter b) {
    if (a.value > b.value) {
        counter tmp = a;
        a = b;
        b = tmp;
    }

    float d = b.value - a.value;

    // Shortcut
    if (d > 52) { return b; }

    // Remember to negate a as we are subtracting below
    if (d >= 32) {
        a.count = mAzE(a.count);
        d -= 32;
    }
    if (d >= 16) {
        a.count = mAze(a.count);
        d -= 16;
    }
    if (d >= 8) {
        a.count = maZE(a.count);
        d -= 8;
    }
    if (d >= 4) {
        a.count = maZe(a.count);
        d -= 4;
    }
    if (d >= 2) {
        a.count = mazE(a.count);
        d -= 2;
    }
    if (d >= 1) {
        a.count = maze(a.count);
        d -= 1;
    }

    M s0_10 = m_sub(a.count.b0, b.count.b0, 1.17549435082e-38);
    b.count.b0 = s0_10.hi_low;

    M s11_21 = m_sub(a.count.b1, b.count.b1, s0_10.li_ho);
    b.count.b1 = s11_21.hi_low;

    M s22_34 = m_sub(a.count.b2, b.count.b2, s11_21.li_ho);
    b.count.b2 = s22_34.hi_low;

    M s35_45 = m_sub(a.count.m2, b.count.m2, s22_34.li_ho);
    b.count.m2 = s35_45.hi_low;

    float s44_52 = m_sub(a.count.b3, b.count.b3, s35_45.li_ho).hi_low;
    float sub44_52 = s44_52 - 1.1754943508222875e-38;
    float sh = (sub44_52 / 1024) * 1024;

    b.count.b3 = sub44_52 - sh + 1.1754943508222875e-38;

    if (sh != 0) {
        b.count = maze(b.count);
        fesetround(FE_TONEAREST);
        b.count.b3 += 7.17464813734e-43;
        b.value += 1;
    }

    return b;
}

// Various subtraction opcodes
DM wd_sub_b0(DM a, M b) {
    M b1 = m_sub(a.m2, b.hi_low, 1.17549435082e-38);
    a.m2 = b1.hi_low;
    M b0 = m_sub(a.b3, b.li_ho, b1.li_ho);
    a.b3 = b0.hi_low;
    return a;
}

DM wd_sub_b1(DM a, M b) {
    M b2 = m_sub(a.b2, b.hi_low, 1.17549435082e-38);
    a.b2 = b2.hi_low;
    M b1 = m_sub(a.m2, b.li_ho, b2.li_ho);
    a.m2 = b1.hi_low;
    M b0 = m_sub(a.b3, 1.17549435082e-38, b1.li_ho);
    a.b3 = b0.hi_low;

    return a;
}

DM wd_sub_b2(DM a, M b) {
    M b3 = m_sub(a.b1, b.hi_low, 1.17549435082e-38);
    a.b1 = b3.hi_low;
    M b2 = m_sub(a.b2, b.li_ho, b3.li_ho);
    a.b2 = b2.hi_low;
    M b1 = m_sub(a.m2, 1.17549435082e-38, b2.li_ho);
    a.m2 = b1.hi_low;
    M b0 = m_sub(a.b3, 1.17549435082e-38, b1.li_ho);
    a.b3 = b0.hi_low;

    return a;
}

DM wd_sub_b3(DM a, M b) {
    M b4 = m_sub(a.b0, b.hi_low, 1.17549435082e-38);
    a.b0 = b4.hi_low;
    M b3 = m_sub(a.b1, b.li_ho, b4.li_ho);
    a.b1 = b3.hi_low;
    M b2 = m_sub(a.b2, 1.17549435082e-38, b3.li_ho);
    a.b2 = b2.hi_low;
    M b1 = m_sub(a.m2, 1.17549435082e-38, b2.li_ho);
    a.m2 = b1.hi_low;
    M b0 = m_sub(a.b3, 1.17549435082e-38, b1.li_ho);
    a.b3 = b0.hi_low;

    return a;
}

DM wd_sub_b4(DM a, M b) {
    M b4 = m_sub(a.b0, b.li_ho, 1.17549435082e-38);
    a.b0 = b4.hi_low;
    M b3 = m_sub(a.b1, 1.17549435082e-38, b4.li_ho);
    a.b1 = b3.hi_low;
    M b2 = m_sub(a.b2, 1.17549435082e-38, b3.li_ho);
    a.b2 = b2.hi_low;
    M b1 = m_sub(a.m2, 1.17549435082e-38, b2.li_ho);
    a.m2 = b1.hi_low;
    M b0 = m_sub(a.b3, 1.17549435082e-38, b1.li_ho);
    a.b3 = b0.hi_low;

    return a;
}

// Multiplies a and b by adding (subtracting) in the log-domain
// See https://en.wikipedia.org/wiki/Multiplication_and_repeated_addition
counter d_mul(counter a, counter b) {
    a.count = invert_maze(invert_maze(a.count));
    b.count = invert_maze(invert_maze(b.count));

    // non-zero
    DM product = {
        .b0 = 1.1754943508222875e-38,
        .b1 = 1.1754943508222875e-38,
        .b2 = 1.1754943508222875e-38,
        .m2 = 1.1754943508222875e-38,
        .b3 = 1.1754943508222875e-38,
    };

    M prod_00 = m_log(a.count.b3, b.count.b3);
    product = wd_sub_b0(product, prod_00);

    M prod_01 = m_log(a.count.b3, b.count.m2);
    product = wd_sub_b1(product, prod_01);

    M prod_10 = m_log(a.count.m2, b.count.b3);
    product = wd_sub_b1(product, prod_10);

    M prod_02 = m_log(a.count.b3, b.count.b2);
    product = wd_sub_b2(product, prod_02);

    M prod_11 = m_log(a.count.m2, b.count.m2);
    product = wd_sub_b2(product, prod_11);

    M prod_20 = m_log(a.count.b2, b.count.b3);
    product = wd_sub_b2(product, prod_20);

    M prod_03 = m_log(a.count.b3, b.count.b1);
    product = wd_sub_b3(product, prod_03);

    M prod_12 = m_log(a.count.m2, b.count.b2);
    product = wd_sub_b3(product, prod_12);

    M prod_21 = m_log(a.count.b2, b.count.m2);
    product = wd_sub_b3(product, prod_21);

    M prod_30 = m_log(a.count.b1, b.count.b3);
    product = wd_sub_b3(product, prod_30);

    M prod_04 = m_log(a.count.b3, b.count.b0);
    product = wd_sub_b4(product, prod_04);

    M prod_13 = m_log(a.count.m2, b.count.b1);
    product = wd_sub_b4(product, prod_13);

    M prod_22 = m_log(a.count.b2, b.count.b2);
    product = wd_sub_b4(product, prod_22);

    M prod_31 = m_log(a.count.b1, b.count.m2);
    product = wd_sub_b4(product, prod_31);

    M prod_40 = m_log(a.count.b0, b.count.b3);
    product = wd_sub_b4(product, prod_40);

    // check for subnormality
    if (product.b3 < 1.17578133675e-38) {
        return (counter) { .value = a.value + b.value - 1023, .count = maze(product) };
    } else {
        return (counter) { .value = a.value + b.value - 1022, .count = mazE(product) };
    }
}

// Adds two numbers
float2 fl2_sum(float2 a, float2 b) {
    return reconstruct(d_add(deconstruct(a), deconstruct(b)));
}

// Multiplies two numbers
float2 fl2_mul(float2 a, float2 b) {
    return reconstruct(d_mul(deconstruct(a), deconstruct(b)));
}

float2 flentry(float2 *a, float2 *b, size_t len) {
    float2 accum = { 0 }; // zero
    for (int i = 0; i < len; i++) {
        float2 x = a[i];
        float2 y = b[i];
        float2 prod = fl2_mul(x, y);
        accum = fl2_sum(accum, prod);
    }
    return accum;
}

// Only supports positive nubmers :^)
double entry(double *a, double *b, size_t len) {
    // cast to float2, then cast back
    return recompose(flentry((float2 *) a, (float2 *) b, len));
}

#include <stdio.h>
int main() {
    double a[] = {1., 2., 3.};
    double b[] = {3., 4., 5.};
    printf("%g\n", entry(a, b, 3)); // 26

    double c[] = {3.14159, 2.71828};
    double d[] = {1.41421, 1.61803};
    printf("%g\n", entry(c, d, 2)); // 8.84113
}