Automatically set higher precision to float in C

Question

I would like to set automatically an higher precision to my delta_prime compared to my typedef t_float in order to avoid cancellation problems so that the user can change t_float if he wants to. Maybe I can try to get the float precision but I do not know how I can do that properly.

In my typedefs.h :

typedef double t_float;

In some code.c :

t_float solve_quadratic(const t_float a, const t_float b_prime, const t_float c)
{
    long double delta_prime = b_prime * b_prime - a * c;

    if (delta_prime < 0.0)
        return INF;

    return (b_prime + sqrt(delta_prime)) / a;
}

Answer 1

You can use templates to build a "type map". Boost's MPL has a "map" for just this, or you can do it yourself:

template <typename Numeric>
struct precision_enlarger;

template <>
struct precision_enlarger<float> { typedef double type; };

template <>
struct precision_enlarger<double> { typedef long double type; };

template <>
struct precision_enlarger<long double> { typedef long double type; };

Then in your code:

typename precision_enlarger<t_float>::type delta_prime = // ...

Answer 2

To avoid large loss of precision and one cannot use higher precision, one can re-arrange the calculation.

const t_float a, const t_float b_prime, const t_float c

// The concern here is when `a*c` is positive, no cancellation when negative.
//long double delta_prime = b_prime * b_prime - a * c;

t_float ac = a*c;
t_float Discriminant;
t_float Discriminant_sr;
if (ac <= 0.0) {
  Discriminant = b_prime * b_prime - ac;
  // or Discriminant_sq = hypot(b_prime, sqrt(-ac));
}
else {
  ac = sqrt(ac);
  // When result approaches 0.0, half the precision loss v. b_prime * b_prime - a*c
  Discriminant = (b_prime - ac) * (b_prime + ac);
  // Note: Discriminant can only be negative via this path
}

// Assume + and - root are equally valid, so use the one that does not cause cancellation.
// b_prime + sqrt(Discriminant)) / a;
Discriminant_sr = sqrt(Discriminant);
t_float Root1;
if (b_prime >= 0.0) {
  Root1 = (b_prime + Discriminant_sr) / a;
else
  Root1 = (b_prime - Discriminant_sr) / a;
return Root1;

// If the other root is needed, it likely may be calculated from _something_ like 
Root2 = c/Root1.