Computations with 80-bit-FPU-precision

Share your advanced PureBasic knowledge/code with the community.
Helle
Enthusiast
Enthusiast
Posts: 178
Joined: Wed Apr 12, 2006 7:59 pm
Location: Germany
Contact:

Computations with 80-bit-FPU-precision

Post by Helle »

Here a little collection of procedures for 80-bit-computations:

Code: Select all

;File "Procfpu80.pb" for including
;"Helle Klaus Helbing, 10.09.2006, PB4.0

Procedure FPUtoDec(Dummyt1)
;-------- Convert a 80-bit-FPU-value to a string with decimal-digits 
 Convert$=""
  !mov [v_Nachkomma],22
  !mov esi,[v_Memory] 

;-------- Clear Memory 128 Byte
  !mov ecx,32
  !xor eax,eax
  !mov ebx,eax
CLMEM:  
  !mov [esi+ebx],eax
  !add ebx,4
  !dec ecx
  !jnz l_clmem 
  
  !mov edi,[esp+4]      ;Read the value for conversion
 
;-------- Signum and Exponent 
  !mov cx,[esi+edi+8]   ;Signum and Exponent            
  !and cx,cx            ;is MSB (Bit15) setting?
  !jns l_w00            ;no, signum is positiv
  !push ecx 
 Convert$=Convert$+"-"
  !pop ecx
  !inc [v_Nachkomma]
  !and cx,7fffh         ;Signum-Bit clear
W00:  
  !sub cx,16383         ;Sub Bias
  !cmp cx,64
  !jl l_egut
  !cmp cx,-64
  !jg l_egut  

;------------------------------------------------
MessageRequester("Error!","Out of range for display!")

End
;------------------------------------------------

;-------- Read the 64-Bit-Mantisse  
EGUT:
  !mov ebx,[esi+edi]
  !mov [v_ML],ebx       ;Mantisse Low-Part
  !mov eax,[esi+edi+4]
  !mov [v_MH],eax       ;Mantisse High-Part

;-------- Test,if Exponent = -16383 (=0 from FPU) 
  !cmp cx,-16383        ;Exponent
  !jne l_w01
  !or eax,eax
  !jnz l_w01  
  !cmp [v_ML],0
  !jnz l_w01
 Convert$="0.0"         ;Result is Null
  !jmp l_w9             ;End
;------------------------------------------------

W01: 
  !and cx,cx
  !jns l_wnn            ;Exponent not negativ

;-------- Exponent <0 
  !pushad 
Convert$=Convert$+"0."
  !popad 
  !cmp cx,-1
  !jne l_w06
  !mov byte[esi+65],5
  !jmp l_w03
W06:
  !neg cx

;-------- For Digit-Counter 
  !push ecx
  !mov edi,ecx
  !bsr edx,ecx
  !mov ecx,edx
  !shr edi,cl
  !dec cl
  !shl edi,cl
  !pop ecx
  !dec edi

  !add [v_Nachkomma],edi 

;------------------------ 
  !dec cx
  !xor edx,edx
W07:
  !shr edx,1
  !shr ebx,1
  !jnc l_w010 
  !or edi,80000000h
W010:
  !shr eax,1
  !jnc l_w09 
  !or ebx,80000000h
W09:
  !dec cx
  !jnz l_w07

  !mov [v_MH],eax
  !mov [v_ML],ebx
  !mov [v_MLL],edx

  !jmp l_w10

;-------- Exponent >=0 
WNN:  
  !xor edx,edx
  !mov edi,edx
W04:
  !shl edi,1
  !shl edx,1
  !adc edi,0  
  !shl eax,1
  !adc edx,0
W05:  
  !shl ebx,1
  !adc eax,0
  !dec cx
  !jns l_w04

  !mov [v_MH],eax
  !mov [v_ML],ebx

  !mov dword[v_Vorkomma],edx
  !mov dword[v_Vorkomma+4],edi
 
  !and eax,eax
  !jns l_w02
  !mov byte[esi+65],5
W02: 
Convert$=Convert$+StrU(Vorkomma,4)+"." 

W03:   
 X4=Len(Convert$)
  !mov eax,[v_X4]
  !sub [v_Nachkomma],eax
