ThinBasic extended numbers in FreeBASIC SDK.

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

Code:

---

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

Code:

---

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

Code:

---

' 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

---

Quote:

---

Originally Posted by Charles Pegge

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

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

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

Charles

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

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

Hi Joshy,

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

Also a few more tests

Charles

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/li...libio/ldouble/

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

Thank you Jack.

Charles

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

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

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

Code:

---

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.

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

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

Code:

---

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

---

Quote:

---

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)

Quote:

---

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.
[edit]
if I change the order of the operands to larger first then it works without casting
that is, instead of

Code:

---

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

---

to

Code:

---

if  tnum >=(tdenm * &hffffL) then

---

it works
[/edit]

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

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

updated file is above.