Exact Double to Decimal Expansion for FreeBASIC

dbl2dec.bas

'' ============================================================
'' Dbl2Dec.bas -- Double float to decimal string expansion
''
'' Copyright (C) 2018 - Jeff Marshall (coder@execulink.com)
''
'' License: 1) Program is provided as-is.  
''          2) There is no warranty or guarantee.
''          3) There are no other conditions.
''
'' Feel free to use, re-use, or abuse :)
''
'' ============================================================

/'
Double float to decimal string expansion.

Every double floating point number has an exact 
decimal value.  Basic steps for expansion are:
	1) decode the sign, exponent, & fraction from
	   DOUBLE floating point data
	2) test for special values and add in
	   the hidden 1 bit (or 0 bit for sub-normals)
	3) depending on the exponent, extract integer
	   and/or fraction part
		a) for small exponents we can use a one or
		   two 64-bit variables, one for integer
		   and one for fraction
		b) for large positive exponents use a
		   BIGNUM integer
		c) for large negative exponents use a
		   BIGNUM fraction

For BIGNUM integer & fraction, three operations are
needed:
	BigLOAD64 - load a 64-bit significand + exponent
	BigIMUL10 - integer multiply by 10
	BigIDIV10 - integer devide by 10.  Also returns
		        the modulo.
'/

'' ------------------------------------
'' BIGNUM stuff
'' ------------------------------------

#define BIGDIGITS 34

'' number of 2^32 digits for the BIGNUM is determined by the fact
'' that double's base-2 exponent can be -1022 to 1023.  The 
'' expansion handles positive & negative exponents separately, so
'' expansion needs at least 1024 bits.  We add another 64 bits 
'' for the significand, of which only 53 bits are significant and
'' leaving us room for the computations.  This gives us a total 
'' of 1088 bits or 34 base 2^32 digits.

''
function BigIsZero( b() as ulongint ) as boolean
	'' return true if all bits are zero
	for i as integer = 0 to BIGDIGITS-1
		if( b(i) <> 0 ) then
			return false
		end if
	next
	return true
end function

sub BigLOAD64( b() as ulongint, f as ulongint, e as longint )

	'' we are using 64-bit data type to hold 32-bit digits.
	'' This makes computations easier since fbc can handle 
	'' 64-bit integers and we don't need to worry about 
	'' overflows when working with 2^32 bit data values in 
	'' 2^64 bit data types.
		
	#define LOMASK &h00000000ffffffffull
	#define HIMASK &hffffffff00000000ull

	'' clear the BIGNUM
	for i as integer = 0 to BIGDIGITS-1
		b(i) = 0
	next

	'' BIGNUM integer
	'' if e >= 0 then this is a BIGNUM integer, and assume
	'' that caller has passed in f and e as follows:
	''        0 <= f < 2^53
	''        0 <= e <= 1023
	''        number = f * 2 ^ e
	
	if( e >= 0 ) then
		
		'' find the starting digit index (equivalent to 
		'' shifting some multiple of 32)
		dim i as integer = e \ 32
		e mod= 32

		'' load the initial value and split 64-bit input into
		'' two 32-bit digits
		b(i+0) = f and LOMASK
		b(i+1) = (f shr 32) and LOMASK

		'' shift to left (multiply by 2) until exponent = 0
		while( e > 0 )

			'' bit shifted out/in - initially zero
			dim c as ulongint = 0

			'' perform SHL 1 for 3 32-bit digits; we never
			'' need to go more than 4 digits because e < 32.
			'' Assume that the caller has allocated enough 
			'' digits that we won't overlow the BIGNUM.

			for j as integer = 0 to 2
				b(i+j) shl= 1
				b(i+j) or= c
				c = b(i+j) 
				c shr= 32
				b(i+j) and= LOMASK
			next
			
			e -= 1
		wend

		return

	end if


	'' BIGNUM fraction

	'' if e < 0, then this is a BIGNUM fraction and we
	'' assume that the fraction has been aligned to
	'' 60 bits on the input.
	''       number = (f / 2^60) * 2^(e + 60)

	'' find the starting digit index (equivalent to 
	'' shifting some multiple of 32)
	dim i as integer = (BIGDIGITS-1) + e \ 32
	e mod= 32

	'' load the initial value
	b(i-0) = (f shr 32) and LOMASK
	b(i-1) = f and LOMASK

	'' shift to right (divide by 2) until e = 0
	while( e < 0 )

		'' bit shifted out/in - initially zero
		dim c as ulongint = 0

		'' perform SHR 1 for 3 32-bit digits; we never
		'' need to go more than 4 digits because e < 32.
		'' Assume that the caller has allocated enough 
		'' digits that we won't underflow the BIGNUM.

		for j as integer = 0 to 2
			b(i-j) or= (c shl 32)
			c = b(i-j) and 1
			b(i-j) shr= 1
		next
		
		e += 1
	wend
		
