Лучшая оптимизированная машиной полиномиальная минимаксная аппроксимация арктангенса на [-1,1]?

Ответ 1

Следующая функция представляет собой полномасштабную реализацию arctan на [0, 1]:

float atan_poly (float a) {
  float s = a * a, u = fmaf(a, -a, 0x1.fde90cp-1f);
  float r1 =               0x1.74dfb6p-9f;
  float r2 = fmaf (r1, u,  0x1.3a1c7cp-8f);
  float r3 = fmaf (r2, s, -0x1.7f24b6p-7f);
  float r4 = fmaf (r3, u, -0x1.eb3900p-7f);
  float r5 = fmaf (r4, s,  0x1.1ab95ap-5f);
  float r6 = fmaf (r5, u,  0x1.80e87cp-5f);
  float r7 = fmaf (r6, s, -0x1.e71aa4p-4f);
  float r8 = fmaf (r7, u, -0x1.b81b44p-3f);
  float r9 = r8 * s;
  float r10 = fmaf (r9, a, a);
  return r10;
}

Следующий тестовый жгут будет прерван, если функция atan_poly не будет точно округлена на [1e-16, 1] и в противном случае напечатает "успех":

int checkit(float f) {
  double d = atan(f);
  float d1 = d, d2 = d;
  if (d1 < d) d2 = nextafterf(d1, 1.0/0.0);
  else d1 = nextafterf(d1, -1.0/0.0);
  float p = atan_poly(f);
  if (p != d1 && p != d2) return 0;
  return 1;
}

int main() {
  for (float f = 1; f > 1e-16; f = nextafterf(f, -1.0/0.0)) {
    if (!checkit(f)) abort();
  }
  printf("success\n");
  exit(0);
}

Проблема с использованием s в каждом умножении состоит в том, что полиномиальные коэффициенты не быстро распадаются. Входы, близкие к 1, приводят к большому количеству отмены почти равных чисел, что означает, что вы пытаетесь найти набор коэффициентов, чтобы накопленное округление в конце вычисления близко аппроксимировало остаточное значение arctan.

Постоянная 0x1.fde90cp-1f - это число, близкое к 1, для которого (arctan(sqrt(x)) - x) / x^3 очень близко к ближайшему плаву. То есть, это константа, которая входит в вычисление u, так что кубический коэффициент почти полностью определяется. (Для этой программы кубический коэффициент должен быть либо -0x1.b81b44p-3f, либо -0x1.b81b42p-3f.)

Альтернативные умножения на s и u влияют на уменьшение эффекта ошибки округления в ri на r{i+2} не более чем на 1/4, так как s*u < 1/4 любое a, Это дает значительную свободу выбора при выборе коэффициентов пятого порядка и далее.

Я нашел коэффициенты с помощью двух программ:

• Одна программа включает в себя множество тестовых точек, записывает систему линейных неравенств и вычисляет оценки коэффициентов из этой системы неравенств. Обратите внимание, что, учитывая a, можно вычислить диапазон r8, который приведет к точно округленному результату. Чтобы получить линейные неравенства, я притворился, что r8 будет вычисляться как многочлен в float s и u в арифметике с вещественным числом; линейные неравенства ограничивали это вещественное число r8 лежать в некотором интервале. Я использовал библиотеку Parma Polyhedra для обработки этих систем ограничений.
• Другая программа случайным образом проверяет наборы коэффициентов в определенных диапазонах, вставляя сначала набор тестовых точек, а затем все float от 1 до 1e-8 в порядке убывания и проверяя, что atan_poly производит точное округление atan((double)x). Если какой-то x не удалось, он распечатал это x и почему он потерпел неудачу.

Чтобы получить коэффициенты, я взломал эту первую программу, чтобы исправить c3, выработать границы на r7 для каждой тестовой точки, а затем получить оценки на коэффициенты более высокого порядка. Затем я взломал его, чтобы исправить c3 и c5 и получить оценки на коэффициенты более высокого порядка. Я сделал это, пока у меня не было всего трех коэффициентов высшего порядка, c13, c15 и c17.

Я увеличил набор тестовых точек во второй программе, пока не прекратил ничего печатать или не распечатывал "успех". Мне потребовалось удивительно мало тестовых точек, чтобы отклонить почти все неправильные полиномы --- я считаю 85 контрольных точек в программе.

Здесь я показываю некоторые из моих работ, выбирая коэффициенты. Чтобы получить достоверно округленный arctan для моего начального набора контрольных точек, предполагая, что r1 через r8 оцениваются в реальной арифметике (и округлены как-то неприятно, но в некотором роде я не могу вспомнить), но r9 и r10 оцениваются в float арифметике, мне нужно:

-0x1.b81b456625f15p-3 <= c3 <= -0x1.b81b416e22329p-3
-0x1.e71d48d9c2ca4p-4 <= c5 <= -0x1.e71783472f5d1p-4
0x1.80e063cb210f9p-5 <= c7 <= 0x1.80ed6efa0a369p-5
0x1.1a3925ea0c5a9p-5 <= c9 <= 0x1.1b3783f148ed8p-5
-0x1.ec6032f293143p-7 <= c11 <= -0x1.e928025d508p-7
-0x1.8c06e851e2255p-7 <= c13 <= -0x1.732b2d4677028p-7
0x1.2aff33d629371p-8 <= c15 <= 0x1.41e9bc01ae472p-8
0x1.1e22f3192fd1dp-9 <= c17 <= 0x1.d851520a087c2p-9

