Machine epsilon estimation in C

Question

I'm new to C and I'm trying to find machine epsilon (1.0 + macheps > 1.0), eta (eta > 0.0) and MAX (MAX < infinity), but my code doesn't work as intended. First of all, macheps is calculated using 80-bit precision. How do I force it to single, double and long double? Secondly, the code doesn't finish calculating the results for double precision at all.

EDIT: Fixed the format error.

/* macheps eta max */

#include <stdio.h>
#include <math.h>
#include <float.h>
#define TYPE long double

int main(void)
{
    TYPE macheps = (TYPE) 1.0;
    TYPE eta = (TYPE) 1.0;
    TYPE maksymilian = (TYPE) 2.0;
    TYPE real_macheps;
    TYPE real_eta;
    TYPE real_maksymilian;

    TYPE something = (TYPE) 1.0 + (TYPE) macheps;

    while ((TYPE) something > (TYPE) 1.0)
    {
        real_macheps = (TYPE) macheps;
        printf("%e ", (TYPE) real_macheps);
        macheps = (TYPE) macheps/(TYPE) 2.0;
        something = (TYPE) 1.0 + (TYPE) macheps;
    }
    printf("%e\n", (TYPE) real_macheps);

    while ((TYPE) eta > (TYPE) 0.0)
    {
        real_eta = (TYPE) eta;
        eta = (TYPE) eta/(TYPE) 2.0;
    }
    printf("%e\n", (TYPE) real_eta);

    while ((TYPE) maksymilian != INFINITY)
    {
        (real_maksymilian) = (TYPE) maksymilian;
        maksymilian = (TYPE) maksymilian*(TYPE) 2.0;
    }
    real_maksymilian = (TYPE) real_maksymilian * (TYPE) (2.0-(TYPE) real_macheps);
    printf("%e\n", (TYPE) real_maksymilian);
}

EDIT2: Shouldn't the above code force precision? What am I missing?

EDIT3: Still doesn't give correct macheps for long double.

Answer 1

You program will invoke undefined behavior. In C, long means long int and not long double.

Answer 2

long is shorthand for long signed int and therefore an integer type. Use long double instead.

Depending on how numbers that cannot be represented are rounded, (1.0 + machepsfloat) > 1.0 might always be true. Check FLT_ROUNDS from <float.h> to determine the rounding behaviour of your implementation.

Answer 3

Forgot to post the corrected code. I had to do a work-around, but it works:

/* macheps eta max */

#include <stdio.h>
#include <math.h>
#include <float.h>
#define TYPE long double /* printf 'e' for float and double, 'Le' for long double */
/* used_type: float 0, double 1, long double 2 */
int main(void)
{
    int used_type = 2;
    TYPE macheps = (TYPE) 1.0;
    TYPE eta = (TYPE) 1.0;
    TYPE maksymilian = (TYPE) 1.0;
    TYPE real_macheps = 0.0;
    TYPE real_eta = 0.0;
    TYPE real_maksymilian = 0.0;
    TYPE something = (TYPE) 1.0 + (TYPE) macheps;

    while ((TYPE) something > (TYPE) 1.0)
    {
        real_macheps = (TYPE) macheps;
        macheps = (TYPE) macheps/(TYPE) 2.0;
        something = (TYPE) 1.0 + (TYPE) macheps;
    }

    while ((TYPE) eta > (TYPE) 0.0)
    {
        real_eta = (TYPE) eta;
        eta = (TYPE) eta/(TYPE) 2.0;
    }

    while ((TYPE) maksymilian != INFINITY)
    {
        (real_maksymilian) = (TYPE) maksymilian;
        maksymilian = (TYPE) maksymilian*(TYPE) 2.0;
    }
    real_maksymilian = (TYPE) real_maksymilian * (TYPE) (2.0-(TYPE) real_macheps);

    if (used_type == 2)
    {
        printf("%Le\n", (TYPE) real_macheps);
        printf("%Le\n", (TYPE) real_eta);
        printf("%Le\n", (TYPE) real_maksymilian);
    }
    else
    {
        printf("%e\n", (TYPE) real_macheps);
        printf("%e\n", (TYPE) real_eta);
        printf("%e\n", (TYPE) real_maksymilian);
    }
}