BCD floating point math

Source-code only - please, don't post questions here.
srvaldez
Posts: 1371
Joined: Sep 25, 2005 21:54

BCD floating point math

Postby srvaldez » Jul 09, 2007 13:02

I translated this routines just for fun, use at your own risk.

Code: Select all

'
' BCD floating-point math package
' Some Assembly Required - BYTE Magazine
' April 1989
' -- Rick Grehan
'
' adapted to FreeBasic by srvaldez

#define NUM_BYTES        32                   'Number of bytes per mantissa
#define NUM_WORDS        NUM_BYTES/2
#define NUM_DWORDS       NUM_BYTES/4
#define NUM_DIGITS       NUM_BYTES*2
#define BIAS             16384

' Error definitions

#define DIVZ_ERR        1                    'Divide by zero
#define EXPO_ERR        2                    'Exponent overflow error
#define EXPU_ERR        3                    'Exponent underflow error

type bcdfloat
    As Ushort sign
    As Ushort exponent
    As Ubyte  mantissa(0 To NUM_BYTES-1)
 
    Declare constructor ( Byref rhs As bcdfloat )
    Declare constructor ( Byref rhs As String ="0" )
    Declare constructor ( Byref rhs As Double )
    Declare constructor ( Byref rhs As Longint )

    Declare operator let( byref rhs As bcdfloat )
    Declare operator let( byref rhs As String )
    Declare operator let( byval rhs As Double )
    Declare operator let( byval rhs As Longint )
    '' implicit step versions
    Declare Operator For ( )
    Declare Operator Step( )
    Declare Operator Next( Byref end_cond As bcdfloat ) As Integer
   
    '' explicit step versions
    Declare Operator For ( Byref step_var As bcdfloat )
    Declare Operator Step( Byref step_var As bcdfloat )
    Declare Operator Next( Byref end_cond As bcdfloat, Byref step_var As bcdfloat ) As Integer
       
    Declare Operator += (byref rhs as bcdfloat)
    Declare Operator -= (byref rhs as bcdfloat)

    Declare operator cast( ) As String
    Declare operator cast( ) As Double
End type

Function str2bcd(Byval value As String , Byref n As bcdfloat) As Integer
  Dim As Integer j,s,d,e,ep,ex,es,i,f,fp,fln
  Dim As String c,f1,f2,f3
  Dim As Ubyte ub
 
  j=1
  s=1
  d=0
  e=0
  ep=0
  ex=0
  es=1
  i=0
  f=0
  fp=0
  f1=""
  f2=""
  f3=""
  value=Ucase(value)
  fln=Len(value)
 
  While j<=fln
    c=Mid(value,j,1)
    If ep=1 Then
      If c=" " Then
        Goto atof1nxtch
      Endif
      If c="-" Then
        es=-es
        c=""
      Endif
      If c="+" Then
        Goto atof1nxtch
      Endif
      If (c="0") And (f3="") Then
        Goto atof1nxtch
      Endif
      If (c>"/") And (c<":") Then 'c is digit between 0 and 9
        f3=f3+c
        ex=10*ex+(Asc(c)-48)
        Goto atof1nxtch
      Endif
    Endif
   
    If c=" " Then
      Goto atof1nxtch
    Endif
    If c="-" Then
      s=-s
      Goto atof1nxtch
    Endif
    If c="+" Then
      Goto atof1nxtch
    Endif
    If c="." Then
      If d=1 Then
        Goto atof1nxtch
      Endif
      d=1
    Endif
    If (c>"/") And (c<":") Then 'c is digit between 0 and 9
      If ((c="0") And (i=0)) Then
        If d=0 Then
          Goto atof1nxtch
        Endif
        If (d=1) And (f=0) Then
          e=e-1
          Goto atof1nxtch
        Endif
      Endif
      If d=0 Then
        f1=f1+c
        i=i+1
      Else
        If (c>"0") Then
          fp=1
        Endif
        f2=f2+c
        f=f+1
      Endif
    Endif
    If c="E" Then
      ep=1
    Endif
    atof1nxtch:
    j=j+1
  Wend
  If fp=0 Then
    f=0
    f2=""
  Endif

  If s=-1 Then s=&h8000 Else s=0
  n.sign=s
  ex=es*ex-1+i+e
  f1=f1+f2
  f1=Mid(f1,1,1)+Right(f1,Len(f1)-1)
  fln=Len(f1)
  If Len(f1)>(NUM_DIGITS+1) Then
    f1=Mid(f1,1,(NUM_DIGITS+1))
  Endif
  While Len(f1)<(NUM_DIGITS+1)
    f1=f1+"0"
  Wend
  ub=Val("&h"+Mid(f1,1,1))
  n.mantissa(0)=ub
  if ub<>0 then fp=1 else fp=0
  For i=1 To NUM_BYTES-1
      ub=Val("&h"+Mid(f1,2*i,2))
      n.mantissa(i)=ub
      if ub<>0 then fp=1
  Next
  if fp then n.exponent=(ex+BIAS+1) Else n.exponent=0
  str2bcd=1