Взятие c3 = -0x1.b81b44p-3, предполагая, что r8 также оценивается в float арифметике:

-0x1.e71df05b5ad56p-4 <= c5 <= -0x1.e7175823ce2a4p-4
0x1.80df529dd8b18p-5 <= c7 <= 0x1.80f00e8da7f58p-5
0x1.1a283503e1a97p-5 <= c9 <= 0x1.1b5ca5beeeefep-5
-0x1.ed2c7cd87f889p-7 <= c11 <= -0x1.e8c17789776cdp-7
-0x1.90759e6defc62p-7 <= c13 <= -0x1.7045e66924732p-7
0x1.27eb51edf324p-8 <= c15 <= 0x1.47cda0bb1f365p-8
0x1.f6c6b51c50b54p-10 <= c17 <= 0x1.003a00ace9a79p-8

Принимая c5 = -0x1.e71aa4p-4, предполагая, что r7 выполняется в float арифметике:

0x1.80e3dcc972cb3p-5 <= c7 <= 0x1.80ed1cf56977fp-5
0x1.1aa005ff6a6f4p-5 <= c9 <= 0x1.1afce9904742p-5
-0x1.ec7cf2464a893p-7 <= c11 <= -0x1.e9d6f7039db61p-7
-0x1.8a2304daefa26p-7 <= c13 <= -0x1.7a2456ddec8b2p-7
0x1.2e7b48f595544p-8 <= c15 <= 0x1.44437896b7049p-8
0x1.396f76c06de2ep-9 <= c17 <= 0x1.e3bedf4ed606dp-9

Принимая c7 = 0x1.80e87cp-5, предполагая, что r6 выполняется в float арифметике:

0x1.1aa86d25bb64fp-5 <= c9 <= 0x1.1aca48cd5caabp-5
-0x1.eb6311f6c29dcp-7 <= c11 <= -0x1.eaedb032dfc0cp-7
-0x1.81438f115cbbp-7 <= c13 <= -0x1.7c9a106629f06p-7
0x1.36d433f81a012p-8 <= c15 <= 0x1.3babb57bb55bap-8
0x1.5cb14e1d4247dp-9 <= c17 <= 0x1.84f1151303aedp-9

Взятие c9 = 0x1.1ab95ap-5, предполагая, что r5 выполняется в float арифметике:

-0x1.eb51a3b03781dp-7 <= c11 <= -0x1.eb21431536e0dp-7
-0x1.7fcd84700f7cfp-7 <= c13 <= -0x1.7ee38ee4beb65p-7
0x1.390fa00abaaabp-8 <= c15 <= 0x1.3b100a7f5d3cep-8
0x1.6ff147e1fdeb4p-9 <= c17 <= 0x1.7ebfed3ab5f9bp-9

Я выбрал точку, близкую к середине диапазона для c11, и произвольно выбрал c13, c15 и c17.

EDIT: теперь я автоматизировал эту процедуру. Следующая функция также является добросовестной реализацией arctan на [0, 1]:

float c5 = 0x1.997a72p-3;
float c7 = -0x1.23176cp-3;
float c9 = 0x1.b523c8p-4;
float c11 = -0x1.358ff8p-4;
float c13 = 0x1.61c5c2p-5;
float c15 = -0x1.0b16e2p-6;
float c17 = 0x1.7b422p-9;

float juffa_poly (float a) {
  float s = a * a;
  float r1 =              c17;
  float r2 = fmaf (r1, s, c15);
  float r3 = fmaf (r2, s, c13);
  float r4 = fmaf (r3, s, c11);
  float r5 = fmaf (r4, s, c9);
  float r6 = fmaf (r5, s, c7);
  float r7 = fmaf (r6, s, c5);
  float r8 = fmaf (r7, s, -0x1.5554dap-2f);
  float r9 = r8 * s;
  float r10 = fmaf (r9, a, a);
  return r10;
}

Мне кажется удивительным, что этот код существует. Для коэффициентов вблизи них вы можете получить привязку на расстоянии между r10 и значением многочлена, оцененного в реальной арифметике порядка нескольких ulps, благодаря медленной сходимости этого многочлена, когда s находится вблизи 1. Я ожидал, что ошибка округления будет вести себя так, как это было бы принципиально "необратимо", просто с помощью коэффициентов настройки.

Ответ 2

Это не ответ на вопрос, но слишком длинный, чтобы вписаться в комментарий:

ваш вопрос об оптимальном выборе коэффициентов C3, C5,..., C17 в полиномиальном приближении к арктангенсу, где вы привязали C1 к 1 и C2, C4,..., C16 к 0.

В названии вашего вопроса говорится, что вы ищете приближения на [-1, 1], и хорошей причиной для приведения четных коэффициентов к 0 является то, что достаточно и необходимо, чтобы приближение было точно нечетной функцией. Код в вашем вопросе "противоречит" названию, применяя полиномиальное приближение только на [0, 1].

Если вы используете алгоритм Ремеза для поиска коэффициентов C2, C3,..., C8 для полиномиальной аппроксимации арктанганта на [0, 1], вы можете получить что-то вроде значений ниже:

#include <stdio.h>
#include <math.h>

float atan_poly (float a)
{
  float r, s;
  s = a;
  //  s = a * a;

  r =             -3.3507930064626076153585890630056286726807491543578e-2;
  r = fmaf (r, s, 1.3859776280052980081098065189344699108643282883702e-1);
  r = fmaf (r, s, -1.8186361916440430105127602496688553126414578766147e-1);
  r = fmaf (r, s, -1.4583047494913656326643327729704639191810926020847e-2);
  r = fmaf (r, s, 2.1335202878219865228365738728594741358740655881373e-1);
  r = fmaf (r, s, -3.6801711826027841250774413728610805847547996647342e-3);
  r = fmaf (r, s, -3.3289852243978319173749528028057608377028846413080e-1);
  r = fmaf (r, s, -1.8631479933914856903459844359251948006605218562283e-5);
  r = fmaf (r, s, 1.2917291732886065585264586294461539492689296797761e-7);

  r = fmaf (r, a, a);
  return r;
}

int main() {
  for (float x = 0.0f; x < 1.0f; x+=0.1f)
    printf("x: %f\n%a\n%a\n\n", x, atan_poly(x), atan(x));
}

Это примерно такая же сложность, как и код в вашем вопросе - количество умножений аналогично. Глядя на этот многочлен, нет никакой причины, в частности, хотеть связать какой-либо коэффициент с 0. Если бы мы хотели аппроксимировать нечетную функцию по сравнению с [-1, 1] без закрепления четных коэффициентов, они автоматически появлялись бы очень маленькими и предметами к поглощению, а затем мы хотели бы привязать их к 0, но для этого приближения по сравнению с [0, 1] они не выполняются, поэтому нам не нужно связывать их.

Это могло быть лучше или хуже, чем нечетный многочлен в вашем вопросе. Оказывается, это хуже (см. Ниже). Это быстрое и грязное приложение LolRemez 0.2 (код, включенный в нижнюю часть вопроса), кажется, однако, достаточно хорош, чтобы поднять вопрос о выборе коэффициентов. В частности, мне было бы любопытно, что произойдет, если вы подвергли коэффициенты в этом ответе тому же шагу оптимизации "самый крутой приличный и симулированный отжиг", который вы применили для получения коэффициентов в своем вопросе.

Итак, чтобы обобщить это замечание, опубликованное как ответ, вы уверены, что ищете оптимальные коэффициенты C3, C5,..., C17? Мне кажется, что вы ищете лучшую последовательность операций с плавающей запятой с одинарной точностью, которые дают точное приближение к арктангенсу и что это приближение не должно быть формой Хорнера градуированного нечетного многочлена.

x: 0.000000
0x0p+0
0x0p+0

x: 0.100000
0x1.983e2cp-4
0x1.983e28938f9ecp-4

x: 0.200000
0x1.94442p-3
0x1.94441ff1e8882p-3

x: 0.300000
0x1.2a73a6p-2
0x1.2a73a71dcec16p-2

x: 0.400000
0x1.85a37ap-2
0x1.85a3770ebe7aep-2

x: 0.500000
0x1.dac67p-2
0x1.dac670561bb5p-2

x: 0.600000
0x1.14b1dcp-1
0x1.14b1ddf627649p-1

x: 0.700000
0x1.38b116p-1
0x1.38b113eaa384ep-1

x: 0.800000
0x1.5977a8p-1
0x1.5977a686e0ffbp-1

x: 0.900000
0x1.773388p-1
0x1.77338c44f8faep-1

Это код, который я связал с LolRemez 0.2, чтобы оптимизировать относительную точность полиномиальной аппроксимации степени-9 арктангенса на [0, 1]:

#include "lol/math/real.h"
#include "lol/math/remez.h"

using lol::real;
using lol::RemezSolver;

real f(real const &y)
{
  return (atan(y) - y) / y;
}

real g(real const &y)
{
  return re (atan(y) / y);
}

int main(int argc, char **argv)
{
  RemezSolver<8, real> solver;
  solver.Run("1e-1000", 1.0, f, g, 50);
  return 0;
}

Ответ 3

Я размышлял над различными идеями, которые я получил в комментариях, а также провел несколько экспериментов, основанных на этой обратной связи. В конце концов я решил, что изысканный эвристический поиск - лучший способ продвижения вперед. Теперь мне удалось уменьшить максимальную ошибку для atanf_poly() до 1.01036 ulps, и только три аргумента превысили мою заявленную цель с ошибкой 1 ulp:

ulp = -1.00829 @ |a| =  9.80738342e-001 0x1.f62356p-1 (3f7b11ab)
ulp = -1.01036 @ |a| =  9.87551928e-001 0x1.f9a068p-1 (3f7cd034)
ulp =  1.00050 @ |a| =  9.99375939e-001 0x1.ffae34p-1 (3f7fd71a)

Основываясь на способе генерации улучшенного приближения, нет никакой гарантии, что это наилучшее приближение; никакого научного прорыва здесь нет. Поскольку ошибка ulp текущего решения еще не полностью сбалансирована, и, поскольку продолжение поиска продолжает обеспечивать лучшие приближения (хотя и с экспоненциально возрастающими временными интервалами), я предполагаю, что вероятность ошибки 1 ulp достижима, но в то же время мы похоже, очень близки к лучшему оптимизированному машиной приближению.

