When long-double calculator loses its precision?

Question

My task was to write a code for a calculator that operates on long double. It turns out that after the certain point the calculator loses precision. For example it computes correctly 9999999 ^ 2, but as far as 99999999 ^ 2 it yields the result, which is 1 bigger than it should be. I've already read something about the float precision, but all questions are concerned with decimals. I know that one can fix it by using the gmp library, but my question is:

Is there any other way to tackle this problem?
How to calculate/show the point to which a long double calculator is accurate?

main.c

int main(int argc, char *argv[])
{
    char str[100];
    char key[] = "exit";

    if (argc > 1)
    {
        switch (argv[1][1])
        {
            case 'V': case 'v':
                setValid(1);
                break;
            case 'E': case 'e':
                setError(1);
                break;
            default:
                setValid(1);
                break;
        }
    }
    else
        setValid(1);

    while (gets(str) != NULL && strcmp(key, str) != 0)
    {
        if (calcFilter(str));
        else
            printf("expression error!\n");
    }

    return 0;
}

evalexpression.c

static long double f1, f2;
static char op;

int isValidExpression(const char *str)
{
    int res;
    char ops[10];

    res = sscanf(str, "%Lf %s %Lf", &f1, ops, &f2);

    if (res == 3 && (ops[0] == '+' || ops[0] == '*' || ops[0] == '-' ||
        ops[0] == '/' || ops[0] == '^'))
    {
        op = ops[0];
        return 1;
    }
    else
        return 0;

long double getExprValue(void)
{
    switch (op)
    {
        case '+':
            return (f1+f2);
        case '*':
            return (f1*f2);
        case '-':
            return (f1-f2);
        case '/':
            return (f1/f2);
        case '^':
            return (pow(f1, f2));
        default:
            perror("something went wrong");
            return 0;
    }
}

calcfilter.c

static int errorsw = 0, validsw = 0;

int calcFilter(const char *str)
{
    if (validsw == 1 && isValidExpression(str))
    {
        printf("%s = %Lf\n", str, getExprValue());
        return 1;
    }
    else if (errorsw == 1 && !isValidExpression(str))
    {
        printf("%s\n", str);
        return 1;
    }
    else if (errorsw)
        return 1;
    else
        return 0;
}

void setError(int mode)
{
    errorsw = mode;
}

void setValid(int mode)
{
    validsw = mode;
}  

Answer 1

Using pow(f1, f2) is counter-productive as that converts the parameters to double before calling pow() and explains OP's woes. double typically has trouble beginning with 16 decimal digts as in OP's case. Recommend powl()

// return (pow(f1, f2));
return (powl(f1, f2));

---

Other notes:

long double can precisely encode numbers with LDBL_DIG decimal digits. That is at least 10 and likely 18+ on your machine.

long double, when printed needs LDBL_DECIMAL_DIG significant digits of output to distinguish that number from all other long double. To see all the significance, avoid using %Lf when trying to determine precision related issues.

printf("%.*Le", LDBL_DECIMAL_DIG - 1, x);   // "%Le"

See also Printf width specifier to maintain precision of floating-point value

---

[edit] When does floating point lose precision?

* / is easy. As long the result is not so small (sub-normal) or overflows, the answer should be expected to be with 0.5 ULP (unit in the last place).

+ - When the result is further away from zero, again lose up to 0.5 ULP. But when there is cancellation like 1.0 - 0.9999999..., nearly all precision can be lost.

z=pow(x,y) can be very imprecise when z is a large number and the function poorly written such as when the math identity is simplistically used: z = exp(log(x)*y). Otherwise a good pow() result will be within 1.0 ULP.

---

@M.M comments that an alternative to using powl() is to #include <tgmath.h> for automatic correct selection of many <math.h> functions.