double_bitbashing.c

#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>

// Analysis of d_osc.c from Pure Data

#define UNITBIT32 1572864.
// #define UNITBIT32 524288.
/* #define UNITBIT32 1048576. */

union fudge
{
    double d;
    uint32_t i[2];
    uint64_t u64;
};

void dump_double(double d)
{
    union fudge f;
    f.d = d;

    printf("d: %f -> 0x%08X%08X\n", d, f.i[1], f.i[0]);
    printf("d: %f -> %" PRIi64 "\n", d, f.u64);

    int bits[64];
    uint64_t ui = f.u64;
    for (int i = 0; i < 64; ++i) {
        bits[64 - i - 1] = ui % 2;
        ui /= 2;
    }
    printf("bits: ");
    for (int i = 0; i < 64; ++i) {
        if (i == 1 || i == 12) printf("|");
        printf("%d", bits[i]);
    }
    printf("\n");
    unsigned sign = (unsigned)(f.u64 >> 63);
    unsigned mantissa = (unsigned)((f.u64 >> 52) & ~(1u << 12));
    uint64_t significand = (uint64_t)(f.u64 & ((1ull << 52) - 1));
    unsigned significand_lo = (unsigned)(significand & 0xFFFFFFFF);
    unsigned significand_hi = (unsigned)(significand >> 32);
    printf("0x%X, 0x%X, 0x%X%08X\n", sign, mantissa, significand_hi, significand_lo);
    printf("========\n");
}

int main(int argc, char *argv[])
{
    // exponent is 2^20, mantissa is 0b10000...
    double d = (double)UNITBIT32;
    // double d = 1.0;

    /* dump_double(1.0); */
    /* dump_double(-1.0); */

    printf("=== UNITBIT32 ===\n");
    dump_double(d);
    
    double phase0 = 0.8;
    printf("=== phase0 ===\n");
    dump_double(phase0);

    double dphase = phase0 + (double)UNITBIT32;
    printf("=== dphase ===\n");
    dump_double(dphase);
    
    union fudge fudge;
    fudge.d = (double)UNITBIT32;
    uint32_t normhipart = fudge.i[1];
    printf("normhipart: 0x%X\n", normhipart);

    double dphase1 = dphase + 0.45;
    printf("=== dphase1 ===\n");
    dump_double(dphase1);

    fudge.d = dphase1;
    fudge.i[1] = normhipart;
    double out1 = fudge.d - UNITBIT32;
    printf("=== out1 ===\n");
    dump_double(out1);

    double dphase2 = dphase1 + 0.45;
    printf("=== dphase2 ===\n");
    dump_double(dphase2);
    fudge.d = dphase2;

    fudge.d = dphase2;
    fudge.i[1] = normhipart;
    double out2 = fudge.d - UNITBIT32;
    printf("=== out2 ===\n");
    dump_double(out2);

    
    
    return 0;
}