Author Topic: Large integers and floats  (Read 391 times)

nidud

  • Member
  • *****
  • Posts: 1386
    • https://github.com/nidud/asmc
Large integers and floats
« on: August 08, 2017, 12:16:04 AM »
There was some problems with the float implementation so I did a small overhaul of the EMU functions. I also added some helper functions for large integers and floats. REAL16 is also added to the assembler.

The new functions are declared in intn.inc and defined in the math directory, and tests added to the test directory.

The real16.asm file is the regress for the Asmc implementation.
Code: [Select]
    ;
    ; v2.25 - REAL16 and float const parameters
    ;
    .486
    .model flat
    .data

    real16 1.e1     ; 0x40024000000000000000000000000000
    real16 1.e2     ; 0x40059000000000000000000000000000
    real16 1.e4     ; 0x400C3880000000000000000000000000
    real16 1.e8     ; 0x40197D78400000000000000000000000
    real16 1.e16    ; 0x40341C37937E08000000000000000000
    real16 1.e32    ; 0x40693B8B5B5056E16B3BE04000000000
    real16 1.e64    ; 0x40D384F03E93FF9F4DAA797ED6E38ED6
    real16 1.e128   ; 0x41A827748F9301D319BF8CDE66D86D62
    real16 1.e256   ; 0x435154FDD7F73BF3BD1BBB77203731FE
    real16 1.e512   ; 0x46A3C633415D4C1D238D98CAB8A978A1
    real16 1.e1024  ; 0x4D4892ECEB0D02EA182ECA1A7A51E317
    real16 1.e2048  ; 0x5A923D1676BB8A7ABBC94E9A519C6536
    real16 1.e4096  ; 0x752588C0A40514412F3592982A7F0095 - F00949524
    real16 1.e8192  ; 0x7FFF0000000000000000000000000000 - Infinity

    REAL16 0.0      ; 0x00000000000000000000000000000000
    REAL16 -0.0     ; 0x80000000000000000000000000000000
    REAL16 1.0      ; 0x3FFF0000000000000000000000000000
    REAL16 -2.0     ; 0xC0000000000000000000000000000000
    ;
    ; 0x7FFEFFFFFFFFFFFFFFFFFFFFFFFFFFFF - max quadruple precision
    ;
    ; 0x7FFB9999999999999999999999999997
    ;
    REAL16 1.189731495357231765085759326628007e4931
    ;
    ; 0x4000921FB54442D18469898CC51701B8
    ;
    REAL16 3.1415926535897932384626433832795028

    .code

p1  proc c a1:real4, a2:real8, a3:real10, a4:real16, a5:oword
    lea eax,a1
    lea eax,a2
    lea eax,a3
    lea eax,a4
    lea eax,a5
    ret
p1  endp

p2  proc stdcall a1:real4, a2:real8, a3:real10, a4:real16, a5:oword
    lea eax,a1
    lea eax,a2
    lea eax,a3
    lea eax,a4
    lea eax,a5
    ret
p2  endp

p3  proc pascal a1:real4, a2:real8, a3:real10, a4:real16, a5:oword
    lea eax,a1
    lea eax,a2
    lea eax,a3
    lea eax,a4
    lea eax,a5
    ret
p3  endp

    p1(1.0, 2.0, 3.0, 4.0, 5.0)
    p2(0.1, 0.2, 0.3, 0.4, 0.5)
    p3(10.1, 10.2, 10.3, 10.4, 10.5)

    end

nidud

  • Member
  • *****
  • Posts: 1386
    • https://github.com/nidud/asmc
Re: Large integers and floats
« Reply #1 on: August 10, 2017, 10:08:41 PM »
Some new changes.

All real numbers are now handled as quadruple precision by the assembler and scaled down accordingly. FPU usage for double (strtod()), and the internal long double implementation for real10 is thereby removed, reducing the binary size by around 10K. This also speed up (somewhat surprisingly) the assembly time for the math test.

The work horse for the ASCII to float conversion is the _atonf() function now able to handle any sized number limited to an INT_MAX exponent. The scaling of numbers is n dwords so the mantissa may be of any length.

Sizing up the maximum digits and significant digits is a bit blurry but it's currently scaled to the mantissa size.
Code: [Select]
; _atonf() - Converts a string to float
;
include intn.inc
include alloc.inc
include ltype.inc

    .code

