Быстрая реализация синуса, косинуса и квадратного корня в С++ (не обязательно должна быть очень точной)

Ответ 1

Самый быстрый способ - предварительно вычислить значения и использовать таблицу, как в этом примере:

Создать таблицу синусоидального поиска на С++

НО, если вы настаиваете на вычислении во время выполнения, вы можете использовать расширение серии Taylor синуса или косинуса...

Подробнее о серии Тейлора... http://en.wikipedia.org/wiki/Taylor_series

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

Также... Не размножайте x ^ n с самого начала каждый раз... например. умножьте x ^ 3 на x еще два раза, а затем на два других, чтобы вычислить показатели.

Ответ 2

Во-первых, серия Taylor НЕ является лучшим/самым быстрым способом реализации синус /cos. Это также не то, как профессиональные библиотеки реализуют эти тригонометрические функции, а знание наилучшей числовой реализации позволяет настроить точность, чтобы повысить скорость более эффективно. Кроме того, эта проблема уже широко обсуждалась в Кару. Вот только один пример.

Во-вторых, большая разница между старым/новым ПК связана с тем, что современная архитектура Intel имеет явный код сборки для вычисления элементарных тригонометрических функций. Сложно превзойти их по скорости исполнения.

Наконец, давайте поговорим о коде на вашем старом ПК. Проверьте реализацию gsl gnu научной библиотеки (или числовые рецепты), и вы увидите, что они в основном используют формулу аппроксимации Чебышева.

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

Для решения этой проблемы, если вам нужны подробности реализации какой-либо специальной функции или числового метода, вам следует взглянуть на код GSL, прежде чем предпринимать какие-либо дальнейшие действия - GSL является СТАНДАРТНОЙ числовой библиотекой.

ОБНОВЛЕНИЕ: вы можете улучшить время выполнения, включив агрессивные флаги оптимизации с плавающей запятой в gcc/icc. Это снизит точность, но кажется, что это именно то, что вы хотите.

ОБНОВЛЕНИЕ 2: Вы можете попытаться создать сетку грубых грехов и использовать процедуру gsl (gsl_interp_cspline_periodic для сплайна с периодическими условиями), чтобы сплайновать эту таблицу (сплайн уменьшит ошибки по сравнению с линейной интерполяцией => вам нужно меньше точек на вашей таблице => меньше кеша не хватает)!

Ответ 3

Здесь гарантированная самая быстрая синусоидальная функция в C++:

double FastSin(double x)
{
    return 0;
}

О, ты хотел лучшей точности, чем | 1.0 |? Хорошо читайте дальше.

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

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

Большинство людей, которые задают этот вопрос, не понимают важности компромисса между производительностью и точностью. В частности, вам придется сделать некоторые выборы в отношении точности вычислений, прежде чем делать что-либо еще. Сколько ошибок вы можете терпеть в результате? 10 ^ -4? 10 ^ -16?

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

Никто не использует только ряды Тейлора для аппроксимации трансцендентных в программном обеспечении. За исключением некоторых весьма специфических случаев, ряды Тейлора обычно медленно приближаются к цели через общие входные диапазоны.

Алгоритмы, которые ваши дедушка и бабушка использовали для эффективного расчета трансцендентальных значений, в совокупности называются CORDIC и были достаточно просты для реализации на аппаратном уровне. Вот хорошо документированная реализация CORDIC на C. Реализации CORDIC, как правило, требуют очень маленькой справочной таблицы, но большинство реализаций даже не требуют наличия аппаратного множителя. Большинство реализаций CORDIC позволяют вам компенсировать производительность за точность, включая ту, которую я связал.

За эти годы было много постепенных улучшений оригинальных алгоритмов CORDIC. Например, в прошлом году некоторые исследователи в Японии опубликовали статью об улучшенном CORDIC с лучшими углами поворота, что сокращает необходимые операции.