End Function

Function bcd2str(Byref n As bcdfloat) As String
    Dim As Integer i, ex, s
    Dim As String v,f
    If n.exponent>0 Then
        ex=(n.exponent And &h7FFF)-BIAS-1
    Else
        ex=0
    End If
    If n.sign Then v="-" Else v=" "
   
    v=v+Hex(n.mantissa(0),1)+"."
    For i=1 To NUM_BYTES-1
        v=v+Hex(n.mantissa(i),2)
    Next
    f=Str(Abs(ex))
    f=String(3-Len(f),"0")+f
    If ex<0 Then v=v+"E-" Else v=v+"E+"
    v=v+f
    bcd2str=v
End Function

Sub RSHIFT_NIBB(Byref mantissa As bcdfloat)
    Asm
        mov esi,[mantissa]
        add esi,4
        mov bl,4
        mov dx,NUM_BYTES
        Xor ah,ah
    0:  mov al,[esi]
        mov cl,bl
        ror ax,cl
        mov cl,bl
        Shr ah,cl
        mov [esi],al
        inc esi
        dec dx
        jnz 0b
    End Asm
End Sub

Sub LSHIFT_NIBB(Byref mantissa As bcdfloat)
    Asm
        mov esi,[mantissa]
        add esi,4
        mov bl,4
        mov dx,NUM_BYTES
        add esi,NUM_BYTES-1
        Xor ah,ah
    0:  mov al,[esi]
        mov cl,bl
        rol ax,cl
        mov cl,bl
        Shl ah,cl
        mov Byte [esi],al
        dec esi
        dec dx
        jnz 0b
    End Asm
End Sub

Function NORM_FAC1(Byref fac1 As bcdfloat) As Integer
    ' normalize the number in fac1
    ' all routines exit through this one.
    ' we leave ax=0 if normalization was
    ' ok...else ax = error code.

    'see if the mantissa is all zeros.
    'if so, set the exponent and sign equal to 0.
    Dim As Integer i,er=0,f=0
    For i=0 To NUM_BYTES-1
        If fac1.mantissa(i)>0 Then f=1
    Next
    If f=0 Then
        fac1.exponent=0
        fac1.sign=0
        Exit Function
    'if the highmost nibble in fac1_man is nonzero,
    'shift the mantissa right 1 nibble and
    'increment the exponent
    Elseif fac1.mantissa(0)>9 Then
        RSHIFT_NIBB(fac1)
        fac1.exponent+=1
    Else
    'now shift fac1_man 1 to the left until a
    'nonzero digit appears in the next-to-highest
    'nibble of fac1_man.  decrement exponent for
    'each shift.
        While fac1.mantissa(0)<1
            LSHIFT_NIBB(fac1)
            fac1.exponent-=1
            If fac1.exponent=0 Then
                NORM_FAC1=EXPU_ERR
                Exit Function
            End If
        Wend
    End If
    'check for overflow/underflow
    If fac1.exponent<0 Then
        NORM_FAC1=EXPO_ERR
    End If