end sub

sub BigIMUL10( b() as ulongint )
	
	'' multiply BIGNUM by 10
	
	'' carry from one digit to the next
	dim as ulongint c = 0

	for i as integer = 0 to BIGDIGITS-1
		
		'' multiply the digit by 10
		b(i) *= 10

		'' add in carry from last operation
		b(i) += c

		'' save the carry for the next digit
		c = b(i) shr 32

		'' make sure digit is < 2^32
		b(i) and= &h00000000ffffffffull
	next

end sub

sub BigIDIV10( b() as ulongint, byref r as ubyte )
	
	'' integer quotient
	dim q as ulongint

	'' modulus 10 (remainder) is never greater than 9
	dim m as ulongint = 0

	'' integer divide each base-2^32 digit by 10
	'' and add the remainder to the next digit
	dim i as integer = BIGDIGITS-1

	'' skip leading zero value digits
	while( i > 0 andalso b(i)=0 )
		i -= 1
	wend
	
	while( i >= 0 )

		'' add in the remainder from last operation
		b(i) += (1ull shl 32) * m

		'' integer division
		q = b(i) \ 10

		'' remainder (carry forward to next operation)
		m = b(i) - q * 10

		'' store result
		b(i) = q

		'' next digit
		i -= 1
	wend

	'' return the remainder (always less than 10)
	r = m

end sub

'' ------------------------------------
'' Dbl2Dec
'' ------------------------------------

