Целочисленный счёт повышенной точности против аппаратного FPU

[Aether]

Добрый день.

В теме: https://forum.shelek.ru/index.php/topic,30591.0.html. Было предложено использовать для хранения результатов вычислений целочисленный формат в 32 бита, а сам счёт производить в формате с плавающей запятой.

Естественно, стало интересно: что быстрее? Преобразовывать числа в double или считать сразу целыми, но с повышенной точностью.

Для тестирования я взял простую формулу расчёта гипотенузы по двум катетам - теорему Пифагора. Итак, функция имеет вид: c = sqrt(a * a + b * b). Как можно заметить, промежуточный аргумент - a * a - имеет разрядность уже не 32, а 64 бита, однако, конечный результат представляется вновь 32-х битным.

Модуль тестирования "a.c"

Код: (C)

//#define DEBUG

#include <locale.h>
#include <stdio.h>
#include <malloc.h>
#include <time.h>

#define TEST_BASE 100000000 //Число вычислений
#define MAX_UTT 2000000000  //4294967295

typedef unsigned int UTT;

int __cdecl _pifagor_fpu(UTT x, UTT y);
int __cdecl _pifagor_hpi(UTT x, UTT y);

int main(void)
{
   #ifdef DEBUG
      FILE* f = fopen("debug.txt", "w");
   #endif

setlocale(LC_ALL, "");
printf("Начало тестирования.\n");

printf("Выделение памяти.\n");
UTT* p_fpu = malloc(TEST_BASE * sizeof(UTT));
if (p_fpu == (void*)0)
   {
   printf("Ошибка выделения памяти. p_fpu");
   return -1;
   }

UTT* p_hpi = malloc(TEST_BASE * sizeof(UTT));
if (p_fpu == (void*)0)
   {
   printf("Ошибка выделения памяти. p_hpi");
   return -1;
   }

printf("Пересчёт значений с использованием FPU.\n");
time_t start_fpu = time(NULL);
UTT* temp = p_fpu;
for (int i = 0; i < TEST_BASE; i++)
   {
   UTT x = i * (MAX_UTT / TEST_BASE);
   UTT y = MAX_UTT - x;
   *temp = _pifagor_fpu(x, y);
   temp++;
   }
time_t end_fpu = time(NULL);

printf("Пересчёт значений с использованием HPI.\n");
time_t start_hpi = time(NULL);
temp = p_hpi;
for (int i = 0; i < TEST_BASE; i++)
   {
   UTT x = i * (MAX_UTT / TEST_BASE);
   UTT y = MAX_UTT - x;
   *temp = _pifagor_hpi(x, y);
   temp++;
   }
time_t end_hpi = time(NULL);

printf("Проверка погрешности.\n");
temp = p_fpu;
UTT* temp_2 = p_hpi;

UTT max_delta = 0;
int max_i = 0;
for (int i = 0; i < TEST_BASE; i++)
   {
   UTT x = i * (MAX_UTT / TEST_BASE);
   UTT y = MAX_UTT - x;
   UTT delta =
   (*temp > *temp_2)?(*temp - *temp_2):(*temp_2 - *temp);
   if (delta > max_delta) {max_delta = delta; max_i = i;}

   #ifdef DEBUG
      fprintf(f, "x: %10u :: y: %10u :: fpu: %10u :: hpi: %10u\n",
              x, y, *temp, *temp_2);
   #endif

   temp++;
   temp_2++;
   }

printf("Максимальная погрешность = %u\n", max_delta);
temp = p_fpu + max_i;
temp_2 = p_hpi + max_i;
printf("На шаге %u >> \n FPU = %u\n HPI = %u\n",
       max_i, *temp, *temp_2);

printf("Время старта FPU      = %u\n", start_fpu);
printf("Время завершения FPU  = %u\n", end_fpu);
time_t d_fpu = end_fpu - start_fpu;
printf("Время работы FPU      = %u\n\n", d_fpu);

printf("Время старта HPI      = %u\n", start_hpi);
printf("Время завершения HPI  = %u\n", end_hpi);
time_t d_hpi = end_hpi - start_hpi;
printf("Время работы HPI      = %u\n\n", d_hpi);

free(p_hpi);
free(p_fpu);

   #ifdef DEBUG
      fclose(f);
   #endif

return 0;
}

Тестируемый модуль "b.asm"

Код: (ASM)

use32
global _pifagor_fpu

_pifagor_fpu:
   push ebp
   mov ebp, esp

 ; Функция: __cdecl Стек:
 ; [ebp + 0]  - значение ebp
 ; [ebp + 4]  - адрес возврата для ret
 ; [ebp + 8]  - параметр х
 ; [ebp + 12] - параметр у

   finit
   fild dword [ebp + 12]
   fmul st0, st0
   fild dword [ebp + 8]
   fmul st0, st0
   fadd st0, st1
   fsqrt
   fist dword [ebp - 4]
   mov eax, dword [ebp - 4]
   leave
   ret

