Реализуем с нуля функцию косинуса на языке C

Я изучил, как реализовать функцию косинуса при помощи нескольких разных подходов. Одна из реализаций почти в три раза быстрее, чем math.h, но придётся смириться с точностью до четырёх знаков после запятой.

Задавались ли вы когда-нибудь вопросом, как в математической библиотеке вашего любимого языка программирования реализованы тригонометрические функции, например, косинус? Это настолько популярная функция, что её можно встретить в каждой математической библиотеке, поэтому реализация должна быть довольно простой, ведь так? Ну уж нет. Почти совершенно точно, что это не так.

Моё исследование началось с того, что мой друг и коллега Стивен Марц работал над ядром операционной системы и я предложил, чтобы он отрисовал на экране функцию косинуса. Я часто использую косинус в качестве «hello, world» для графических приложений. Возникла проблема: его ядро не задействовало стандартную библиотеку C (а значит, прощай math.h!), а целевой платформой являлась архитектура RISC-V (а значит, никаких подобий команды fcos Intel!).

Так началось моё долгое приключение.

Стоит сказать, что я не математик и не гуру программирования систем. На самом деле, за десять лет изучения в колледже математики и computer science я каким-то образом так и не прошёл курс тригонометрии. Поэтому я расскажу о своих трудностях изучения и реализации косинуса на C. Мои цели были следующими:

• Достаточно простая реализация, чтобы даже я мог понять код.
• Достаточно точная реализация, чтобы ни один смертный не заметил погрешности.
• Достаточно быстрая, чтобы её можно было использовать в видеоигре.

Мы изучим множество способов вычисления косинуса и несколько оптимизаций. Так как я не мастер языка C, то буду избегать сложных трюков и микрооптимизаций (но если вам они известны, то напишите мне!). На самом деле, я часто буду делать выбор в пользу более читаемого кода вместо чуть более производительного. В процессе работы мы будем проводить бенчмарки, чтобы понимать последствия своих действий. Весь код и бенчмарки можно найти в репозитории GitHub.

▍ Почему мне важен косинус?

За всё время кодинга косинус я использовал только в одной ситуации: игры, игры, игры.

function move(obj) {
    obj.x += obj.speed * Math.sin(obj.rotation);
    obj.y += obj.speed * Math.cos(obj.rotation);
}

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

Визуально это выглядит так:


Возможно, не стоит смотреть слишком долго…

▍ Способы вычисления косинуса

Когда я решил забраться в эту потустороннюю кроличью нору, то обнаружил способ аппроксимации косинуса при помощи ряда Тейлора.


Похоже, это стандартный способ вычисления косинуса в математических уровнях, по крайней мере, в высокоуровневых. Чем больше членов, тем точнее будет аппроксимация.

Мне в голову пришла мысль об использовании таблицы поиска. Это массив предварительно вычисленных значений, в котором можно искать ближайшее значение косинуса для введённых данных. Во времена, когда вычислительные ресурсы были ограниченными, такой способ использовался довольно часто. Мне не удалось найти на GitHub проектов, в которых для тригонометрических функций используется таблица, но уверен, что они существуют.

Также в моих поисках часто встречался термин CORDIC. Это итеративный способ вычисления функций наподобие синуса и косинуса при помощи только сложения, вычитания, битового сдвига и небольшой таблицы поиска. Его часто реализовывали аппаратно, начиная с конца 1950-х, или в ПО, которое работало в маломощных CPU или микроконтроллерах, например, применяемых в калькуляторах. Этот способ был довольно популярен в прошлом, он использовался в математическом сопроцессоре Intel 8087, калькуляторах TI-85 и HP-48G. Однако я не смог найти ссылок на использование его сегодня. Подробнее о нём можно прочитать в статье Википедии или в исходной научной статье с описанием способа; также можно изучить его реализацию на C. Я не буду сравнивать свои способы с ним, но мне немного любопытно, насколько хорошо он себя проявляет. Вот иллюстрация из исходной статьи:


Также есть команда fcos процессоров Intel. Эта команда вычисляет косинус по значению с плавающей запятой в радианах в интервале от -2^63 до 2^63, сохраняемому в регистре FPU, и заменяет это значение результатом. Не знаю, используется ли по-прежнему в современных процессорах Intel способ CORDIC, а если да, то реализован ли он аппаратно или программно. Дизассемблировав несколько программ, использующих косинус, я не смог найти ни одной, задействующей команду fcos. Хотя она быстра, была задокументирована неточность этой команды, которая наиболее заметна, когда аргумент приближается к числам, кратным pi/2. Подробности можно прочитать в неофициальной документации или в официальной справке Intel.