Лучшее качество нового приближения является результатом изысканного процесса поиска. Я заметил, что все наибольшие ошибки ulp в полиноме близки к единице, например, в [0.75,1.0], являются консервативными. Это позволяет быстро сканировать интересные наборы коэффициентов, максимальная погрешность которых меньше некоторой границы, скажем, 1,08 ulps. Затем я могу подробно и подробно проанализировать все наборы коэффициентов в эвристически выбранном гиперконусе, закрепленном в этой точке. Этот второй шаг ищет минимальную ошибку ulp в качестве основной цели и максимальный процент правильно округленных результатов в качестве вторичной цели. Используя этот двухэтапный процесс на всех четырех ядрах моего процессора, я смог значительно ускорить процесс поиска: до сих пор я смог проверить коэффициенты коэффициентов 2 ²¹.

В зависимости от диапазона каждого коэффициента по всем "близким" решениям я теперь оцениваю, что общее полезное пространство поиска для этой задачи аппроксимации составляет >= 2 ²⁴ множители коэффициентов, а не более оптимистичное число 2 ²⁰ Я выбросил раньше. Это похоже на возможную проблему для решения для кого-то, кто либо очень терпелив, либо имеет множество вычислительных лошадей в их распоряжении.

Мой обновленный код выглядит следующим образом:

// max ulp err = 1.01036
float atanf_poly (float a)
{
    float r, s;
    s = a * a;
    r =              0x1.7ed22cp-9f;
    r = fmaf (r, s, -0x1.0c2c2ep-6f);
    r = fmaf (r, s,  0x1.61fdf6p-5f);
    r = fmaf (r, s, -0x1.3556b4p-4f);
    r = fmaf (r, s,  0x1.b4e12ep-4f);
    r = fmaf (r, s, -0x1.230ae0p-3f);
    r = fmaf (r, s,  0x1.9978eep-3f);
    r = fmaf (r, s, -0x1.5554dap-2f);
    r = r * s;
    r = fmaf (r, a, a);
    return r;
}

// max ulp err = 1.51871
float my_atanf (float a)
{
    float r, t;
    t = fabsf (a);
    r = t;
    if (t > 1.0f) {
        r = 1.0f / r;
    }
    r = atanf_poly (r);
    if (t > 1.0f) {
        r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, -r); // pi/2 - r
    }
    r = copysignf (r, a);
    return r;
}

Обновить (после повторения вопроса через два с половиной года)

Используя T. Myklebust проект публикации в качестве отправной точки, я нашел арктангенциальное приближение на [-1,1], которое имеет наименьший ошибки, чтобы иметь максимальную ошибку 0,94528 ulp.

/* Based on: Tor Myklebust, "Computing accurate Horner form approximations 
   to special functions in finite precision arithmetic", arXiv:1508.03211,
   August 2015. maximum ulp err = 0.94528
*/
float atanf_poly (float a)
{
    float r, s;
    s = a * a;                        
    r =              0x1.6d2086p-9f;  //  2.78569828e-3
    r = fmaf (r, s, -0x1.03f2ecp-6f); // -1.58660226e-2
    r = fmaf (r, s,  0x1.5beebap-5f); //  4.24722321e-2
    r = fmaf (r, s, -0x1.33194ep-4f); // -7.49753043e-2
    r = fmaf (r, s,  0x1.b403a8p-4f); //  1.06448799e-1
    r = fmaf (r, s, -0x1.22f5c2p-3f); // -1.42070308e-1
    r = fmaf (r, s,  0x1.997748p-3f); //  1.99934542e-1
    r = fmaf (r, s, -0x1.5554d8p-2f); // -3.33331466e-1
    r = r * s;
    r = fmaf (r, a, a);
    return r;
}

Ответ 4

Это не ответ, но расширенный комментарий.

Последние процессоры Intel и некоторые будущие процессоры AMD имеют AVX2. В Linux найдите флаг avx2 в /proc/cpuinfo, чтобы узнать, поддерживает ли ваш процессор эти файлы.

AVX2 - это расширение, которое позволяет нам строить и вычислять с использованием 256-битных векторов - например, восемь чисел с одной точностью или четыре числа с двойной точностью - вместо просто скаляров. Он включает поддержку FMA3, что означает плавное многократное добавление для таких векторов. Проще говоря, AVX2 позволяет нам оценивать восемь полиномов параллельно, почти в то же время, когда мы оцениваем один, используя скалярные операции.

Функция error8() анализирует один набор коэффициентов с использованием предопределенных значений x, сравнивая их с предварительно вычисленными значениями atan(x) и возвращает ошибку в ULP (ниже и выше желаемого результата отдельно), а также количество результатов, которые точно соответствуют требуемому значению с плавающей запятой. Они не нужны для простого тестирования того, лучше ли набор коэффициентов, чем самый известный в настоящий момент набор, но допускают различные стратегии, по которым будут проверяться коэффициенты. (В принципе, максимальная ошибка в ULP формирует поверхность, и мы пытаемся найти самую низкую точку на этой поверхности, зная, что "высота" поверхности в каждой точке позволяет нам сделать обоснованные догадки о том, в каком направлении идти - - как изменить коэффициенты.)