global _pifagor_hpi

_pifagor_hpi:
   push ebx
   push ecx
   push edx
   push esi
   push edi
   push ebp
   mov ebp, esp
   
 ; Функция: __cdecl Стек:
 ; [ebp + 0]  - значение ebp
 ; [ebp + 4]  - значение edi
 ; [ebp + 8]  - значение esi
 ; [ebp + 12] - значение edx
 ; [ebp + 16] - значение ecx
 ; [ebp + 20] - значение ebx
 ; [ebp + 24] - адрес возврата для ret
 ; [ebp + 28] - параметр x
 ; [ebp + 32] - параметр y

   mov eax, dword [ebp + 32]
   mul eax
   mov esi, eax
   mov edi, edx

   mov eax, dword [ebp + 28]
   mul eax
   add esi, eax
   adc edi, edx

   xor ebx, ebx
   xor ecx, ecx
   xor ebp, ebp
   dec ecx

sqrt_loop:
   mov eax, ebx
   add eax, ecx
   rcr eax, 1
   cmp eax, ebp
   je sqrt_out
   mov ebp, eax
   mul eax
   cmp edx, edi
   je sqrt_esi
   jb sqrt_low
sqrt_high:
   mov ecx, ebp
   jmp sqrt_loop  
sqrt_low:
   mov ebx, ebp
   jmp sqrt_loop
sqrt_esi:
   cmp eax, esi
   je sqrt_out
   jb sqrt_low
   jmp sqrt_high
sqrt_out:
   mov eax, ebp
   pop ebp
   add sp, 20
   ret

Сборка "make.bat"

Код:

@nasm b.asm -f elf32
@tcc b.o a.c -o a.exe
@a.exe
@PAUSE

Результат

Начало тестирования.
Выделение памяти.
Пересчёт значений с использованием FPU.
Пересчёт значений с использованием HPI.
Проверка погрешности.
Максимальная погрешность = 1
На шаге 2237 >>
 FPU = 1999955261
 HPI = 1999955260
Время старта FPU      = 1469983664
Время завершения FPU  = 1469983671
Время работы FPU      = 7

Время старта HPI      = 1469983671
Время завершения HPI  = 1469983717
Время работы HPI      = 46

То есть на процессоре Intel n3700 с аппаратным FPU всё таки выгоднее его использовать, и чем более сложные вычисления, тем это более актуально. Прошу высказать Ваше мнение, указать на ошибки, возможно, предложить меры по ускорению процессов и их более правильному написанию.

[RXL]

1. Почему использован 32-битный код, а не 64-битный?
2. Зачем столько пушей в _pifagor_hpi?
3. Зачем finit в _pifagor_fpu?
4. Зачем enter, leave и адресация по ebp, если можно сразу по esp? Кстати, в _pifagor_fpu вместо enter использован эквивалент: push ebp / mov ebp, esp. А почему тогда в _pifagor_hpi использован еще один подход.
5. В _pifagor_fpu обращение к стеку по незарезервированному месту ([ebp - 4]). Переключение задачи между fist и последующим mov перезатрет значение.

По моему, подчистив можно серьезно ускорить. И нужен другой метод измерения: time не учитывает переключение задач.

[Aether]

1. Почему использован 32-битный код, а не 64-битный?
2. Зачем столько пушей в _pifagor_hpi?
3. Зачем finit в _pifagor_fpu?
4. Зачем enter, leave и адресация по ebp, если можно сразу по esp? Кстати, в _pifagor_fpu вместо enter использован эквивалент: push ebp / mov ebp, esp. А почему тогда в _pifagor_hpi использован еще один подход.
5. В _pifagor_fpu обращение к стеку по незарезервированному месту ([ebp - 4]). Переключение задачи между fist и последующим mov перезатрет значение.

По моему, подчистив можно серьезно ускорить. И нужен другой метод измерения: time не учитывает переключение задач.

1. У меня в практике ещё много 32-х битных аппаратов. И опыта с 64 битами мало, возможно, стоит попробовать реализовать, но после разбора с 32 битами.
2. Перестраховался: мало ли компилятор С перед вызовом процедуры что-то будет держать в регистрах? Он обязан сделать что-то, вроде, pushad или нет? Кстати, сопутствующее, если изменить передачу параметров на __fastcall, то можно ли будет уйти от стека вообще?
3. Без finit вычисления были очень замедлены, и через 10 вычислений начиналась ошибка. Первым делом решил переинициализировать - помогло, но причина самому интересна.
4. Не знал, как события развернутся. Сначала думал регистров не хватит, придётся использовать стек, а потом оказалось всё в порядке, хотя была бы третья или четвёртая степень, например, для сплайна, тут уже не развернёшься. enter я подзабыл, можно записать enter 4, 0, если правильно помню первый параметр, как раз аренда под локальные переменные. В hpi ebp играет роль сначала указателя, а потом переменной. leave сначала поместит в esp ebp, что недопустимо. Поэтому в конце программы сначала извлекается истинный ebp, затем сдвигаю esp, чтобы удалить то, что было в пушах.
5. Сорри.