_atonf proc uses esi edi ebx number:ptr, string:ptr, endptr:ptr, expptr:ptr, bits, expbits, n

    local buffer, flags, exponent, radix
    local digits, sigdig, maxdig, maxsig

    mov radix,10    ; assume desimal
    mov eax,bits    ; 53-bit --> 17
    mov ecx,eax     ; 64-bit --> 22
    mov edx,eax     ;113-bit --> 39
    shr eax,2
    shr edx,4
    shr ecx,5
    add eax,edx
    add eax,ecx
    shr ecx,2
    sub eax,ecx
    mov maxdig,eax
    add eax,10
    mov maxsig,eax
    inc eax
    mov buffer,alloca(eax)
    mov edi,number
    mov ecx,n
    xor eax,eax
    rep stosd
    mov sigdig,eax
    mov exponent,eax
    mov flags,_ST_ISZERO
    mov esi,string

    .repeat

        .repeat
            lodsb
            .break(1) .if !al
        .until !(_ltype[eax+1] & _SPACE)
        dec esi

        xor ecx,ecx
        .if al == '+'
            inc esi
            or  ecx,_ST_SIGN
        .endif
        .if al == '-'
            inc esi
            or  ecx,_ST_SIGN or _ST_NEGNUM
        .endif

        lodsb
        .break .if !al

        or  al,0x20
        .if al == 'n'
            mov ax,[esi]
            or  ax,0x2020
            .if ax == 'na'
                add esi,2
                or ecx,_ST_ISNAN
                movzx eax,byte ptr [esi]
                .if al == '('
                    lea edx,[esi+1]
                    movzx eax,byte ptr [edx]
                    .while al == '_' || _ltype[eax+1] & _DIGIT or _UPPER or _LOWER
                        inc edx
                        mov al,[edx]
                    .endw
                    .if al == ')'
                        lea esi,[edx+1]
                    .endif
                .endif
            .else
                dec esi
                or ecx,_ST_INVALID
            .endif
            mov flags,ecx
            .break
        .endif

        .if al == 'i'
            mov ax,[esi]
            or  ax,0x2020
            .if ax == 'fn'
                add esi,2
                or ecx,_ST_ISINF
            .else
                dec esi
                or ecx,_ST_INVALID
            .endif
            mov flags,ecx
            .break
        .endif

        mov ah,[esi]
        or  ah,0x20
        .if ax == 'x0'
            or ecx,_ST_ISHEX
            add esi,2
            mov edx,bits
            shl edx,2
            mov maxdig,edx
        .endif
        dec esi

        mov edi,buffer
        xor ebx,ebx
        xor edx,edx
        xor eax,eax
        ;
        ; Parse the mantissa
        ;
        .while 1
            lodsb
            .break .if !al
            .if al == '.'
                .break .if ecx & _ST_DOT
                or ecx,_ST_DOT
            .else
                .if ecx & _ST_ISHEX
                    .break .if !(_ltype[eax+1] & _HEX)
                .else
                    .break .if !(_ltype[eax+1] & _DIGIT)
                .endif
                .if ecx & _ST_DOT
                    inc sigdig
                .endif
                or ecx,_ST_DIGITS
                or edx,eax
                .continue .if edx == '0' ; if a significant digit
                .if ebx < maxsig
                    stosb
                .endif
                inc ebx
            .endif
        .endw
        mov byte ptr [edi],0
        mov digits,ebx
        ;
        ; Parse the optional exponent
        ;
        xor edx,edx
        .if ecx & _ST_DIGITS

            xor edi,edi ; exponent

            or  al,0x20
            mov ah,'e'
            .if ecx & _ST_ISHEX
                mov ah,'p'
            .endif

            .if al == ah

                mov al,[esi]
                lea edx,[esi-1]
                .if al == '+'
                    inc esi
                .endif
                .if al == '-'
                    inc esi
                    or  ecx,_ST_NEGEXP
                .endif
                and ecx,not _ST_DIGITS

                .while 1
                    movzx eax,byte ptr [esi]
                    .break .if !(_ltype[eax+1] & _DIGIT)
                    .if edi < 100000000 ; else overflow
                        lea ebx,[edi*8]
                        lea edi,[edi*2+ebx-'0']
                        add edi,eax
                    .endif
                    or  ecx,_ST_DIGITS
                    inc esi
                .endw
                .if ecx & _ST_NEGEXP
                    neg edi
                .endif
                .if !(ecx & _ST_DIGITS)
                    mov esi,edx
                .endif
            .else
                dec esi ; digits found, but no e or E
            .endif

            mov edx,edi
            mov eax,sigdig
            .if ecx & _ST_ISHEX
                shl eax,2
            .endif
            sub edx,eax
            mov ebx,digits
            mov eax,maxdig
            .if ebx > eax
                add edx,ebx
                mov ebx,eax
                .if ecx & _ST_ISHEX
                    shl eax,2
                .endif
                sub edx,eax
            .endif
            mov eax,buffer
            .while 1
                .break .ifs ebx <= 0
                .break .if byte ptr [eax+ebx-1] != '0'
                inc edx
                .if ecx & _ST_ISHEX
                    add edx,3
                .endif
                dec ebx
            .endw
            mov digits,ebx
        .else
            mov esi,string
        .endif
        .if ecx & _ST_ISHEX
            mov radix,16
        .endif
        mov flags,ecx
        mov exponent,edx
        mov eax,digits
        add eax,buffer
        mov byte ptr [eax],0
        ;
        ; convert string to binary
        ;
        _atond(number, buffer, radix, n)
        ;
        ; get bit-count of number
        ;
        .if _bsrnd(number, n)
            ;
            ; shift number to bit-size
            ;
            ; - 0x0001 --> 0x8000
            ;
            mov edi,bits
            mov ebx,edi
            sub ebx,eax

            .if eax > edi
                sub eax,edi
                ;
                ; or 0x10000 --> 0x1000
                ;
                _shrnd(number, eax, n)
            .elseif ebx
                _shlnd(number, ebx, n)
            .endif
            ;
            ; create exponent bias and mask
            ;
            mov ecx,expbits
            mov eax,1
            shl eax,cl
            mov ecx,eax     ; sign bit
            dec eax         ; mask
            mov edx,eax
            shr edx,1       ; bias
            add edi,edx     ; bits + bias
            sub edi,ebx     ; - shift count
            dec edi
            and edi,eax     ; remove sign bit
            .if flags & _ST_NEGNUM
                or edi,ecx
            .endif
            mov ebx,ecx
            xor edx,edx     ; get dword-id from bits
            mov eax,bits
            mov ecx,32
            div ecx
            mov ecx,ebx
            shl ecx,1
            dec ecx
            .if edx
                ;
                ; then shift exponent and sign to high word
                ;
                mov edx,ecx
                mov ecx,31
                bsf ebx,ebx
                sub ecx,ebx
                shl edi,cl
                shl edx,cl
                mov ecx,edx
            .endif

            lea ebx,[eax*4]
            add ebx,number
            not ecx
            and [ebx],ecx
            or  [ebx],edi
        .else
            or flags,_ST_ISZERO
        .endif

    .until 1

    mov eax,expptr
    .if eax
        mov ecx,exponent
        mov [eax],ecx
    .endif
    mov eax,endptr
    .if eax
        mov [eax],esi
    .endif
    mov eax,flags
    ret