End Function

Sub fpadd_aux(Byref fac1 As bcdfloat, Byref fac2 As bcdfloat)
    Asm
        mov edi,[fac1]
        add edi,4 '@fac1.mantissa
        mov esi,[fac2]
        add esi,4 '@fac2.mantissa
        mov ebx,NUM_BYTES-1
        clc
    0:  mov al,[esi+ebx]
        adc al,[edi+ebx]
        daa
        mov [edi+ebx],al
        dec ebx
        jns 0b
    End Asm
    NORM_FAC1(fac1)
End Sub

Sub fpsub_aux(Byref fac1 As bcdfloat, Byref fac2 As bcdfloat)
    Asm
        mov        edi,[fac1]
        add edi,4   '@fac1.mantissa
        mov        esi,[fac2]
        add esi,4   '@fac2.mantissa
        mov        ebx,NUM_BYTES-1
        clc
fad5:        mov        al,[edi+ebx]
        sbb        al,[esi+ebx]
        das
        mov        [edi+ebx],al
        dec        ebx
        jns        fad5
'if there's a borrow, flip sign of fac1
        jnc        fad6
        mov esi,[fac1]  '@fac1.sign
        Xor        word Ptr [esi],&h8000
'adjust
        add        esi,4   '@fac1.mantissa
        mov        ebx,NUM_BYTES-1
        stc
        Xor        ah,ah
fad5a:        mov        al,ah
        sbb        al,[esi+ebx]
        das
        mov        [esi+ebx],al
        dec        ebx
        jns        fad5a
fad6:
    End Asm
    NORM_FAC1(fac1)
End Sub

Function fpadd(Byref x As bcdfloat, Byref y As bcdfloat) As bcdfloat
   
    Dim As bcdfloat fac1,fac2,fac3
    Dim As Integer t
   
    fac1=x
    fac2=y
    If fac1.exponent<fac2.exponent Then 'swap
        fac3=fac1
        fac1=fac2
        fac2=fac3
    End If
   
    'The difference between the two
    'exponents indicate how many times
    'we have to multiply the mantissa
    'of FAC2 by 10 (i.e., shift it right 1 place).
    'If we have to shift more times than
    'we have digits, the result is already in FAC1.
    t=fac1.exponent-fac2.exponent
    If t>0 And t<(NUM_DIGITS+8) Then 'shift
        While t
            RSHIFT_NIBB(fac2)
            t-=1
        Wend
    End If
    'See if the signs of the two numbers
    'are the same.  If so, add; if not, subtract.
    If fac1.sign=fac2.sign Then 'add
        fpadd_aux(fac1,fac2)
    Else
        fpsub_aux(fac1,fac2)
    End If
    fpadd=fac1
End Function

Function fpsub(Byref x As bcdfloat, Byref y As bcdfloat) As bcdfloat
    Dim As bcdfloat fac1,fac2
    fac1=x
    fac2=y
    Asm
        lea edi,[fac2]
        Xor        word Ptr [edi],&h8000
    End Asm
    fac1=fpadd(fac1,fac2)
    fpsub=fac1
End Function

