Скачать презентацию Магия машинных чисел Август 2008 С А Майданов Скачать презентацию Магия машинных чисел Август 2008 С А Майданов

a4d728b07e04acd327ee616212975f2b.ppt

  • Количество слайдов: 52

Магия машинных чисел Август 2008 С. А. Майданов Руководитель подразделения Intel Numerics в России Магия машинных чисел Август 2008 С. А. Майданов Руководитель подразделения Intel Numerics в России

Обзор семинара • Вместо предисловия • Незнание не освобождает от ответственности • Арифметика, доступная Обзор семинара • Вместо предисловия • Незнание не освобождает от ответственности • Арифметика, доступная на компьютере • IEEE-754 арифметика • Что должен знать каждый • Numerics – численные методы

Вместо предисловия • Выберите корректный термин, применимый к вещественным вычислениям на компьютере: 1. 2. Вместо предисловия • Выберите корректный термин, применимый к вещественным вычислениям на компьютере: 1. 2. 3. 4. 5. Арифметика с плавающей точкой Арифметика с плавающей запятой Плавающая арифметика Вещественная арифметика Компьютерная арифметика

Вместо предисловия • Численное оценивание производной функции в точке • Точность аппроксимации Вместо предисловия • Численное оценивание производной функции в точке • Точность аппроксимации

Вместо предисловия • Посмотрим на вывод программы • В чем же дело? Вместо предисловия • Посмотрим на вывод программы • В чем же дело?

Незнание не освобождает от ответственности • Два примера пренебрежения проблемами компьютерной арифметики • Сбой Незнание не освобождает от ответственности • Два примера пренебрежения проблемами компьютерной арифметики • Сбой системы наведения боеголовки ракеты Patriot – – 25 февраля 1991 года во время боевых действий в Персидском заливе американская ракета Patriot не смогла перехватить иракскую ракету, в результате чего та поразила одну из казарм армии США. Погибло 28 человек. Из официального заключения: «Результатом сбоя системы наведения послужила программная ошибка – неточное вычисление времени, прошедшего с момента инициализации системы. Небольшая погрешность представления 1/10 секунды в двоичной арифметике с плавающей запятой накопилась до 0. 34 секунды после 100 часов нахождения ракеты на боевом дежурстве, что не позволило надежно рассчитать траекторию перехвата. . . » • Взрыв ракетоносителя Ariane 5 – – 4 июля 1996 года Европейское космическое агентство запустило ракетоноситель Ariane 5, который взорвался вскоре после старта. Это был первый старт после почти 10 -летней подготовки миссии, которая обошлась Агентству порядка 7 миллиардов долларов. Сама ракета и груз, находившийся на ее борту оцениваются порядка 500 миллионов долларов. Из официального заключения: «Результатом крушения ракетоносителя послужила программная ошибка в одном из вспомогательных устройств. 64 битное число с плавающей запятой, используемое для представления горизонтальной компоненты скорости ракеты относительно стратовой платформы преобразовывалось к 16 -битному знаковому целому. Результирующее число оказалось больше 32767. . . » http: //www. ima. umn. edu/~arnold/455. f 96/disasters. html

Незнание не освобождает от ответственности • О чем это шутят и по сей день? Незнание не освобождает от ответственности • О чем это шутят и по сей день? • “United we stay, divided we fall” • “Life is like a box of chocolates” Ошибка в инструкции FDIV процессора Intel® Pentium Pro, стоившая компании 475 млн. долларов

Арифметика, доступная на компьютере • 0, 1, 2, 3, . . . • Беззнаковые Арифметика, доступная на компьютере • 0, 1, 2, 3, . . . • Беззнаковые целые с операциями +, -, *, / • … -3, -2, -1, 0, +1, +2, +3, … • Знаковые целые с операциями +, -, *, /

Арифметика, доступная на компьютере • 0, 1, 2, 3, . . . • Беззнаковые Арифметика, доступная на компьютере • 0, 1, 2, 3, . . . • Беззнаковые целые с операциями +, -, *, / • … -3, -2, -1, 0, +1, +2, +3, … • Знаковые целые с операциями +, -, *, / • Сложение, вычитание, умножение точные • Если не возникает переполнения • Граница переполнения определяется размером типа • Операция деления a/b возвращает частное • Остаток r может быть вычислен – r = a-b*c – r=a%b – Что эффективнее?

