Спустя всего несколько месяцев после исторической прогулки по Луне Нила Армстронга, ученик массачусетской школы Lexington High School Джим Сторер написал первую версию игры Lunar Landing. К 1973 году она стала самой популярной компьютерной игрой с большим отрывом от остальных. В этой простой текстовой игре вы управляете аппаратом для посадки на Луну, стремясь максимально плавно приземлиться. Всё движение происходит вертикально, каждые десять симулируемых секунд игрок решает, сколько топлива нужно сжечь.
Недавно я исследовал график оптимального сжигания топлива для наиболее мягкого приземления с максимумом оставшегося топлива. К моему удивлению, теоретически наилучшая стратегия не сработала. Игра ошибочно думает, что аппарат не коснулся поверхности, хотя на самом деле это так. Углубившись в анализ, я был потрясён сложной физикой и вычислениями в игре. В конечном итоге я обнаружил баг: отсутствующее деление пополам; похоже, почти за 55 лет никто не заметил этого.
▍ Приземление с максимумом топлива
Чтобы потратить при посадке минимум топлива, нужно спуститься за самое короткое время. Изначально мы максимизируем скорость тем, что не включаем двигатель, затем в последнюю возможную секунду жмём на полную мощность, снижая скорость до нуля как раз в момент касания поверхности. Игроки в Kerbal Space Program называют это «суицидальным сжиганием», потому что подобрать время правильно очень сложно и никаких возможностей для ошибки нет.
Путём проб и ошибок, а также двоичного поиска (проводимого вручную) можно найти график, позволяющий приземлиться на грани возможного. В течение 70 секунд мы не тратим ни капли топлива, потом 164,31426784 фунтов/с в течение следующих 10 секунд, а затем максимальные 200 фунтов/с:
Идеальным игра считает приземление при скорости менее 1 миля/ч, но мы приземляемся со скоростью более чем 3,5 мили/ч и получаем сообщение, что «можно и получше». Но если сжечь ещё на 0,00000001 больше фунтов/с, то мы полностью промахнёмся мимо поверхности, направившись вверх со скоростью 114 миль/ч:
Как мы скакнули от жёсткой посадки к взлёту без мягкой посадки между ними?
▍ Симуляция физики: один умный парнишка
Изучая игры, я ожидаю увидеть простое интегрирование Эйлера, которое часто встречается даже в современных проектах. Это процесс вычисления сил в начале каждого шага времени, далее по формуле F=ma вычисляется ускорение, потом ускорение принимается постоянным в течение шага времени 
c:
Так как в течение шага времени масса меняется, будет меняться и ускорение, поэтому то, что мы приняли его за константу — это только аппроксимация. Но к моему удивлению, Джим использовал точное решение по формуле Циолковского с разложением логарифма в ряд Тейлора. Также он использовал алгебру для упрощения вычислений и снижения величины погрешности округления. Очень впечатляющее достижение для ученика старшей школы в 1969 году! Я связался с ним и получил такой ответ:
«В то время я хорошо разбирался в математическом анализе и был знаком с такими концепциями, как ряд Тейлора; но помню я и то, что в выводе уравнений мне помогал отец, работавший физиком».
Именно благодаря формуле Циолковского становится оптимальным «суицидальное сжигание», а пять первых использованных Джимом членов ряда Тейлора, аргумент которых максимум равен 0,1212, делает формулу точной до шести знаков после запятой. Так что источник проблемы, который мы ищем, находится не в ней.
▍ Когда мы касаемся поверхности, допущения накрываются медным тазом
Формула Циолковского работает, пока мы не приземляемся. В общем случае коллизии между твёрдыми объектами — очень сложная часть движков динамики, от их реализации зависит, будет ли движок просто хорошим или отличным. Об этом я узнал, поучаствовав в разработке Open Dynamics Engine в качестве постдока в Массачусетском технологическом институте.
И в этом состоит самая большая трудность для игры Lunar Landing Game. Представьте, что аппарат снижается в начале состоящего из десяти секунд хода, но поднимается в его конце. При этом недостаточно просто проверить, будет ли он над поверхностью в обеих точках. Между началом и концом хода он мог опуститься ниже поверхности. Когда это происходит, программа должна откатиться назад и изучить более ранний момент.
Наиболее очевидный кандидат — это самая нижняя точка траектории, в которой скорость равна нулю. Оказалось, что для формулы Циолковского для этой нижней точки нет выражения в замкнутом виде, в котором бы задействовались только базовые математические функции. [Для этого нужна W-функция Ламберта, обратная 
] Поэтому вместо этого можно воспользоваться любимым приёмом физиков и взять только первые несколько членов ряда Тейлора. Если взять только первые два члена логарифма, то задачу упрощается до квадратичного уравнения, и можно воспользоваться старой доброй формулой корней квадратного уравнения из старшей школы. И аппроксимация будет достаточно точной для хода длиной десять секунд, с погрешностью около 0,1%.
Но Джим поступил иначе. В его формуле квадратный корень есть в знаменателе, но не в числителе. Кроме того, её погрешность в тридцать раз больше.
▍ Как узнать, что ты пробил дно
Какое же решение он мог придумать? Я долго смотрел на это, пытаясь придумать какой-то другой способ аппроксимации нижней точки траектории, в котором использовалось бы только сложение, вычитание, умножение, деление и квадратный корень. Если взять только первый член логарифма, то аппроксимация будет хуже, чем квадратичная, но не потребуется вычислять квадратный корень. Если взять третий член, то нам придётся решать кубическое уравнение, что в общем случае требует гораздо больше кода; похоже, в нашем случае нет никакого особого вида уравнения, имеющего простое решение. [И вряд ли Джим его придумал, потому что это достаточно сложная тема для ученика старших классов в 1969 году. Кроме того, альтернативная форма была мало известна, к тому же погрешности от приведения гораздо меньше, чем погрешности от взятия только двух членов ряда Тейлора.] Существует множество аппроксимаций 
, но в тех, где не используются ряды Тейлора, задействованы сложные функции наподобие 
, которые трудно обратить.
Но потом я внимательнее присмотрелся к его квадратному корню. Он имеет следующий вид:
Что ужасно похоже на квадратичную функцию, где деление на 4 находится внутри квадратного корня. Они должны быть как-то связаны. Но почему оно в знаменателе? Он нашёл квадратный многочлен в 
? Нет, потому что t может очень близко к нулю, поэтому его формулу нужно будет аппроксимировать в широких пределах очень больших значений, а квадратный многочлен плохо для этого подходит. Может, он вывел квадратичную аппроксимацию логарифма, а затем подставил 
, вычислил T и выполнил обратную подстановку? Изучив этот вопрос, я заново открыл для себя альтернативный вид формулы корней квадратного уравнения, в которой квадратный корень находится внизу! И в самом деле, он соответствует формуле в коде Джима.
А у этого альтернативного вида есть ключевое преимущество. Обнаружив, что мы столкнулись с землёй, нужно вернуться назад и найти время столкновения. И мы делаем то же самое, усекая ряд Тейлора, чтобы получить квадратное уравнение. В этом случае в исходном виде возникает деление на ноль, когда старший член равен нулю. Такое происходит, когда сила тяги ракетного двигателя точно уравновешивает силу гравитации. И, вероятно, такое часто случалось у многих игроков, зависавших над поверхностью и плавно снижавшихся. При этом они даже не обязаны быть абсолютно равны. Если сила тяги близка к силе гравитации, то можно получить критическую потерю значащих цифр в числителе, после чего крошечный знаменатель резко увеличит погрешность. Альтернативный вид гораздо более стабилен численно и даже работает, когда член второй степени равен нулю, то есть когда вместо квадратного у нас получается линейное уравнение. Меня потрясло, что Джим как-то самостоятельно вывел этот вид или где-то о нём узнал. Не думаю, что его можно найти в обычных учебниках не по числовым вычислениям, и вряд ли он часто применяется в физике.
▍ Давайте ещё раз проверим вывод
Но если эта формула эквивалентна, то почему погрешность аппроксимации в тридцать раз больше? Самостоятельно выведя формулу, мы получим:
Что идентично коду Джима… если не считать отсутствующую двойку в знаменателе внутри квадратного корня! Вероятно, это просто ошибка, возникшая или при выводе формулы, или при написании кода на компьютере. В конце концов, алгебраическая компьютерная система MACSYMA появилась лишь за год до этого и была недоступна в старшей школе, поэтому ему приходилось делать все вычисления ручкой на бумаге.
Из-за этого код постоянно занижает время до самой нижней точки. Джим компенсирует это двумя способами: прибавляет 0,05 с, а затем заново выполняет вычисление из этой новой, более близкой точки. И это объясняет, почему игра пропускает время посадки: первое вычисление производится, когда аппарат ещё над поверхностью и всё ещё снижается, а второе — когда он достиг низа и снова поднимается, что занимает меньше 0,05 с.
Что произойдёт, если добавить потерявшуюся двойку и избавиться от 0,05? Теперь наилучший показатель, которого можно добиться при помощи «суицидального сжигания», будет таким:
Скорость снизилась до 1,66 миль/ч, сместившись почти на три четвёртых в сторону идеального приземления при 1 миле/ч. Она неидеальна, потому что мы по-прежнему используем только первые два члена ряда Тейлора. Кроме того, когда мы решили, что самая низкая точка находится под поверхностью, нам нужно будет найти, когда мы впервые коснулись поверхности, для чего требуется схожая аппроксимация. Помогла бы ещё одна итерация, однако после устранения бага мы начали завышать время, так что нам не нужно возвращаться назад во времени, что может означать необходимость поиска другого решения квадратного уравнения. Это можно упростить, использовав только один член ряда Тейлора, как это сделано в методе Ньютона. Затем можно остановиться, когда величина скорости окажется меньше некого порогового значения, и использовать получившуюся высоту, чтобы понять, приземлились ли мы. Но для этого требуется больше работы, а в игру и так уже интересно играть, так что зачем усложнять код?
Кроме того, возможно даже совершить мягкую посадку, достаточно лишь завершить 14-й код с низкой высотой и скоростью, а затем задействовать низкую тягу на 15-м ходе, приземлившись примерно через 150 секунд. Мы не можем только достичь теоретического «суицидального сжигания» с полной тягой при приземлении, занимающего примерно 148 секунд
▍ ЦУП, мы идём на снижение с включённым двигателем
Если говорить в целом, то это очень впечатляющий результат для 18-летнего ученика старшей школы, писавшего код в 1969 году на PDP-8. Это было задолго до того, как в старшей школе начали преподавать Computer Science. В те времена были плохо известны такие аспекты числовых вычислений, как итеративное улучшение приблизительного результата при помощи метода Ньютона и проблема критической потери значащих цифр. Я узнал о них, только когда начал учиться на докторскую степень по робототехнике.
Ещё меня удивляет то, что, насколько я знаю, баг присутствовал в игре почти 55 лет, но никто его ранее не заметил. Вероятно, дело в том, что даже с багом игра была увлекательной, она и оставалась сложной, и позволяла совершить мягкую посадку. Если поставить цель не просто выиграть, а найти оптимальную стратегию, то неизбежно обнаружишь небольшие расхождения. Подозреваю, что все остальные просто радовались игре и получали от неё удовольствие.
Автор: ru_vds
