ThinBasic extended numbers in FreeBASIC SDK.

[D.J.Peters]

hello Charles
i opened a new post for Extended numbers in FreeBASIC SDK
(how to call a user script function from TB_Module? isn't the right place)

why do you make a string from DOUBLE (this done by FreeBASIC) we need EXT

i changed this

operator EXT.cast as string
dim as any ptr p = @value
dim as double d = any
asm
mov eax,[p]
fldt [eax]
fstp qword ptr [d]
end asm
return float_to_ascii(d)
end operator
to
operator EXT.cast as string
dim as any ptr p = @value
return float_to_ascii(p)
end operator]
now you must increase the max decimal places from 17 to ?
(is it 23 with Extended numbers i don't know)

and take a look to this line "print c,dV,fV"

Joshy

' Extended Numbers for FreeBASIC
' by Joshy and Charles

' assignment (Ext,Double,Single,integer,short,byte)
' casting (String,Double,Single,integer,short,byte)
' Math (+,-,*,/)
' Ext=SQRT(Ext,Double,Single,integer,short,byte)
type NumberFormat
dp as long ' DECIMAL PLACES
sn as long ' SCIENTIFIC NOTATION
trz as long ' STRIP END ZEROS
sdp as long ' ALLOW START WITH DECIMAL POINT
end type
dim shared as NumberFormat num

' EXTENDED TYPE
type EXT
as byte Value(9)
declare constructor
declare constructor(byref l as EXT)
declare constructor(byval l as double)
declare operator let(byval r as double)
declare operator cast as string
declare operator cast as double
end type
declare function float_to_ascii(pa as ext ptr) as string

' constructors
constructor EXT
' dummy constructor
end constructor

constructor EXT(byref l as EXT)
dim as any ptr d =@value(0),s=@l.value(0)
asm
mov eax,[s]
fldt [eax]
mov eax,[d]
fstpt [eax]
end asm
end constructor

constructor EXT(byval l as double)
dim as any ptr p =@value(0)
asm
fld qword ptr [l]
mov eax,[p]
fstpt [eax]
end asm
end constructor

' asignment
operator EXT.let(byval l as double)
dim as any ptr p =@value(0)
asm
fld qword ptr [l]
mov eax,[p]
fstpt [eax]
end asm
end operator

operator EXT.cast as string
dim as any ptr p =@value(0)
return float_to_ascii(p)
end operator

operator EXT.cast as double
dim as any ptr p =@value(0)
dim as double d
asm
mov eax,[p]
fldt [eax]
fstp qword ptr [d]
end asm
return d
end operator

Operator + (ByRef l As EXT, ByRef r As EXT) As EXT
dim as EXT t
dim as any ptr p = @t.value(0)
asm
mov eax,[l]
fldt [eax]
mov eax,[r]
fldt [eax]
faddp
mov eax,[p]
fstpt [eax]
end asm
Return t
End Operator

Operator - (ByRef l As EXT, ByRef r As EXT) As EXT
dim as EXT t
dim as any ptr p = @t.value(0)
asm
mov eax,[l]
fldt [eax]
mov eax,[r]
fldt [eax]
fsubp
mov eax,[p]
fstpt [eax]
end asm
Return t
End Operator


Operator / (ByRef l As EXT, ByRef r As EXT) As EXT
dim as EXT t
dim as any ptr p = @t.value(0)
asm
mov eax,[l]
fldt [eax]
mov eax,[r]
fldt [eax]
fdivp
mov eax,[p]
fstpt [eax]
end asm
Return t
End Operator

Operator * (ByRef l As EXT, ByRef r As EXT) As EXT
dim as EXT t
dim as any ptr p = @t.value(0)
asm
mov eax,[l]
fldt [eax]
mov eax,[r]
fldt [eax]
fmulp
mov eax,[p]
fstpt [eax]
end asm
Return t
End Operator

declare function SQREXT overload(byref r as EXT) as EXT
declare function SQREXT overload(byval r as double) as EXT


function SQREXT(byref r as EXT) as EXT
dim as EXT t
dim as any ptr p = @t.value(0)
asm
mov eax,[r]
fldt [eax]
fsqrt
mov eax,[p]
fstpt [eax]
end asm
Return t
end function

function SQREXT(byval r as double) as EXT
dim as EXT t
dim as any ptr p = @t.value(0)
asm
fld qword ptr [r]
fsqrt
mov eax,[p]
fstpt [eax]
end asm
Return t
end function

function SQRT(byval r as integer) as EXT
dim as EXT t
dim as any ptr p = @t.value(0)
asm
fild dword ptr [r]
fsqrt
mov eax,[p]
fstpt [eax]
end asm
Return t
end function


function float_to_ascii(pa as ext ptr) as string
static as ubyte s(23), t(23) 'BUFFERS
static as ubyte bcd(11)
dim as long esize,tempdw,dp,sn,snv,b,nzero,oldcw,truncw
static as zstring ptr ps,pt

dim as long num_trz,num_dp,num_sn,num_sdp

num_trz = num.trz
num_dp = num.dp
num_sn = num.sn
num_sdp = num.sdp

ps=varptr(s(0))
pt=varptr(t(0))

asm
mov eax,0
mov [snv],eax
mov [sn],eax
mov [nzero],eax
mov ecx,[pa]
mov eax,[ecx+8]
and eax,&h7fff
cmp eax,0
jnz xfa2 'exit EXCLUDE NON ZERO
cmp dword ptr [num_trz],0
jz xfa2 'exit ZERO STRIPPER FLAG
or eax,[ecx+4]
jnz xfa2 'exit
test byte ptr [ecx+9],&h80
jz xfa1 'exit
end asm
*pt="-0"
asm
jmp fadonez 'NEGATIVE ZERO
xfa1:
end asm

*pt="0"
asm
jmp fadonez 'POSITIVE ZERO
xfa2:
cmp eax,&h7fff
jnz xfa5n 'exit
mov dword ptr [nzero],1
mov eax,[ecx+4]
and eax,0x7fffffff 'exclude bit 63
or eax,[ecx]
cmp eax,0
jz xfa4 'exit
test byte ptr [ecx+7],&h40
jnz xfa3 'exit
end asm

*pt="#sNAN"
asm
jmp fadonez 'SIGNALLING NAN

xfa3:
end asm

*pt="#qNAN"
asm
jmp fadonez 'QUIET NAN
xfa4:
'NEGATIVE / POSITIVE INFINITY
test byte ptr [ecx+9],&h80
jz xfa5 'exit
end asm

*pt="#-INF"
asm
jmp fadonez 'NEGATIVE INFINITY
xfa5:
end asm

*pt="#INF"
asm
jmp fadonez 'POSITIVE INFINITY
xfa5n:
mov eax,[pa]
fldt [eax]
'CHECK FOR ZERO
mov dword ptr [esize],0
fldz
fcomip st(0),st(1)
jz xfa6 'exit
mov dword ptr [nzero],1

fldlg2 ' log10(2)
fld st(1) ' copy Src
fabs ' insures a positive value
fyl2x ' ->[log2(Src)]*[log10(2)] = log10(Src)
fstcw [oldcw] ' get current control word
fwait
mov ax,[oldcw]
or ax,&hc00 ' code it for truncating
mov [truncw],ax
fldcw [truncw] ' insure rounding code of FPU to truncating

fist dword ptr [esize] ' store characteristic of logarithm
fldcw [oldcw] ' load back the former control word

ftst ' test logarithm for its sign
fstsw ax ' get result
fwait
sahf ' transfer to CPU flags
sbb dword ptr [esize],0 ' decrement esize if log is negative
fstp st(0) ' get rid of the logarithm

xfa6:
' DECIMAL PLACES LIMIT
mov eax,[num_dp]
cmp eax,17
jle xfa7 'exit
mov eax,17 'LIMIT DECIMAL PLACES
xfa7:
mov [dp],eax
' IS SCIENTIFIC NOTATION ALWAYS REQUIRED
cmp byte ptr [num_sn],0
jnz ENotation
' VERY LARGE NUMBERS
mov ecx,[esize]
cmp ecx,18
jl xfa8 'exit
jmp ENotation
xfa8:
'SMALL NUMBERS
cmp dword ptr [esize],0
jge xfa9 'exit
mov ecx,[dp]
mov edx,[esize]
neg edx
cmp edx,4
jg ENotation 'LIMIT FOR SIMPLE FORMAT
mov eax,ecx
mov [dp],ecx
jmp PowerAdjust
xfa9:
' NUMBERS NOT REQUIRING SCIENTIFIC NOTATION
mov eax,[dp]
mov ecx,eax 'DECIMAL PLACES
add ecx,[esize] 'INTEGER DIGITS
sub ecx,17
jle xfa10 'exit
' TOO MANY DIGITS? (ecx contains excess digits)
sub eax,ecx 'REDUCE MULTIPLIER PLACES IF NECESSARY
sub [dp],ecx 'REDUCE DECIMAL PLACES ALSO
xfa10:
jmp PowerAdjust
ENotation:
mov ecx,[esize]
mov [snv],ecx
mov eax,[dp]
sub eax,ecx
mov dword ptr [sn],1 'SCIENTIFIC NOTATION FLAG

PowerAdjust:
mov [tempdw],eax 'ADJUSTED MULTIPLIER
' Multiply the number by the power of 10
mov eax,tempdw
cmp eax,0
jz xfa11 'exit
fild dword ptr [tempdw]
fldl2t
fmulp st(1),st(0) '->log2(10)*exponent
fld st(0)
frndint 'get the characteristic of the log
fxch st(1)
fsub st(0),st(1) 'get only the fractional part but keep the characteristic
f2xm1 '->2^(fractional part)-1
fld1
faddp st(1),st(0) 'add 1 back
fscale 're-adjust the exponent part of the REAL number
fstp st(1) 'get rid of the characteristic of the log
fmulp st(1),st(0) '->16-digit integer
xfa11:
fbstp [bcd] 'SAVE AS PACKED BINARY CODED DECIMAL
' EXPAND DIGITS
lea edx,[bcd]
lea ecx,
push ebx
mov bl,10
rfa12:
dec bl
jl xfa12 'exit
mov ah,[edx]
inc edx
mov al,ah
and al,15
add al,48
mov [ecx],al
inc ecx
mov al,ah
shr al,4
and al,15
add al,48
mov [ecx],al
inc ecx
jmp rfa12 'repeat

xfa12:
mov byte ptr [ecx],0
' COPY FORMATTED
lea edx,[t]
' NEGATIVE SIGN NEEDED?
cmp al,48
jz xfa13 'exit
mov byte ptr [edx],45
inc edx
xfa13:
lea ebx,
add ebx,18
mov cl,19
mov ah,0
mov ch,[dp]
dec ch
rfa18:
dec cl
jl xfa18 'exit
' INSERT DECIMAL POINT
cmp cl,ch
jnz xfa15 'exit
' PLACE LEADING ZERO
cmp dword ptr [num_sdp],0
jnz xfa14 'exit SDP FLAG TO INHIBIT
cmp ah,0
jnz xfa14 'exit
mov byte ptr [edx],48
inc edx
xfa14:
mov byte ptr [edx],46
inc edx
mov ah,1 'STOP STRIPPING ZEROS

xfa15:
mov al,[ebx]
dec ebx
cmp ah,0
jnz faok

xfa16:
cmp al,48
jz xfa17 'exit STRIP LEADING ZEROS
mov ah,1 'INHIBIT FUTURE STRIPPING
faok:
mov [edx],al
inc edx

xfa17:
jmp rfa18 'repeat

xfa18:
pop ebx
fadone:
' REMOVE ENDING ZEROS
cmp dword ptr [num_trz],0
jz xfa21 'exit
lea ecx,[t] 'BASE ADDRESS OF NUMBER STRING
rfa20:
dec edx
cmp edx,ecx
jle xit1 'exit LEAVE FIRST CHARACTER ALONE
mov al,[edx]
cmp al,46
jnz xfa19 'exit
dec edx
jmp xit1 'STRIP DOT AND EXIT
xfa19:
cmp al,48
jnz xfa20 'exit ONLY LOOK AT RIGHT HAND ZEROS
jmp rfa20 'repeat 'STRIP ZERO AND CONTINUE
xfa20:
xit1:
inc edx
xfa21:
' ENSURE AT LEAST ONE DIGIT
lea ecx,[t]
mov al,[ecx]
cmp al,45
jnz xfa22 'exit
inc ecx
xfa22:
cmp ecx,edx
jnz xfa23 'exit
mov byte ptr [edx],48
inc edx
xfa23:
' CHECK FOR SCIENTIC NOTATION
cmp dword ptr [sn],0
jz xfa27 'exit
cmp dword ptr [nzero],0
jz xfa27 'exit
mov eax,[snv]
cmp eax,0
jz xfa27 'exit E VALUE ZERO SO OMIT
mov byte ptr [edx],69
inc edx 'E'
mov cl,43
cmp eax,0
jge xfa24 'exit
neg eax
mov cl,45
xfa24:
mov [edx],cl
inc edx 'SIGN
mov cl,100
div cl
push eax
and eax,&hff
mov cl,10
div cl
cmp ax,0
jz xfa26 'exit
or eax,&h3030 'TO ASCII THOUSANDS AND HUNDREDS
cmp al,48
jz xfa25 'exit
mov [edx],al
inc edx
xfa25:
mov [edx],ah
inc edx
xfa26:
pop eax
shr eax,8
div cl
or eax,&h3030 'TO ASCII TENS AND UNITS
mov [edx],aX
add edx,2
xfa27:
mov byte ptr [edx],0 'APPEND NULL TERMINATOR
fadonez:
end asm
function=rtrim(*pt)
end function



dim as EXT a = 1
dim as EXT b = 3.0
dim as EXT c = a+b*2
dim as EXT d = a/b

' WITH NUMBER FORMAT CONTROL
num.dp =17 ' DECIMAL PLACES
num.trz= 0 ' STRIP TRAILING ZEROS
num.sn = 0 ' SCIENTIFIC NOTATION BY DEFAULT
num.sdp= 0 ' INHIBIT ZERO BEFORE DECIMAL POINT
print a,b,c,d
print

num.dp =6 ' DECIMAL PLACES
num.trz=1 ' STRIP TRAILING ZEROS
num.sn =1 ' SCIENTIFIC NOTATION BY DEFAULT
num.sdp=0 ' INHIBIT ZERO BEFORE DECIMAL POINT
print a,b,c,d
print

num.dp =23 ' DECIMAL PLACES
num.trz=0 ' STRIP TRAILING ZEROS
num.sn =1 ' SCIENTIFIC NOTATION BY DEFAULT
num.sdp=0 ' INHIBIT ZERO BEFORE DECIMAL POINT
print a,b,c,d
print

d=16:a=SQRT(d):b=SQRT(16.0):c=SQRT(16)
print a,b,c

' double,single,int,short,byte casting
c=1.23456789012345678
dim as double dV=c
dim as single fV=c
dim as integer iV=c
dim as short sV=c
dim as byte bV=c

num.sn =1 ' SCIENTIFIC NOTATION BY DEFAULT
num.sdp=1 ' INHIBIT ZERO BEFORE DECIMAL POINT
print "Ext"," ","Double"," ","Single"
print c,dV,fV
print "Integer","Short","Byte"
print iV,sV,bV

sleep

[jack]

Hi Zak,

Yes there is something strange going on. With the Oxygen version, I get an accurate 0.333333333333333. But with large powers, say 1e+200 / 3 the final 2 digits are out: 3.33333333333333320E+19.

I'll keep on experimenting.

Could you try running some big thirds and see what results you get with your code?

Charles

I get a string of 18 3's, I had trouble with fscale function way back it was not very accurate, but on newer procceessors it should work OK, when I sober up I will do some testing.

here's my code, like I said I got the NaN's wrong

[Charles Pegge]

Many thanks Joshy,

We can see that handling all these different types takes a lot of code!

Extended numbers have a 64 bit fraction which imposes a natural limit of around 18 significant digits. FBSTP stores 18 digits in packed BCD form, taking up 9 bytes. The final byte carries the sign bit.

For many operations the last one or two digits may be inaccurate. I am investigating what situations cause this.

Charles

[Charles Pegge]

Thanks for your code Jack. I will study it later. My bed time :)