Арифметика, доступная на компьютере • Можно ли представить дроби, имея в арсенале лишь целочисленную Арифметика, доступная на компьютере • Можно ли представить дроби, имея в арсенале лишь целочисленную арифметику? 1. Да 2. Нет 3. Да, предварительно помолясь 4. Что такое дробь?

Арифметика, доступная на компьютере • Можно ли представить дроби, имея в арсенале лишь целочисленную Арифметика, доступная на компьютере • Можно ли представить дроби, имея в арсенале лишь целочисленную арифметику? 1. Да 2. Нет 3. Да, предварительно помолясь 4. Что такое дробь? Представлять дробь парой целых (числитель и знаменатель) • r 1=a/b, r 2=c/d • r 1+r 2=(a*d+c*b)/(b*d) • r 1*r 2=(a*c)/(b*d) • r 1/r 2=(a*d)/(b*c) Непрактично

Арифметика, доступная на компьютере • Можно ли представить дроби, имея в арсенале лишь целочисленную Арифметика, доступная на компьютере • Можно ли представить дроби, имея в арсенале лишь целочисленную арифметику? 1. Да 2. Нет 3. Да, предварительно помолясь 4. Что такое дробь? Арифметика с фиксированной запятой • Ввести шкалирующий множитель 2 -n (или с другим показателем, например 10) • r 1=2 -n*a, r 2=2 -n*b • r 1+r 2=2 -n*(a+b) • r 1*r 2=2 -n*((a*b)>>n) • r 1/r 2=2 -n*((a/b)<

Арифметика, доступная на компьютере • Можно ли представить дроби, имея в арсенале лишь целочисленную Арифметика, доступная на компьютере • Можно ли представить дроби, имея в арсенале лишь целочисленную арифметику? 1. Да 2. Нет 3. Да, предварительно помолясь 4. Что такое дробь? Арифметика с плавающей запятой • Число представляется парой целых (шкалирующий множитель и мантисса) • r 1=2 n*a, r 2=2 m*b Основа современных компьютерных систем с поддержкой плавающей запятой

Арифметика, доступная на компьютере • Истоки арифметики с плавающей запятой прослеживаются вплоть до 1800 Арифметика, доступная на компьютере • Истоки арифметики с плавающей запятой прослеживаются вплоть до 1800 г. до н. э. (Вавилон) • Применяли арифметические операции с плавающей запятой по основанию 60, но не имели обозначения для показателя. Нужный показатель подразумевался лицом, производившим вычисления. • По-крайней мере в одном случае обнаружен ошибочный ответ при выполнении сложения из-за неправильного выравнивания операндов

Арифметика, доступная на компьютере • IEEE-754 – это 1. Прородитель стандарта IEEE-1394 Fire. Wire Арифметика, доступная на компьютере • IEEE-754 – это 1. Прородитель стандарта IEEE-1394 Fire. Wire 2. Результат проекта KCS, не только разделивший мир компьютерной арифметики на несколько лагерей, но и объединивший многих крупных игроков на рынке процессоров • Kahan-Сooner-Stone 3. Рингтон Интела в России по мотивам песни бременских музыкантов (E—E-E повторяется 754 раза) 4. Кодовое название следующего поколения i. Phone от Apple по цене $754, который будет поддерживать стандарт 3 E 5. Ни одно из вышеперечисленного Стандарт IEEE-754 (Арифметика с плавающей запятой) утвержден в 1985 г. в результате титанических усилий проф. Кахана, его единомышленников и учеников, Интел и других ведущих производителей микропроцессоров

Арифметика, доступная на компьютере • Хаос, царивший в вычислениях с плавающей запятой до IEEE-754 Арифметика, доступная на компьютере • Хаос, царивший в вычислениях с плавающей запятой до IEEE-754 • X*1. 0 -> переполнение • X*1. 0 -> сокращение 4 младших битов • X-Y=0, когда X<>Y • X==Y и X-Y==0. 0 дают разные ответы • Невозможно определить ситуацию 0/0 или переполнение Проф. Кахан принял активное участие в разработке со-процессора i 8087, полностью удовлетворяющего спецификациям IEEE-754 За выдающийся вклад по стандартизации вычислений с плавающей запятой проф. Кахан был награжден премией Тюринга в 1989 г.

Арифметика, доступная на компьютере • Слева направо • Проф. Вильям Кахан (William Kahan), Stanford Арифметика, доступная на компьютере • Слева направо • Проф. Вильям Кахан (William Kahan), Stanford University • Проф. Джером Кунен (Jerome Coonen), U. C. Berkley (ученик Кахана) • Питер Танг (Peter Tang), Intel (ученик Кахана)

IEEE-754 арифметика • Число с плавающей запятой • X=(-1)s*2 k*b – s {0, 1} IEEE-754 арифметика • Число с плавающей запятой • X=(-1)s*2 k*b – s {0, 1} – Emin<=k<=Emax – 1<=b<2

IEEE-754 арифметика • Формат представления чисел с плавающей запятой IEEE-754 арифметика • Формат представления чисел с плавающей запятой

IEEE-754 арифметика • Формат представления чисел с плавающей запятой (SP, DP) • Поле экспоненты IEEE-754 арифметика • Формат представления чисел с плавающей запятой (SP, DP) • Поле экспоненты e=k+bias • Поле мантиссы f - дробная часть

IEEE-754 арифметика • Формат представления чисел с плавающей запятой (SP, DP) • Emin-1 зарезервировано IEEE-754 арифметика • Формат представления чисел с плавающей запятой (SP, DP) • Emin-1 зарезервировано для представления – +/-0 – Денормализованных чисел • Emax+1 зарезервировано для представления – +/- – Na. N (Not-a-Number) • 0. 0 f • Поля s=0, e=0, f=0 • -1. 0 f • Поля s=1, e=127, f=0 • 1. 25 f • Поля s=0, e=127, f=0. 25 (2 -23*010… 0 b=2 -2) 23

IEEE-754 арифметика • 0. 1 f представляется как 1. s=0, e=124, f=0. 6 2. IEEE-754 арифметика • 0. 1 f представляется как 1. s=0, e=124, f=0. 6 2. s=0, e=123, f=1. 6 3. s=0, e=-4, f=0. 6 4. s=0, e=123, f=0. 599999904632568359375 5. s=0, e=124, f=1. 599999904632568359375 0. 1 f представляется в машине числом 1. 599999904632568359375/16= = 0. 0999999940395355224609375

IEEE-754 арифметика • Бесконечности • • + : s=0, e=255, f=0 – : s=1, IEEE-754 арифметика • Бесконечности • • + : s=0, e=255, f=0 – : s=1, e=255. f=0 1. 0/+0. 0 = + Na. N • • s=0, e=255, f<>0 0. 0/0. 0 = Na. N

Что должен знать каждый • IEEE-754 совместимая аппаратура реализует операции +, -, *, /, Что должен знать каждый • IEEE-754 совместимая аппаратура реализует операции +, -, *, /, sqrt с точностью до последнего бита (корректно округленные) • Относительная ошибка округления не превышает • = 2 -24 ~= 6*10 -8 для SP • = 2 -53 ~= 1. 1*10 -16 для DP • Более чем достаточно для подавляющего числа задач • Нужно ли программисту прикладывать дополнительные усилия при такой точной аппаратуре?

Что должен знать каждый float s; int i; s=0. 0 f; _controlfp(_PC_24, _MCW_PC); for( Что должен знать каждый float s; int i; s=0. 0 f; _controlfp(_PC_24, _MCW_PC); for( i=0; i<(1<<25); i++ ) { s += 1. 0 f; } printf(“s=%. 0 fn”, s); • Каков будет ответ машины? 1. 224=16777216 2. 225=33554432 3. +Inf 4. Na. N Со временем разность порядков s и 1. 0 f становится настолько велика, что добавление 1. 0 f к s не меняет суммы

Что должен знать каждый • Более реалистичный пример – Метод Монте-Карло • Суммируем n Что должен знать каждый • Более реалистичный пример – Метод Монте-Карло • Суммируем n случайных чисел u, равномерно распределенных на интервале (0, 1) • По центральной предельной теореме s ~ N(n/2, n/12) • Возьмем n=226 • • Ответ машины 16777216 • Pr(s<=16777216) = Ожидаемый ответ 33554432 3 с вероятностью 99. 7% Со временем разность порядков s и u увеличивается, что приводит к постепенной потери младших битов u при суммировании. Начиная с некоторого момента, она становится настолько велика, что добавление u к s не меняет суммы.

Что должен знать каждый • Вычисление корней квадратного уравнения • • a = 1. Что должен знать каждый • Вычисление корней квадратного уравнения • • a = 1. 222919 b = 3. 340000 c = 2. 280527 Дискриминант d = 8. 1*10 -7>0

Что должен знать каждый • Вычисление корней квадратного уравнения • • • a = Что должен знать каждый • Вычисление корней квадратного уравнения • • • a = 1. 222919 b = 3. 340000 c = 2. 280527 Дискриминант d = 8. 1*10 -7>0 Машинное представление • • = 1. 222918987 = 3. 339999914 = 2. 280527115 Дискриминант d = -2. 1*10 -7<0 Неточность машинного представления параметров задачи может качественно изменить решение

Что должен знать каждый • Вычисление корней квадратного уравнения float a b c r Что должен знать каждый • Вычисление корней квадратного уравнения float a b c r 1, r 2, a, b, c, d, s; = 0. 6623223423957825 f; = 3. 3399999141693115 f; = 4. 2107892036437988 f; d = b*b - 4. 0 f*a*c; s = sqrtf(d); r 1 = (-b + s)/a*0. 5 f; r 2 = (-b - s)/a*0. 5 f; printf("%. 16 f, %. 16 fn", r 1, r 2); • Ожидаемый ответ • r 1=-2. 520984188, r 2=-2. 521877422

Что должен знать каждый • Вычисление корней квадратного уравнения float a b c r Что должен знать каждый • Вычисление корней квадратного уравнения float a b c r 1, r 2, a, b, c, d, s; = 0. 6623223423957825 f; = 3. 3399999141693115 f; = 4. 2107892036437988 f; d = b*b - 4. 0 f*a*c; s = sqrtf(d); r 1 = (-b + s)/a*0. 5 f; r 2 = (-b - s)/a*0. 5 f; printf("%. 16 f, %. 16 fn", r 1, r 2); • Ожидаемый ответ • • r 1=-2. 520984188, r 2=-2. 521877422 Вычисленный • r 1=-2. 520693540, r 2=-2. 522167921

Что должен знать каждый • Теорема: Пусть и ненулевые числа одного знака, тогда справедливо Что должен знать каждый • Теорема: Пусть и ненулевые числа одного знака, тогда справедливо неравенство для относительной погрешности: • (< > – < >)≤| + |/| - |*max( (< >), (< >)) • В случае SP max( (< >), (< >)) ≤ ~ 6*10 -8 • В примере детерминант оказывается близок к нулю • • • = b*b; = 4*a*c ~ ~ 11. 16 - ~ 3. 5*10 -7 | + |/| - | ~ 22. 3/ 3. 5*10 -7 ~ 6*108 (< > – < >)≤ 1

Что должен знать каждый • Способы побороть катастрофическую потерю точности 1. 2. 3. 4. Что должен знать каждый • Способы побороть катастрофическую потерю точности 1. 2. 3. 4. 5. Перечитать IEEE-754 и помолиться Перечитать IEEE-754 и застрелиться Авось, сойдет и так Ваш вариант ответа

Что должен знать каждый • Способы побороть катастрофическую потерю точности • Повысить точность формата Что должен знать каждый • Способы побороть катастрофическую потерю точности • Повысить точность формата SP->DP->QP

Что должен знать каждый • Способы побороть катастрофическую потерю точности • • Повысить точность Что должен знать каждый • Способы побороть катастрофическую потерю точности • • Повысить точность формата SP->DP->QP Эмулировать расширенную мантиссу

Что должен знать каждый • Способы побороть катастрофическую потерю точности • • • Повысить Что должен знать каждый • Способы побороть катастрофическую потерю точности • • • Повысить точность формата SP->DP->QP Эмулировать расширенную мантиссу Реорганизовать выражение (формулы приведения, разложения, . . . )

Что должен знать каждый • Вычисление • • • Способ 1: x*x-y*y Способ 2: Что должен знать каждый • Вычисление • • • Способ 1: x*x-y*y Способ 2: (x-y)*(x+y) Какой из способов даст более точный ответ? 1. Способ 1 2. Способ 2 3. Оба способа равноценны

Что должен знать каждый • Способ 1: x*x-y*y • • = x*x; = y*y Что должен знать каждый • Способ 1: x*x-y*y • • = x*x; = y*y | + |/| - | - катастрофическая потеря точности, если x близко к y (по модулю)

Что должен знать каждый • Способ 1: x*x-y*y • • • = x*x; = Что должен знать каждый • Способ 1: x*x-y*y • • • = x*x; = y*y | + |/| - | - катастрофическая потеря точности, если x близко к y (по модулю) Способ 2: (x-y)*(x+y) • Операции –, +, * привнесут ошибки 1, 2 и 3 соответственно – • (x-y)*(1+ 1)*(x+y)*(1+ 2)*(1+ 3), где i не превысят по модулю Относительная погрешность – – Не превысит (1+ )3 Первый порядок аппроксимации ~3 Запись выражения в тождественной форме зачастую позволяет избежать проблемы катастрофической потери точности

Что должен знать каждый • Вычисление корней квадратного уравнения • Случай b*b>>4*|a*c| • Катастрофическая Что должен знать каждый • Вычисление корней квадратного уравнения • Случай b*b>>4*|a*c| • Катастрофическая потеря точности • • • b>0: в r 1 при вычислении b<0: в r 2 при вычислении Как победить проблему? • • Повысить точность формата Обратиться к специалисту Помолиться еще раз Посмотреть следующий слайд

Что должен знать каждый • Вычисление корней квадратного уравнения • Случай b*b>>4*|a*c| • Катастрофическая Что должен знать каждый • Вычисление корней квадратного уравнения • Случай b*b>>4*|a*c| • Катастрофическая потеря точности • • • b>0: в r 1 при вычислении b<0: в r 2 при вычислении Как победить проблему? • • b>0: Домножить числитель и знаменатель r 1 на b<0: Домножить числитель и знаменатель r 2 на Запись выражения в тождественной форме зачастую позволяет избежать проблемы катастрофической потери точности

Что должен знать каждый • Если не удается найти «хорошее» тождественное преобразование. . . Что должен знать каждый • Если не удается найти «хорошее» тождественное преобразование. . . • А всегда ли нужно прибегать к повышению точности формата или эмулированию расширеной мантиссы, чтобы результат IEEE-754 операции был точным? • Некоторые полезные свойства IEEE-754 арифметики • Пусть x=2 k*b, 1<=b<2, m>0. Тогда <-2 k+m> вытеснит (занулит) m младших битов x • • =x*2 k • Пусть ½<=x/y<=2, тогда =x-y Пусть по крайнтей мере p/2 младших битов мантисы чисел x, y нулевые, тогда =x*y В некоторых случаях IEEE-754 операции точны

Что должен знать каждый • Эмулирование расширенной мантиссы • Точное сложение/вычитание – z=<x+y> – Что должен знать каждый • Эмулирование расширенной мантиссы • Точное сложение/вычитание – z= – z+zz=x+y • Теорема Деккера • • Для любых IEEE-754 чисел x и y, где |x|>=|y|, коррекционный член zz может быть вычислен как zz: => Теорема Мюллера-Кнута • Для любых IEEE-754 чисел x и y коррекционный член zz может быть вычислен как w: =, v: =, z 1: =, z 2: =, zz: =

Что должен знать каждый • Эмулирование расширенной мантиссы • Точное умножение – z=<x*y> – Что должен знать каждый • Эмулирование расширенной мантиссы • Точное умножение – z= – z+zz=x*y • Алгоритм Деккера • Получить из x и y числа xh, yh, где ~p/2 младших битов нулевые • • • xl : = , yl : = p : = xh*yh, q : = xh*yl+xl*yh z : = p+q, zz : = ((p-z)+q)+xl*yl

Что должен знать каждый • Вернемся к корням квадратного уравнения • Как избежать катастрофической Что должен знать каждый • Вернемся к корням квадратного уравнения • Как избежать катастрофической потери точности в b*b-4*a*c • b=bh+bl, a=ah+al, c=ch+cl • b 2 h = bh*bh, b 2 l = 2*bh*bl+bl*bl ach = ah*ch, acl = ah*cl+al*ch+al*cl d 1 = (b 2 h-4*ach) d 2 = (b 2 l-4*acl) • d = d 1+d 2

Что должен знать каждый • Формула Герона (площадь треугольника) Что должен знать каждый • Формула Герона (площадь треугольника)

Что должен знать каждый • Формула Герона (площадь треугольника) • Когда возможна катастрофическая потеря Что должен знать каждый • Формула Герона (площадь треугольника) • Когда возможна катастрофическая потеря точности?

Что должен знать каждый • Формула Герона (площадь треугольника) • Когда возможна катастрофическая потеря Что должен знать каждый • Формула Герона (площадь треугольника) • Когда возможна катастрофическая потеря точности? • Один из углов очень острый – Например, a ~= b+c • Как избежать катастрофической потери точности? c b a

Что должен знать каждый • Формула Герона (площадь треугольника) • Когда возможна катастрофическая потеря Что должен знать каждый • Формула Герона (площадь треугольника) • Когда возможна катастрофическая потеря точности? • Один из углов очень острый – Например, a ~= b+c • Как избежать катастрофической потери точности? • Положим a>b>c • P 1: =a+(b+c) P 2: =c-(a-b) P 3: =c+(a-b) P 4: =a+(b-c) S: =0. 25*sqrt(P 1*P 2*P 3*P 4) c b a

Что должен знать каждый • Численное оценивание производной функции в точке • Точность аппроксимации Что должен знать каждый • Численное оценивание производной функции в точке • Точность аппроксимации Помните об особенностях IEEE-754 арифметики

Numerics – Численные методы • Разработчики Intel • 80% Россия, 20% США • Эксперты Numerics – Численные методы • Разработчики Intel • 80% Россия, 20% США • Эксперты в областях • Арифметика с плавающей запятой • Стандарты математических библиотек • Статистический анализ и моделирование – В приложениях биоинженерии и финансовой математики • Микроархитектурная оптимизация и анализ производительности • Поставщик математических компонентов в программные продукты Интел • Intel® compilers: LIBM, SVML (Short Vector Math Library) • Intel® Math Kernel Library: VML (Vector Math Library), VSL (Vector Statistical Library) • Intel® Integrated Performance Primitives: ipp. VM (Vector Math) • Intel® Summary Statistics Library • Консультационные проекты по численным и статистическим проблемам • Научно-исследовательские проекты с университетами США и России Миссия команды Развивать и расширять экосистему вычислительных методов путем распространения эксперизы по вопросам вычислений с плавающей точкой и статистического анализа нашим клиентам, а также повышения конкурентноспособности платформ и программных продуктов Интел посредством разработки математически ориетированных решений

Numerics – Численные методы QP Method/Parameter Inaccuracies Dominate DP LA SP HA Apps not Numerics – Численные методы QP Method/Parameter Inaccuracies Dominate DP LA SP HA Apps not very demanding for DP accuracy EP SP LA Apps not very demanding for accuracy SP EP Media apps Meets accuracy requirements of most apps DP HA Rigorous accuracy specs DP CR Bitwise Identical Results Позиционирование точностей математических библиотек acy cur Ac & • DP LA is sufficient for majority of apps on i • DP HA is used sometimes to meet (at times cis e artificial) accuracy specs for customer’s Pr benchmarks/acceptance tests ter • DP CR and/or higher precision (Quad Precision) rea may address bitwise compatible results issue in G certain customer apps • SP HA, SP LA, DP EP are targeted to apps where math function inaccuracies don’t dominate method/parameter inaccuracies, e. g. Monte Carlo simulations • SP EP is targeted to class of media/graphics apps

Numerics – Численные методы • Кто нам нужен • Целеустремленный, инициативный, самостоятельный. . . Numerics – Численные методы • Кто нам нужен • Целеустремленный, инициативный, самостоятельный. . . • Математик и программист. . . • Боготворящий численный и функциональный анализ • Способный к самостоятельному обучению • Готовый к открытию неизведанного. . . в области численных методов • Знающий, как написать работу и представить ее на конференции • Командный игрок