flat assembler
Message board for the users of flat assembler.

Index > Main > 10^256bit .fraction

Thread Post new topic Reply to topic

Joined: 18 Jul 2009
Posts: 549
edemko 24 Aug 2010, 14:50
How can we get 10^.256 bit fraction?
Have a look at a <f2xey> procedure.
I put 2 HLT instros to easily debug.
The reasons i'm comming so complicated way described in same proc too.
To easily manipulate with numbers, use OllyDbg's built-in convertors or this tool i've long-long time working on:

format pe gui 4.0
include 'win32ax.inc'

section '' code executable import readable writeable
  library kernel32,'kernel32.dll',\
  include 'api\kernel32.inc'
  include 'api\user32.inc'

;tf f1,$3fff-0000,$8000000000000000 ;2^0*1
;tf f2,$bfff-0064,$8000000000000000 ;2^-64*1
macro tf name*, exponent*, significand*{
  label name tbyte at $
  dq significand
  dw exponent

;include '../macro.inc'
;include '../proc.inc'

;buf80bytes db 80 dup 0
f dt -65535.0
;tf f,$403e,$ffffffffffffffff

entry $
        mov     eax,f
        call    f2xey

  Prepares floating value for text formating.  This is an intermediate procedure.
  eax         -> tbyte float location
  ecx:ebx:eax <- decimal mantissa given as a float, absolute value
  edx         <- bit 31..31: sign of decimal mantissa, 1='-', 0='+'
                 bit 30..30: sign of decimal exponent, 1='-', 0='+'
                 bit 00..15: decimal exponent given as an integer, absolute value
  flags.cf    <- 0 if not zero
                 1 if zero, ecx.31 equals zero's sign, ebx:eax=0, edx=?
  flags       <- ?
proc f2xey
        ;decom`pose number
        stc                          ;ieee stupido mode ON
        call    fxtract_
        jnc     .not_zero
        ret     0

        ;10^x = 2^unbisased * normalized      =
        ;     = 2^unb * 2^log(2;norm)         =
        ;     = 2^(unb + log(2;norm))        ->
        ;   x = log(10;2^(unb + log(2;norm))) =
        ;     = log(10;2)*(unb + log(2;norm))
        ;log(10;2)       = 3FFD 9A209A84 FBCFF799
        ;                =     .4D104D42'7DE7FBCC'8xxxxxxx
        ;unb             = -16'445..+16'384
        ;log(2;norm)     = [0..1)
        ;min log(2;norm) = log(2;$3fff'80000000'00000001)
        ;                = 3FC0 B8AA3B29 5C17F0BB ->
        ;               -> ie $3fff-$3fc0=$3f=63
        ;Performing "unb + log(2;norm)", fpu loses a lot of fractional bits.
        ;That's why we are to use 128 bits for fraction.
        ;Mention: for .fraction only, hence 3FFD 9A209A84 FBCFF799 =
        ;=9A209A84 FBCFF799 shr 1 but not 2
        push    esi edi ebp

        ;get log(2;normalized) into ST0
        push    word $3fff ebx eax
        fld     tbyte[esp]
        fyl2x                        ;got = [0..1)

        ;unnormalize ST0 into BP:.EBX:EAX:ESI:EDI as common whole.fraction
        fstp    tbyte[esp]
        mov     ebp,ecx              ;keep value sign(ebp.31) and unb(ebp.00..15)
        pop     eax ebx cx           ;operate with registers
        neg     cx                   ;getting shift count to get actual value
        jnc     .idle                ;log(2;norm)=0
        add     cx,$3fff-1           ;-1..-62
.idle:  xor     esi,esi              ;complete fraction
        xor     edi,edi              ;complete fraction
        cmp     cl,32
        jb      .b32                 ;small shift
        xchg    eax,esi
        xchg    ebx,eax
        sub     cl,32                ;denormalized partially
.b32:   shrd    edi,esi,cl
        shrd    esi,eax,cl
        shrd    eax,ebx,cl
        shr     ebx,cl               ;denormalized fully - actual log(2; normalized) got

        ;add unbiased in BP to log(2;normalized) in .EBX:EAX:ESI:EDI
        test    bp,bp
        jns     .unbiased_added      ;whole.fraction :)^
        neg     bp                   ;how would you do with -1+0.3? -> -(1-0.3)=-0.7
        bts     ebp,30               ;as floats stored positive, they need a bit saying "-" or "+"
        neg     ebx                  ;subtract
        sbb     ebp,0                ;and chain; we could also use "sbb bp,0", nah
        neg     eax
        sbb     ebx,0
        sbb     ebp,0
        neg     esi
        sbb     eax,0
        sbb     ebx,0
        sbb     ebp,0
        neg     edi
        sbb     esi,0
        sbb     eax,0
        sbb     ebx,0
        sbb     ebp,0                ;negative whole.fraction got

        ;multiply BP:.EBX:EAX:ESI:EDI by log(10;2)
        push    ebx eax esi edi
        mov     esi,esp              ;<_128mul128>.m2
        mov     ecx,$4D10'4D42
        mov     ebx,$7DE7'FBCC
        mov     eax,$8000'0000
        push    ecx ebx eax 0
        mov     edi,esp              ;<_128mul128>.m1
        stdcall _128mul128,edi,esi
        movzx   esi,bp               ;we've multiplied .fractions, now whole.*log(10;2)
        xor     bp,bp                ;clean whole. part
        mul     esi
        add     [edi+5*4],eax
        adc     [edi+6*4],edx
        adc     dword[edi+7*4],0
        adc     ebp,0
        mov     eax,ebx
        mul     esi
        add     [edi+6*4],eax
        adc     [edi+7*4],edx
        adc     ebp,0
        mov     eax,ecx
        mul     esi
        add     [edi+7*4],eax
        adc     ebp,edx              ;got whole.fraction

        ;getting 10^.fraction = (10^-1..10^+1) = (0.1..10)
        std                          ;searching for the highest non-zero dword...
        add     edi,7*4              ;assume we've found it
        xor     eax,eax
        mov     ecx,8                ;256bits=8dwords
        repe    scasd
        jne     .fraction_nz         ;found, else .fraction zero
        mov     ecx,$3fff            ;exponent=0 ie 2^0
        mov     ebx,$8000'0000       ;normalized 1.0
        xor     eax,eax              ;too
        mov     edx,ebp              ;edx.31=input sign, edx.30=10's integral power sign, edx.00.15=10's integral power
        add     esp,8*4              ;get rid of 256bit .fraction
        jmp     .tired_edemko        ;bye

        mov     ch,cl
        xor     esi,esi
        mov     ebx,[edi+4]
        test    cl,cl
        jz      .normalize_fraction
        mov     eax,[edi]
        dec     cl
        jz      .normalize_fraction
        mov     esi,[edi-4]
        neg     ch
        add     ch,7
        shl     ch,5
        bsr     edx,ebx
        neg     edx
        add     edx,31
        add     ch,dl
        mov     cl,dl
        shld    ebx,eax,cl
        shld    eax,esi,cl
        movzx   ecx,ch
        neg     ecx
        add     ecx,$3fff-1
        shl     ecx,1
        bt      ebp,30
        rcr     cx,1
        push    ecx ebx eax
        fld     tbyte[esp]
        call    power10

        pop     ebp edi esi
        ret     0

;  address of tbyte float -> eax
;  flags.cf               -> ignored if qword[eax]=0
;                            else it defines how to process a $0000'xxxxxxxx'xxxxxxxx case
;                            if cf=0: no stupido IEEE processing will be done
;                                     this is not recommended if no exotic wanted
;                                     $0000'8xxxxxxx'xxxxxxxx ->fxtract-> cx = -16'383
;                                     $0000'4xxxxxxx'xxxxxxxx ->fxtract-> cx = -16'384
;                            if cf=1: increase exponent by 1 like fpu's FXTRACT does
;                                     $0000'8xxxxxxx'xxxxxxxx ->fxtract-> cx = -16'382
;                                     $0000'4xxxxxxx'xxxxxxxx ->fxtract-> cx = -16'383
;  ebx:eax                <- normalized mantissa
;  ecx                    <- bit 31   : sign of mantissa
;                            bit 0..15: unbiased integer exponent within -16'446(MINUS stupido.ieee.flags.cf)..+16'384
;                            bit rest : 0
;  flags.cf               <- 1 if mantissa was zero
;  ------------------------------------------------------
;  possible usage:
;    2^unbiased * normalized = 10^x
;    x = log(10;2^unbiased * normalized)
;      = log(10;2^unbiased) + log(10;normalized)
;      = log(10;2)*unbiased + log(10;2)*log(2;normalized)
;      = log(10;2)*(unbiased + log(2;normalized))
;      = Frank Airline in my dreams, thank you
;  possible usage:
;    2^unbiased = 10^x
;    x=log(10;2^unbiased)=log(10;2)*unbiased
;  possible usage:
;    10^x = 2^(unbiased + log(2;normalized))
;    x    = log(10;2) * (unbiased + log(2;normalized))
;  possible usage:
;    boundless division and multiplication
proc fxtract_
        push    edx
        pushf                    ;save ieee.stupido
        movzx   ecx,word[eax+8]
        shl     cx,1             ;drop sign bit
        rcr     ecx,1            ;into the top
        sub     cx,$3fff         ;to easily get unbiased exponent
        mov     ebx,[eax+4]      ;mantissa
        mov     eax,[eax]        ;mantissa
        mov     edx,eax
        or      edx,ebx
        jnz     .nz              ;else zero found
        or      byte[esp],1      ;flags.cf=1
.done:  popf
        pop     edx
        ret     0                ;fasmme
.nz:    btr     dword[esp],0     ;flags.cf=0, get ieee.stupido
        jnc     .no_ieee         ;no ieee.stupido
        cmp     cx,$c001         ;have we met a stupido condition(-16'383)?
        jne     .no_ieee         ;no
        inc     ecx              ;do stupido
        test    ebx,ebx
        js      .done            ;there is nothing to normalize in mantissa, return -16'383
        test    ebx,ebx
        jnz     .ebx_nz          ;end normalization
        xchg    ebx,eax
        sub     cx,32            ;partial normalization done
.ebx_nz:bsr     edx,ebx          ;get times to normalize
        neg     edx
        add     edx,31
        sub     cx,dx            ;_balance_
        xchg    cl,dl
        shld    ebx,eax,cl
        shl     eax,cl
        mov     cl,dl
        jmp     .done            ;2010_08_21

;  [m2:m1] = [m2]*[m1].
;  flags   <-?
;  low 128 bits stored at <m1> then high 128 bits stored at <m2>
;  there will be fun when buffers overlap :)
proc _128mul128; m1,m2
        push    ebp
        mov     ebp,esp
        push    edi
        mov     edi,[ebp+2*4]           ;m1
        push    esi
        mov     esi,[ebp+3*4]           ;m2
        push    edx ecx ebx eax         ;we'll use stepping _64mul64
                                        ;[m1]*[m2] =  [m1.1]*[m2.1]          + ;00I
                                        ;          + ([m1.1]*[m2.2]) shl 064 + ;0II
                                        ;          + ([m1.2]*[m2.1]) shl 064 + ;III
                                        ;          + ([m1.2]*[m2.2]) shl 128   ;0IV
        mov     edx,[esi+4]
        mov     ecx,[esi]
        mov     ebx,[edi+4]
        mov     eax,[edi]
        call    _64mul64
        push    0 0 0 0 edx ecx ebx eax ;00I
        mov     edx,[esi+12]
        mov     ecx,[esi+8]
        mov     ebx,[edi+4]
        mov     eax,[edi]
        call    _64mul64
        add     [ebp-12*4],eax
        adc     [ebp-11*4],ebx
        adc     [ebp-10*4],ecx
        adc     [ebp-9*4],edx
        adc     dword[ebp-8*4],0        ;0II
        mov     edx,[esi+4]
        mov     ecx,[esi]
        mov     ebx,[edi+12]
        mov     eax,[edi+8]
        call    _64mul64
        add     [ebp-12*4],eax
        adc     [ebp-11*4],ebx
        adc     [ebp-10*4],ecx
        adc     [ebp-9*4],edx
        adc     dword[ebp-8*4],0        ;III
        mov     edx,[esi+12]
        mov     ecx,[esi+8]
        mov     ebx,[edi+12]
        mov     eax,[edi+8]
        call    _64mul64
        add     [ebp-10*4],eax
        adc     [ebp-9*4],ebx
        adc     [ebp-8*4],ecx
        adc     [ebp-7*4],edx           ;0IV
        pop     dword[edi] dword[edi+4] dword[edi+8] dword[edi+12] dword[esi] dword[esi+4] dword[esi+8] dword[esi+12]

        pop     eax ebx ecx edx
        pop     esi edi
        ret     8                       ;2010_08_20


;you may skip it though it's interesting
; $ffffffff(edx)ffffffff(ecx) * $ffffffff(ebx)ffffffff(eax)      =
;=$ffffffff'ffffffff * $1'00000000'00000000 - $ffffffff'ffffffff =
;=$ffffffff'ffffffff'00000000'00000000      - $ffffffff'ffffffff =
;=$ffffffff'ffffffff'00000000'00000000  - 1 - $ffffffff'fffffffe =
;=$ffffffff'fffffffe'ffffffff'ffffffff      - $ffffffff'fffffffe =
;=$ffffffff'fffffffe'00000000'00000001, 128 bits :b

;edx:ecx:ebx:eax = edx:ecx*ebx:eax = ecx*ebx:eax+(edx*ebx:eax)shl 32
;flags           = ?
proc _64mul64
        push    eax ebx edx
        call    _64mul32
        xchg    eax,[esp+8]
        xchg    ebx,[esp+4]
        xchg    ecx,[esp]
        call    _64mul32
        mov     edx,ecx
        mov     ecx,ebx
        mov     ebx,eax
        pop     eax
        add     ecx,eax
        adc     edx,0
        pop     eax
        add     ebx,eax
        adc     ecx,0
        adc     edx,0
        pop     eax
        ret     0

;ecx:ebx:eax = ecx*ebx:eax = ecx*ebx:0+ecx*0:eax
;edx         = ?, ecx*0:eax shr 32 actually
;flags       = ?
proc _64mul32
        xchg    eax,ebx
        mul     ecx
        xchg    ebx,eax
        xchg    ecx,edx
        mul     edx
        add     ebx,edx
        adc     ecx,0
        ret     0

;power   -> st0
;st0     <- 10^power
;fpu.tos <- returned
;ensure st5..st7 are free before calling this proc
proc power10
        fldl2t         ;log(2;10)
        fmulp          ;x

                       ;2^x=2^x.fraction*2^x.whole=2^x.fraction shl x.whole
        fld     st0    ;x
        frndint        ;x.whole
        fsub    st1,st0;x.fraction
        f2xm1          ;2^x.fraction -1
        faddp          ;2^x.fraction
        fscale         ;shl x.whole
        fstp    st1    ;10^x
        ret     0      ;-)
Post 24 Aug 2010, 14:50
View user's profile Send private message Reply with quote

Joined: 21 Jul 2003
Posts: 3911
Location: vpcmipstrm
bitRAKE 25 Aug 2010, 02:21
Check this out:


fld [t0_256]
f2xm1 ; this is okay because -1 < 0.256 * log2(10) < 1    
...works on my system.

Edit: actually, http://fornax.sk/mate/Download/dosprog/tasm/show87.asm has the display routine to convert FP to ASCII.
Post 25 Aug 2010, 02:21
View user's profile Send private message Visit poster's website Reply with quote

Joined: 18 Jul 2009
Posts: 549
edemko 25 Aug 2010, 08:56
Their formula seems unseen for me :)
To be more precise i got a .fraction which occupies 256 bits = 8 dwords of memory <> 0.256 :)
I've not seen a tool which could decode $403e'ffffffff'ffffffff into 18446744073709551615 decimal exactly.
<f2xey> has 2 HLTs inserted to stop OllyDbg i use, and see 10^.fraction_64_bits_size.
Of course, there is an error as rest of fraction stays lazy.
See pictures i like to speak in so:
Tomasz seems not to use floats often but i do, here is a sample procedure. Add smth like to fasm please.
;  Recognize hex-string into a tbyte float.
;  src         = string
;  options     = bits 00..007: dummy char not causing errors
;  flags.cf    <-error
;  ecx:ebx:eax <- result
proc bsa_h2f; src, options
        pushf                    ;error=1
        push    esi
        mov     esi,[esp+3*4]    ;src
        test    esi,esi
        jz      .done            ;src=0
        movzx   ebx,byte[esi-4]
        lea     esi,[esi+ebx-1]  ;switch LeastSignificantByte
        shl     ebx,8
        jz      .done            ;src.body_len=0 | src.body_len>255
        mov     bl,[esp+4*4]     ;options(dummy)
        xor     ecx,ecx
        push    ecx ecx ecx      ;result
.load:  lodsb
        cmp     al,bl
        je      .continue        ;dummy
        call    ASCII
        sub     al,'0'
        jz      .register        ;0*smth=0
        cmp     ecx,80-4
        ja      .pop             ;over(80bits)flow
        cmp     al,10
        jb      .input           ;00..10
        cmp     al,'F'-'0'
        ja      .pop             ;not hex
        cmp     al,'A'-'0'
        jb      .pop             ;not hex
        sub     al,7             ;10..15
.input: movzx   eax,al
        cmp     cl,64
        jb      .32?             ;we are in 64..79
        sub     cl,64
        shl     eax,cl
        add     cl,64
        or      [esp+2*4],eax    ;update result.64..79
        jmp     .register
.32?:   cmp     cl,32
        jb      .32              ;we are in 32..63
        sub     cl,32
        shl     eax,cl
        add     cl,32
        or      [esp+1*4],eax    ;update result.32..63, .., i've got used counting from 1
        jmp     .register
.32:    shl     eax,cl           ;we are in 0..31
        or      [esp],eax
        add     ecx,4
        dec     bh
        jnz     .load            ;till the end of the sun
        test    ecx,ecx
        jz      .pop             ;dummy string
        dec     byte[esp+(3+1)*4];drop flags.cf
.pop:   pop     eax ebx ecx
.done:  pop     esi
        ret     2*4              ;2010_08_09

Finally Tomasz you form dword, qword or tbyte due the directive used.
Post 25 Aug 2010, 08:56
View user's profile Send private message Reply with quote

Joined: 21 Jul 2003
Posts: 3911
Location: vpcmipstrm
bitRAKE 25 Aug 2010, 11:46
I'd use integer to convert 256 bit number. For a fraction, multiply by largest power of 10 which fits into your integer size (word, dword, qword). For DWORD we are talking about 10^9. Then you have first nine digits of fraction. iterate to get more digits. This works for fraction of any bit size.
Post 25 Aug 2010, 11:46
View user's profile Send private message Visit poster's website Reply with quote

Joined: 18 Jul 2009
Posts: 549
edemko 25 Aug 2010, 13:27
   10^.fraction_256                                =
 = wait, let us limit it to .fraction128           =
 = 10^.fraction[127..064] * 10^.fraction[063..000] =
 = happy                                           =
 = multiply by 10^19 or 10^20 to get 19 digits of fraction

 10^19 and 10^20, suppose those are 10^1 and 10^2 ->
 1.0 * 10^1 = 10
 0.1 * 10^2 = 10
 =thank you
Post 25 Aug 2010, 13:27
View user's profile Send private message Reply with quote

Joined: 18 Jul 2009
Posts: 549
edemko 29 Aug 2010, 07:21
edit #3: wrong approach, found correct formula, anyway addition used, this way loss of precision occurs, eh, stopping

edit №2:
;getting power of ten we use a formula: 10^x=2^y, y=log(2;10)*x
;finally when fpu executes F2XM1 we must add 1 where rounding occurs
;hence a new formula, assuming that .fraction takes 128 bits:
;      1_F2XM1*2^1_FSCALE * 2_F2XM1*2^2_FSCALE         =
;->    removing M1 factor for very member              ->
;->    1_F2XM1*2^(1_FSCALE+1) * 2_F2XM1*2^(2_FSCALE+1) =
;=     1_F2XM1*2_F2XM1*2^(1_FSCALE+2_FSCALE+2)         eg
;eg    (2-1)*2^0 * (2-1)*2^0 =
;=     1*2^0     * 1*2^0     ->
;->    1*1*2^(0+0+2)         =
;=     4                     check
;check 2*2^0 * 2*2^0         =
;=     4                     =
;=     :)
Post 29 Aug 2010, 07:21
View user's profile Send private message Reply with quote
Display posts from previous:
Post new topic Reply to topic

Jump to:  

< Last Thread | Next Thread >
Forum Rules:
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
You cannot attach files in this forum
You can download files in this forum

Copyright © 1999-2024, Tomasz Grysztar. Also on GitHub, YouTube.

Website powered by rwasa.