Test of chopping off bits from mantissa of golden ratio and see how much the golden point set degrades

auuu.cpp

#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstring>

void print_mantissa(double value) {
    uint64_t v;
    memcpy(&v, &value, sizeof(double));

    for(int i=0;i<52;i++){
        int bit = (v >> (51-i) & 1);
        printf("%d", bit);
    }
}

double drand48() {
    // Dilbert random number generator
    double prng = 5.0;

    return 1.0 / prng;
}

double golden_ratio(int precision_loss) {
    double au = (1.0 + sqrt(5.0)) / 2.0 - 1.0;

    // The mantissa is the last 52 bits, so we can start clearing bits from the end
    uint64_t value;
    memcpy(&value, &au, sizeof(double)); // Using memcpy to not run into UB
    uint64_t mask = ~((1ull<<precision_loss)-1); // FIXME: There's an off-by-one here right now

    value = value & mask;

    memcpy(&au, &value, sizeof(double));

    return au;
}

void golden_set(double points[][2], const unsigned int N, int precision_loss) {
    // set the initial first coordinate
    double x = drand48();
    double min = x;
    unsigned int idx = 0;
    // set the first coordinates
    for(unsigned int i = 0; i < N; ++i) {
        points[i][1] = x;
        // keep the minimum
        if(x < min)
            min = x, idx = i;
        // increment the coordinate
        x += golden_ratio(precision_loss);
        if(x >= 1) --x;
    }
    // find the first Fibonacci >= N
    unsigned int f = 1, fp = 1, parity = 0;
    while(f + fp < N) {
        unsigned int tmp = f; f += fp, fp = tmp;
        ++parity;
    }
    // set the increment and decrement
    unsigned int inc = fp, dec = f;
    if(parity & 1)
        inc = f, dec = fp;
    // permute the first coordinates
    points[0][0] = points[idx][1];
    for(unsigned int i = 1; i < N; ++i) {
        if(idx < dec) {
            idx += inc;
            if(idx >= N) idx -= dec;
        } else
            idx -= dec;
        points[i][0] = points[idx][1];
    }
    // set the initial second coordinate
    double y = drand48();
    // set the second coordinates
    for(unsigned int i = 0; i < N; ++i) {
        points[i][1] = y;
        // increment the coordinate
        y += golden_ratio(precision_loss);
        if(y >= 1) --y;
    }
}


int main() {
    unsigned int N = 20;
    auto points = new double[N][2];
    int precision_loss = 46;

    for( int i=0;i<51;i++) {
        double au = golden_ratio(i);
        printf("%d -- %1.60f -- ", i, au);

        print_mantissa(au);
        printf("\n");
    }


    golden_set(points, N, precision_loss);
 
    for (int i=0;i<N;i++){
        printf("%f %f\n", points[i][0], points[i][1]);
    }
    delete[] points;

    return 0;
}