Надо сказать, что любой человек в здравом уме и памяти воспользуется math.h языка C. И для повседневного использования, вероятно, так и нужно поступать. Но помните, что мы не можем использовать ничего из стандартной библиотеки. К тому же это было бы неинтересно. Я сравнил точность функции косинуса из math.h на своём компьютере с Wolfram Alpha. Выяснилось, что math.h точна до 15 знаков после запятой, а это больше, чем мне когда-либо понадобится.

Чтобы понять, как косинус вычисляет стандартная библиотека, я взглянул на исходники множества реализаций стандартной библиотеки C. Glibc, Newlibc, Musl и так далее. И хотя, похоже, во всех них использовался ряд Тейлора, мне сложно было во всех них разобраться. Все они сильно отличались друг от друга, часто использовали множество компактных функций, применяли магические числа, имели таблицы предварительно вычисленных значений и активно задействовали трюки с битами. Разработчики потратили много времени на то, чтобы сделать их быстрыми.

Вот скриншот того, как я изучал соответствующий код Musl. От cos() к __cos(), а потом к __rem_pio2().


Уф…

▍ Подготовка бенчмарка

В процессе реализации различных способов вычисления косинуса я буду сравнивать их по двум параметрам: времени выполнения и точности. Для проверки времени выполнения каждая функция вычисляется 100 миллионов раз с интервалом входных значений, время исполнения замеряется при помощи функции clock файла time.h. Для проверки точности я беру разность результатов моей функции и результатов math.h для всего интервала входных значений, и возвращаю наихудшее значение. Например, точность 0.0002 означает, что в наихудшем случае моя реализация для одного входного значения среди большого интервала значений на 0.0002 отличалась от результатов math.h.

Вот соответствующий код:

// Замеряем время выполнения в секундах. Чем меньше число, тем лучше.
double runtime(double (*func)(double))
{
    clock_t start = clock();
    for (int i = 0; i < 100000000; i++)
    {
        volatile double c = func(i / 10000.0);
        (void)c;
    }
    return (clock() - start) / (double)CLOCKS_PER_SEC;
}

// Находим наихудший случай точности по сравнению с math.h. Чем меньше, тем лучше.
double accuracy(double (*func)(double))
{
    double w = -1;
    double start = 0;
    double stop = CONST_2PI;
    double step = 0.0000001;
    for (double i = start; i < stop; i += step)
    {
        double c = absd(func(i) - cos(i));
        if (c > w)
        {
            w = c;
        }
    }
    return w;
}

Бенчмарк компилировался при помощи clang 11.0.3 и запускался на 13-дюймовом MacBook Pro 2018 с CPU i7 на 2,7 ГГц и 16 ГБ ОЗУ.

Весь код бенчмарков можно найти в репозитории GitHub. Благодарю доктора Марца за то, что он переписал его, повысив удобство интерфейса.

▍ Первая попытка: наивный ряд Тейлора с литералами

Сначала я попробовал трансляцию ряда Тейлора с четырьмя членами:

double cos_taylor_literal_4terms_naive(double x) {
    return 1 - ((x*x)/(2)) + ((x*x*x*x)/(24)) - ((x*x*x*x*x*x)/(720)) + ((x*x*x*x*x*x*x*x)/(40320));
}

Он казался на удивление точным, когда я начал тестировать его со входными значениями наподобие 0.1 и 0.235. Мой энтузиазм угас, когда я сравнил его график с графиком math.h:


Фиолетовая линия — это math.h, а зелёная — моя функция. Между -pi и +pi она выглядит вполне неплохо, но потом начинаются резкие отклонения.

Возможно, поможет увеличение количества членов.

double cos_taylor_literal_6terms_naive(double x) {
    return 1 - ((x*x)/(2)) + ((x*x*x*x)/(24)) - ((x*x*x*x*x*x)/(720)) + ((x*x*x*x*x*x*x*x)/(40320)) - ((x*x*x*x*x*x*x*x*x*x)/(3628800)) + ((x*x*x*x*x*x*x*x*x*x*x*x)/(479001600)); 
}

Попробуем нарисовать график…


Гораздо лучше. Но всё равно при приближении к 2pi начинаются сильные погрешности.

▍ Вторая попытка: улучшенный ряд Тейлора с литералами

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

Первым делом я попробовал уменьшить интервал. Чем больше входное значение, тем менее точен этот способ. Так как косинус повторяется каждые 2pi, нам достаточно просто выполнить x = x % (2*pi);. Однако в C оператор деления с остатком не работает с числами с плавающей запятой, поэтому мы создали собственный.

#define modd(x, y) ((x) - (int)((x) / (y)) * (y))

double cos_taylor_literal_6terms_2pi(double x)
{
    x = modd(x, CONST_2PI);
    return 1 - ((x * x) / (2)) + ((x * x * x * x) / (24)) - ((x * x * x * x * x * x) / (720)) + ((x * x * x * x * x * x * x * x) / (40320)) - ((x * x * x * x * x * x * x * x * x * x) / (3628800)) + ((x * x * x * x * x * x * x * x * x * x * x * x) / (479001600));
}

