deadbraincoral's stats

def entry(grid, x, y):
    """
    Fill the region at (x, y) in image with value true.
    """
    if grid[x][y]:
        return

    # Fill this cell and any adjacent true cells.
    grid[x][y] = True

    if x > 0:
        entry(grid, x - 1, y)
    if x < len(grid) - 1:
        entry(grid, x + 1, y)
    if y > 0:
        entry(grid, x, y - 1)
    if y < len(grid[x]) - 1:
        entry(grid, x, y + 1)

#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
}

def entry(list): return sorted(list)