Charles

[jack]

Charles, my code is messy and ugly but maybe you can use some portions after cleaning it up. :)

[Charles Pegge]

Hi Jack,

Your code looks fine too me and very comprehensive. It shows that producing a full implementation of the Extended type for FreeBasic takes a lot of code. You have over 100k of source code there. It is a valuable resource. Many thanks!

One technique I am developing in Oxygen is to implement a type called FPU which can be used to cover any type which can be passed on the FPU stack. This radically reduces all the combinatorial overloading that is necessary in FB.

Charles

[Charles Pegge]

Hi Joshy,

I've added Ascii to Float, complementing Float to Ascii.

Also a few more tests

Charles

[jack]

Charles here's some code in C to do io of long double, specifically have a look at ioldouble.c which is in the public domain http://www.mastodon.biz/~orc/Code/libc/libc-current/libio/ldouble/

[jack]

here's ioldouble.c with the unused code removed, makes it a bit easier to follow.

[Charles Pegge]

Thank you Jack.

Charles

[jack]

hello Charles
for some time I have been thinking about translating some long double io routines that were placed in the public domain by Stephen L Moshier the author of the cephes math library, I can't find the file anymore by using Google but the name of the file was ioldouble.c, anyway here's the translation so far, not complete and not perfect.

here's the C file and the FreeBASIC translation

