Трюки с float: быстрое вычисление логарифмов. (комментарии) (9 стр)

_declspec(naked) float _vectorcall ln1(float x)
{
  static const float ct[10] =           // Таблица констант
  {
    0.00830399152f,                     // a3
    0.0903585330f,                      // b3
    1.35235262f,                        // a1
    1.31841099f,                        // b1
    0.377887428f,                       // a2
    0.667115986f,                       // b2
    1.09942472f,                        // a0
    0.762063146f,                       // b0
    1.0f,                               // 1
    0.693147181f                        // ln 2
  };
  _asm
  {
    mov edx,-126                        // В edx - смещение для коррекции порядка нормализованных float
    vmovd eax,xmm0                      // Прочитать в eax двоичное представление числа x
    cmp eax,0x7F800000                  // Сравнить число x со значением Inf
    jge ln_exit                         // Если x=+Inf или +NaN, вернуть аргумент x без изменений
    ror eax,23                          // Циклически сдвинуть число так, чтобы порядок оказался в al
    movzx ecx, al                       // Получить в ecx порядок числа со смещением +127
    jecxz ln_denorm                     // Перейти, если x=0 или x денормализованное (т.е. близко к 0)
      ln_calc:                          // Точка входа для вычисления логарифма после нормализации числа 
    cmp eax,0x570A3E00                  // Сравнить мантиссу с 1.34
    setb al                             // al=1, если мантисса меньше 1.34; иначе al=0
    sub cl,al                           // ecx=k+126, где k - добавка к двоичному логарифму (целая часть)
    or al,126                           // Изменить порядок x так, чтобы привести к диапазону 0.67<=x<1.34
    add ecx,edx                         // ecx=k
    ror eax,9                           // Получить в eax двоичное представление числа и проверить его знак
    jnc ln_cont                         // Перейти, если x>0, т.е. аргумент корректный. Продолжить, если x<0
      ln_nan:                           // Точка входа для формирования результата NaN
    sbb eax,eax                         // eax=-1 (значение NaN) в случае x<0
      ln_inf:                           // Точка входа для формирования результата -Inf в случае x=0
    or eax,0xFF800000                   // Сформировать в eax значение -Inf, если x=0; или NaN, если x<0
    vmovd xmm0,eax                      // Перенести результат из eax в xmm0
      ln_exit:                          // В xmm0 готов результат для нештатного значения x
    ret                                 // Вернуться
      ln_denorm:                        // Точка входа для обработки денормализованного x, в т.ч. x=0
    btr eax,8                           // Прочитать в cf и сбросить знаковый бит числа
    bsr ecx,eax                         // Найти номер старшего единичного бита; установить zf, если x=0
    jz ln_inf                           // Если x=0, перейти для формирования результата -Inf
    jc ln_nan                           // Если x<0 (или -NaN), перейти для формирования результата NaN
    lea edx,[edx+2*ecx-31]              // Обновить смещение порядка для денормализованного числа
    neg ecx                             // 32-ecx - это необходимое количество сдвигов мантиссы влево
    shl eax,cl                          // Сформировать мантиссу нормализованного числа
    jmp ln_calc                         // Число нормализовано. Перейти для вычисления логарифма
      ln_cont:                          // Продолжаем вычисления в штатном случае, когда x>0
    vmovd xmm0, eax                     // Перенести из eax в xmm0 приведённое значение x
    mov eax, offset ct                  // В eax адрес таблицы констант
    vsubss xmm0, xmm0, [eax+32]         // xmm0 = x-1 = t
    vmovups xmm1,[eax]                  // xmm1 = b1 : a1 | b3 : a3
    vshufps xmm0,xmm0,xmm0,0            // xmm0 = t : t | t : t
    vcvtsi2ss xmm4,xmm4,ecx             // xmm4 = k
    vmulps xmm2,xmm0,xmm0               // xmm2 = t^2 : t^2 | t^2 : t^2
    vfmadd213ps xmm1, xmm0, [eax+16]    // xmm1 = b1*t+b0 : a1*t+a0 | b3*t+b2 : a3*t+a2
    vmovhlps xmm3,xmm3,xmm1             // xmm3 = b1*t+b0 : a1*t+a0
    vfmadd213ps xmm1, xmm2, xmm3        // xmm1 = b3*t^3+b2*t^2+b1*t+b0 : a3*t^3+a2*t^2+a1*t+a0 = Q : P
    vmovshdup xmm2,xmm1                 // xmm2 = Q
    vdivss xmm1,xmm1,xmm2               // xmm1 = P/Q
    vfmadd213ss xmm0, xmm1, xmm4        // xmm1 = P/Q*t+k - готовый двоичный логарифм
    vmulss xmm0, xmm0, [eax+36]         // Перевести двоичный логарифм в натуральный
    ret                                 // Вернуться
  }
}