_atonf endp

    end

nidud

  • Member
  • *****
  • Posts: 1386
    • https://github.com/nidud/asmc
Re: Large integers and floats
« Reply #2 on: August 22, 2017, 10:11:54 AM »
Some more __float128 stuff.

The mantissa in a quad is 113 bit where only 112 is stored. To apply math this has to be expanded to 128, in RDX:RAX in this case, where the high bit and 15 low bits are shifted out and the rest stored in the end.

Code: [Select]
    mov rax,[r8] ; low 48-bit
    shl rax,16 ;
    mov rdx,[r8+6] ; high 64-bit
    mov cx,[r8+14]
    and cx,0x7FFF ; bit 112 set if not zero
    neg cx
    rcr rdx,1 ; 113 bits..
    rcr rax,1

Here's a test version for adding and subtracting two REAL16 numbers. This is similar to the REAL10 version using a 64-bit mantissa in EDX:EAX.

Code: [Select]
include math.inc

option win64:rsp nosave noauto

.code

_subfq proc result:ptr, a:ptr, b:ptr

    mov r9d,0x80000000
    jmp _lk_add

_subfq endp

_addfq proc result:ptr, a:ptr, b:ptr

    xor r9d,r9d

_addfq endp

_lk_add proc private uses rsi rdi rbx result:ptr, a:ptr, b:ptr, sign:dword
    ;
    ; quad float [rcx] = quad float [rdx] +- quad float [r8]
    ;
    mov r10,rcx ; save dest.

    mov rbx,[r8]
    shl rbx,16
    mov rdi,[r8+6]
    mov si,[r8+14]
    and si,0x7FFF
    neg si
    mov si,[r8+14]
    rcr rdi,1
    rcr rbx,1
    shl esi,16
    mov r8,rdx
    mov rdx,[r8+6]
    mov rax,[r8]
    shl rax,16
    mov si,[r8+14]
    and si,0x7FFF
    neg si
    mov si,[r8+14]
    rcr rdx,1
    rcr rax,1

    .repeat             ; create frame -- no loop

        add si,1            ; add 1 to exponent
        jc  er_NaN_A        ; quit if NaN
        jo  er_NaN_A        ; ...
        add esi,0xFFFF      ; readjust low exponent and inc high word
        jc  er_NaN_B        ; quit if NaN
        jo  er_NaN_B        ; ...
        sub esi,0x10000     ; readjust high exponent
        xor esi,r9d         ; flip sign if subtract

        .if !rax && !rdx    ; A is 0 ?
            shl si,1        ; place sign in carry
            .ifz
                shr esi,16
                mov rax,rbx     ; return B
                mov rdx,rdi
                shl esi,1
                or  rbx,rdi     ; check for 0
                or  bx,si
                .ifnz           ; if not zero
                    shr esi,1   ; -> restore sign bit
                .endif
                .break
            .endif
            rcr si,1        ; put back the sign
        .endif
        .if !rbx && !rdi    ; B is 0 ?
            .break .if !(esi & 0x7FFF0000)
        .endif

        mov ecx,esi         ; exponent and sign of A into ECX
        rol esi,16          ; shift to top
        sar esi,16          ; duplicate sign
        sar ecx,16          ; ...
        and esi,0x80007FFF  ; isolate signs and exponent
        and ecx,0x80007FFF  ; ...
        mov r9d,ecx         ; assume A < B
        rol esi,16          ; rotate signs to bottom
        rol ecx,16          ; ...
        add cx,si           ; calc sign of result
        rol esi,16          ; rotate signs to top
        rol ecx,16          ; ...
        sub cx,si           ; calculate difference in exponents
        .ifnz               ; if different
            .ifc                ; if B < A
                mov r9d,esi     ; get larger exponent for result
                neg cx          ; negate the shift count
                xchg rbx,rax    ; flip operands
                xchg rdi,rdx
            .endif
            .if cx > 128        ; if shift count too big
                mov esi,r9d
                shl esi,1       ; get sign
                rcr si,1        ; merge with exponent
                mov rax,rbx     ; get result
                mov rdx,rdi
                .break
            .endif
        .endif

        mov esi,r9d
        mov ch,0            ; zero extend B
        or  ecx,ecx         ; get bit 0 of sign word - value is 0 if
                            ; both operands have same sign, 1 if not
        .ifs                ; if signs are different
            mov ch,-1       ; - set high part to ones
            neg rdi         ; - negate the fraction of B
            neg rbx
            sbb rdi,0
            xor esi,0x80000000 ; - flip sign
        .endif

        xor r11,r11         ; get a zero for sticky bits
        .if cl              ; if shifting required
            .if cl >= 64        ; if shift count >= 64
                .if rax         ; check low order qword for 1 bits
                    inc r11     ; r11=1 if RAX non zero
                .endif
                .if cl == 128   ; if shift count is 128
                    or  r11,rdx ; get rest of sticky bits from high part
                    xor rdx,rdx ; zero high part
                .endif
                mov rax,rdx     ; shift right 64
                xor rdx,rdx
            .endif
            xor  r8,r8
            shrd r8,rax,cl      ; get the extra sticky bits
            or   r11,r8         ; save them
            shrd rax,rdx,cl     ; align the fractions
            shr  rdx,cl
        .endif

        add rax,rbx
        adc rdx,rdi
        adc ch,0
        .ifs
            .if cl == 128
                xor r8b,r8b
                mov r9,0x7FFFFFFFFFFFFFFF
                .if r11 & r9
                    inc r8b     ; make single sticky bit
                .endif
                shr r8b,1
                adc rax,0       ; round up fraction if required
                adc rdx,0
            .endif
            neg rdx
            neg rax
            sbb rdx,0
            mov ch,0
            xor esi,0x80000000
        .endif

        .if rax || rdx || ch
            .if !si
                add esi,esi
                rcr si,1
                .break
            .endif
            ;
            ; if top bits are 0
            ;
            .if !ch
                rol r11,1   ; set carry from last sticky bit
                rol r11,1
                .repeat
                    dec si  ; decrement exponent
                    .ifz
                        add esi,esi
                        rcr si,1
                        .break(1)
                    .endif
                    adc rax,rax ; shift fraction left one bit
                    adc rdx,rdx
                .untilb         ; until carry
            .endif
            inc si
            cmp si,0x7FFF
            je  overflow
            stc                 ; set carry
            rcr rdx,1           ; shift fraction right 1
            rcr rax,1
            bt  eax,15
            .ifc                ; if guard bit is on
                shl r11,1       ; get top sticky bit
                .ifz            ; if no more sticky bits
                    ror rax,1   ; set carry with bottom bit of DX
                    rol rax,1
                .endif
                adc rax,0       ; round up fraction if required
                adc rdx,0
                .ifc            ; if we got a carry
                    rcr rdx,1   ; shift fraction right 1
                    rcr rax,1
                    inc si      ; increment exponent
                    cmp si,0x7FFF
                    je  overflow
                .endif
            .endif
        .else
            xor esi,esi
        .endif
        add esi,esi
        rcr si,1
        .break

      ; A is a NaN or infinity

      er_NaN_A:

        dec si
        add esi,0x10000
        .ifnc
            .break .ifno
        .endif
        sub esi,0x10000
        mov r11,0x8000000000000000
        .if !rax && !rbx
            .if rdx == rdi && rdx == r11
                shr rdx,2
                or  si,-1 ; -NaN - FFFF40 ?
                .break
            .endif
        .endif
        .if rdi == rdx
            cmp rbx,rax
        .endif
        jna return_B
        .break

      ; B is a NaN or infinity

      er_NaN_B:

        sub esi,0x10000
        .if !rbx
            mov rdx,0x8000000000000000
            .if rdx == rdi
                xor esi,r9d
            .endif
        .endif
      return_B:
        mov rdx,rdi
        mov rax,rbx
        shr esi,16
        .break

      overflow:
        mov si,0x7FFF
        xor rax,rax
        xor rdx,rdx

    .until 1

    shl rax,1           ; shift bits back
    rcl rdx,1
    shr rax,16          ; 16 low bits
    mov [r10],rax
    mov [r10+6],rdx
    mov [r10+14],si     ; exponent and sign
    mov rax,r10         ; return result
    ret

