Комментарии 23

А нельзя было по порядку определить примерное число десяток, на которые следует умножить/разделить число (порядок — целая часть логарифма, а значит, зависимость будет линейная)? Это, по идее, позволило бы уменьшить число циклов (за счет применения возведения 10 в нужную степень: цикл останется, но меньший по размеру). А то в double десятичный порядок может достигать 308, за вычетом (для отрицательных порядков — с прибавлением) 15-16 десятичных цифр, получается довольно большое количество итераций.
Т.е., условно, мы сперва "сдвигаем" число так, чтобы оно получалось примерно целым, затем используем описанный алгоритм.
Тестировать алгоритм, имеющий цикл, зависящий от порядка, т.е., по сути, сложность O(порядок), только на одном на конкретном непредельном значении с небольшим порядком, не совсем правильно.

интересный подход, понравилось удаление периодических 9-ок в конце, но есть нюанс — в 64-битном режиме на x86 используется SSE+ для арифметики с плавающей точкой. В принципе и в 32-битном режиме тоже можно использовать SSE. Все использованные инструкции в коде имеют аналоги в SSE (FBSTP, увы нет) и позволяют напрямую адресовать регистры, вместо стека. FXTRACT можно в целочисленных регистрах сделать, и вместо пошагового умножения/деления на 10 сделать «пристрелку» по степени 2 — умножить на константу log10 (2), получить степень 10ки, возвести 10 в эту степень по таблице или быстрым способом и скорректировать число. деления вообще заметно медленнее, вместо деления на 10 быстрее умножать на 1/10, но теряется точность.
я не заметил ваш комментарий, но выбрал именно такой вариант, внизу кусочек кода который работает таким образом.
Денормализованные числа, а также всякие nan, inf алгоритм не поддерживает? Почему выбрано 16 знаков после запятой? Ведь точное представление double (64-бита) может иметь гораздо больше знаков: например 0.1 из-за непредставимости будет иметь для double точное значение 1.000000000000000055511151231257827021181583404541015625e-1 = 0x3FB999999999999A (ну или 9.999999999999999167332731531132594682276248931884765625e-2 = 0x3FB9999999999999 если стоит округление к -inf), при этом там как раз 16 нулей/девяток, но после них идет еще хвост, который хочется выкинуть, но по хорошему нельзя, ведь он бывает важен. Или цель была получить как можно более скоростной алгоритм?
Откуда столько знаков? Под мантиссу отводится 53 бита если перевести 53 единицы в десятичное число получится 9007199254740991, как раз 16 цифр. Внутренне представление FPU 80 бит, но и там столько десятичных цифр не набирается
Откуда столько знаков? Под мантиссу отводится 53 бита если перевести 53 единицы в десятичное число получится 9007199254740991, как раз 16 цифр.

Причем тут перевод 53 единиц в десятичное число?
Возьмем пример проще — абсолютно точно представимое число даже для 32-х битного float = 2**-45 (0x29000000). А теперь возведем 2 в степень -45 и посчитаем сколько это: 2.8421709430404007434844970703125e-14, как видите знаков сильно больше 16-ти, и все точные, т.к. число представимо точно. Вообще 32-х битные float могут иметь до 105 точных знаков 2**-149 = 1.40129846432481707092372958328991613128026194187651577175706828388979108268586060148663818836212158203125e-45, а уж 64-битный double так вообще 752 (2**-1074), а уж если брать 80-битный long double, то и вовсе 11497! Не буду копировать это число, текст будет слишком большим :)

Вы исходите из изначально ложного предположения, что в float хранится точное число. Но это не так. Вы знаете только первые 53 бит мантиссы, но не знаете значение 54-го бита. Значению (0x29000000), как и любому другому, будет соответствовать целый диапазон десятичных чисел, а не какое-то одно число. Поэтому логично из этого диапазона выбирать число, имеющее минимальное количество значащих цифр.

Вы знаете только первые 53 бит мантиссы, но не знаете значение 54-го бита.

Согласен, но это нисколько не противоречит моему примеру выше. В примере 2**-45 вполне представимое точно число (т.е. оно только одно соответствует этой записи, а не какой-нибудь диапазон чисел), использующее вообще только 1 бит мантиссы (о чем и говорит запись в сыром виде — 0x29000000), тем не менее имеющее более 16 десятичных знаков при точном переводе из двоичного вида. На то она и «плавающая точка», что имеет «плавающую точность», зависящую от диапазона переводимых чисел — от 8 до 105 знаков (для 32-х битных float). Если не верите мне, смотрите кучу примеров в сети, например randomascii.wordpress.com/2012/03/08/float-precisionfrom-zero-to-100-digits-2

Утверждение, что число 2⁻⁴⁵ может быть представимо без округления в двоичном формате — верное, с этим никто не спорит. Да, в принципе, для любого двоичного значения существует десятичное, которое может быть переведено в двоичное без округления. Вот только практической пользы в этом мало.

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