Используются четыре таблицы с заранее рассчитанными таблицами: known_x для аргументов known_f для корректно округленных результатов с одной точностью, known_a для точного значения с двойной точностью (я просто надеюсь, что библиотека atan() достаточно точен для этого, но на него не следует полагаться, не проверяя!) и known_m для масштабирования разности двойной точности для ULP. Учитывая желаемый диапазон аргументов, функция precalculate() будет предварительно вычислять их с помощью функции библиотеки atan(). (Он также полагается на форматы с плавающей точкой IEEE-754, а порядок байтов с плавающей точкой и целочисленным байтом - один и тот же, но это верно для процессоров, на которых выполняется этот код.)

Обратите внимание, что массивы known_x, known_f и known_a могут храниться в двоичных файлах; содержание known_m тривиально получено из known_a. Использование библиотеки atan() без проверки это не очень хорошая идея, но поскольку мои результаты соответствуют результатам njuffa, я не стал искать лучшую ссылку atan().

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

#define _POSIX_C_SOURCE 200809L
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <immintrin.h>
#include <math.h>
#include <errno.h>

/** poly8() - Compute eight polynomials in parallel.
 * @x - the arguments
 * @c - the coefficients.
 *
 * The first coefficients are for degree 17, the second
 * for degree 15, and so on, down to degree 3.
 *
 * The compiler should vectorize the expression using vfmaddXXXps
 * given an AVX2-capable CPU; for example, Intel Haswell,
 * Broadwell, Haswell E, Broadwell E, Skylake, or Cannonlake;
 * or AMD Excavator CPUs. Tested on Intel Core i5-4200U.
 *
 * Using GCC-4.8.2 and
 *     gcc -O2 -march=core-avx2 -mtune=generic
 * this code produces assembly (AT&T syntax)
 *     vmulps       %ymm0, %ymm0, %ymm2
 *     vmovaps      (%rdi), %ymm1
 *     vmovaps      %ymm0, %ymm3
 *     vfmadd213ps  32(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  64(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  96(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  128(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  160(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  192(%rdi), %ymm2, %ymm1
 *     vfmadd213ps  224(%rdi), %ymm2, %ymm1
 *     vmulps       %ymm2, %ymm1, %ymm0
 *     vfmadd132ps  %ymm3, %ymm3, %ymm0
 *     ret
 * if you omit the 'static inline'.
*/
static inline __v8sf poly8(const __v8sf x, const __v8sf *const c)
{
    const __v8sf xx = x * x;
    return (((((((c[0]*xx + c[1])*xx + c[2])*xx + c[3])*xx + c[4])*xx + c[5])*xx + c[6])*xx + c[7])*xx*x + x;
}

/** error8() - Calculate maximum error in ULPs
 * @x  - the arguments
 * @co - { C17, C15, C13, C11, C9, C7, C5, C3 }
 * @f  - the correctly rounded results in single precision
 * @a  - the expected results in double precision
 * @m  - 16777216.0 raised to the same power of two as @a normalized
 * @n  - number of vectors to test
 * @max_under - pointer to store the maximum underflow (negative, in ULPs) to
 * @max_over  - pointer to store the maximum overflow (positive, in ULPs) to
 * Returns the number of correctly rounded float results, 0..8*n.
*/
size_t error8(const __v8sf *const x, const float *const co,
              const __v8sf *const f, const __v4df *const a, const __v4df *const m,
              const size_t n,
              float *const max_under, float *const max_over)
{
    const __v8sf c[8] = { { co[0], co[0], co[0], co[0], co[0], co[0], co[0], co[0] },
                          { co[1], co[1], co[1], co[1], co[1], co[1], co[1], co[1] },
                          { co[2], co[2], co[2], co[2], co[2], co[2], co[2], co[2] },
                          { co[3], co[3], co[3], co[3], co[3], co[3], co[3], co[3] },
                          { co[4], co[4], co[4], co[4], co[4], co[4], co[4], co[4] },
                          { co[5], co[5], co[5], co[5], co[5], co[5], co[5], co[5] },
                          { co[6], co[6], co[6], co[6], co[6], co[6], co[6], co[6] },
                          { co[7], co[7], co[7], co[7], co[7], co[7], co[7], co[7] } };
    __v4df min = { 0.0, 0.0, 0.0, 0.0 };
    __v4df max = { 0.0, 0.0, 0.0, 0.0 };
    __v8si eqs = { 0, 0, 0, 0, 0, 0, 0, 0 };
    size_t i;

    for (i = 0; i < n; i++) {
        const __v8sf v = poly8(x[i], c);
        const __v4df d0 = { v[0], v[1], v[2], v[3] };
        const __v4df d1 = { v[4], v[5], v[6], v[7] };
        const __v4df err0 = (d0 - a[2*i+0]) * m[2*i+0];
        const __v4df err1 = (d1 - a[2*i+1]) * m[2*i+1];
        eqs -= (__v8si)_mm256_cmp_ps(v, f[i], _CMP_EQ_OQ);
        min = _mm256_min_pd(min, err0);
        max = _mm256_max_pd(max, err1);
        min = _mm256_min_pd(min, err1);
        max = _mm256_max_pd(max, err0);
    }

    if (max_under) {
        if (min[0] > min[1]) min[0] = min[1];
        if (min[0] > min[2]) min[0] = min[2];
        if (min[0] > min[3]) min[0] = min[3];
        *max_under = min[0];
    }

    if (max_over) {
        if (max[0] < max[1]) max[0] = max[1];
        if (max[0] < max[2]) max[0] = max[2];
        if (max[0] < max[3]) max[0] = max[3];
        *max_over = max[0];
    }

    return (size_t)((unsigned int)eqs[0])
         + (size_t)((unsigned int)eqs[1])
         + (size_t)((unsigned int)eqs[2])
         + (size_t)((unsigned int)eqs[3])
         + (size_t)((unsigned int)eqs[4])
         + (size_t)((unsigned int)eqs[5])
         + (size_t)((unsigned int)eqs[6])
         + (size_t)((unsigned int)eqs[7]);
}

