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

_declspec(naked) float _vectorcall ln1(float x)
{
static const float ct[8] =          // Таблица констант
  {
    1.0f,                           // a0 = 1
    0.693147182f,                   // b0 = ln 2
    0.343713552f,                   // a2
    0.606786430f,                   // b2
    1.23005533f,                    // a1
    1.19918299f,                    // b1
    0.00755311921f,                 // a3
    0.0821869969f                   // b3
  };
  _asm
  {
    vmovd eax,xmm0                  // Прочитать в eax двоичное представление числа x
    mov edx,-126                    // В edx - смещение для коррекции порядка нормализованных float
    cmp eax,0x7F800000              // Сравнить число x со значением Inf
    jge ln_exit                     // Если x=+Inf или +NaN, вернуть аргумент x без изменений
    rol eax,9                       // Циклически сдвинуть число так, чтобы порядок оказался в al
    movzx ecx,al                    // Получить в ecx порядок числа со смещением +127
    jecxz ln_denorm                 // Перейти, если x=0 или x денормализованное (т.е. близко к 0)
      ln_calc:                      // Точка входа для вычисления логарифма после нормализации числа 
    cmp eax,0x55555600              // Сравнить мантиссу с 4/3
    setb al                         // al=1, если мантисса меньше 4/3; иначе al=0
    sub cl,al                       // ecx=k+126, где k - добавка к двоичному логарифму (целая часть)
    or al,126                       // Изменить порядок x так, чтобы привести к диапазону 2/3<=x<4/3
    add ecx,edx                     // ecx=k
    ror eax,9                       // Получить в eax двоичное представление числа и проверить его знак
    jc ln_nan                       // Перейти, если x<0 (аргумент некорректный); продолжить для x>0
    mov edx,offset ct               // В edx - адрес таблицы констант
    vmovd xmm0,eax                  // Перенести из eax в xmm0 приведённое значение x
    vmovups xmm1,[edx]              // xmm1 = b2 : a2 | b0 : a0
    vsubss xmm0,xmm0,xmm1           // xmm0 = x-1 = t
    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
    vfmadd231ps xmm1,xmm0,[edx+16]  // xmm1 = b3*t+b2 : a3*t+a2 | b1*t+b0 : a1*t+a0
    vmovhlps xmm3,xmm3,xmm1         // xmm3 = b3*t+b2 : a3*t+a2
    vfmadd231ps 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,[edx+4]        // Перевести двоичный логарифм в натуральный
    ret                             // Вернуться
      ln_cont:                      // Продолжаем в штатном режиме, когда денормализованное x>0
    lea edx,[edx+2*ecx-31]          // Обновить смещение порядка для денормализованного числа
    neg ecx                         // 32-ecx - это необходимое количество сдвигов мантиссы влево
    shl eax,cl                      // Сформировать мантиссу нормализованного числа
    jmp ln_calc                     // Число нормализовано. Перейти для вычисления логарифма
      ln_denorm:                    // Точка входа для обработки денормализованного x, в т.ч. x=0
    btr eax,8                       // Прочитать в cf и сбросить знаковый бит числа
    bsr ecx,eax                     // Найти номер старшего единичного бита; установить zf, если x=0
    ja ln_cont                      // Перейти, если денормализованное x>0
    jz ln_inf                       // Если x=0, перейти для формирования результата -Inf
      ln_nan:                       // Точка входа для формирования результата NaN, если x<0 (или -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                             // Вернуться
  }
}