06-floating1 (1238881), страница 2

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 2)

В частности следует➢ Избегать сравнения на равенство➢ Быть очень аккуратным с ошибками сложения➢ Учитывать конечный размер плавающих чисел➢ Иметь в виду, что операции над числами не всегда возвращают числа•Далее каждый из этих пунктов будет разобран детальноИзбегайте сравнивать на равенство•Следующий код не слишком сложен и сравнение должно выполняться, но увыd1 = 10.0;d2 = sqrt(d1);d3 = d2 * d2;if (d1 == d3) {// сюда мы можем и не попасть (зависит от округления)}•Правильный способ сравнивать: в пределах некоего if (fabs(d1 - d3) < tolerance) {•Хорошая ли идея выбирать tolerance == ulp(result)?Аккуратнее с ошибками сложения•В следующем примере мы пытаемся вычислить как можно более точнуюпроизводную, измельчая шаг до предела double диапазонаdouble h, cosval;for (i = 1; i < 20; ++i) {h = pow(10.0, -i);cosval = (sin(1.0 + h) - sin(1.0)) / h;printf("%d:\t%.15lf\n", i, cosval); // всё лучше и лучше?}cosval = cos(1.0);printf("True result: %.15lf\n", cosval);•Результаты при этом могут быть неожиданны (показать demo)Помните о конечности диапазона•Даже числа одинарной точности предоставляют гигантские диапазоны.

Ноконечные•Это заметно при сложении с очень большими числамиf = 16777216.0f; // 2^24nextf = f + 1.0f; // побитово не отличается от f•И при сложении с очень малымиfone = 1.0f;feps = 0.00000005f;fenext = fone + feps; // побитово не отличается от fone•И тут встаёт вопрос: а с насколько маленькими можно складывать?Ваш результат это не всегда число•Следующий код позволяет получить бесконечность за конечное времяdouble d = 1.79e+308;double infd = 2.0 * d;printf("d: %le\tinfd: %le\n", d, infd);•Дальше над результатом будут работать другие правила (напримерумножение на 0 даст не 0, а NaN)Сложение внутри диапазона•По определению FLT_EPSILON (и её аналог DBL_EPSILON) это минимальнаяконстанта, такая, что1.0f + FLT_EPSILON != 1.0f•Обычно FLT_EPSILON это число около 1e-5•Эта константа указывает сколько порядков внутри диапазона.

То есть числапорядка 1e-22 безопасно складывать с числами где-то от 1e-27 и до 1e-17.Но это не совсем линейно. Проведите эксперименты самостоятельно!•Сам диапазон задаётся константами FLT_MIN (DBL_MIN) и FLT_MAX(DBL_MAX) где имеется в виду минимум и максимум по модулюЛогарифмы для больших чисел•Сложение важно, потому что часто для больших чисел им можно заменить200!умножение. Например пусть хочется точно подсчитать190!∙10!•Определим функцию*double log_fact(int n) {double sum = 0.0;for (int i = 2; i <= n; ++i) sum += log((double)i);return sum;}•Теперь можно подсчитать непосредственноx = exp(log_fact(200) – log_fact(190) – log_fact(10));*функция log_fact несколько наивна. На практике lgamma даёт лучшую точностьОбсуждение•Допустим вы хотите найти корень уравнения 2 ∗ sin − 5 + 7 = 0•Вы знаете, что он лежит где-то в диапазоне от -3 до 3•Как вы его найдёте?Обсуждение•Допустим вы хотите найти корень уравнения 2 ∗ sin − 5 + 7 = 0•Вы знаете, что он лежит где-то в диапазоне от -3 до 3•Как вы его найдёте?•В данном случае нам повезло: −3 > 0 и 3 < 0•Для решения можно воспользоваться дихотомией: на каждом шаге делитьотрезок пополам и если там значение совпадает по знаку с левым, то братьправый интервал и наоборот•Это очень похоже на бинарный поиск в сортированном массивеProblem DH: дихотомия уравнений•Дано уравнение 2 ∗ sin − 5 + 7 = 0•Найдите делением пополам корень, лежащий в интервале от -3 до 3•Попробуйте теперь найти один из семи корней, лежщих в интервале от -13.5до 13•Творческая задача: сможете ли вы написать программу, которая находит всесемь? Как вы сделаете чтобы она работала в общем случае?•Проверить численный ответ можно здесь:https://www.wolframalpha.com/input/?i=x%5E2*sin(x)+-+5x+%2B+7+%3D+0Алгоритм SC – метод Риддерса•Метод Риддерса обычно сходится несколько быстрее, чем дихотомия(математические подробности в )typedef double (*func_t)(double x);double fsgn(double x) { return signbit(x) ? -1.0 : 1.0; }double secant(func_t f, double xleft, double xright) {assert(fsgn(f(xleft)) != fsgn(f(xright)));// далее в цикле:// xmid = (xleft + xright) / 2.0;// fl = f(xleft); fr = f(xright); fm = f(xmid);// xnew = xmid + (xmid-xleft) * fsgn(fl - fr) * fm / sqrt(fm * fm – fl * fr);// заменяем xleft = xnew или xright = xnew в зависимости от знака f(xnew)// проверяем условие выхода из цикла fabs(f(xnew)) < precisionreturn xnew;}Problem EС – исследование сходимости•Замерьте количество итераций методом Риддерса и дихотомией дляуравнения 2 ∗ sin − 5 + 7 = 0•Попробуйте разные начальные интервалы•Попробуйте float и double precision•Подтверждают ли ваши результаты теоретическое превосходство метода?Обсуждение•Рассмотрим уравнение 2 + − 0.827185 = 0•У него два действительных корня, но довольно сложно выбрать двазначения, в которых функция принимала бы разные знаки (попробуйте!)•Что делать в этом случае?Алгоритм N – метод Ньютонаstruct func_deriv { double func; double der; };// возвращает значение функции и производной в точке xtypedef struct func_deriv (*fder_t) (double x);double newton(fder_t f, double x) {// реализуйте самостоятельно +1 = −return x;}( )′ ( )Problem EN: решение методом Ньютона•Даны два уравнения 2 + − 0.827185 = 0 2 ∗ sin − 5 + 7 = 0•Решите оба используя алгоритм N•Были ли у вас какие-нибудь сложности со вторым?Проблемы со сходимостью•У метода Ньютона есть проблемы со сходимостью (картинки из )Фракталы•Несмотря на трудности которые создают проблемы со сходимостью, они жепорождают фрактальную структуру•Например рассмотрим в комплексных числах уравнение3 − 1 = 0•Для него есть такие 0 для которых метод Ньютона сходится и такие, длякоторых нет•Из-за нестабильного поведения около локальных максимумов, областьсходимости образует самоподобную кривую, то есть собственно фрактал•Рисунки таких фракталов на комплексной плоскости могут быть крайнекрасивыПодробнее: https://en.wikipedia.org/wiki/Newton_fractalПримеры для более сложных функцийИсточник: https://vk.com/fracgenProblem RI – инверсный корень•Метод Ньютона часто используется, чтобы вычислять функции•Действительно, пусть дано и надо найти =•Это всё равно, что решить уравнение =•Напишите функцию121− =0double inv_sqrt(double a);•Не используйте при реализации стандартную функцию sqrt и деление,воспользуйтесь методом Ньютона•Как вы будете тестировать вашу функцию?Problem AN – улучшение приближений•Возьмите быстрые приближения из Problem AQ и улучшите каждое из нихдополнительным шагом метода Ньютона•Существенно ли это улучшит точность по сравнению с библиотечнымифункциями?Магический инверсный корень•В качестве приближённого решения предыдущей проблемы будет работатьследующая магическая процедураfloat magic_inv_sqrt (float y) {double x2 = 0.5f * y;long i = *(long *) &y;i = 0x5f3759df - (i >> 1); // Magic!y = *(float *) &i;y = y * (1.5f - (x2 * y * y)); // one additional Newton stepreturn y;}•Оцените разницу в точности с вашим решением для Problem RI•Догадываетесь ли вы как это работает?Подробнее: https://en.wikipedia.org/wiki/Fast_inverse_square_rootЛитература•С11 ISO/IEC – "Information technology – Programminglanguages – C", 2011• ANSI/IEEE Std 754 – "Standard for Binary Floating-PointArithmetic", 1985•& Brian W.

Kernighan, Dennis Ritchie – The C programminglanguage, 1988• W. Press, S. Teukolsky – Numerical Recipes in C, 2ndedition, 1992• John D. Cook – Five Tips for Floating Point Programming,2014• James F. Blinn – Floating point tricks, 1997.

Характеристики

Список файлов лекций