/** precalculate() - Allocate and precalculate tables for error8().
 * @x0   - First argument to precalculate
 * @x1   - Last argument to precalculate
 * @xptr - Pointer to a __v8sf pointer for the arguments
 * @fptr - Pointer to a __v8sf pointer for the correctly rounded results
 * @aptr - Pointer to a __v4df pointer for the comparison results
 * @mptr - Pointer to a __v4df pointer for the difference multipliers
 * Returns the vector count if successful,
 * 0 with errno set otherwise.
*/
size_t precalculate(const float x0, const float x1,
                 __v8sf **const xptr, __v8sf **const fptr,
                 __v4df **const aptr, __v4df **const mptr)
{
    const size_t align = 64;
    unsigned int i0, i1;
    size_t       n, i, sbytes, dbytes;
    __v8sf      *x = NULL;
    __v8sf      *f = NULL;
    __v4df      *a = NULL;
    __v4df      *m = NULL;

    if (!xptr || !fptr || !aptr || !mptr) {
        errno = EINVAL;
        return (size_t)0;
    }

    memcpy(&i0, &x0, sizeof i0);
    memcpy(&i1, &x1, sizeof i1);

    i0 ^= (i0 & 0x80000000U) ? 0xFFFFFFFFU : 0x80000000U;
    i1 ^= (i1 & 0x80000000U) ? 0xFFFFFFFFU : 0x80000000U;

    if (i1 > i0)
        n = (((size_t)i1 - (size_t)i0) | (size_t)7) + (size_t)1;
    else
    if (i0 > i1)
        n = (((size_t)i0 - (size_t)i1) | (size_t)7) + (size_t)1;
    else {
        errno = EINVAL;
        return (size_t)0;
    }

    sbytes = n * sizeof (float);
    if (sbytes % align)
        sbytes += align - (sbytes % align);

    dbytes = n * sizeof (double);
    if (dbytes % align)
        dbytes += align - (dbytes % align);

    if (posix_memalign((void **)&x, align, sbytes)) {
        errno = ENOMEM;
        return (size_t)0;
    }
    if (posix_memalign((void **)&f, align, sbytes)) {
        free(x);
        errno = ENOMEM;
        return (size_t)0;
    }
    if (posix_memalign((void **)&a, align, dbytes)) {
        free(f);
        free(x);
        errno = ENOMEM;
        return (size_t)0;
    }
    if (posix_memalign((void **)&m, align, dbytes)) {
        free(a);
        free(f);
        free(x);
        errno = ENOMEM;
        return (size_t)0;
    }

    if (x1 > x0) {
        float *const xp = (float *)x;
        float        curr = x0;

        for (i = 0; i < n; i++) {
            xp[i] = curr;
            curr = nextafterf(curr, HUGE_VALF);
        }

        i = n;
        while (i-->0 && xp[i] > x1)
            xp[i] = x1;
    } else {
        float *const xp = (float *)x;
        float        curr = x0;

        for (i = 0; i < n; i++) {
            xp[i] = curr;
            curr = nextafterf(curr, -HUGE_VALF);
        }

        i = n;
        while (i-->0 && xp[i] < x1)
            xp[i] = x1;
    }

    {
        const float *const xp = (const float *)x;
        float *const       fp = (float *)f;
        double *const      ap = (double *)a;
        double *const      mp = (double *)m;

        for (i = 0; i < n; i++) {
            const float curr = xp[i];
            int         temp;

            fp[i] = atanf(curr);
            ap[i] = atan((double)curr);

            (void)frexp(ap[i], &temp);
            mp[i] = ldexp(16777216.0, temp);
        }
    }

    *xptr = x;
    *fptr = f;
    *aptr = a;
    *mptr = m;

    errno = 0;
    return n/8;
}

static int parse_range(const char *const str, float *const range)
{
    float fmin, fmax;
    char  dummy;

    if (sscanf(str, " %f %f %c",   &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %f:%f %c",   &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %f,%f %c",   &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %f/%f %c",   &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %ff %ff %c", &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %ff:%ff %c", &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %ff,%ff %c", &fmin, &fmax, &dummy) == 2 ||
        sscanf(str, " %ff/%ff %c", &fmin, &fmax, &dummy) == 2) {
        if (range) {
            range[0] = fmin;
            range[1] = fmax;
        }
        return 0;
    }

    if (sscanf(str, " %f %c",  &fmin, &dummy) == 1 ||
        sscanf(str, " %ff %c", &fmin, &dummy) == 1) {
        if (range) {
            range[0] = fmin;
            range[1] = fmin;
        }
        return 0;
    }

    return errno = ENOENT;
}