W10:
  !mov byte[esi+1],5    ;Value for start (0.5 for 2^-1)
  
;-------- Conversion mantisses to decimal with unpacked BCD´s
  !mov edx,1
W1:
  !xor ah,ah
  !mov edi,edx
W2:
  !mov ecx,4
  !mov al,[esi+edi]
  !mov bl,al
  !add al,ah 
  !xor ah,ah
W3:  
  !add al,bl
  !aaa
  !dec ecx
  !jnz l_w3
  !mov [esi+edi+1],al
  !dec edi
  !jnz l_w2
  !mov [esi+edi+1],ah   ;write the (2^-n)-value
;-------- Test, if add to value 
  !mov eax,[v_MH]       ;Mantisse High-Part
  !shl eax,1

  !mov ebx,[v_ML]
  !shl ebx,1
  !adc eax,0

  !mov ecx,[v_MLL]
  !shl ecx,1
  !adc ebx,0
 
  !mov [v_MLL],ecx
  !mov [v_ML],ebx
  !mov [v_MH],eax 
 
  !and eax,eax
  !jns l_w5

  !xor ah,ah
  !mov edi,63
W4:  
  !mov bl,[esi+edi]
  !mov al,[esi+edi+64]
  !add al,bl 
  !aaa
  !mov [esi+edi+64],al
  !add [esi+edi+63],ah
  !xor ah,ah
  !dec edi
  !jnz l_w4
W5:
  !inc edx
  !cmp edx,64
  !jb l_w1
   
;-------- Write digits for display to string
AUSGABE: 
  !mov edi,1
W8:
  !mov al,[esi+edi+64] 
  !mov [v_X1],al
 Convert$=Convert$+Str(X1)
  !inc edi
  !cmp edi,[v_Nachkomma]
  !jb l_w8

W9:   
EndProcedure 
;-------- 