function Dbl2Dec( byval fp as double ) as string

	dim as ulongint i = any  '' fp converted to 64 bits
	dim as ulongint b = any '' biased exponent 0 to &h7ff
	dim as byte     s = any '' sign 0=positive or 1=negative
	dim as longint  e = any '' exponent -1022 to 1023 for normal numbers
	dim as ulongint f = any '' fraction part
	dim as ulongint n = any '' integer part
	
	i = *cast( ulongint ptr, @fp )

	'' decimal result
	dim as string d

	'' extract sign; sign 0=positive or 1=negative
	s = i shr 63

	'' extract biased exponent
	b = (i and &h7ff0000000000000ull) shr 52

	'' extract mantissa
	f = (i and &h000fffffffffffffull)

	'' NAN? or INF?
	if( b = &h7ffull ) then
		if( f <> 0 ) then
			d = "NAN"
		else
			d = "1.#INF"
		end if
		if( s ) then
			d = "-" & d
		end if
		return d
	end if

	'' subnormal? or signed-zero?
	if( b = &h000 ) then
		if( f <> 0 ) then
			'' subnormal
			e = 1-1023
		else
			'' signed zero
			e = 0
		end if
	else
		e = b - 1023
		'' add the hidden 1 bit - making fraction 53 bits
		f or= &h0010000000000000ull			  
	end if

	'' we now have:
	''    s = sign     (1 bit)
	''    e = exponent (11 bits)
	''    n = integer  (1 bit) 
	''    f = fraction (52 bits)
	''  decimal = s * 2^e * f, -1022 <= e <= 1023

	'' move decimal place to the left by one and add one
	'' to the exponent to include the integer part in the
	'' fraction, so now we have:
	''    n = integer  (0 bits )
	''    f = fraction (53 bits)
	e += 1

	'' exponent determines what method to use for
	'' decimal expansion.  There is overlap in the ranges at
	'' which each method will work for a given exponent.
	'' Only one method actually gets executed.

	if( e >= 0 and e <= 53 ) then
		'' expand INTEGER + FRACTION 
		'' to decimal using two 64 bit variables

		'' extract integer part and fraction part
		n = f shr (53 - e)
		f = (f shl e) and &h001fffffffffffffull

		'' expand integer
		if( n ) then
			while( n )
				d = chr(48 + (n mod 10)) & d
				n \= 10
			wend
		else
			d &= "0"
		end if

		'' need decimal if there is a fraction
		if( f ) then
			d &= "."

			'' expand 53 bit binary fraction to decimal
			while( f )
				f *= 10
				d &= chr( asc("0") + (f shr 53) )
				f and= &h001fffffffffffffull
			wend
		end if

	elseif( e >= 53 andalso e <= 64 ) then
		'' expand INTEGER 
		'' to decimal using one 64-bit variable
		
		'' convert      n  = f/(2^53) * 2^(e)
		''      to      n  = f'/(2^64)
		''  where       f' = f * 2^(e-53)

		n = f shl (e-53)

		'' expand integer
		if( n ) then
			while( n )
				d = chr(48 + (n mod 10)) & d
				n \= 10
			wend
		else
			d &= "0"
		end if

	elseif( e >= -7 andalso e <= 0 ) then
		'' expand FRACTION
		'' to decimal using one 64-bit variable
	
		'' need to reserve 4 bits for fraction expansion 
		'' convert      n = f/(2^53) * 2^(e)
		''      to      n = f'/(2^60)
		''   where     f' = f * 2^(e+7)
		f = f shl (e+7)

		d &= "0."

		'' expand 60 bit binary fraction to decimal
		while( f )
			f *= 10
			d &= chr( asc("0") + (f shr 60) )
			f and= &h0fffffffffffffffull
		wend

	elseif( e > 53 ) then
		'' expand INTEGER
		'' to decimal using BIGNUM integer

		'' place decimal point to the right
		e -= 53

		'' load a BIGNUM with the integer
		dim b( 0 to BIGDIGITS-1 ) as ulongint
		BigLOAD64( b(), f, e )

		'' expand big integer to decimal
		while( not BigIsZero(b()) )
			dim r as ubyte = 0
			BigIDIV10( b(), r )
			d = chr(48 + r) & d
		wend

	elseif( e < 0 ) then
		'' expand FRACTION
		'' to decimal using BIGNUM fraction

		'' shift fraction to align to 60 bits; we only
		'' need top 4 bits available to expand the decimal
		f shl= 7

		'' load a BIGNUM with the integer
		dim b( 0 to BIGDIGITS-1 ) as ulongint
		BigLOAD64( b(), f, e )

		d &= "0."

		'' expand binary fraction to decimal, next
		'' digit is in the top 4 bits
		while( not BigIsZero( b() ) )
			BigIMUL10( b() )
			d &= chr( asc("0") + (b(BIGDIGITS-1) shr 28) )
			b(BIGDIGITS-1) and= &h0fffffffull
		wend

	end if

	'' add the sign
	if( s ) then
		d = "-" & d
	end if

	'' the expanded decimal as string is returned
	return d

end function

dbl2dec_test.bas

#include "dbl2dec.bas"

sub show( byval fp as double )
	print left( fp & space(25), 25 ) & " = ";Dbl2Dec( fp )
end sub

show( 0.6804801726248115 )
show( 0.5 )
show( -0.5 )
show( 0.89 )
show( 1234.5 )
show( 1.5e+23 )
show( 1D+50 )
show( 1D+51 )
show( 2^64 )
show( 1d-15 )

/' OUTPUT:

0.6804801726248115        = 0.68048017262481153011322021484375
0.5                       = 0.5
-0.5                      = -0.5
0.89                      = 0.89000000000000001332267629550187848508358001708984375
1234.5                    = 1234.5
1.5e+023                  = 150000000000000004194304
1e+050                    = 100000000000000007629769841091887003294964970946560
1e+051                    = 999999999999999993220948674361627976461708441944064
1.844674407370955e+019    = 18446744073709551616
1e-015                    = 0.00000000000000100000000000000007770539987666107923830718560119501514549256171449087560176849365234375

'/