The flip side of this question is figuring out how many decimal digits it takes to uniquely identify a float. Again, we aren’t concerned here with converting the exact value of the float to a decimal (we’ll get to that), but merely having enough digits to uniquely identify a particular float.

И далее выводится, что для одинарной точности "enough digits" — это 9 для худшего случая, для некоторых может быть меньше (для 0.5, например, хватает 2).


Теоретически, конечно, может потребоваться полное представление, но это только если число из строки читается в формат более высокой разрядности, чем тот, в котором строка была выведена — тогда могут, действительно, быть неоднозначности (ближайшее к 0.1 число в одинарной точности не то же самое, что в двойной). Но если выводить тупо 16 десятичных цифр, то это проявляется при чтении в четверную точность или выше, что крайне редко на практике.

Не хватает сравнения с популярными реализациями C/C++/Rust и внезапного обнаружения того факта, что они работают быстрее.

Ну да. Вычислить десятичный логарифм одной FPU операцией — порядка 50 тактов. Автор же зачем-то для этого дела приспособил цикл, где одное только деление — около 15 тактов, а суммарно будет порядка 25–30 на один порядок: это почти 10000 тактов, если на входе будет 1e308, а если это extended (80-bit) формат, т.е. до 10^4932, так вообще под 150000 тактов.


Вычислять порядок

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

Ну так считаете логарифм, берёте от него целую часть, потенцируете, делите одно на другое — и опа, мантисса уже в нужном виде.

Вроде такого?
Ноутбук с Intel® Core(TM) i5-8265U CPU @ 1.60GHz
Python 3.7 (тот язык что "тормозной")


>>> timeit.timeit('f"{n:.16E}"',"n=1234.567890", number=100000000)
86.52470250002807

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

так Вы потратили время не на решение, а на тестирование. А оно везде простое.
test:proc main;
dcl x float(53), s char(*) var;

x=1234.567890;
put data(time);

do to 100_000_000; s=x; end;

put data(time,s);
end test;

TIME= 19:17:44 TIME= 19:18:20 S= 1.23456788999999E+003
процессор i5-2310M 2.5 GHz
популярная реализация была придумана кем то или существовала изначально от сотворения мира?
Популярные реализации по определению существовали до написания данной статьи. Поэтому, предлагая им замену, претендующую на «максимальную эффективность», стоило проверить, действительно ли эффективность максимальная.
Скажите а к себе вы предъявляете те же требования? Вы всегда пишите максимально эффективный код либо не пишите вообще?
Я никаких требования не предъявлял, их предъявил сам автор:
Причем хотелось сделать это максимально эффективно, в полной мере используя аппаратные возможности обработки чисел.
При этом полный отказ от использования команд FPU вряд ли позволит увеличить скорость требуемого преобразования.
Вот это надо было хотя бы попытаться доказать.
Похожую вещь делал давно в институте, даже исходник нашел )) Может кому пригодится — там ввод и вывод числел с плавающей запятой:
Код
.386
model use16 small
.stack 100h

.data
bufstr db 82 dup(?)
dopvar db 12 dup(?)
newstr db 13,10,'$'
rconst_0 dw 2
rconst_1 dd 1000000
rconst_2 dw 10

.code

pstep_2 proc

fabs
fyl2x
fld1
fscale
fld st(1)
fld st
fstcw word ptr dopvar
fwait
or word ptr dopvar,0c00h
fldcw word ptr dopvar
frndint
fsub
fidiv rconst_0
f2xm1
fld1
fadd
fmul st(1),st
fmul
fxch
ffree st
fincstp
ret
pstep_2 endp

write_0 proc
xor si,si
cmp ax,0
je w1c_x
mov bx,ax
and ah,80h
jz w1c
neg bx
mov ah,2
mov dl,'-'
int 21h

w1c: mov cx,10000

w1c0: mov ax,bx
xor dx,dx
div cx
mov bx,dx
cmp ax,0
jne w1c3
cmp si,0
je w1c1
w1c3: or si,1
mov dl,al
mov ah,2
add dl,30h
int 21h

w1c1: xor dx,dx
mov ax,cx
mov cx,10
div cx
mov cx,ax
cmp cx,0
jne w1c0

w1c_x: cmp si,0
jne w1c_y
mov dl,30h
mov ah,2
int 21h

w1c_y: lea dx,newstr
mov ah,9
int 21h
ret
write_0 endp

write_real proc

fxtract ;st-мантисса,st1-показатель
fxch ;меняем местами
fldl2t ;st-log[2](10),st1-p,st2-m
fdivp ;st-dp(десятичный показатель),st1-m
fstcw word ptr dopvar
fwait

or word ptr dopvar,0c00h
fldcw word ptr dopvar

fist word ptr dopvar+10 ;сохр. целую часть
fisub word ptr dopvar+10 ;выделим дробную часть
fabs ;ее модуль
fldl2t ;st-log,st1-modp,st2-m
fmul ;st-binp,st1-m
fld st ;st-binp,st1-binp,st2-m
frndint ;st-mbinp,st1-binp,st2-m
fxch
fsub st,st(1)
fidiv rconst_0 ;st=st/2,st1-mbinp,st2-m
f2xm1 ;st=2^st-1
fld1
fadd ;st=st+1
fld st
fmul ;st=st^2
fscale ;st=st*2^st1,st1-m

