Быстрая функция вычисления логарифма с одиночной точностью

// Функция быстрого вычисления натурального логарифма с двойной точностью
// Для повышения производительности алгоритм не поддерживает денормализованные числа (они считаются нулём)
double ln(double x)
{
  static const double ct[10] =          // Таблица констант
  {
    1.414213562373095049,               // 2^0.5
    2.885390081777926811,               // b0
    0.9617966939259898017,              // b1
    0.5770780163454805578,              // b2
    0.4121985858489949842,              // b3
    0.3205985341276662424,              // b4
    0.2623343917068792287,              // b5
    0.2209121151421818382,              // b6
    0.2136683673423465570,              // b7
    0.6931471805599453094               // ln 2
  };
  _asm
  {
    movsd xmm0, [x]                     // xmm0 = x
    vmovsd xmm0, [x]                    // xmm0 = x
    or edx, -1                          // edx = -1
    vpextrw eax, xmm0, 3                // Прочитать в ax старшее слово числа, содержащее знаковый бит и порядок
    sahf                                // Проверить знаковый бит числа
    js ln_break                         // Перейти для обработки ошибки, если x <= -0
    mov ecx, 0x00003FF0                 // Выставить в ch страший байт порядка и подготовить маску в cl
    or cl, al                           // Выставить младший байт порядка x
    shr eax, 4                          // Получить в eax порядок числа плюс смещение 1023
    jnz ln_ok                           // Перейти, если x>0 и не является денормализованным
    mov dl, 240                         // Подготовить в dx старшее слово значения -Inf
      ln_break :                        // Здесь формируется результат в случае ошибки
    vmovd xmm0, edx                     // Поместить edx в младшую часть xmm0
    vpsllq xmm0, xmm0, 48               // Сдвинуть слово 0xFFFF или 0xFFF0 в старшую часть числа
    jmp ln_end                          // Перейти в конец функции
      ln_ok :                           // Продолжаем с допустимым значением x > 0
    vcvtsi2sd xmm2, xmm2, edx           // xmm2 = -1
    vpinsrw xmm0, xmm0, ecx, 3          // Заменить порядок числа x. Теперь 1 <= x < 2
    mov edx, offset ct                  // В edx адрес таблицы констант
    sub eax, 1023                       // Вычислить в eax порядок без смещения; это будет добавка к логарифму
    vcomisd xmm0, [edx]                 // Сравнить x и sqrt(2)
    jc ln_x_done                        // Перейти, если 1<=x<sqrt(2), т.е. не требуется коррекция порядка
    and ecx, -17                        // Вычесть из порядка 1, чтобы стало sqrt(2)/2 <= x < 1
    inc eax                             // Добавка к логарифму соотвтетственно увеличивается на 1
    vpinsrw xmm0, xmm0, ecx, 3          // Обновить значение xmm0
      ln_x_done :                       // В xmm0 приведённое значение x, такое что sqrt(2)/2 <= x < sqrt(2)
    vcvtsi2sd xmm3, xmm3, eax           // В xmm3 добавка к двоичному логарифму
    vsubsd xmm1, xmm0, xmm2             // xmm1 = x+1
    vaddsd xmm0, xmm0, xmm2             // xmm0 = x-1
    shr ecx, 11                         // ecx = 7 - инициализируем счётчик цикла
    vdivsd xmm0, xmm0, xmm1             // xmm0 = (x-1)/(x+1) = u
    vmovsd xmm1, [edx+64]               // Инициализировать сумму младшим коэффициентом: xmm1 = b7
    vmulsd xmm2, xmm0, xmm0             // xmm2 = u^2
      ln_loop :                         // Цикл вычисления полинома (ecx пробегает от 7 до 1)
    vfmadd213sd xmm1, xmm2, [edx+ecx*8] // Обновляем сумму ряда (схема Горнера)
    loop ln_loop                        // Учесть все 8 коэффициентов
    vfmadd213sd xmm0, xmm1, xmm3        // Получить в xmm0 двоичный логарифм и прибавить к нему порядок
    vmulsd xmm0, xmm0, [edx+72]         // Перевести двоичный логарифм в натуральный
      ln_end :                          // Перенос результата в стек сопроцессора в соответствии с соглашениями среды
    vmovsd [x], xmm0                    // Сохранить результат на месте переменной x
    fld [x]                             // Затем загрузить его в стек FPU
  }
}