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

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

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

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

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

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

    double e = (eps + eps1) / 2.0;
    while (eps1 > e && e > eps) {
        assert(1. + eps <= 1. && 1. + eps1 > 1.);
        sum = 1.0 + e;
        if (sum > 1.0) {
            eps1 = e;
        } else {
            eps = e;
        }
        e = (eps + eps1) / 2.0;
    }

    assert(1. + eps <= 1. && 1. + eps1 > 1.);
    printf("machEps=%.16G\n", eps);

    return 0;
}