fwait

mov al,byte ptr dopvar+11
and al,80h
jz mw1

fdivr st,st(2) ;st-десятичн. мантисса(показат.<0)
fwait
jmp mw2

mw1: fmul st,st(2) ;st-десятичн. мантисса(показат.>0)
mw2: fimul rconst_1 ;st=st*1000000
fbstp tbyte ptr dopvar
ffree st
fincstp
ffree st
fincstp
fwait

mov ah,byte ptr dopvar+9
and ah,80h
jz mw5
mov ah,02h
mov dl,'-'
int 21h

mw5: mov bl,byte ptr dopvar+3
mov cl,bl
and cl,15
mov dl,bl
and dl,240
shr dl,4
cmp dl,0
je mw7
add dl,30h
mov ah,2
int 21h
mw7: mov dl,cl
add dl,30h
mov ah,2
int 21h

mov dl,'.'
mov ah,2
int 21h

mov si,2

mw4: mov bl,byte ptr dopvar+si
mov cl,bl
and cl,15
mov dl,bl
and dl,240
shr dl,4
add dl,30h
mov ah,2
int 21h
mov dl,cl
add dl,30h
mov ah,2
int 21h
cmp si,0
je mw3
dec si
jmp mw4

mw3: mov dl,'e'
mov ah,2
int 21h

mov ax,word ptr dopvar+10
call write_0
ret
write_real endp

read_real proc
r20: readln
xor ah,ah
mov cl,byte ptr bufstr+1
add cl,2
xor ch,ch
mov si,cx
and byte ptr bufstr[si],0
mov si,2
and byte ptr bufstr,0

r23: cmp si,cx
je r20
mov al,byte ptr bufstr[si]
cmp al,'-'
je r22
cmp al,'.'
je r27
cmp al,30h
jb r24
cmp al,40h
jb r21
r24: inc si
jmp r23

r22: or byte ptr bufstr,1
inc si
jmp r23

r27: fild rconst_2
fldz
jmp r28

r21: sub al,30h
mov word ptr bufstr+1,ax
fild rconst_2
fild word ptr bufstr+1


r25: inc si
mov al,byte ptr bufstr[si]
cmp al,'.'
je r28
sub al,30h
cmp al,9
ja r2_x
mov word ptr bufstr+1,ax

fmul st,st(1)
fild word ptr bufstr+1
fadd
jmp r25

r28: inc si
mov al,byte ptr bufstr[si]
sub al,30h
cmp al,9
ja r2_x
mov word ptr bufstr+1,ax
fild word ptr bufstr+1
fdiv st,st(2)
fadd
fld st(1)
fmul st(2),st
ffree st
fincstp
jmp r28

r2_x: and byte ptr bufstr,1
jz r2_y
fchs
r2_y: fxch
ffree st
fincstp
add ax,30h
and al,95
cmp al,'E'
jne r2_z
inc si
mov al,byte ptr bufstr[si]
and byte ptr bufstr,0
cmp al,'-'
je rx0
cmp al,'+'
je rx1
sub al,30h
cmp al,9
ja r2_z
jmp rx2
rx0: or byte ptr bufstr,1
rx1: inc si
mov al,byte ptr bufstr[si]
sub al,30h
cmp al,9
ja r2_z
rx2: inc si
xor bh,bh
mov bl,byte ptr bufstr[si]
sub bl,30h
cmp bl,9
ja rx3
mul rconst_2
add ax,bx
inc si
mov bl,byte ptr bufstr[si]
sub bl,30h
cmp bl,9
ja rx3
mul rconst_2
add ax,bx
rx3: and byte ptr bufstr,1
jz rx4
neg ax
rx4: mov word ptr bufstr+1,ax
fild word ptr bufstr+1
fild rconst_2
call pstep_2
fmul
r2_z: lea dx,newstr
mov ah,9
int 21h
ret
read_real endp

main proc
mov ax,@data
mov ds,ax

call read_real
call write_real

mov ax,4c00h
int 21h

main endp
end main

что вы скажите про этот кусочек кода?

movd qword ptr[esp - xmmword], xmm0
fld qword ptr[esp - xmmword]
fld st(0)
fxtract
fldl2t
fst st(1)
fdivr st(0),st(2)
frndint
; fldcw [esp - word]

fist dword ptr[esp - dword]
mov edx,[esp - dword]
mov dword ptr[esp - dword], 10h
fisubr dword ptr[esp - dword]

fmulp st(1),st(0)
fst st(1)
frndint
fsub st(1),st(0)
fld1
fscale
fstp st(1)
fmulp st(2),st(0)
f2xm1
fld1
faddp st(1),st(0)
fmulp st(1),st(0)
frndint

fbstp tbyte ptr[esp - xmmword * 2]
Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.