;------------------ "Math-Library" -------------- 
;-------- Addition ---------- D3=D1+D2 ----------
Procedure Addt(Dummyt1,Dummyt2,Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !mov edi,[esp+4]
  !fld tword[esi+edi]            
  !mov edi,[esp+8]  
  !fld tword[esi+edi]  
  !faddp st1,st0
  !mov edi,[esp+12]
  !fstp tword[esi+edi]  ;Result
EndProcedure 
;----------------------------

;-------- Cosinus ----------- D3=Cos(D1(rad)) ---
Procedure Cost(Dummyt1,Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !mov edi,[esp+4]  
  !fld tword[esi+edi]  
  !fcos
  !mov edi,[esp+8]  
  !fstp tword[esi+edi]  ;Result
EndProcedure 
;----------------------------

;-------- Cotangens --------- D3=Cot(D1(rad)) ---
Procedure Cott(Dummyt1,Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !mov edi,[esp+4]  
  !fld tword[esi+edi]  
  !fptan
  !fdivrp st1,st0
  !mov edi,[esp+8] 
  !fstp tword[esi+edi]  ;Result
EndProcedure 
;----------------------------

;-------- Division ---------- D3=D1/D2 ----------
Procedure Divt(Dummyt1,Dummyt2,Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !mov edi,[esp+4]
  !fld tword[esi+edi]            
  !mov edi,[esp+8]  
  !fld tword[esi+edi]  
  !fdivp st1,st0
  !mov edi,[esp+12]
  !fstp tword[esi+edi]  ;Result
EndProcedure 
;----------------------------

;-------- Ln2 --------------- D3=Ln2 ------------
Procedure Ln2t(Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !fldln2 
  !mov edi,[esp+4] 
  !fstp tword[esi+edi]  ;Result
 
EndProcedure
;----------------------------

;-------- Log2 -------------- D3=Log2 -----------
Procedure Log2t(Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !fldlg2 
  !mov edi,[esp+4] 
  !fstp tword[esi+edi]  ;Result
 
EndProcedure
;----------------------------

;-------- Multiplikation ---- D3=D1*D2 ----------
Procedure Mult(Dummyt1,Dummyt2,Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !mov edi,[esp+4] 
  !fld tword[esi+edi]            
  !mov edi,[esp+8]   
  !fld tword[esi+edi]  
  !fmulp st1,st0
  !mov edi,[esp+12]
  !fstp tword[esi+edi]  ;Result
EndProcedure 
;----------------------------

;-------- Pi ---------------- D3=Pi -------------
Procedure Pit(Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !fldpi
  !mov edi,[esp+4] 
  !fstp tword[esi+edi]  ;Result
EndProcedure 
;----------------------------

;-------- Pow --------------- D3=D1^D2 ----------
Procedure Powt(Dummyt1,Dummyt2,Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !mov edi,[esp+8]
  !fld tword[esi+edi]            
  !mov edi,[esp+4]  
  !fld tword[esi+edi]
  !fyl2x
  !fld st0
  !frndint
  !fsub st1,st0
  !fxch st1
  !f2xm1
  !fld1
  !faddp st1,st0
  !fscale
  !fstp st1 
  !mov edi,[esp+12]
  !fstp tword[esi+edi]  ;Result
EndProcedure 
;----------------------------

;-------- Quadratwurzel ----- D3=D1^0.5 ---------
Procedure Sqrt(Dummyt1,Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !mov edi,[esp+4]  
  !fld tword[esi+edi]  
  !fsqrt
  !mov edi,[esp+8] 
  !fstp tword[esi+edi]  ;Result
EndProcedure 
;----------------------------

;-------- Sinus ------------- D3=Sin(D1(rad)) ---
Procedure Sint(Dummyt1,Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !mov edi,[esp+4]  
  !fld tword[esi+edi]  
  !fsin
  !mov edi,[esp+8] 
  !fstp tword[esi+edi]  ;Result
EndProcedure 
;----------------------------

;-------- Subtraktion ------- D3=D1-D2 ----------
Procedure Subt(Dummyt1,Dummyt2,Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !mov edi,[esp+4]
  !fld tword[esi+edi]            
  !mov edi,[esp+8]  
  !fld tword[esi+edi]  
  !fsubp st1,st0
  !mov edi,[esp+12]
  !fstp tword[esi+edi]  ;Result
EndProcedure 
;----------------------------

;-------- Tangens ----------- D3=Tan(D1(rad)) ---
Procedure Tant(Dummyt1,Dummyt3)
  !mov esi,[v_Memory] 
  !fninit                       
  !mov edi,[esp+4]  
  !fld tword[esi+edi]  
  !fptan
  !fdivp st1,st0
  !mov edi,[esp+8] 
  !fstp tword[esi+edi]  ;Result
;----------------------------
EndProcedure
For tests:

Code: Select all

;File "testfpu80.pb"
;"Helle Klaus Helbing, 10.09.2006, PB4.0

;-------- For the Procedures in Procfpu80
Global X1.b
Global MH.l             ;Mantisse High
Global ML.l             ;Mantisse Low
Global MLL.l            ;Mantisse Low Low (for Exponent <0)
Global Memory.l
Global Dummyt1.l        ;Parameters for 80-Bit-Operations
Global Dummyt2.l
Global Dummyt3.l
Global X4.l
Global Nachkomma.l      ;limit digit-counter
Global Vorkomma.q
Global Convert$

;-------- For this Tests
Global ZD1.l=128        ;Pointer in Memory, For Alignment 16 Bytes (need
Global ZD2.l=144        ;10 Bytes)
Global ZD3.l=160        


XIncludeFile "Procfpu80.pb"

;-------- Programm-Beginn
Memory = AllocateMemory(1000)     ;128 Bytes for Procedure FPUtoDec + variables
  !mov esi,[v_Memory]

;-------- Initialisation Variables with values (if needed)  
  !fninit

  !mov edi,[v_ZD1]      ;pointer ZD1 to memory
  !fld [D1]             ;read D1
  !fstp tword[esi+edi]  ;write the value of variable D1 to memory

  !mov edi,[v_ZD2]
  !fld [D2]
  !fstp tword[esi+edi]

;-------- Tests
;Addt(ZD1,ZD2,ZD3)
;Cost(ZD1,ZD3)
;Cott(ZD1,ZD3)
;Divt(ZD1,ZD2,ZD3)
;Ln2t(ZD3)
;Log2t(ZD3)
;Mult(ZD1,ZD2,ZD3)
;Pit(ZD3)
Powt(ZD1,ZD2,ZD3)
;Sint(ZD1,ZD3)
;Sqrt(ZD1,ZD3)
;Subt(ZD1,ZD2,ZD3)
;Tant(ZD1,ZD3)

;-------- Display result
FPUtoDec(ZD3)           ;Convert result to String (Convert$) for display
MessageRequester("80-bit computation", Convert$)


End

;-------- write a value to the 80-Bit-Variables
!section '.data' Data readable writeable

!D1 dt 123.456789
!D2 dt 6.789
An example:

Code: Select all

;Example for 80-bit-computation : Earth-Surface in m²
;"Helle Klaus Helbing, 10.09.2006, PB4.0

;-------- For the Procedures in Procfpu80
Global X1.b
Global MH.l
Global ML.l
Global MLL.l
Global Memory.l
Global Dummyt1.l
Global Dummyt2.l
Global Dummyt3.l
Global X4.l
Global Nachkomma.l
Global Vorkomma.q
Global Convert$

;-------- For this Example
Global ZPI.l=128
Global ZV4.l=144
Global ZE1.l=160        
Global ZE2.l=176
Global ZRE.l=192
Global ZSE.l=208

XIncludeFile "Procfpu80.pb"

;-------- 
Memory = AllocateMemory(1000)
  !mov esi,[v_Memory]

;-------- Example : Earth-Surface (4*Pi*R*R) in m²
  !fninit

  !mov edi,[v_ZV4]      ;value 4.0
  !fld [V4]
  !fstp tword[esi+edi]

  !mov edi,[v_ZRE]      ;Earth-Radius
  !fld [RE]
  !fstp tword[esi+edi]

;--------
Pit(ZPI)                ;loading the value of Pi
Mult(ZV4,ZPI,ZE1)       ;write 4*Pi to E1
Mult(ZRE,ZRE,ZE2)       ;write Re*Re to E2
Mult(ZE1,ZE2,ZSE)       ;mul E1,E2 is result 

FPUtoDec(ZSE)
MessageRequester("80-bit computation","Earth-Surface in m² :  " + Convert$)

End

;-------- If variables have a start-value
!section '.data' Data readable writeable

!V4 dt 4.0              ;V4 = value 4.0
!RE dt 6371007.176      ;Earth-Radius in meters (GRS 80-Ellipsoid)
Greetings
Helle
jack
Addict
Addict
Posts: 1358
Joined: Fri Apr 25, 2003 11:10 pm

Post by jack »

some suggestions.
1. define a structure of 16 bytes and use that instead of AllocateMemory, much friendlier.
2. !finit is not the fastest way to set the fpu to extended mode, do something like this:

Code: Select all

;define at start of main
Global ctrlwrd.w = 4927;&B0001001100111111

;then instead of finit

! fldcw [v_ctrlwrd] ; this guarantees extended precision
of course you will have to make sure you leave the fpu stack empty, otherwise you get some unexpected results.
3. instead of using esp to access the arguments to a function use p.v_Dummy1 for example or p.p_Dummy1 if it's a pointer.
4. using a structure to hold the extended values, return the address in eax, that way you can nest function calls.

example

Code: Select all

Structure ext
  t.l[4]
EndStructure

Global ctrlwrd.w = 4927;&B0001001100111111 

;-------- Addition ---------- D3=D1+D2 ---------- 
Procedure Addt(*Dummyt1.ext,*Dummyt2.ext,*Dummyt3.ext) 
  !fldcw [v_ctrlwrd] ; this guarantees extended precision
  !mov edi,[p.p_Dummyt1]
  !fld tword[edi]            
  !mov edi,[p.p_Dummyt2] 
  !fld tword[edi]  
  !faddp st1,st0 
  !mov eax,[p.p_Dummyt3] ;use eax to hold pointer to result
  !fstp tword[eax]  ;Result 
  ProcedureReturn ;return eax
EndProcedure 

Define.ext x,y,z

;example
Addt(x,Addt(x,y,z),z)
SoulReaper
Enthusiast
Enthusiast
Posts: 372
Joined: Sun Apr 03, 2005 2:14 am
Location: England

Post by SoulReaper »

Still Blown away @helle Nice assembly code/functions :) :lol: :wink:
Post Reply