Он показывает себя лучше для значений в пределах 2pi, но всё равно очень неточен при значениях, близких к 2pi. Можно улучшить ситуацию, потому что значения косинуса эквивалентны при каждом числе, кратном pi, меняется только знак. Мы можем выполнить mod на 2pi, вычесть pi, если значение больше pi, вычислить ряд Тейлора, а затем добавить нужный знак. То есть на самом деле мы вычисляем косинус только от 0 до pi.

double cos_taylor_literal_6terms_pi(double x)
{
    x = modd(x, CONST_2PI);
    char sign = 1;
    if (x > CONST_PI)
    {
        x -= CONST_PI;
        sign = -1;
    }
    return sign * (1 - ((x * x) / (2)) + ((x * x * x * x) / (24)) - ((x * x * x * x * x * x) / (720)) + ((x * x * x * x * x * x * x * x) / (40320)) - ((x * x * x * x * x * x * x * x * x * x) / (3628800)) + ((x * x * x * x * x * x * x * x * x * x * x * x) / (479001600)));
}

Это помогло существенно улучшить точность при близких к 2pi входных значениях.

Следующая оптимизация включает в себя устранение некоторых избыточных вычислений. Можно заметить, что в коде часто используется x*x. Я всего лишь уменьшил количество умножений, добавив double xx = x * x;.

double cos_taylor_literal_6terms(double x)
{
    x = modd(x, CONST_2PI);
    char sign = 1;
    if (x > CONST_PI)
    {
        x -= CONST_PI;
        sign = -1;
    }
    double xx = x * x;

    return sign * (1 - ((xx) / (2)) + ((xx * xx) / (24)) - ((xx * xx * xx) / (720)) + ((xx * xx * xx * xx) / (40320)) - ((xx * xx * xx * xx * xx) / (3628800)) + ((xx * xx * xx * xx * xx * xx) / (479001600)));
}

Производительность существенно улучшилась! Также я попробовал double xxxx = xx * xx;, но не увидел особой разницы, поэтому двинулся дальше.

Я всё ещё не знал точно, сколько членов использовать, поэтому попробовал до десяти членов, чтобы проверить, как это улучшит точность:

double cos_taylor_literal_10terms(double x)
{
    x = modd(x, CONST_2PI);
    char sign = 1;
    if (x > CONST_PI)
    {
        x -= CONST_PI;
        sign = -1;
    }
    double xx = x * x;

    return sign * (1 - ((xx) / (2)) + ((xx * xx) / (24)) - ((xx * xx * xx) / (720)) + ((xx * xx * xx * xx) / (40320)) - ((xx * xx * xx * xx * xx) / (3628800)) + ((xx * xx * xx * xx * xx * xx) / (479001600)) - ((xx * xx * xx * xx * xx * xx * xx) / (87178291200)) + ((xx * xx * xx * xx * xx * xx * xx * xx) / (20922789888000)) - ((xx * xx * xx * xx * xx * xx * xx * xx * xx) / (6402373705728000)) + ((xx * xx * xx * xx * xx * xx * xx * xx * xx * xx) / (2432902008176640000)));
}

После этого график с десятью членами начал полностью накладываться на math.h. Прогресс! При сравнении наихудшего случая точности с math.h четыре члена были ужасны. Шесть членов отличались на 0.0001, а это больше точности, чем мне нужно. Десять членов отклонялись всего на 0.00000000007. Ура!

Однако увеличение количества членов приводит к резкому росту времени исполнения. Судя по бенчмарку, наивное решение с четырьмя членами занимает всего 0,4 секунды, с шестью членами — 0,94 секунды, а с десятью — 1,46 секунды. А math.h требуется всего 1,04 секунды при большей точности.

Добавление новых членов при такой методике показалось мне довольно нелепым, поэтому я двинулся дальше.

▍ Третья попытка: ряд Тейлора со скользящим произведением

Когда я показал свои результаты доктору Марцу, он провернул какую-то алгебраическую магию и отправил мне свою модифицированную версию. В его методике устраняется большой объём избыточных вычислений благодаря сохранению результатов, а также появилась возможность указания количества нужных членов. В некоторых областях применения может быть удобно иметь параметр, управляющий степенью точности/скоростью.

double cos_taylor_running_yterms(double x, int y)
{
    int div = (int)(x / CONST_PI);
    x = x - (div * CONST_PI);
    char sign = 1;
    if (div % 2 != 0)
        sign = -1;

    double result = 1.0;
    double inter = 1.0;
    double num = x * x;
    for (int i = 1; i <= y; i++)
    {
        double comp = 2.0 * i;
        double den = comp * (comp - 1.0);
        inter *= num / den;
        if (i % 2 == 0)
            result += inter;
        else
            result -= inter;
    }
    return sign * result;
}