Function fpmul(Byref x As bcdfloat, Byref y As bcdfloat) As bcdfloat
    Dim As bcdfloat fac1,fac2,fac3
    Dim As Integer count, ex, er
    Dim As Ubyte low_nibble
   
    fac1=x
    fac2=y
    'check exponents.  if either is zero,
    'the result is zero
    If fac1.exponent=0 Or fac2.exponent=0 Then 'result is zero...clear fac1.
        fac1.sign=0
        fac1.exponent=0
        For count=0 To NUM_BYTES-1
            fac1.mantissa(count)=0
        Next
        NORM_FAC1(fac1)
        fpmul=fac1
        Exit Function
    Else
        'now determine exponent of result.
        'as you do...watch for overflow.
        ex=fac2.exponent-BIAS+fac1.exponent-1
        If ex<0 Then
            er=EXPO_ERR
            Exit Function
        End If
        fac1.exponent=ex
        'determine the sign of the product
        fac1.sign=fac1.sign Xor fac2.sign
        'copy fac1 mantissa to fac3 and clear fac1's mantissa
        For count=0 To NUM_BYTES-1
            fac3.mantissa(count)=fac1.mantissa(count)
            fac1.mantissa(count)=0
        Next
        'for number of digits in the floating
        'point number, add fac2 to fac1 number
        'of times given by lowest nibble in
        'fac3.  then shift fac1 and fac3
        'right 1 digit and repeat.
        count=NUM_DIGITS-1
        While count
            low_nibble=fac3.mantissa(NUM_BYTES-1) And &h0f
            While low_nibble
                low_nibble-=1
                If low_nibble>=0 Then
                    Asm
                        lea        edi,[fac1]+4
                        lea        esi,[fac2]+4
                        mov        ebx,NUM_BYTES-1
                        clc
                    0:  mov        al,[esi+ebx]
                        adc        al,[edi+ebx]
                        daa
                        mov        [edi+ebx],al
                        dec        ebx
                        jns        0b
                    End asm
                End If
            Wend
            count-=1
            If count Then
                RSHIFT_NIBB(fac1)
                RSHIFT_NIBB(fac3)
            End If
        Wend
    End If   
    NORM_FAC1(fac1)
    fpmul=fac1
End Function

Function fpdiv(Byref x As bcdfloat, Byref y As bcdfloat) As bcdfloat
    Dim As bcdfloat fac1,fac2,fac3
    Dim As Integer count, er
   
    fac1=x
    fac2=y
    If fac2.exponent=0 Then ' if fac2 = 0, return
        ' a divide-by-zero error and
        ' bail out.
        er=DIVZ_ERR
        Exit Function
    Elseif fac1.exponent=0 Then 'fact1=0, just return
        er=0
        fpdiv=fac1
        Exit Function
    Else
        'calculate the exponent of the quotient.
        fac1.exponent=fac1.exponent-fac2.exponent+BIAS
        'calculate the sign of quotient.
        fac1.sign=fac1.sign Xor fac2.sign
        ' clear fac3 - this will hold the mantissa of the result

        For count=0 To NUM_BYTES-1
            fac3.mantissa(count)=0
        Next
        'set up for looped-subtract
        'subtract fac2 from fac1
        Asm
            mov        ecx,NUM_DIGITS
    fpd2a:  lea        edi,[fac1]+4
            lea        esi,[fac2]+4
            mov        ebx,NUM_BYTES-1
            clc
        0:        mov        al,[edi+ebx]
            sbb        al,[esi+ebx]
            das
            mov [edi+ebx],al
            dec ebx
            jns 0b
            'did the subtraction fail? if so,
            'add fac2 back into fac1
            jnc fpd5
            lea        edi,[fac1]+4
            lea        esi,[fac2]+4
            mov        ebx,NUM_BYTES-1
            clc
        0:  mov al,[esi+ebx]
            adc al,[edi+ebx]
            daa
            mov [edi+ebx],al
            dec ebx
            jns 0b
            jmp fpd6
            'if subtract succeeded, add 1 to
            'lowmost nibble of fac3 and try another
            'subtract.
    fpd5:   lea edi,[fac3]+4+NUM_BYTES-1
            inc Byte Ptr [edi]
            jmp        fpd2a
            'shift fac3 and fac1 left 1 nibble
    fpd6:        dec        ecx                        'done?
            jz        fpd7
            push ecx
            lea esi,[fac3]+4
            mov bl,4
            mov dx,NUM_BYTES
            add esi,NUM_BYTES-1
            Xor ah,ah
    0:      mov al,[esi]
            mov cl,bl
            rol ax,cl
            mov cl,bl
            Shl ah,cl
            mov Byte [esi],al
            dec esi
            dec dx
            jnz 0b
            lea esi,[fac1]+4
            mov bl,4
            mov dx,NUM_BYTES
            add esi,NUM_BYTES-1
            Xor ah,ah
    0:      mov al,[esi]
            mov cl,bl
            rol ax,cl
            mov cl,bl
            Shl ah,cl
            mov Byte [esi],al
            dec esi
            dec dx
            jnz 0b
            'loop and repeat for number of digits
            pop ecx
            jmp        fpd2a
    fpd7:       
        End asm
    End If
        'move fac3_man into fac1_man
        For count=0 To NUM_BYTES-1
            fac1.mantissa(count)=fac3.mantissa(count)
        Next
        'exit through normalization routine.       
    NORM_FAC1(fac1)
    fpdiv=fac1