static int fix_range(float *const f)
{
    if (f && f[0] > f[1]) {
        const float tmp = f[0];
        f[0] = f[1];
        f[1] = tmp;
    }
    return f && isfinite(f[0]) && isfinite(f[1]) && (f[1] >= f[0]);
}

static const char *f2s(char *const buffer, const size_t size, const float value, const char *const invalid)
{
    char  format[32];
    float parsed;
    int   decimals, length;

    for (decimals = 0; decimals <= 16; decimals++) {
        length = snprintf(format, sizeof format, "%%.%df", decimals);
        if (length < 1 || length >= (int)sizeof format)
            break;

        length = snprintf(buffer, size, format, value);
        if (length < 1 || length >= (int)size)
            break;

        if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value)
            return buffer;

        decimals++;
    }

    for (decimals = 0; decimals <= 16; decimals++) {
        length = snprintf(format, sizeof format, "%%.%dg", decimals);
        if (length < 1 || length >= (int)sizeof format)
            break;

        length = snprintf(buffer, size, format, value);
        if (length < 1 || length >= (int)size)
            break;

        if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value)
            return buffer;

        decimals++;
    }

    length = snprintf(buffer, size, "%a", value);
    if (length < 1 || length >= (int)size)
        return invalid;

    if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value)
        return buffer;

    return invalid;
}

int main(int argc, char *argv[])
{
    float xrange[2] = { 0.75f, 1.00f };
    float c17range[2], c15range[2], c13range[2], c11range[2];
    float c9range[2], c7range[2], c5range[2], c3range[2];
    float c[8];

    __v8sf *known_x;
    __v8sf *known_f;
    __v4df *known_a;
    __v4df *known_m;
    size_t  known_n;

    if (argc != 10 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) {
        fprintf(stderr, "\n");
        fprintf(stderr, "Usage: %s [ -h | --help ]\n", argv[0]);
        fprintf(stderr, "       %s C17 C15 C13 C11 C9 C7 C5 C3 x\n", argv[0]);
        fprintf(stderr, "\n");
        fprintf(stderr, "Each of the coefficients can be a constant or a range,\n");
        fprintf(stderr, "for example 0.25 or 0.75:1. x must be a non-empty range.\n");
        fprintf(stderr, "\n");
        return EXIT_FAILURE;
    }

    if (parse_range(argv[1], c17range) || !fix_range(c17range)) {
        fprintf(stderr, "%s: Invalid C17 range or constant.\n", argv[1]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[2], c15range) || !fix_range(c15range)) {
        fprintf(stderr, "%s: Invalid C15 range or constant.\n", argv[2]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[3], c13range) || !fix_range(c13range)) {
        fprintf(stderr, "%s: Invalid C13 range or constant.\n", argv[3]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[4], c11range) || !fix_range(c11range)) {
        fprintf(stderr, "%s: Invalid C11 range or constant.\n", argv[4]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[5], c9range) || !fix_range(c9range)) {
        fprintf(stderr, "%s: Invalid C9 range or constant.\n", argv[5]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[6], c7range) || !fix_range(c7range)) {
        fprintf(stderr, "%s: Invalid C7 range or constant.\n", argv[6]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[7], c5range) || !fix_range(c5range)) {
        fprintf(stderr, "%s: Invalid C5 range or constant.\n", argv[7]);
        return EXIT_FAILURE;
    }
    if (parse_range(argv[8], c3range) || !fix_range(c3range)) {
        fprintf(stderr, "%s: Invalid C3 range or constant.\n", argv[8]);
        return EXIT_FAILURE;
    }

    if (parse_range(argv[9], xrange) || xrange[0] == xrange[1] ||
        !isfinite(xrange[0]) || !isfinite(xrange[1])) {
        fprintf(stderr, "%s: Invalid x range.\n", argv[9]);
        return EXIT_FAILURE;
    }

    known_n = precalculate(xrange[0], xrange[1], &known_x, &known_f, &known_a, &known_m);
    if (!known_n) {
        if (errno == ENOMEM)
            fprintf(stderr, "Not enough memory for precalculated tables.\n");
        else
            fprintf(stderr, "Invalid (empty) x range.\n");
        return EXIT_FAILURE;
    }

    fprintf(stderr, "Precalculated %lu arctangents to compare to.\n", 8UL * (unsigned long)known_n);
    fprintf(stderr, "\nC17 C15 C13 C11 C9 C7 C5 C3 max-ulps-under max-ulps-above correctly-rounded percentage cycles\n");
    fflush(stderr);

    {
        const double  percent = 12.5 / (double)known_n;
        size_t        rounded;
        char          c17buffer[64], c15buffer[64], c13buffer[64], c11buffer[64];
        char          c9buffer[64], c7buffer[64], c5buffer[64], c3buffer[64];
        char          minbuffer[64], maxbuffer[64];
        float         minulps, maxulps;
        unsigned long tsc_start, tsc_stop;

        for (c[0] = c17range[0]; c[0] <= c17range[1]; c[0] = nextafterf(c[0], HUGE_VALF))
        for (c[1] = c15range[0]; c[1] <= c15range[1]; c[1] = nextafterf(c[1], HUGE_VALF))
        for (c[2] = c13range[0]; c[2] <= c13range[1]; c[2] = nextafterf(c[2], HUGE_VALF))
        for (c[3] = c11range[0]; c[3] <= c11range[1]; c[3] = nextafterf(c[3], HUGE_VALF))
        for (c[4] = c9range[0]; c[4] <= c9range[1]; c[4] = nextafterf(c[4], HUGE_VALF))
        for (c[5] = c7range[0]; c[5] <= c7range[1]; c[5] = nextafterf(c[5], HUGE_VALF))
        for (c[6] = c5range[0]; c[6] <= c5range[1]; c[6] = nextafterf(c[6], HUGE_VALF))
        for (c[7] = c3range[0]; c[7] <= c3range[1]; c[7] = nextafterf(c[7], HUGE_VALF)) {
            tsc_start = __builtin_ia32_rdtsc();
            rounded = error8(known_x, c, known_f, known_a, known_m, known_n, &minulps, &maxulps);
            tsc_stop = __builtin_ia32_rdtsc();
            printf("%-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %lu %.3f %lu\n",
                   f2s(c17buffer, sizeof c17buffer, c[0], "?"),
                   f2s(c15buffer, sizeof c15buffer, c[1], "?"),
                   f2s(c13buffer, sizeof c13buffer, c[2], "?"),
                   f2s(c11buffer, sizeof c11buffer, c[3], "?"),
                   f2s(c9buffer, sizeof c9buffer, c[4], "?"),
                   f2s(c7buffer, sizeof c7buffer, c[5], "?"),
                   f2s(c5buffer, sizeof c5buffer, c[6], "?"),
                   f2s(c3buffer, sizeof c3buffer, c[7], "?"),
                   f2s(minbuffer, sizeof minbuffer, minulps, "?"),
                   f2s(maxbuffer, sizeof maxbuffer, maxulps, "?"),
                   rounded, (double)rounded * percent,
                   (unsigned long)(tsc_stop - tsc_start));
            fflush(stdout);

        }
    }

    return EXIT_SUCCESS;
}

