2013-12-09 3 views
0

Я хотел бы установить автоматически более высокую точность моего delta_prime по сравнению с моим typedef t_float, чтобы избежать проблем с отменой, чтобы пользователь мог изменить t_float, если захочет. Может быть, я могу попытаться получить точность поплавка, но я не знаю, как я могу это сделать правильно.Автоматически устанавливать более высокую точность для float в C

В моей typedefs.h:

typedef double t_float; 

В некоторых 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; 
} 
+0

См. Http://math.stackexchange.com/questions/187242/quadratic-equation-error. – lhf

+0

Это не отвечает на мой вопрос. И в моем случае я хочу избежать b^2 = ac отмены ... – coincoin

+0

Рассматривали ли вы, если что-то вроде [MPFR] (http://www.mpfr.org/) является жизнеспособной альтернативой? –

ответ

1

Вы можете использовать шаблоны для создания "типа карты". MPL поднимать торг имеет «карту» именно для этого, или вы можете сделать это самостоятельно:

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

Затем в коде:

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

Спасибо, я сохраню это для кода на C++, но я пишу на C. – coincoin

0

Чтобы избежать большой потери точности и один не может использовать более высокую точность, можно переупорядочить расчет.

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. 
Смежные вопросы