End Function

Declare Function cbcdfloat Overload ( Byref lhs As bcdfloat ) As bcdfloat
Declare Function cbcdfloat ( Byref lhs As String )  As bcdfloat

Function cbcdfloat ( Byref lhs As bcdfloat ) As bcdfloat
    Return lhs
End Function

Function cbcdfloat ( Byref lhs As String ) As bcdfloat
    Dim As bcdfloat retval
    str2bcd (  lhs, retval )
    Return retval
End Function

Operator + ( Byref lhs As bcdfloat, Byref rhs As bcdfloat ) As bcdfloat
    Dim As bcdfloat retval
    retval = fpadd(lhs, rhs )
    Return retval
End Operator

'=============================================================
Operator bcdfloat.+= ( byref rhs as bcdfloat )
    this = fpadd(this, rhs)
End Operator
'=============================================================

Operator - ( Byref lhs As bcdfloat, Byref rhs As bcdfloat ) As bcdfloat
    Dim As bcdfloat retval
    retval = fpsub( lhs, rhs )
    Return retval
End Operator

'=============================================================
Operator bcdfloat.-= ( byref rhs as bcdfloat )
    this = fpsub(this, rhs)
End Operator
'=============================================================

Operator - ( Byref lhs As bcdfloat ) As bcdfloat
    dim as bcdfloat z=lhs
    z.sign = z.sign xor &H8000
    operator = z
End Operator

Operator * ( Byref lhs As bcdfloat, Byref rhs As bcdfloat ) As bcdfloat
    Dim As bcdfloat retval
    retval = fpmul(lhs, rhs )
    Return retval
End Operator

Operator / ( Byref lhs As bcdfloat, Byref rhs As bcdfloat ) As bcdfloat
    Dim As bcdfloat retval
    retval = fpdiv( lhs, rhs )
    Return retval
End Operator

'=============================================================

function compare ( Byref lhs As bcdfloat, Byref rhs As bcdfloat ) As Integer
    Dim As integer relop,i
    Dim as bcdfloat comp
    comp = fpsub(lhs,rhs)
    relop=comp.exponent
    for i=0 to NUM_BYTES-1
        relop=relop + comp.mantissa(i)
    next
    if relop>0 then
        if comp.sign=&h8000 then
            relop=-1
        else
            relop=1
        end if
    end if
    Return relop
End function

Operator = ( Byref lhs As bcdfloat, Byref rhs As bcdfloat ) As Integer
    Dim As integer relop
    relop=compare(lhs,rhs)
    Return relop=0
End Operator

Operator < ( Byref lhs As bcdfloat, Byref rhs As bcdfloat ) As Integer
    Dim As integer relop
    relop=compare(lhs,rhs)
    Return relop<0