[Charles Pegge]

Thanks Jack,

Long quad looks interesting. This could be useful in physics and anything requiring super-high precision in the real world. I think this would also translate quite easily into Oxygen with few alterations to the original C.

I noticed there are are many short integers involved. I wonder if we could take advantage of 32 bits and considerably improve performance.

Charles

[jack]

using 32-bit is likely to speed things up, but would take considerable work.
I still need to go over my translation and clean things up.
one line of code that caused me to almost give up was this line in the sub edivm


if( (tdenm * 0xffffL) < tnum )

where tdenm is an unsigned short and tnum is unsigned long
looks inocent enough but I had to cast tdenm to ulong to get it to work properly.

[Charles Pegge]

Oxygen would have trouble with this too since it does not currently cast integer literals to be long or short. Therefore the left side expression is treated as 16 bit and a 16 bit unsigned comparison is made between left side and right side.

What happens if you use a dword variable instead of &hffff literal? You may be able to avoid casting.

Charles

[jack]

Charles, here's the output of the following test program with the cast to ulong in place


dim as ushort ip(NI-1)
dim as string s
dim as integer i
? "Pi = ";etoasc(@epi(0), 20)
ediv(@eten(0),@epi(0),@ip(0))
? "Pi/10 = ";etoasc(@ip(0), 20)
? "1.0 = ";etoasc(@eone(0), 20)
emul(@epi(0), @epi(0), @ip(0))
? "Pi^2 = ";etoasc(@ip(0), 20)
for i=4 to 20
? etoasc(@epi(0), i)
next
Print "Press RETURN to end ";
sleep



