Комментарии 23
А нельзя было по порядку определить примерное число десяток, на которые следует умножить/разделить число (порядок — целая часть логарифма, а значит, зависимость будет линейная)? Это, по идее, позволило бы уменьшить число циклов (за счет применения возведения 10 в нужную степень: цикл останется, но меньший по размеру). А то в double десятичный порядок может достигать 308, за вычетом (для отрицательных порядков — с прибавлением) 15-16 десятичных цифр, получается довольно большое количество итераций.
Т.е., условно, мы сперва "сдвигаем" число так, чтобы оно получалось примерно целым, затем используем описанный алгоритм.
Тестировать алгоритм, имеющий цикл, зависящий от порядка, т.е., по сути, сложность O(порядок), только на одном на конкретном непредельном значении с небольшим порядком, не совсем правильно.
округление в FBSTP управляется флагом CW.RC
Откуда столько знаков? Под мантиссу отводится 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 десятичных цифр, то это проявляется при чтении в четверную точность или выше, что крайне редко на практике.
Ну да. Вычислить десятичный логарифм одной 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 у автора, но потратил я на это "решение" времени меньше чем писал этот комментарий.
Причем хотелось сделать это максимально эффективно, в полной мере используя аппаратные возможности обработки чисел.
При этом полный отказ от использования команд 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
Перевод числа в строку с помощью FPU