Ради бенчмарков я не использовал эту версию со вторым параметром. Вместо этого я дублировал функцию и прописал в цикле постоянные значения (например, 6, 10 и 16).

При добавлении новых членов точность возрастает крайне несущественно. При 16 членах наихудший случай точности отличается на 0.0000000000000009, однако бенчмарк времени исполнения занимает 2,57 секунды. Это ни в коем случае не медленно… если не сравнивать с math.h.

▍ Четвёртая попытка: таблица поиска

Мне хотелось попробовать способ с таблицей поиска. Смысл заключался в предварительном вычислении множества значений и записи их в массив. Таблицы поиска существовали задолго до компьютеров, поэтому методика не нова. В данном случае я надеялся, что увеличив объём используемой памяти, я получу большой рост скорости и обеспечу при этом достаточную степень точности.

Для генерации таблицы поиска доктор Марц написал скрипт на Python, генерирующий файл заголовка C, содержащий массив, в котором каждый элемент — это значение косинуса, вычисленное с помощью math.h. Очень умно!

from math import cos, pi

def main(f, PRECISION, NAME):
    f.write("double %s[] = {\n" % NAME)
    j = 0
    p = 0.0
    while True:
        f.write("{:.20f}, ".format(cos(p)))
        j += 1
        p += PRECISION
        if p > 2*pi:
            break
    f.write("1.0 };\n")
    f.write("const int %s_size = %d;\n" % (NAME, j+1))


if __name__ == '__main__':
    main(open("costable_1.h", "w"), 1.0, "costable_1")
    main(open("costable_0_1.h", "w"), 0.1, "costable_0_1")
    main(open("costable_0_01.h", "w"), 0.01, "costable_0_01")
    main(open("costable_0_001.h", "w"), 0.001, "costable_0_001")
    main(open("costable_0_0001.h", "w"), 0.0001, "costable_0_0001")

Мы хотели протестировать таблицы с разной точностью, поэтому сгенерировали таблицы с 8 значениями, 64 значениями, 630 значениями, 6285 значениями и 62833 значениями. Расплачиваться за это пришлось ценой увеличения размера исполняемого файла. Влияние таблиц 1.0 и 0.1 незаметно, однако остальные таблицы увеличили размер исполняемого файла примерно на 5 КБ, 50 КБ и 500 КБ.

double absd(double a) { *((unsigned long *)&a) &= ~(1UL << 63); return a; }

double cos_table_0_01(double x)
{
    x = absd(x);
    x = modd(x, CONST_2PI);
    return costable_0_01[(int)(x * 100 + 0.5)];
}

Таблицы обеспечивают хороший баланс точности. Наихудший случай точности для самой маленькой таблицы составляет 0.49, то есть пользоваться ею нельзя. Но с каждым увеличением размера таблицы мы получаем на один знак после запятой больше точности: 0.049, 0.0049, 0.00049 и 0.000049. Тест времени исполнения для каждой таблицы равен примерно 0.38. А это быстро! (На самом деле, нам удалось снизить время исполнения примерно до 0.33, но код стал некрасивым.)

▍ Пятая попытка: таблица поиска с LERP

Да, таблица поиска — это отлично, но мы можем улучшить показатели для значений, которых в ней нет. Представляю вашему вниманию линейную интерполяцию. Звучит круто, но на самом деле это просто вычисление взвешенного среднего между двумя значениями. Если входного значения нет в таблице, мы вычисляем аппроксимацию на основании того, какой элемент таблицы ближе. Код:

#define lerp(w, v1, v2) ((1.0 - (w)) * (v1) + (w) * (v2))

double cos_table_0_01_LERP(double x)
{
    x = absd(x);
    x = modd(x, CONST_2PI);
    double i = x * 100.0;
    int index = (int)i;
    return lerp(i - index,        /* вес      */
        costable_0_01[index],     /* нижнее значение */
        costable_0_01[index + 1]  /* верхнее значение */
        );
}

▍ Окончательное сравнение

Вот таблица сравнения наших функций при наихудшем случае точности по сравнению с math.h (чем меньше, тем лучше!):


А вот сравнение времени исполнения наших функций при вычислении 100 000 000 значений (чем меньше, тем лучше!):

▍ Вывод

Так что же мы порекомендуем использовать? По возможности math.h. Ни одна из этих функций не является особо медленной, а большинство из них достаточно точно. Но когда в следующий раз я буду делать игру, активно использующую тригонометрические функции, то воспользуюсь таблицей 0_001.

До встречи в новой кроличьей норе!