What it does is to calculate in a hardware-independent manner a range of values relating to the Floating Point processor.
Here is an example: Suppose you are writing an iterative looping routine (perhaps to calculate the roots of a mathematical function) and you need a criterion for deciding when the solution has been determined to the maximum possible accuracy i.e. least possible error.
Perhaps you have arbitrarily decided to stop looping when the difference between successive solution approximations is say 1E-9 (10 to the power of -9). But the hardware may be capable of better accuracy. This is where my program could be useful. The value of 'epsilon' that it generates will be a better criterion than the arbitrarily chosen 1E-9.
There is no need to actually know the value of epsilon, merely assign it to a double precision variable and compare it with the loop criterion.
Hopefully this program will work correctly on all platforms supported by PureBasic. The only part which is platform specific is the procedure StrE() and this is used only to display the values of some of the FP attributes. But as I said in the previous paragraph, you do not need to know the values, just use them.
The program is translated from C++ from the third edition of the Book "Numerical Recipes" within the chapter "Less-Numerical Algorithms".
Note that some of the program's results relate to the radix of the floating point unit being used. For example 'Exponent Digits' which has a value of 11 for a typical Pentium PC implies an exponent range of about 2^11 (= 2048) for this binary FP unit, not 10^11. Incidently, some exponent values are often reserved to indicate special cases such as infinity, and therefore not used for normal numbers.
Code: Select all
; Machar
; Determines parameters affecting floating-point arithmetic
; These values are machine specific
; From: Numerical Recipes - Less-Numerical Algorithms
EnableExplicit
Structure MacharType
; Double fields are: eps, epsneg, xmin, xmax
ibeta.i ; Radix in which numbers are represented (2, 16, 10)
it.i ; Mantissa digits (in base ibeta)
machep.i; Exponent of most negative power of ibeta
; that when added to 1.0 results in <>1.0
eps.d ; Floating point precision = ibeta^machep
negep.i ; Exponent of most negative power of ibeta
; that when subtracted from 1.0 results in <>1.0
epsneg.d; Floating point precision = ibeta^negep
iexp.i ; Exponent bits (including sign or bias)
minexp.i; Exponent of most negative power of ibeta
; consistent with no leading mantissa zeros
xmin.d ; Magnitude of smallest usable FP value (ibeta^minexp)
maxexp.i; Exponent of most negative power of ibeta that overflows
xmax.d ; Magnitude of largest usable FP value
; (1-epsneg)*ibeta^maxexp
irnd.i ; How rounding and underflow are handled (0..5)
; 0 or 3 Result truncated, not rounded (bad)
; 1 or 4 Result rounded, but not to IEEE standard
; 2 or 5 Result rounded to IEEE standard (good)
; 0,1,2 Value<xmin underflows abruptly to zero
; 3,4,5 Value<xmin underflows to IEEE standard
; (exponent frozen, mantissa has leading zeroes)
ngrd.i ; Number of guard digits when truncating a mantissa product
EndStructure
Procedure.s StrE(fp.d, length=16) ; AKJ 17-Nov-08
; Convert FP number to fixed length scientific format
; The length is the number of significant mantissa digits
Protected mnt.d, exp ; mantissa, unbiased exponent [radix 10]
Protected fp$ ; The result
Protected test.d ; Test whether F2XM1 argument is in range
If length<0 Or length>20: length = 20: EndIf ; Sanity check
; Special case
If fp=0.0
ProcedureReturn "+"+StrD(0.0, length-1)+" E+0"
EndIf
; Extract sign
If fp<0.0: fp = -fp: fp$ = "-": Else: fp$ = "+": EndIf
; Calculate radix 10 mantissa and exponent
EnableASM
FNINIT ; Round to nearest, all exceptions masked, 64-bit precision
FLD qword [p.v_fp]
FXTRACT ; Exponent (radix 2) -> ST1, mantissa -> ST0
FXCH ; Exponent (radix 2) -> ST0, mantissa -> ST1
FLD ST0 ; Push ST0 (duplicate exp2)
FLDLG2 ; Push Log10(2) = 0.301
FMULP ; Approx exponent (radix 10) -> ST0
FRNDINT ; Round exponent to integer
; Must round (not truncate) for F2XM1 to work
FIST exp ; Approx exponent (radix 10), may be wrong by 1
; Calculate mantissa correction = 2^(exp2-exp*log2(10))
FLDL2T ; Push Log2(10) = 3.322
FMULP ; exp2*log2(10)
FSUBP ; exp2-exp*log2(10) Hopefully in the range -1.0..+1.0
FST test ; Magnitude upper limit is 0.5*Log2(10)=1.661
If Abs(test)<1 ; If already in range -1.0..+1.0
F2XM1 ; (2^x)-1
FLD1 ; +1
FADDP ; Mantissa correction
Else ; Use 2^x = Square_of(2^(x/2))
FLD1 ; 1
FLD ST0 ; Duplicate
FADDP ; 2
FDIVP ; x/2
F2XM1 ; (2^(x/2))-1
FLD1 ; +1
FADDP ; [Square root of mantissa correction]
FLD ST0 ; Duplicate
FMULP ; Square to give mantissa correction
EndIf
FMULP ; Mantissa
FSTP mnt
DisableASM
; Normalise mantissa to be >=1 and <10 (i.e. 1 digit before DP)
; In theory at most one adjustment will ever be needed
While mnt<1.0: mnt * 10.0: exp - 1: Wend
While mnt>=10.0: mnt / 10.0: exp + 1: Wend
; Build result
fp$ + StrD(mnt, length-1)+" E"
If exp>=0: fp$ + "+": EndIf
ProcedureReturn fp$ + Str(exp)
EndProcedure
Procedure Machar()
; Int i,itemp,iz,j,k,mx,nxres
; Doub a,b,beta,betah,betain,one,t,temp,temp1,tempa,two,y,z,zero
Protected i, itemp, iz, j, k, mx, nxres
Protected a.d, b.d, beta.d, betah.d, betain.d, t.d
Protected temp.d, temp1.d, tempa.d
Protected zero.d, one.d, two.d, y.d, z.d
Static m.MacharType
one = 1.0
two = one+one
zero = one-one
; Determine ibeta and beta by the method of M.Malcolm
a = one
Repeat
a + a
temp = a+one
temp1 = temp-a
Until temp1-one<>zero
b = one
Repeat
b + b
temp = a+b
itemp = Int(temp-a) ; !!!
Until itemp<>0
m\ibeta = itemp
beta = m\ibeta ; Convert to double
; Determine it
m\it = 0
b = one
Repeat
m\it + 1
b * beta
temp = b+one
temp1 = temp-b
Until temp1-one<>zero
; Determine irnd (adjusted later)
m\irnd = 0
betah = beta/two
temp = a+betah
If temp-a<>zero: m\irnd = 1: EndIf
tempa = a+beta
temp = tempa+betah
If m\irnd=0 And temp-tempa<>zero: m\irnd = 2: EndIf
; Determine negep and epsneg
m\negep = m\it+3
betain = one/beta
a = one
For i=1 To m\negep: a * betain: Next i
b = a
Repeat
temp = one-a
If temp-one<>zero: Break: EndIf
a * beta
m\negep - 1
ForEver
m\negep = -m\negep
m\epsneg = a
; Determine machep and eps
m\machep = -m\it-3
a = b
Repeat
temp = one+a
If temp-one<>zero: Break: EndIf
a * beta
m\machep + 1
ForEver
m\eps = a
; Determine ngrd
m\ngrd = 0
temp = one+m\eps
If m\irnd=0 And temp*one-one<>zero: m\ngrd=1: EndIf
; Determine iexp
i = 0
k = 1
z = betain
t = one+m\eps
nxres = 0
Repeat ; Loop until underflow occurs, then exit
y = z
z = y*y
a = z*one ; Check for underflow
temp = z*t
If a+a=zero Or Abs(z)>=y: Break: EndIf
temp1 = temp*betain
If temp1*beta=z: Break: EndIf
i + 1
k + k
ForEver
If m\ibeta<>10
m\iexp = i+1
mx = k+k
Else ; For decimal FP units only
m\iexp = 2
iz = m\ibeta
While k>=iz
iz * m\ibeta
m\iexp + 1
Wend
mx = iz+iz-1
EndIf
; Determine minexp and xmin
Repeat ; Loop until underflow occurs, then exit
m\xmin = y
y * betain
a = y*one ; Check for underflow
temp = y*t
If a+a<>zero And Abs(y)<m\xmin
k + 1
temp1 = temp*betain
If temp1*beta=y And temp<>y
nxres = 3
m\xmin = y
Break
EndIf
Else
Break
EndIf
ForEver
; Determine maxexp
m\minexp = -k
If mx<=k+k-3 And m\ibeta<>10
mx + mx
m\iexp + 1
EndIf
m\maxexp = mx+m\minexp
m\irnd + nxres ; Adjust irnd to reflect partial underflow
If m\irnd>=2: m\maxexp - 2: EndIf ; Adjust if meets IEEE standard
i = m\maxexp+m\minexp
; Adjust for FP units with implicit leading bit in mantissa
; and those with radix point at extreme right of mantissa
If m\ibeta=2 And i=0: m\maxexp - 1: EndIf ; Was && !i !!!
If i>20: m\maxexp - 1: EndIf
If a<>y: m\maxexp - 2: EndIf
; Determine xmax
m\xmax = one-m\epsneg
If m\xmax*one<>m\xmax: m\xmax=one-beta*m\epsneg: EndIf
m\xmax / (m\xmin*beta*beta*beta)
i = m\maxexp+m\minexp+3
For j=1 To i
If m\ibeta=2
m\xmax + m\xmax
Else
m\xmax * beta
EndIf
Next j
ProcedureReturn @m ; Return structure
EndProcedure
Define *mach.MacharType
*mach = Machar() ; Calculate machine dependent values
With *mach
Debug "Quantity: Calculated Numeric Value"
Debug "Radix: "+Str(\ibeta)
Debug "Mantissa digits: "+Str(\it)
Debug "Rounding style: "+Str(\irnd)
Debug "Guard digits: "+Str(\ngrd)
Debug "Epsilon (precision): "+StrE(\eps)
Debug "Neg epsilon (precision): "+StrE(\epsneg)
Debug "Epsilon power: "+Str(\machep)
Debug "Neg epsilon power: "+Str(\negep)
Debug "Exponent digits: "+Str(\iexp)
Debug "Min exponent: "+Str(\minexp)
Debug "Max exponent: "+Str(\maxexp)
Debug "Minimum: "+StrE(\xmin)
Debug "Maximum: "+StrE(\xmax)
EndWith
End