//
// Find the machine epsilon for the type "float"
//
#include <stdio.h>
#include <assert.h>

int main() {
    float eps1 = 2.0f;
    float eps = 1.0f;
    int exponent = 0;       // eps == 2^exponent
    float sum = 1.0f + eps;
    while (sum > 1.0f) {
        eps1 = eps;
        eps /= 2.0f;
        --exponent;
        sum = 1.0f + eps;
    }

    sum = 1.0f + eps;
    assert(sum <= 1.0f);

    sum = 1.0f + eps1;
    assert(sum > 1.0f);

    printf("eps=%.8G  (2 power %d)\n", eps, exponent);

    // Find the maximal epsilon, using the
    // bisection method.
    // We have
    //     1.0f + eps == 1.0f
    //     1.0f + eps1 > 1.0f
    // So the machEps is between eps and eps1
    //     eps <= e <= eps1

    float e = (eps + eps1) / 2.0f;
    while (eps1 > e && e > eps) {
        sum = 1.0f + e;
        if (sum > 1.0f) {
            eps1 = e;
        } else {
            eps = e;
        }
        e = (eps + eps1) / 2.0f;
    }

    printf("machEps=%.8G\n", eps);

    return 0;
}