Pi = 3.14159265358979323846e0
Pi/10 = 3.14159265358979323846e-1
1.0 = 1.00000000000000000000e0
Pi^2 = 9.86960440108935861883e0
3.1416e0
3.14159e0
3.141593e0
3.1415927e0
3.14159265e0
3.141592654e0
3.1415926536e0
3.14159265359e0
3.141592653590e0
3.1415926535898e0
3.14159265358979e0
3.141592653589793e0
3.1415926535897932e0
3.14159265358979324e0
3.141592653589793238e0
3.1415926535897932385e0
3.14159265358979323846e0
Press RETURN to end

if I take out the cast statement afore mentioned you get this (using a long variable instead of &hffffL makes no difference)


Pi = 2.68427263874999999985e15
Pi/10 = 2.19895614566399999988e19
1.0 = 1.00000000000000000000e0
Pi^2 = 1.16411769059045078135e6
2.6843e15
2.68427e15
2.684273e15
2.6842726e15
2.68427264e15
2.684272639e15
2.6842726387e15
2.68427263875e15
2.684272638750e15
2.6842726387500e15
2.68427263875000e15
2.684272638750000e15
2.6842726387500000e15
2.68427263875000000e15
2.684272638750000000e15
2.6842726387499999999e15
2.68427263874999999985e15
Press RETURN to end

the L suffix does not appear to be important but the cast is.
it took me quite some time to find the cause of the problem.

if I change the order of the operands to larger first then it works without casting
that is, instead of


if (cast(ulong, tdenm) * &hffffL) < tnum then

to


if tnum >=(tdenm * &hffffL) then

it works

[Charles Pegge]

Jack,

Definitely a bug in FreeBasic since the Left side expression should be upgraded from 16 bits to 32 bits when the long variable is encountered.


I tested this situation with Oxygen and found to my satisfaction that it treated the condition correctly.

However I soon lost my complacency when I realised that the type casting is determined only by the left side. Both Left and right sides should be considered to ensure that the higher type is used on both sides of the comparator.

Then there is the issue of comparing Longs with Dwords...

Anyway I'm grateful that you have raised the issue.

Charles

[jack]

I updated the translation, added string and number conversions, it's complete enough to update my real10 package.

updated file is above.