Когда я пишу эту статью, то чувствую себя довольно глупо. На самом деле, это история с моралью «Прежде, чем действовать, изучи вопрос и понимай, в чём заключается твоя цель, потому что иначе потеряешь много времени».
Я продолжаю работать над проектом PSRayTracing. Как ни стараюсь я положить его на полку, время от времени слышу о чём-то «новом» и задаюсь вопросом: «а можно ли засунуть это в мой трассировщик лучей, чтобы выжать из него ещё немного скорости?». На этот раз такой темой стали аппроксимации Паде. Моя цель заключалась в обеспечении более быстрых (и точных) тригонометрических аппроксимаций.
Увы, это не помогло... однако я обнаружил нечто иное, позволившее существенно ускорить мой трассировщик!
Более быстрая тригонометрия
Тригонометрические функции часто применяются во всех графических приложениях. Но с точки зрения времени вычислений они могут быть немного затратными. Конечно, точность — это здорово, но нас обычно заботит скорость. Поэтому если мы можем найти «достаточно хорошую» аппроксимацию, которая быстрее реальной функции, то обычно вполне приемлемо её использовать.
Для текстурирования объектов в PSRayTracing (основанном на информации из книг Ray Tracing in One Weekend редакции приблизительно 2020 года) используется функция std::asin(). При профилировании некоторых сцен я заметил, что к этой функции выполняется существенное количество вызовов, поэтому решил, что неплохо было бы найти её аппроксимацию.
В конечном итоге я написал собственную аппроксимацию на основе рядов Тейлора. Она быстрее, но есть у неё и изъян: когда x на входе меньше -0.8 или больше 0.8, у неё начинаются сильные отклонения. Поэтому для обеспечения корректности мне пришлось при выходе за эти границы возвращаться к использованию std::asin().
Код на C++ выглядит так:
double _asin_approx_private(const double x) { // Здесь используется аппроксимация на основе ряда Тейлора. // См. http://mathforum.org/library/drmath/view/54137.html // // В случае, когда x=[-1, -0.8) или (0.8, 1.0], к сожалению, возникают большие погрешности // по сравнению с arcsin, поэтому в этом случае мы используем реальную функцию if ((x < -0.8) || (x > 0.8)) { return std::asin(x); } // Аппроксимация рядом Тейлора constexpr double a = 0.5; constexpr double b = a * 0.75; constexpr double c = b * (5.0 / 6.0); constexpr double d = c * (7.0 / 8.0); const double aa = (x * x * x) / 3.0; const double bb = (x * x * x * x * x) / 5.0; const double cc = (x * x * x * x * x * x * x) / 7.0; const double dd = (x * x * x * x * x * x * x * x * x) / 9.0; return x + (a * aa) + (b * bb) + (c * cc) + (d * dd); }
Путём проб и ошибок я нашёл ряд Тейлора четвёртого порядка, обеспечивающий наибольшую производительность на моём «железе». Он был осязаемо быстрее (на +5%), поэтому я оставил его и перешёл к следующей оптимизации.
Аппроксимации Паде
Не помню, откуда я о них узнал, у меня не сохранилось никаких воспоминаний. Подробнее об этих аппроксимациях можно прочитать в статье Википедии. Если вкратце, то это математический инструмент, который может помочь с аппроксимацией функции. Для его вычисления необходимо начать с ряда Тейлора (или Маклорена). Хоть PSRayTracing в основном пишется на C++, для упрощения будет использоваться Python; когда это будет важно, мы вернёмся к нашему любимому компилируемому языку.
Для описанной выше аппроксимации арксинуса на основе ряда Тейлора код на Python будет таким:
def taylor_fourth_order(x: float) -> float: return x + (x**3)/6 + (3*x**5)/40 + (5*x**7)/112
Выполнив преобразования, мы получим то, что называется аппроксимацией Паде [3/4]:
def asin_pade_3_4(x): a1 = -367.0 / 714.0 b1 = -81.0 / 119.0 b2 = 183.0 / 4760.0 n = 1.0 + (a1 * x**2) d = 1.0 + (b1 * x**2) + (b2 * x**4) return x * (n / d)
Я также покажу аппроксимации на основе ряда Тейлора пятого порядка, то есть аппроксимацию Паде [5/4]:
def asin_pade_5_4(x): a1 = -1709.0 / 2196.0 a2 = 69049.0 / 922320.0 b1 = -2075.0 / 2196.0 b2 = 1075.0 / 6832.0 n = 1.0 + (a1 * x**2) + (a2 * x**4) d = 1.0 + (b1 * x**2) + (b2 * x**4) return x * (n / d)
Составив график всех трёх, получим следующее:
Ого, это уже выглядит намного лучше. Меньше погрешностей! Это не так легко разглядеть, поэтому давайте увеличим правую часть функций:
Погрешность не исчезла полностью, но стала намного меньше. Вместо того, чтобы откатываться к встроенному методу asin(), давайте вытянем из рукава ещё одну карту: применение обратных тригонометрических функций/преобразований с половинным углом. Взгляните:
Может, они кажутся немного непонятными, зато позволяют нам сделать хитрый ход. Когда |x| больше значения, мы можем «телепортироваться» с края функции ближе к центру arcsin(), выполнить вычисления, а затем вернуться назад и использовать результат там. Вот, как выглядит наша новая аппроксимация asin(x) на Python:
def asin_pade_3_4_half_angle_correction(x: float) -> float: abs_x = abs(x) # Если мы вышли за пределы диапазона, то можно использовать преобразования с половинным углом, чтобы учесть погрешность if abs_x > 0.85: small_x = math.sqrt(0.5 * (1.0 - abs_x)) r = (math.pi / 2) - (2.0 * asin_pade_3_4(small_x)) return -r if x < 0 else r # В пределах диапазона мы можем, как обычно, использовать аппроксимацию 3/4 return asin_pade_3_4(x)
С корректировкой края выглядят примерно так:
Возможно, это не так легко разглядеть, но пунктирные линии — это график с методикой коррекции на основе преобразования с половинным углом; он очень близок к прямой y=0. Небольшую погрешность можно разглядеть, только приблизив график.
Можно использовать и дополнительную оптимизацию: применить (и адаптировать) аппроксимацию Паде [1/2] внутри тела if, потому что small_x всегда будет меньше квадратного корня 0.075 (равного ~0.27). Аппроксимация Паде [1/2] для asin() вычисляется гораздо быстрее, но только для малых значений x. Чтобы оптимизировать код ещё больше, её даже можно встроить в функцию:
def asin_pade_1_2(x): b1 = -1.0 / 6.0 d = 1.0 + (b1 * x**2) return (x / d) # ... def asin_pade_3_4_half_angle_correction(x: float) -> float: abs_x = abs(x) # Если мы вышли за пределы диапазона, то можно использовать преобразования с половинным углом, чтобы учесть погрешность, а также "малую" аппроксимацию Паде if abs_x > 0.85: z = (1.0 - abs_x) / 2 b1 = -1.0 / 6.0 d = 1.0 + (b1 * z) small_pade = math.sqrt(z) / d r = (math.pi / 2) - (2.0 * small_pade) return -r if x < 0 else r # В пределах диапазона мы можем, как обычно, использовать аппроксимацию 3/4 return asin_pade_3_4(x)
Результат всё равно выглядел, как график выше, поэтому я не думаю, что необходимо добавлять ещё одну. В коде на C++ это выглядит так:
constexpr double HalfPi = 1.5707963267948966; inline double asin_pade_3_4(const double x) { constexpr double a1 = -367.0 / 714.0; constexpr double b1 = -81.0 / 119.0; constexpr double b2 = 183.0 / 4760.0; const double x2 = x * x; const double n = 1.0 + (a1 * x2); const double d = 1.0 + (b1 * x2) + (b2 * x2 * x2); return x * (n / d); } double asin_pade_3_4_half_angle_correction(const double x) { const double abs_x = std::abs(x); if (abs_x <= 0.85) { return asin_pade_3_4(x); } else { // Края кривой Паде const double z = 0.5 * (1.0 - abs_x); constexpr double b1 = -1.0 / 6.0; const double d = 1.0 + (b1 * z); const double pade_result = std::sqrt(z) / d; const double r = HalfPi - (2.0 * pade_result); return std::copysign(r, x); } }
По сравнению с исходным способом аппроксимации он более сложен, но обладает и преимуществами:
В большом диапазоне значений [-0.85, 0.85] будут использоваться упрощённые вычисления
В пограничных случаях он использует более быстрые вычисления, чем
std::asin()Меньше погрешность
Давайте подставим этот код в PSRayTracing и посмотрим, повысилась ли скорость. Мы используем стандартную сцену (финальный рендер из книги 2):
Измерения
asin() используется в модели земного шара. Все изображения обычно выглядят одинаково (за исключением небольшого шума). Для тестирования мы будем выполнять рендеринг изображения 1080p с 250 сэмплами на пиксель и задействуем несколько ядер. Тестирование выполнялось на M4 Mac Mini (в macOS Tahoe с компилятором GCC15 и флагом оптимизации -O3). Выполнив по несколько прогонов, возьмём медиану:
С std::asin() для рендеринга сцены потребовалось примерно 111 секунд:
ben@Benjamins-Mac-mini build_gcc15 % ./PSRayTracing -j 4 -n 250 -s 1920x1080 -o render_std_asin.png Scene: book2::final_scene Render size: 1920x1080 Samples per pixel: 250 Max number of ray bounces: 50 Number of render threads: 4 Copy per thread: on Saving to: render_std_asin.png Seed: `ASDF` Rendering: [==================================================] 100% 111s Render took 110.891 seconds
Старой аппроксимации asin() потребовалось примерно 105 секунд (ускорение на ~5%):
ben@Benjamins-Mac-mini build_gcc15 % ./PSRayTracing -j 4 -n 250 -s 1920x1080 - o render_asin_taylor.png ... Render took 104.674 seconds
А это новое решение с аппроксимацией Паде выполнялось... то же количество времени:
ben@Benjamins-Mac-mini build_gcc_15 % ./PSRayTracing -j 4 -n 250 -s 1920x1080 -o render_asin_pade.png ... Render took 104.87 seconds
Новый метод asin() лучше (с точки зрения корректности), но различий в производительности практически нет.
Промпт LLM
Невозможно спорить, что ИИ-инструментарий сильно изменился за последние три с лишним года. На CppCon 2023 Андрей Александреску выступил с замечательным докладом «Роботы хотят забрать у вас работу: исследование генеративного ИИ для C++», рекомендую его послушать. После всей описанной выше работы и размышлений я решил спросить у LLM (в данном случае Gemini): «Какую быструю аппроксимацию asin(x) можно использовать в C++?». Он выдал мне следующее:
constexpr double HalfPi = 3.1415926535897932385 / 2.0; double fast_asin_cg(const double x) { // Исходные коэффициенты Minimax constexpr double a0 = 1.5707288; constexpr double a1 = -0.2121144; constexpr double a2 = 0.0742610; constexpr double a3 = -0.0187293; // Убираем знак const double abs_x = fabs(x); // Вычисляем многочлен по схеме Горнера double p = a3 * abs_x + a2; p = p * abs_x + a1; p = p * abs_x + a0; // Применяем sqrt и смещение на pi/2 const double x_diff = sqrt(1.0 - abs_x); const double result = HalfPi - (x_diff * p); // Возвращаем знак return copysign(result, x); }
Слишком просто, чтобы быть правдой. В качестве источника ИИ указал на документацию по asin в Nvidia Cg Toolkit; этот продукт не обновлялся с 2012 года. Эталонная реализация содержала вышеприведённый код в переписанном виде. Можно с лёгкостью портировать его в Python/Jupyter, чтобы проверить:
def asin_cg(x: float) -> float: ''' Fast branchless asin(x) approximation. Based on Abramowitz and Stegun formula 4.4.45 ''' # https://developer.download.nvidia.com/cg/asin.html # https://personal.math.ubc.ca/~cbm/aands/page\_81.htm # Исходные коэффициенты Minimax Абрамовица и Стиган a0 = 1.5707288 a1 = -0.2121144 a2 = 0.0742610 a3 = -0.0187293 abs_x = abs(x) # Вычисляем многочлен по схеме Горнера p = a3 * abs_x + a2 p = p * abs_x + a1 p = p * abs_x + a0 result = (math.pi / 2) - math.sqrt(1.0 - abs_x) * p # Нативным образом восстанавливаем знак return math.copysign(result, x)
Я поверить не мог, настолько всё это оказалось чистым и изящным. Реализация, погрешность и вывод. Взгляните сами:
Кривая идеально накладывается на функцию arcsin(). А погрешность практически отсутствует. Однако настоящей проверкой будет сам трассировщик лучей:
ben@Benjamins-Mac-mini build_gcc_15_new_asin_cg % ./PSRayTracing -j 4 -n ... Render took 101.462 seconds
Ничего себе! Это существенно быстрее, чем все остальные методики. А при сравнении рендера с результатом использования std::asin() они оказались неразличимы. Мы нашли более качественную реализацию asin().
Ещё измерения
Это заставило меня выполнить бенчмаркинг этой реализации на нескольких чипах и операционных системах.
Intel i7-10750H, Ubuntu 24.04 (с GCC 14 и clang 19):
./test_gcc_O3 "ASDF" "100" "10000000" std::asin() time: 29197.9 ms asin_cg() time: 19839.8 ms Verification sums: std::asin(): -34549.5 asin_cg(): -34551.1 Difference: 1.60886 Error: -0.00465669 % Speedup: 1.47169x faster ./test_clang_O3 "ASDF" "100" "10000000" std::asin() time: 29520.7 ms asin_cg() time: 19044.3 ms ... Speedup: 1.55011x faster
Intel i7-10750H, Windows 11 (с MSVC 2022):
C:\Users\Benjamin\Projects\PSRayTracing\experiments\asin_cg_approx>test_msvc_O2.exe ASDF 100 10000000 std::asin() time: 12458.1 ms asin_cg() time: 6562.1 ms ... Speedup: 1.8985x faster
Apple M4, macOS Tahoe (с GCC 15 в Homebrew и clang 17):
./test_gcc_O3 "ASDF" "100" "10000000" std::asin() time: 10469 ms asin_cg() time: 10251 ms ... Speedup: 1.02126x faster ./test_clang_O3 "ASDF" "100" "10000000" std::asin() time: 12650 ms asin_cg() time: 12073.2 ms ... Speedup: 1.04777x faster
Во всех случаях лидером оказывалась CG-аппроксимация asin(). На чипе Intel он быстрее на крайне малую долю. Мне было бы любопытно протестировать её на системе x86_64 с процессором AMD, но это упражнение я оставлю читателям. Предполагаю, что там эта аппроксимация проявит себя столь же замечательно. На чипе Apple M4 разрыв не такой уж большой, но всё равно заметный (и воспроизводимый). Любое различие больше 2% важно; я имею в виду старый доклад Николаса Ормрода.
Выводы
Думаю, изначально я пошёл по пути рядов Тейлора, потому что пробовал их с sin() и cos(), а потом предположил, что их можно применить и к другим тригонометрическим функциям. Я не догадался сразу посмотреть, не решил ли кто-то уже мою задачу нахождения более быстрого арксинуса в компьютерной графике.
Хуже всего то, что решение существовало ещё до появления LLM. Сейчас у меня не получается воспроизвести запрос, но я ввёл в поисковый движок что-то вроде "fast c++ asin approximation cg". Первым результатом была ссылка на страницу документации Nvidia Cg Toolkit. Я обнаружил её только несколько дней назад.
Меня удивило, что мне никто об этом не сказал. Я даже указал ускоренный asin() в README, как моё достижение, но никто и не подумал меня исправить... Я знаю, что этот проект (и статьи) известны в кругах разработчиков на C++ и создателей компьютерной графики. Даже более опытные и умные люди ни слова мне не сказали.
Этот потрясающий фрагмент кода ждал меня в документации к мёртвому ПО, которое, в свою очередь, позаимствовало исходную формулу из учебника по математике из 60-х. Ещё меня расстроило то, что при попытке поисков не обнаружилось никаких бенчмарков. Надеюсь, теперь эта информация станет известнее.
Думаю, моя главная проблема заключается в том, что я ни разу не подумал остановиться, задуматься о своих истинных целях и поискать, не достиг ли их уже кто-то ещё. Вот, что мне дал этот опыт.
Ну и ещё красивые графики.