Код компилируется с использованием GCC-4.8.2 в Linux, но может потребоваться изменить для других компиляторов и/или ОС. (Я был бы счастлив включить/принять исправления, исправляющие их, но у меня просто нет Windows или ICC, поэтому я мог проверить.)

Чтобы скомпилировать это, я рекомендую

gcc -Wall -O3 -fomit-frame-pointer -march=native -mtune=native example.c -lm -o example

Запустите без аргументов, чтобы увидеть использование; или

./example 0x1.7ed24ap-9f -0x1.0c2c12p-6f  0x1.61fdd2p-5f -0x1.3556b0p-4f 0x1.b4e138p-4f -0x1.230ae2p-3f  0x1.9978eep-3f -0x1.5554dap-2f 0.75:1

чтобы проверить, что он сообщает для набора коэффициентов njuffa, по сравнению со стандартной библиотекой C library atan(), со всеми возможными x из [0.75, 1].

Вместо фиксированного коэффициента вы также можете использовать min:max для определения диапазона сканирования (сканирование всех уникальных значений с плавающей запятой с одной точностью). Проверяется каждая возможная комбинация коэффициентов.

Поскольку я предпочитаю десятичную нотацию, но мне нужно сохранить значения точными, я использую функцию f2s() для отображения значений с плавающей запятой. Это простая вспомогательная функция грубой силы, которая использует кратчайшее форматирование, которое дает одно и то же значение при анализе обратно на float.

Например,

./example 0x1.7ed248p-9f:0x1.7ed24cp-9f -0x1.0c2c10p-6f:-0x1.0c2c14p-6f 0x1.61fdd0p-5f:0x1.61fdd4p-5f -0x1.3556aep-4f:-0x1.3556b2p-4f 0x1.b4e136p-4f:0x1.b4e13ap-4f -0x1.230ae0p-3f:-0x1.230ae4p-3f 0x1.9978ecp-3f:0x1.9978f0p-3f -0x1.5554d8p-2f:-0x1.5554dcp-2f 0.75:1

вычисляет все комбинации коэффициентов 6561 (3 ⁸) ± 1 ULP вокруг набора njuffa для x в [0.75, 1]. (Действительно, это показывает, что уменьшение C ₁₇ на 1 ULP до 0x1.7ed248p-9f дает точные результаты.)

(Этот запуск занял 90 секунд на Core i5-4200U на частоте 2,6 ГГц - почти в моей оценке 30 коэффициентов в секунду на ГГц на ядро. Хотя этот код не имеет резьбы, ключевыми функциями являются поточно- безопасный, поэтому нарезание резьбы не должно быть слишком сложным. Это Core i5-4200U - это ноутбук, и он становится очень жарким даже при напряжении только одного ядра, поэтому я не беспокоился.)

(Я считаю, что вышеуказанный код находится в общественном достоянии или CC0-лицензирован, если преданность общественному достоянию невозможна. На самом деле, я не уверен, что он достаточно креативен, чтобы быть полностью защищенным авторским правом. использовать его где угодно, как угодно, если вы не обвиняете меня, если он ломается.)

Вопросы? Улучшения? Разрешения для исправления Linux/GCC'ism приветствуются!