Если у вас есть аппаратные множители (и вы почти наверняка), или если вы не можете позволить себе справочную таблицу, как требует CORDIC, вы всегда можете использовать полином Чебышева, чтобы сделать то же самое. Полиномы Чебышева требуют умножения, но это редко является проблемой на современном оборудовании. Нам нравятся полиномы Чебышева, потому что они имеют максимально предсказуемые максимальные ошибки для данного приближения. Максимум последнего слагаемого в многочлене Чебышёва по всему входному диапазону ограничивает ошибку в результате. И эта ошибка уменьшается по мере увеличения числа членов. Вот один пример полинома Чебышева, дающий синус-аппроксимацию в огромном диапазоне, игнорирующий естественную симметрию функции синуса и просто решающий проблему аппроксимации, добавляя в нее больше коэффициентов.

Нам также нравятся полиномы Чебышева, потому что ошибка в приближении равномерно распределена по диапазону выходных данных. Если вы пишете аудио-плагины или выполняете цифровую обработку сигналов, полиномы Чебышева дают вам дешевый и предсказуемый эффект сглаживания "бесплатно".

Если вы хотите найти свои собственные полиномиальные коэффициенты Чебышева в определенном диапазоне, многие математические библиотеки называют процесс поиска этих коэффициентов "Чебышевская подгонка" или что-то в этом роде.

Квадратные корни, как и сейчас, обычно рассчитываются по некоторому варианту алгоритма Ньютона-Рафсона, обычно с фиксированным числом итераций. Обычно, когда кто-то запускает "удивительно новый" алгоритм для создания квадратных корней, это просто замаскированный Ньютон-Рафсон.

Полиномы Ньютона-Рафсона, CORDIC и Чебышева позволяют вам поменять скорость на точность, поэтому ответ может быть настолько неточным, насколько вы хотите.

Наконец, когда вы закончили все свои модные тесты и микрооптимизацию, убедитесь, что ваша "быстрая" версия на самом деле быстрее, чем версия библиотеки. Вот типичная реализация библиотеки fsin() , ограниченной в домене от -pi/4 до pi/4. И чертовски медленно.

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

ТЛ: д-р; перейдите в Google "синус-аппроксимация" или "косинус-аппроксимация" или "квадратный корень приближения" или "теория аппроксимации."

Ответ 4

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

Число с плавающей точкой, определенное IEEE-754, использует некоторый определенный бит, представляющий время описания нескольких основанных 2. Некоторые биты представляют собой базовое значение.

float squareRoot(float x)
{
  unsigned int i = *(unsigned int*) &x;

  // adjust bias
  i  += 127 << 23;
  // approximation of square root
  i >>= 1;

  return *(float*) &i;
}

Это постоянное время вычисления корня квадрата

Ответ 5

Основываясь на идее http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648 и некотором ручном переписывании для улучшения производительности в микро-бенчмарке, я в итоге получил следующую косинусную реализацию, которая используется в физическом моделировании HPC, который является узким местом из-за повторных вызовов cos на большом пространстве чисел. Он достаточно точный и намного быстрее, чем справочная таблица, в особенности деление не требуется.

template<typename T>
inline T cos(T x) noexcept
{
    constexpr T tp = 1./(2.*M_PI);
    x *= tp;
    x -= T(.25) + std::floor(x + T(.25));
    x *= T(16.) * (std::abs(x) - T(.5));
    #if EXTRA_PRECISION
    x += T(.225) * x * (std::abs(x) - T(1.));
    #endif
    return x;
}

Компилятор Intel по крайней мере также достаточно умен для векторизации этой функции при использовании в цикле.

Если определено EXTRA_PRECISION, максимальная ошибка составляет около 0,00109 для диапазона от -π до π, при условии, что T является double как это обычно определяется в большинстве реализаций C++. В противном случае максимальная ошибка составляет около 0,056 для того же диапазона.

Ответ 6

Для x86 аппаратные команды FP с квадратным корнем выполняются быстро (sqrtss is sqrt Scalar Single-precision). Одинарная точность быстрее, чем двойная точность, поэтому используйте float вместо double для кода, где вы можете позволить себе использовать меньшую точность.

Для 32-битного кода вам обычно нужны параметры компилятора, чтобы заставить его выполнять математику FP с инструкциями SSE, а не x87. (например, -mfpmath=sse)