_lk_add endp

    end

nidud

  • Member
  • *****
  • Posts: 1386
    • https://github.com/nidud/asmc
Re: Large integers and floats
« Reply #3 on: October 22, 2017, 01:29:26 AM »
Added formatted output for quad float to the STDIO functions.
'L' is assigned to REAL10 and 'll' to REAL16.

Code: [Select]
include stdio.inc

.data

n equ <3.141592653589793238462643383279502884197169399375105820974945>

d REAL8  n
l REAL10 n
q REAL16 n

.code

main proc

    printf("%g\n",     d)
    printf("%Lg\n",    l)
    printf("%llg\n\n", q)

    printf("%f\n",     d)
    printf("%Lf\n",    l)
    printf("%llf\n\n", q)

    printf("%e\n",     d)
    printf("%Le\n",    l)
    printf("%lle\n\n", q)

    printf("%.34f\n",  d)
    printf("%.34Lf\n", l)
    printf("%.34llf\n",q)
    ret

main endp

    end

Code: [Select]
3.14159
3.14159
3.14159

3.141593
3.141593
3.141593

3.141593e+000
3.141593e+000
3.141593e+000

3.1415926535897931160000000000000000
3.1415926535897932385128000000000000
3.1415926535897932384626433832795020

nidud

  • Member
  • *****
  • Posts: 1386
    • https://github.com/nidud/asmc
Re: Large integers and floats
« Reply #4 on: November 19, 2017, 04:07:02 AM »
Quad float is now also added to the 64-bit STDIO functions.

Note that REAL8 are past in register and REAL10/REAL16 are past as address (in the same reg) regardless of attributes (&/addr) in format arguments.

test case:
Code: [Select]
- format %g, %Lg, %llg
3.14159
3.14159
3.14159
- format %f, %Lf, %llf
3.145933
3.145933
3.145933
- format %e, %Le, %lle
3.145933e+000
3.145933e+000
3.145933e+000
- format %.34f, %.34Lf, %.34llf
3.1459226535897895633000000000000000
3.1459226535897932385128000000000000
3.1459226535897932384626433832795020