End Operator

Operator > ( Byref lhs As bcdfloat, Byref rhs As bcdfloat ) As Integer
    Dim As integer relop
    relop=compare(lhs,rhs)
    Return relop>0
End Operator

Operator <= ( Byref lhs As bcdfloat, Byref rhs As bcdfloat ) As Integer
    Dim As integer relop
    relop=compare(lhs,rhs)
    Return relop<=0
End Operator

Operator >= ( Byref lhs As bcdfloat, Byref rhs As bcdfloat ) As Integer
    Dim As integer relop
    relop=compare(lhs,rhs)
    Return relop>=0
End Operator

'=============================================================

Operator bcdfloat.cast ( ) As String
    Operator = bcd2str( this )
End Operator

Operator bcdfloat.cast ( ) As Double
    Operator = val(bcd2str( this ))
End Operator

constructor bcdfloat ( Byref rhs As bcdfloat )
   this = rhs
end constructor

constructor bcdfloat ( Byref rhs As String )
    str2bcd ( rhs, this )
end constructor

constructor bcdfloat ( Byref rhs As Double )
    str2bcd ( str(rhs), this )
end constructor

constructor bcdfloat ( Byref rhs As Longint )
    str2bcd ( str(rhs), this )
end constructor

operator bcdfloat.let ( Byref rhs As bcdfloat )
    Dim As integer i
    this.sign=rhs.sign
    this.exponent=rhs.exponent
    for i=0 to NUM_BYTES-1
        this.mantissa(i)=rhs.mantissa(i)
    next
end operator

operator bcdfloat.let ( Byref rhs As String )
    str2bcd( rhs, this )
end operator

operator bcdfloat.let ( Byval rhs As Double )
    str2bcd( str(rhs), this )
end operator

operator bcdfloat.let ( Byval rhs As Longint )
    str2bcd( str(rhs), this )
end operator

'=========================================================
'' For Next for extended type
''
'' implicit step versions
''
'' In this example, we interpret implicit step
'' to mean 1
Operator bcdfloat.for( )
End Operator
 
Operator bcdfloat.step( )
        this += 1 'this = this+1 '
End Operator
 
Operator bcdfloat.next( Byref end_cond As bcdfloat ) As Integer
        Return this <= end_cond
End Operator
 
 
'' explicit step versions
''
Operator bcdfloat.for( Byref step_var As bcdfloat )
End Operator
 
Operator bcdfloat.step( Byref step_var As bcdfloat )
        this += step_var 'this = this + step_var '   
End Operator
 
Operator bcdfloat.next( Byref end_cond As bcdfloat, Byref step_var As bcdfloat ) As Integer
        If step_var < 0 Then
                Return this >= end_cond
        Else
                Return this <= end_cond
        End If
End Operator

'=============================================================

Function even(Byref x As bcdfloat) As bcdfloat
    Dim As bcdfloat mxs, term, sigma, last_sigma

    mxs = -x * x
    term = 1.0
    sigma = 1.0
    Dim As Longint n = 0
    Do
        n = n + 2
        last_sigma = sigma
        term = term * mxs / (n * (n-1))
        sigma = sigma + term
    Loop Until last_sigma = sigma
    Return sigma
End Function

'=====================================================================
dim as double et=timer
Dim As bcdfloat last_x, x, y
#if 0
x = 1.57
Do
    Print "2*x = ";x * 2
    last_x = x
    y = even(x)
    x = x + y
Loop Until last_x = x
Print timer-et;" seconds"
print "       3.14159265358979323846264338327950288419716939937510582097494459230781640"

'=====================================================================
Print "Normal exit encountered."
#endif

for x=1 to 10
    ? using("###");x;"===";
    ? x/9
next
Print "press Return to End ";
Sleep

Return to “Tips and Tricks”

Who is online

Users browsing this forum: No registered users and 0 guests