Чтобы получить функции C sqrt() или sqrtf() для встраивания как только sqrtsd или sqrtss, , вам нужно скомпилировать с помощью -fno-math-errno. Наличие математических функций errno на NaN обычно считается ошибкой дизайна, но стандарт требует этого. Без этой опции gcc вставляет ее, а затем имеет ветку compare +, чтобы увидеть, был ли результат NaN, и если это так вызывает библиотечную функцию, поэтому он может установить errno. Если ваша программа не проверяет errno после математических функций, нет опасности при использовании -fno-math-errno.

Вам не нужны какие-либо "небезопасные" части -ffast-math, чтобы получить sqrt и некоторые другие функции для встроенного улучшения или вообще, но -ffast-math может иметь большое значение (например, разрешить компилятору авто-векторизация в тех случаях, когда это изменяет результат, потому что математика FP не является ассоциативной.

например. с компиляцией gcc6.3 float foo(float a){ return sqrtf(a); }

foo:    # with -O3 -fno-math-errno.
    sqrtss  xmm0, xmm0
    ret

foo:   # with just -O3
    pxor    xmm2, xmm2   # clang just checks for NaN, instead of comparing against zero.
    sqrtss  xmm1, xmm0
    ucomiss xmm2, xmm0
    ja      .L8          # take the slow path if 0.0 > a
    movaps  xmm0, xmm1
    ret

.L8:                     # errno-setting path
    sub     rsp, 24
    movss   DWORD PTR [rsp+12], xmm1   # store the sqrtss result because the x86-64 SysV ABI has no call-preserved xmm regs.
    call    sqrtf                      # call sqrtf just to set errno
    movss   xmm1, DWORD PTR [rsp+12]
    add     rsp, 24
    movaps  xmm0, xmm1    # extra mov because gcc reloaded into the wrong register.
    ret

Код gcc для случая NaN кажется слишком сложным; он даже не использует возвращаемое значение sqrtf! В любом случае, это тот беспорядок, который вы фактически получаете без -fno-math-errno, для каждого сайта вызова sqrtf() в вашей программе. В основном это просто код-раздувание, и ни один из блоков .L8 никогда не будет выполняться при принятии sqrt числa >= 0.0, но на ускоренном пути все еще есть несколько дополнительных инструкций.

Если вы знаете, что ваш ввод в sqrt отличен от нуля, вы можете использовать быструю, но очень приблизительную ответную инструкцию sqrt, rsqrtps (или rsqrtss для скалярной версии), Одна итерация Newton-Raphson доводит ее до почти той же точности, что и аппаратная одноточечная инструкция sqrt, но не совсем.

sqrt(x) = x * 1/sqrt(x), для x!=0, поэтому вы можете вычислить sqrt с rsqrt и умножить. Это быстро, даже на P4 (это все еще актуально в 2013 году)?

На P4 может потребоваться использовать rsqrt + итерацию Newton, чтобы заменить один sqrt, даже если вам не нужно что-либо делить на него.

См. также ответ, который я написал недавно об обработке нулей при расчете sqrt(x) как x*rsqrt(x) с помощью итерации Newton. Я включил некоторое обсуждение ошибки округления, если вы хотите преобразовать значение FP в целое число и ссылки на другие соответствующие вопросы.

P4:

• sqrtss: 23c латентность, а не конвейерная
• sqrtsd: задержка 38 с, не конвейерная
• fsqrt (x87): 43c латентность, а не конвейерная

• rsqrtss/mulss: задержка 4c + 6c. Возможно, один на 3c пропускную способность, поскольку они, очевидно, не нуждаются в одном и том же исполнительном блоке (mmx против fp).

• SIMD-упакованные версии несколько медленнее

Skylake:

• sqrtss/sqrtps: 12c латентность, по одной на 3 c пропускную способность
• sqrtsd/sqrtpd: латентность 15-16 с, пропускная способность 4-6 с
• fsqrt (x87): 14-21cc латентность, по одной на пропускную способность 4-7c
• rsqrtss/mulss: задержка 4c + 4c. Один на 1 пропускную способность.
• Именные версии SIMD 128b имеют одинаковую скорость. 256b - немного более высокая латентность, почти половина пропускной способности. Версия rsqrtss имеет полную производительность для векторов 256b.

С помощью Newton Iteration версия rsqrt не намного быстрее, чем раньше.

Числа из Экспериментальные испытания Agner Fog. Обратитесь к руководству по микрочипу, чтобы понять, что делает код быстрым или медленным. Также см. Ссылки в x86 теги wiki.

IDK, как лучше всего рассчитать sin/cos. Я читал, что аппаратное обеспечение fsin/fcos (и только немного медленнее fsincos, которое делает оба сразу) - это не самый быстрый способ, но IDK, что есть.

Ответ 7

QT имеет быстрые реализации синуса (qFastSin) и косинуса (qFastCos), которые используют справочную таблицу с интерполяцией. Я использую его в своем коде, и они работают быстрее, чем std: sin/cos и достаточно точны для того, что мне нужно (ошибка ~ = 0.01%, я думаю):

https://code.woboq.org/qt5/qtbase/src/corelib/kernel/qmath.h.html#_Z8qFastSind

#define QT_SINE_TABLE_SIZE 256


inline qreal qFastSin(qreal x)
{
   int si = int(x * (0.5 * QT_SINE_TABLE_SIZE / M_PI)); // Would be more accurate with qRound, but slower.
   qreal d = x - si * (2.0 * M_PI / QT_SINE_TABLE_SIZE);
   int ci = si + QT_SINE_TABLE_SIZE / 4;
   si &= QT_SINE_TABLE_SIZE - 1;
   ci &= QT_SINE_TABLE_SIZE - 1;
   return qt_sine_table[si] + (qt_sine_table[ci] - 0.5 * qt_sine_table[si] * d) * d;
}

inline qreal qFastCos(qreal x)
{
   int ci = int(x * (0.5 * QT_SINE_TABLE_SIZE / M_PI)); // Would be more accurate with qRound, but slower.
   qreal d = x - ci * (2.0 * M_PI / QT_SINE_TABLE_SIZE);
   int si = ci + QT_SINE_TABLE_SIZE / 4;
   si &= QT_SINE_TABLE_SIZE - 1;
   ci &= QT_SINE_TABLE_SIZE - 1;
   return qt_sine_table[si] - (qt_sine_table[ci] + 0.5 * qt_sine_table[si] * d) * d;
}

LUT и лицензию можно найти здесь: https://code.woboq.org/qt5/qtbase/src/corelib/kernel/qmath.cpp.html#qt_sine_table

Ответ 8

Я использую следующий код для вычисления тригонометрических функций с четной точностью. Константа N определяет количество требуемых бит точности (например, N = 26 даст точность точности). В зависимости от требуемой точности, предварительно вычислительное хранилище может быть небольшим и входить в кеш. Это требует только операций сложения и умножения, а также очень легко векторизовать.

Алгоритм предварительно вычисляет значения sin и cos для 0.5 ^ i, я = 1,..., N. Затем мы можем объединить эти предварительно вычислимые значения, вычислить sin и cos для любого угла до разрешения 0,5 ^ N

template <class QuadReal_t>
QuadReal_t sin(const QuadReal_t a){
  const int N=128;
  static std::vector<QuadReal_t> theta;
  static std::vector<QuadReal_t> sinval;
  static std::vector<QuadReal_t> cosval;
  if(theta.size()==0){
    #pragma omp critical (QUAD_SIN)
    if(theta.size()==0){
      theta.resize(N);
      sinval.resize(N);
      cosval.resize(N);

      QuadReal_t t=1.0;
      for(int i=0;i<N;i++){
        theta[i]=t;
        t=t*0.5;
      }

      sinval[N-1]=theta[N-1];
      cosval[N-1]=1.0-sinval[N-1]*sinval[N-1]/2;
      for(int i=N-2;i>=0;i--){
        sinval[i]=2.0*sinval[i+1]*cosval[i+1];
        cosval[i]=sqrt(1.0-sinval[i]*sinval[i]);
      }
    }
  }

  QuadReal_t t=(a<0.0?-a:a);
  QuadReal_t sval=0.0;
  QuadReal_t cval=1.0;
  for(int i=0;i<N;i++){
    while(theta[i]<=t){
      QuadReal_t sval_=sval*cosval[i]+cval*sinval[i];
      QuadReal_t cval_=cval*cosval[i]-sval*sinval[i];
      sval=sval_;
      cval=cval_;
      t=t-theta[i];
    }
  }
  return (a<0.0?-sval:sval);
}

Ответ 9

Это реализация метода серии Тейлора, ранее приведенного в akellehe answer.

unsigned int Math::SIN_LOOP = 15;
unsigned int Math::COS_LOOP = 15;

// sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...
template <class T>
T Math::sin(T x)
{
    T Sum       = 0;
    T Power     = x;
    T Sign      = 1;
    const T x2  = x * x;
    T Fact      = 1.0;
    for (unsigned int i=1; i<SIN_LOOP; i+=2)
    {
        Sum     += Sign * Power / Fact;
        Power   *= x2;
        Fact    *= (i + 1) * (i + 2);
        Sign    *= -1.0;
    }
    return Sum;
}

// cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + ...
template <class T>
T Math::cos(T x)
{
    T Sum       = x;
    T Power     = x;
    T Sign      = 1.0;
    const T x2  = x * x;
    T Fact      = 1.0;
    for (unsigned int i=3; i<COS_LOOP; i+=2)
    {
        Power   *= x2;
        Fact    *= i * (i - 1);
        Sign    *= -1.0;
        Sum     += Sign * Power / Fact;
    }
    return Sum;
}

Ответ 10

Разделяя мой код, это полином 6-й степени, ничего особенного, но переставленный, чтобы избежать pow. В Core i7 это в 2,3 раза медленнее, чем в стандартной реализации, хотя и немного быстрее для диапазона [0..2 * PI]. Для старого процессора это может быть альтернативой стандартному sin/cos.

/*
    On [-1000..+1000] range with 0.001 step average error is: +/- 0.000011, max error: +/- 0.000060
    On [-100..+100] range with 0.001 step average error is:   +/- 0.000009, max error: +/- 0.000034
    On [-10..+10] range with 0.001 step average error is:     +/- 0.000009, max error: +/- 0.000030
    Error distribution ensures there no discontinuity.
*/

const double PI          = 3.141592653589793;
const double HALF_PI     = 1.570796326794897;
const double DOUBLE_PI   = 6.283185307179586;
const double SIN_CURVE_A = 0.0415896;
const double SIN_CURVE_B = 0.00129810625032;

double cos1(double x) {
    if (x < 0) {
        int q = -x / DOUBLE_PI;
        q += 1;
        double y = q * DOUBLE_PI;
        x = -(x - y);
    }
    if (x >= DOUBLE_PI) {
        int q = x / DOUBLE_PI;
        double y = q * DOUBLE_PI;
        x = x - y;
    }
    int s = 1;
    if (x >= PI) {
        s = -1;
        x -= PI;
    }
    if (x > HALF_PI) {
        x = PI - x;
        s = -s;
    }
    double z = x * x;
    double r = z * (z * (SIN_CURVE_A - SIN_CURVE_B * z) - 0.5) + 1.0;
    if (r > 1.0) r = r - 2.0;
    if (s > 0) return r;
    else return -r;
}

double sin1(double x) {
    return cos1(x - HALF_PI);
}

Ответ 11

Просто используйте FPU со встроенным x86 для приложений Wintel. По сообщениям, прямая функция ЦП ЦП по-прежнему опережает любые другие алгоритмы по скорости. Мой пользовательский код библиотеки x86 Math предназначен для стандартного MSVC++ 2005 года и более поздних версий. Вам нужны отдельные версии float/double, если вы хотите большей точности, о которой я говорил. Иногда стратегия компилятора "__inline" выходит из строя, поэтому, чтобы быть в безопасности, вы можете удалить ее. Имея опыт, вы можете переключаться на макросы, чтобы каждый раз полностью избегать вызова функции.

extern __inline float  __fastcall fs_sin(float x);
extern __inline double __fastcall fs_Sin(double x);
extern __inline float  __fastcall fs_cos(float x);
extern __inline double __fastcall fs_Cos(double x);
extern __inline float  __fastcall fs_atan(float x);
extern __inline double __fastcall fs_Atan(double x);
extern __inline float  __fastcall fs_sqrt(float x);
extern __inline double __fastcall fs_Sqrt(double x);
extern __inline float  __fastcall fs_log(float x);
extern __inline double __fastcall fs_Log(double x);

extern __inline float __fastcall fs_sqrt(float x) { __asm {
FLD x  ;// Load/Push input value
FSQRT
}}

extern __inline double __fastcall fs_Sqrt(double x) { __asm {
FLD x  ;// Load/Push input value
FSQRT
}}

extern __inline float __fastcall fs_sin(float x) { __asm {
FLD x  ;// Load/Push input value
FSIN
}}

extern __inline double __fastcall fs_Sin(double x) { __asm {
FLD x  ;// Load/Push input value
FSIN
}}    

extern __inline float __fastcall fs_cos(float x) { __asm {
FLD x  ;// Load/Push input value
FCOS
}}

extern __inline double __fastcall fs_Cos(double x) { __asm {
FLD x  ;// Load/Push input value
FCOS
}}

extern __inline float __fastcall fs_tan(float x) { __asm {
FLD x  ;// Load/Push input value
FPTAN
}}

extern __inline double __fastcall fs_Tan(double x) { __asm {
FLD x  ;// Load/Push input value
FPTAN
}}

extern __inline float __fastcall fs_log(float x) { __asm {
FLDLN2
FLD x
FYL2X
FSTP ST(1) ;// Pop1, Pop2 occurs on return
}}

extern __inline double __fastcall fs_Log(double x) { __asm {
FLDLN2
FLD x
FYL2X
FSTP ST(1) ;// Pop1, Pop2 occurs on return
}}

Ответ 12

Более 100000000 тестов, milianw ответ на 2 раза медленнее, чем реализация std:: cos. Однако вы можете выполнить его быстрее, выполнив следующие шаги:

- > использовать float

- > не использовать пол, а static_cast

- > не использовать абс, но тернарный условный

- > использовать #define константу для деления

- > использовать макрос, чтобы избежать вызова функции

// 1 / (2 * PI)
#define FPII 0.159154943091895
//PI / 2
#define PI2 1.570796326794896619

#define _cos(x)         x *= FPII;\
                        x -= .25f + static_cast<int>(x + .25f) - 1;\
                        x *= 16.f * ((x >= 0 ? x : -x) - .5f);
#define _sin(x)         x -= PI2; _cos(x);

Более 100000000 вызов для std:: cos и _cos (x), std:: cos run на ~ 14s vs ~ 3s для _cos (x) (немного больше для _sin (x))

Ответ 13

Итак, позвольте мне перефразировать, что эта идея исходит из аппроксимации функций косинуса и синуса на интервале [-pi/4, +pi/4] с ограниченной ошибкой с использованием алгоритма Ремеза. Затем с использованием остатка с плавающей точкой с уменьшенным диапазоном и LUT для выходных значений cos & sine целочисленного отношения, аппроксимация может быть перемещена к любому угловому аргументу.

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

void sincos_fast(float x, float *pS, float *pC){
     float cosOff4LUT[] = { 0x1.000000p+00,  0x1.6A09E6p-01,  0x0.000000p+00, -0x1.6A09E6p-01, -0x1.000000p+00, -0x1.6A09E6p-01,  0x0.000000p+00,  0x1.6A09E6p-01 };

    int     m, ms, mc;
    float   xI, xR, xR2;
    float   c, s, cy, sy;

    // Cody & Waite range reduction Algorithm, [-pi/4, pi/4]
    xI  = floorf(x * 0x1.45F306p+00 + 0.5);
    xR  = (x - xI * 0x1.920000p-01) - xI*0x1.FB5444p-13;
    m   = (int) xI;
    xR2 = xR*xR;

    // Find cosine & sine index for angle offsets indices
    mc = (  m  ) & 0x7;     // two complement permits upper modulus for negative numbers =P
    ms = (m + 6) & 0x7;     // two complement permits upper modulus for negative numbers =P, note phase correction for sine.

    // Find cosine & sine
    cy = cosOff4LUT[mc];     // Load angle offset neighborhood cosine value 
    sy = cosOff4LUT[ms];     // Load angle offset neighborhood sine value 

    c = 0xf.ff79fp-4 + xR2 * (-0x7.e58e9p-4);               // TOL = 1.2786e-4
    // c = 0xf.ffffdp-4 + xR2 * (-0x7.ffebep-4 + xR2 * 0xa.956a9p-8);  // TOL = 1.7882e-7

     s = xR * (0xf.ffbf7p-4 + x2 * (-0x2.a41d0cp-4));   // TOL = 4.835251e-6
    // s = xR * (0xf.fffffp-4 + xR2 * (-0x2.aaa65cp-4 + xR2 * 0x2.1ea25p-8));  // TOL = 1.1841e-8

     *pC = c*cy - s*sy;     
    *pS = c*sy + s*cy;
}

float sqrt_fast(float x){
    union {float f; int i; } X, Y;
    float ScOff;
    uint8_t e;

    X.f = x;
    e = (X.i >> 23);           // f.SFPbits.e;

    if(x <= 0) return(0.0f);

    ScOff = ((e & 1) != 0) ? 1.0f : 0x1.6a09e6p0;  // NOTE: If exp=EVEN, b/c (exp-127) a (EVEN - ODD) := ODD; but a (ODD - ODD) := EVEN!!

    e = ((e + 127) >> 1);                            // NOTE: If exp=ODD,  b/c (exp-127) then flr((exp-127)/2)
    X.i = (X.i & ((1uL << 23) - 1)) | (0x7F << 23);  // Mask mantissa, force exponent to zero.
    Y.i = (((uint32_t) e) << 23);

    // Error grows with square root of the exponent. Unfortunately no work around like inverse square root... :(
    // Y.f *= ScOff * (0x9.5f61ap-4 + X.f*(0x6.a09e68p-4));        // Error = +-1.78e-2 * 2^(flr(log2(x)/2))
    // Y.f *= ScOff * (0x7.2181d8p-4 + X.f*(0xa.05406p-4 + X.f*(-0x1.23a14cp-4)));      // Error = +-7.64e-5 * 2^(flr(log2(x)/2))
    // Y.f *= ScOff * (0x5.f10e7p-4 + X.f*(0xc.8f2p-4 +X.f*(-0x2.e41a4cp-4 + X.f*(0x6.441e6p-8))));     // Error =  8.21e-5 * 2^(flr(log2(x)/2))
    // Y.f *= ScOff * (0x5.32eb88p-4 + X.f*(0xe.abbf5p-4 + X.f*(-0x5.18ee2p-4 + X.f*(0x1.655efp-4 + X.f*(-0x2.b11518p-8)))));   // Error = +-9.92e-6 * 2^(flr(log2(x)/2))
    // Y.f *= ScOff * (0x4.adde5p-4 + X.f*(0x1.08448cp0 + X.f*(-0x7.ae1248p-4 + X.f*(0x3.2cf7a8p-4 + X.f*(-0xc.5c1e2p-8 + X.f*(0x1.4b6dp-8))))));   // Error = +-1.38e-6 * 2^(flr(log2(x)/2))
    // Y.f *= ScOff * (0x4.4a17fp-4 + X.f*(0x1.22d44p0 + X.f*(-0xa.972e8p-4 + X.f*(0x5.dd53fp-4 + X.f*(-0x2.273c08p-4 + X.f*(0x7.466cb8p-8 + X.f*(-0xa.ac00ep-12)))))));    // Error = +-2.9e-7 * 2^(flr(log2(x)/2))
    Y.f *= ScOff * (0x3.fbb3e8p-4 + X.f*(0x1.3b2a3cp0 + X.f*(-0xd.cbb39p-4 + X.f*(0x9.9444ep-4 + X.f*(-0x4.b5ea38p-4 + X.f*(0x1.802f9ep-4 + X.f*(-0x4.6f0adp-8 + X.f*(0x5.c24a28p-12 ))))))));   // Error = +-2.7e-6 * 2^(flr(log2(x)/2))

    return(Y.f);
}