Какой метод измерения лучше?

Добавлено через 32 минуты и 19 секунд:
__fastcall в tcc работает, но с особенностями.

_pifagor_fpu

Код: (ASM)

_pifagor_fpu:

 ; Функция: __fastcall
 ; ecx - параметр х
 ; edx - параметр у

   push ecx
   push edx
   finit
   fild dword [esp + 4]
   fmul st0, st0
   fild dword [esp]
   fmul st0, st0
   fadd st0, st1
   fsqrt
   fist dword [esp + 4]
   pop eax
   pop eax  
   ret

Удаление finit, также приводит к ошибкам.

[RXL]

fist оставляет значение в стеке FPU: используй fistp. Аналогично с fadd: используй faddp.

Выделять/очищать стек лучше через sub/add esp. В отличии от push не пишет в память. Если конечно не нужно инициализировать.

Время с точностью до такта можно померить посредством rdtsc. Но прерывание другими задачами оно не учитывает.
Можно посчитать чисто процессорное время:
 для винды (точность 100 нс): https://msdn.microsoft.com/ru-ru/library/windows/desktop/ms683223%28v=vs.85%29.aspx
 для posix: (точность 1 мкс): http://pubs.opengroup.org/onlinepubs/009695399/functions/times.html

[Sla]

ППЦ, над х-нёй задумывайтесь

Есть исполнительный механизм, есть его точность . С такой точностью, должны делаться вычисления!

Я не Смотрел в код!!!! я обращаю внимание на ОСНОВЫ!!!

[Aether]

Спасибо за рекомендации, та же база испытаний улучшила результат с (7...8 - 45...48) до (3...4 - 42...44). Полагаю, всё же сильная просадка в HPI идёт при извлечении корня - алгоритм последовательного приближения довольно медленный. Ряд Тейлора тоже не вариант - деление тяжёлая операция. Задумываюсь, что, возможно, можно сделать как-то табличный пересчёт, потратить, скажем 65Кб памяти в угоду скорости?

Есть исполнительный механизм, есть его точность . С такой точностью, должны делаться вычисления!
я обращаю внимание на ОСНОВЫ!!!

Sla, приведу пример про то, как малая погрешность приводит к значительному отклонению. Вам знаком пузырьковый уровень для строительных работ - такой кусочек алюминиевого профиля длиной от полуметра до двух? Так вот, его качество в лучшем случае даёт ошибку в 1 мм на 1 м, согласен, на 1 м это не так уж и много. Теперь представим, что дом имеет длину в 100 м. На этой дистанции результат будет уже 100 мм, плюс погрешность перестановок, как Вам перекос дома в полтора метра? Естественно, для решения этой задачи используется другой инструмент. Задача математической модели: решать одинаково хорошо широкий спектр задач. А это лишь небольшие соображения по поводу частных вопросов.

[RXL]

А чем FPU не устраивает? Мантиса 64 бита для временного значения.

[Aether]

А чем FPU не устраивает? Мантиса 64 бита для временного значения.

Просто хотел взглянуть на картину в цифрах - практика лучший критерий. Теперь есть лучшее понимание, что и где стоит применять.

Ещё два вопроса:
1. Что происходит с регистрами FPU во время, например, переключения задач? В том числе и с регистрами управления? Какие конфликты возможны.
2. У него есть особые состояния число NaN и т.д. Приходилось ли сталкиваться практически с их разрешением, в каких задачах их использование стоит предусматривать?

[RXL]

1. При переключении задач контекст FPU аппаратно не сохраняется. Видимо по старой архитектуре 286 и 386, где FPU был сопроцессором. Сохранение/восстановление контекста реализуется операционкой. Т.ч. можно считать, что FPU в полном распоряжении твоей задачи.

2. Операции с NaN приведут к исключению. Как можно получить NaN. Например, корень из отрицательного значения.

[darkelf]

1. При переключении задач контекст FPU аппаратно не сохраняется. Видимо по старой архитектуре 286 и 386, где FPU был сопроцессором. Сохранение/восстановление контекста реализуется операционкой. Т.ч. можно считать, что FPU в полном распоряжении твоей задачи.

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

В итоге - с FPU могут работать параллельно множество задач, но сохранение/восстановление его регистров делают при помощи исключений: запоминается процесс текущий владелец FPU, при попытке использования сопроцессора другим процессом происходит исключение и тогда выполняется сохранение контекста исполнения FPU в описателе текущего процесса-владельца, загрузка значений из для нового процесса - владельца FPU и запоминание, что теперь уже этот процесс владеет FPU. В итоге получается, что если владелец FPU не менялся, то и лишних сохранений/восстановлений не происходит. Но как и было сказано - да, для задачи создаётся иллюзия полного владения ресурсами FPU.