sqrt (1.0 - pow (1.0,2)) возвращает -nan

Question

Я нашел интересную проблему с плавающей точкой. Я должен вычислить несколько квадратных корней в моем коде, и выражение выглядит так:

sqrt(1.0 - pow(pos,2)) 

, где позиция идет от -1,0 до 1,0 в цикле. -1.0 отлично подходит для pow, но когда pos = 1.0, я получаю -nan. Делая некоторые тесты, с помощью GCC 4.4.5 и КМК 12.0, выход

1.0 - pow(pos,2) = -1.33226763e-15 

и

1.0 - pow(1.0,2) = 0 

или

poss = 1.0 
1.0 - pow(poss,2) = 0 

Где ясно первый один будет давать проблемы, будучи отрицательным. Кто-нибудь знает, почему pow возвращает число меньше 0? Полный код обижая ниже:

int main() { 
    double n_max = 10; 
    double a = -1.0; 
    double b = 1.0; 
    int divisions = int(5 * n_max); 
    assert (!(b == a)); 

    double interval = b - a; 
    double delta_theta = interval/divisions; 
    double delta_thetaover2 = delta_theta/2.0; 
    double pos = a; 
    //for (int i = 0; i < divisions - 1; i++) { 
    for (int i = 0; i < divisions+1; i++) { 

    cout<<sqrt(1.0 - pow(pos, 2)) <<setw(20)<<pos<<endl; 

    if(isnan(sqrt(1.0 - pow(pos, 2)))){ 
     cout<<"Danger Will Robinson!"<<endl; 
     cout<< sqrt(1.0 - pow(pos,2))<<endl; 
     cout<<"pos "<<setprecision(9)<<pos<<endl; 
     cout<<"pow(pos,2) "<<setprecision(9)<<pow(pos, 2)<<endl; 
     cout<<"delta_theta "<<delta_theta<<endl; 
     cout<<"1 - pow "<< 1.0 - pow(pos,2)<<endl; 
     double poss = 1.0; 
     cout<<"1- poss "<<1.0 - pow(poss,2)<<endl; 


    } 

    pos += delta_theta; 

} 

return 0; 
} 

Answer 1

Проблема в том, что вычисления с плавающей запятой не являются точными и что 1 - 1^2 может давать небольшие отрицательные результаты, что приводит к недействительному вычислению sqrt.

Рассмотрим укупорки ваш результат:

double x = 1. - pow(pos, 2.); 
result = sqrt(x < 0 ? 0 : x); 

или

result = sqrt(abs(x) < 1e-12 ? 0 : x); 

Answer 2

Когда вы держите Инкрементирование поз в цикле, ошибки округления и накапливаются в вашем случае окончательного значения> 1,0. Вместо этого вычислите pos путем умножения на каждый раунд, чтобы получить минимальное количество ошибок округления.

Answer 3

setprecision(9) будет вызывать округление. Используйте отладчик, чтобы узнать, что это за значение. Помимо этого, по крайней мере, установите точность за пределами возможного размера используемого вами типа.

Answer 4

Вы почти всегда ошибок округления при расчете с двойниками, так как двойной тип имеет только 15 значащих десятичных цифр (52 бит) и много десятичных чисел не конвертируются в двоичные числа с плавающей запятой без округления. Стандарт IEEE содержит много усилий, чтобы сохранить эти ошибки на низком уровне, но в принципе это не всегда удается. Для подробного ознакомления см. this document

В вашем случае вы должны рассчитать pos на каждый цикл и округлить до 14 или менее цифр. Это должно дать вам чистый 0 для sqrt.

Вы можете известковы позы внутри цикла, как

pos = round(a + interval * i/divisions, 14); 

с раундом определяется как

double round(double r, int digits) 
{ 
    double multiplier = pow(digits,10); 
    return floor(r*multiplier + 0.5)/multiplier; 
}