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.

It is more robust and portable to forbid usage outside the minimum standard C requirements. Ex: if you have an int variable that may overflow, watch it so that it doesn't overflow beyond the minimal required precison from the standard even if the machine allows it. This ensure an architecture/platform independent behavior and is easier to maintain and debug. If you are performing scientific computing and want as much precision as possible, in my experience, the algorithms are typically more severe a bottleneck than machine precicion. — Sebastien, Oct 20 '13 at 14:44
Because that's the point of the exercise. To confirm that they are the same.=] — user2126752, Oct 20 '13 at 16:59

Oswald · Answer 1 · 2013-10-20T14:59:02.467

0

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.

edited Oct 20 '13 at 14:59

answered Oct 20 '13 at 14:42

Oswald

31,254
3
43
68

I made a second edit of the code. What am I missing? Why isn't the precision force by the selected TYPE? – user2126752 Oct 20 '13 at 15:24

score 0 · Accepted Answer · answered Oct 20 '13 at 14:45

0

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

answered Oct 20 '13 at 14:45

haccks

104,019
25
176
264

Thanks, fixed. Now it does finish, but single and double still give wrong machine epsilon (80-bit). How do I force 32-bit and 64-bit? And long double still give wrong answers to each one of them. – user2126752 Oct 20 '13 at 14:55
Editted it again, but it still doesn't give correct macheps for long double. – user2126752 Oct 20 '13 at 16:12
Your format specifier for `long double` is seem to be wrong. Try to change it by `%Le`. – haccks Oct 20 '13 at 16:24
The problem is with the values used for calculations, not the printed output. – user2126752 Oct 20 '13 at 17:00
If you will try to print a value with a format specifier other than that is defined for that value then it will invoke undefined behavior. – haccks Oct 20 '13 at 17:06

score 0 · Answer 3 · answered Oct 20 '13 at 20:39

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);
    }
}

Machine epsilon estimation in C

3 Answers3