decfloat.pbi
Code: Select all
; functions that take number of digits instead of number of dwords as optional parameter
; fp2str
; fp2str_exp
; fp2str_fix
; fpeps
; pi_brent_salamin
; the rest take number of dwords as optional parameter
EnableExplicit
CompilerIf Not Defined(number_of_digits,#PB_Constant)
#number_of_digits = 80
CompilerEndIf
#num_digits = #number_of_digits
CompilerIf ((#num_digits / 8)*8)<#num_digits
#num_dwords = #num_digits / 8 +1
CompilerElse
#num_dwords = #num_digits / 8
CompilerEndIf
#bias = 1073741824 ;2 ^ 30
Structure decfloat
sign.l
exponent.l
mantissa.l[#num_dwords+1]
EndStructure
; error definitions
#divz_err = 1 ;divide by zero
#expo_err = 2 ;exponent overflow error
Procedure.s String(n.l, s.s)
Define.s sn
sn=Space(n)
ReplaceString(sn," ",s,#PB_String_InPlace)
ProcedureReturn sn
EndProcedure
CompilerIf Not Defined(maximum_Pi_Dec_precision,#PB_Constant)
#maximum_Pi_Dec_precision = 1024
CompilerEndIf
Declare pi_brent_salamin (*pi_bs.decfloat, digits_in.l= #number_of_digits)
Global.decfloat Pi_dec
pi_brent_salamin (Pi_dec, #maximum_Pi_Dec_precision)
Procedure copydec(*n.decfloat, *result.decfloat)
Define.l i
*result\sign=*n\sign
*result\exponent=*n\exponent
For i=0 To #num_dwords
*result\mantissa[i]=*n\mantissa[i]
Next
EndProcedure
Procedure.s fp2str_exp (*n.decfloat, places.l=#number_of_digits-2)
Define.l i, ex
Define.s v, f, ts
If *n\exponent > 0
ex = (*n\exponent & $7fffffff) - #bias - 1
Else
ex = 0
EndIf
If *n\sign
v = "-"
Else
v = " "
EndIf
ts = Str(*n\mantissa[0])
ts = Trim(ts)
If Len(ts) < 8
ts = ts + String(8 - Len(ts), "0")
EndIf
v = v + Left(ts, 1) + "." + Mid(ts, 2)
For i = 1 To #num_dwords - 1
ts = Str(*n\mantissa[i])
ts = Trim(ts)
If Len(ts) < 8
ts = String(8 - Len(ts), "0") + ts
EndIf
v = v + ts
Next
v = Left(v, places + 3)
f = Trim(Str(Abs(ex)))
;f = String(9 - Len(f), "0") + f
If ex < 0
v = v + "e-"
Else
v = v + "e+"
EndIf
v = v + f
ProcedureReturn v
EndProcedure
;long part of num
Procedure fptrunc (*result.decfloat, *num.decfloat)
Define.decfloat ip
Define.l ex, ex2, j, k
ex = (*num\exponent & $7fffffff) - #bias
If ex < 1
copydec(ip, *result) ;result = ip
ProcedureReturn
EndIf
If ex >= (#num_digits)
copydec(*num, *result) ;result = num
ProcedureReturn
EndIf
ex2 = ex / 8
k = ex2
j = ex % 8
While ex2 > 0
ex2 = ex2 - 1
ip\mantissa[ex2]=*num\mantissa[ex2]
Wend
If j = 1
ip\mantissa[k]= 10000000 * (*num\mantissa[k] / 10000000)
ElseIf j = 2
ip\mantissa[k]= 1000000 * (*num\mantissa[k] / 1000000)
ElseIf j = 3
ip\mantissa[k]= 100000 * (*num\mantissa[k] / 100000)
ElseIf j = 4
ip\mantissa[k]= 10000 * (*num\mantissa[k] / 10000)
ElseIf j = 5
ip\mantissa[k]= 1000 * (*num\mantissa[k] / 1000)
ElseIf j = 6
ip\mantissa[k]= 100 * (*num\mantissa[k] / 100)
ElseIf j = 7
ip\mantissa[k]= 10 * (*num\mantissa[k] / 10)
ElseIf j = 8
ip\mantissa[k]= *num\mantissa[k]
EndIf
ip\exponent = ex + #bias
ip\sign = *num\sign
copydec(ip, *result) ;result = ip
EndProcedure
Procedure.l fptrunc_is_odd (*num.decfloat)
Define.l ex, j, k
ex = (*num\exponent & $7fffffff) - #bias
If ex < 1
ProcedureReturn 0
EndIf
If ex >= (#num_digits)
PrintN( "error in function fptrunc_is_odd")
ProcedureReturn 99999999
EndIf
k = ex / 8
j = ex % 8
If j = 1
ProcedureReturn (*num\mantissa[k] / 10000000) & 1
ElseIf j = 2
ProcedureReturn (*num\mantissa[k] / 1000000) & 1
ElseIf j = 3
ProcedureReturn (*num\mantissa[k] / 100000) & 1
ElseIf j = 4
ProcedureReturn (*num\mantissa[k] / 10000) & 1
ElseIf j = 5
ProcedureReturn (*num\mantissa[k] / 1000) & 1
ElseIf j = 6
ProcedureReturn (*num\mantissa[k] / 100) & 1
ElseIf j = 7
ProcedureReturn (*num\mantissa[k] / 10) & 1
ElseIf j = 8
ProcedureReturn *num\mantissa[k] & 1
EndIf
ProcedureReturn 0
EndProcedure
Procedure.d fp2dbl (*n.decfloat)
Define.l ex
Define.s v, f, ts
If *n\exponent > 0
ex = (*n\exponent & $7fffffff) - #bias - 1
Else
ex = 0
EndIf
If *n\sign
v = "-"
Else
v = " "
EndIf
ts = Trim(Str(*n\mantissa[0]))
If Len(ts) < 8
ts = ts + String(8 - Len(ts), "0")
EndIf
v = v + Left(ts, 1) + "." + Mid(ts, 2)
ts = Trim(Str(*n\mantissa[1]))
If Len(ts) < 8
ts = String(8 - Len(ts), "0") + ts
EndIf
v = v + ts
f = Str(Abs(ex))
f = String(5 - Len(f), "0") + f
If ex < 0
v = v + "e-"
Else
v = v + "e+"
EndIf
v = v + f
ProcedureReturn ValD(v)
EndProcedure
Procedure.s fp2str_fix (*n.decfloat, places.l=#number_of_digits-2)
Define.l i, ex
Define.s v, ts, s
If *n\exponent > 0
ex = (*n\exponent & $7fffffff) - #bias - 1
Else
ex = 0
EndIf
If *n\sign
s = "-"
Else
s = " "
EndIf
ts = Trim(Str(*n\mantissa[0]))
If Len(ts) < 8
ts = ts + String(8 - Len(ts), "0")
EndIf
v = ts
For i = 1 To #num_dwords - 1
ts = Trim(Str(*n\mantissa[i]))
If Len(ts) < 8
ts = String(8 - Len(ts), "0") + ts
EndIf
v = v + ts
Next
If places < #num_digits
v = Left(v, places)
EndIf
If ex = 0
v = Left(v, 1) + "." + Mid(v, 2)
ElseIf ex < 0
v = "0." + String(Abs(ex) - 1, "0") + v
ElseIf ex > 0
v = Left(v, ex + 1) + "." + Mid(v, ex + 2)
EndIf
ProcedureReturn s + v
EndProcedure
Declare fpadd (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
Declare fpsub (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
Declare si2fp (*result.decfloat, m.q, dwords_in.l=#num_dwords)
Procedure.s fp2str (*n.decfloat, places_in.l=#number_of_digits)
Define.l ex, e, sgn, places=places_in
Define.decfloat z
If places>#number_of_digits-1
places=#number_of_digits-2
EndIf
If *n\exponent <> 0
ex = (*n\exponent & $7fffffff) - #bias - 1
Else
ex = 0
EndIf
If *n\exponent=0
ProcedureReturn " 0"
EndIf
e=ex
If e<0
e=-e
EndIf
si2fp (z, 5)
sgn=*n\sign
If (ex > (-5)) And (e < places)
z\exponent = (-(places+e) + #bias + 1)
If sgn=0
fpadd(z, z, *n)
Else
fpsub(z, *n, z)
EndIf
ProcedureReturn fp2str_fix(z, places)
Else
z\exponent = (-(places+e+1) + #bias + 1)
If sgn=0
fpadd(z, z, *n)
Else
fpsub(z, *n, z)
EndIf
ProcedureReturn fp2str_exp(z, places)
EndIf
EndProcedure
; in case the above doesn't work right use this instead
; Procedure.s fp2str (*n.decfloat, places.l=#number_of_digits-2)
; Define.l ex
; If *n\exponent <> 0
; ex = (*n\exponent & $7fffffff) - #bias - 1
; Else
; ex = 0
; EndIf
; If (Abs(ex) < places) And (places<(#num_digits-1))
; ProcedureReturn fp2str_fix(*n, places)
; Else
; ProcedureReturn fp2str_exp(*n, places)
; EndIf
; EndProcedure
Procedure.s fp2strT (*n.decfloat, places_in.l=#number_of_digits)
Define.l ex, e, sgn, places=places_in, k
Define.decfloat z
Define.s s, s2
If places>#number_of_digits-1
places=#number_of_digits-2
EndIf
If *n\exponent <> 0
ex = (*n\exponent & $7fffffff) - #bias - 1
Else
ex = 0
EndIf
If *n\exponent=0
ProcedureReturn " 0"
EndIf
e=ex
If e<0
e=-e
EndIf
si2fp (z, 5)
sgn=*n\sign
If (ex > (-5)) And (e < places)
z\exponent = (-(places+e) + #bias + 1)
If sgn=0
fpadd(z, z, *n)
Else
fpsub(z, *n, z)
EndIf
s=fp2str_fix(z, places)
s=RTrim(s, "0")
If Right(s,1)="."
s=Left(s,Len(s)-1)
EndIf
ProcedureReturn s
Else
z\exponent = (-(places+e+1) + #bias + 1)
If sgn=0
fpadd(z, z, *n)
Else
fpsub(z, *n, z)
EndIf
s=fp2str_exp(z, places)
k=FindString(s, "e")
s2=Left(s, k-1)
s2=RTrim(s2, "0")
If Right(s2,1)="."
s2=Left(s2,Len(s2)-1)
EndIf
s=s2+Right(s,Len(s)-k+1)
ProcedureReturn s
EndIf
EndProcedure
Procedure str2fp (*result.decfloat, value_in.s)
Define.l j, s, d, e, ep, ulng
Define.l ex, es, i, f, fp, fln
Define.s c, f1, f2, f3, ts, value
Define.decfloat n
value=value_in
j = 1
s = 1
d = 0
e = 0
ep = 0
ex = 0
es = 1
i = 0
f = 0
fp = 0
f1 = ""
f2 = ""
f3 = ""
value = UCase(value)
fln = Len(value)
While j <= fln
c = Mid(value, j, 1)
If ep = 1
If c = " "
j = j + 1
Goto skip_while
EndIf
If c = "-"
es = -es
c = ""
EndIf
If c = "+"
j = j + 1
Goto skip_while
EndIf
If (c = "0") And (f3 = "")
j = j + 1
Goto skip_while
EndIf
If (c > "/") And (c < ":") ;c is digit between 0 and 9
f3 = f3 + c
ex = 10 * ex + (Asc(c) - 48)
j = j + 1
Goto skip_while
EndIf
EndIf
If c = " "
j = j + 1
Goto skip_while
EndIf
If c = "-"
s = -s
j = j + 1
Goto skip_while
EndIf
If c = "+"
j = j + 1
Goto skip_while
EndIf
If c = "."
If d = 1
j = j + 1
Goto skip_while
EndIf
d = 1
EndIf
If (c > "/") And (c < ":") ;c is digit between 0 and 9
If ((c = "0") And (i = 0))
If d = 0
j = j + 1
Goto skip_while
EndIf
If (d = 1) And (f = 0)
e = e - 1
j = j + 1
Goto skip_while
EndIf
EndIf
If d = 0
f1 = f1 + c
i = i + 1
Else
If (c > "0")
fp = 1
EndIf
f2 = f2 + c
f = f + 1
EndIf
EndIf
If c = "E" Or c = "D"
ep = 1
EndIf
j = j + 1
skip_while:
Wend
If fp = 0
f = 0
f2 = ""
EndIf
If s = -1
s = $8000
Else
s = 0
EndIf
n\sign = s
ex = es * ex - 1 + i + e
f1 = f1 + f2
f1 = Mid(f1, 1, 1) + Right(f1, Len(f1) - 1)
fln = Len(f1)
If Len(f1) > (#num_digits + 1 + 8)
f1 = Mid(f1, 1, (#num_digits + 1 + 8))
EndIf
While Len(f1) < (#num_digits + 1 + 8)
f1 = f1 + "0"
Wend
j = 1
For i = 0 To #num_dwords
ts = Mid(f1, j, 8)
ulng = Val(ts)
n\mantissa[i]= ulng
If ulng <> 0
fp = 1
EndIf
j = j + 8
Next
If fp
n\exponent = (ex + #bias + 1)
Else
n\exponent = 0
EndIf
copydec(n, *result) ;result = n
EndProcedure
Procedure si2fp (*result.decfloat, m.q, dwords_in.l=#num_dwords)
Define.l dwords, i, ndw
Define.q n
Define.decfloat fac1
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
n = Abs(m)
ndw = n
If n > 9999999999999999
str2fp (fac1, Str(m))
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
For i = 1 To dwords
fac1\mantissa[i] = 0
Next
If m = 0
fac1\exponent = 0
fac1\sign = 0
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
fac1\exponent = #bias
If n < 100000000
If n < 10
fac1\mantissa[0] = n * 10000000
fac1\exponent = fac1\exponent + 1
ElseIf n < 100
fac1\mantissa[0] = n * 1000000
fac1\exponent = fac1\exponent + 2
ElseIf n < 1000
fac1\mantissa[0] = n * 100000
fac1\exponent = fac1\exponent + 3
ElseIf n < 10000
fac1\mantissa[0] = n * 10000
fac1\exponent = fac1\exponent + 4
ElseIf n < 100000
fac1\mantissa[0] = n * 1000
fac1\exponent = fac1\exponent + 5
ElseIf n < 1000000
fac1\mantissa[0] = n * 100
fac1\exponent = fac1\exponent + 6
ElseIf n < 10000000
fac1\mantissa[0] = n * 10
fac1\exponent = fac1\exponent + 7
ElseIf n < 100000000
fac1\mantissa[0] = ndw
fac1\exponent = fac1\exponent + 8
EndIf
EndIf
If n > 99999999
fac1\exponent = fac1\exponent + 8
If n < 1000000000
fac1\mantissa[0] = n / 10
fac1\mantissa[1] = (n % 10) * 10000000
fac1\exponent = fac1\exponent + 1
ElseIf n < 100000000000
fac1\mantissa[0] = n / 100
fac1\mantissa[1] = (n % 100) * 1000000
fac1\exponent = fac1\exponent + 2
ElseIf n < 1000000000000
fac1\mantissa[0] = n / 1000
fac1\mantissa[1] = (n % 1000) * 100000
fac1\exponent = fac1\exponent + 3
ElseIf n < 10000000000000
fac1\mantissa[0] = n / 10000
fac1\mantissa[1] = (n % 10000) * 10000
fac1\exponent = fac1\exponent + 4
ElseIf n < 100000000000000
fac1\mantissa[0] = n / 100000
fac1\mantissa[1] = (n % 100000) * 1000
fac1\exponent = fac1\exponent + 5
ElseIf n < 1000000000000000
fac1\mantissa[0] = n / 1000000
fac1\mantissa[1] = (n % 1000000) * 100
fac1\exponent = fac1\exponent + 6
ElseIf n < 10000000000000000
fac1\mantissa[0] = n / 10000000
fac1\mantissa[1] = (n % 10000000) * 10
fac1\exponent = fac1\exponent + 7
ElseIf n < 100000000000000000
fac1\mantissa[0] = n / 100000000
fac1\mantissa[1] = n % 100000000
fac1\exponent = fac1\exponent + 8
EndIf
EndIf
If m < 0
fac1\sign = $8000
Else
fac1\sign = 0
EndIf
copydec(fac1, *result) ;result = fac1
EndProcedure
Procedure rshift_1 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = dwords To 1 Step -1
v1 = *mantissa\mantissa[i] / 10
v2 = *mantissa\mantissa[i - 1] % 10
v2 = v2 * 10000000 + v1
*mantissa\mantissa[i] = v2
Next
*mantissa\mantissa[0] = *mantissa\mantissa[0] / 10
EndProcedure
Procedure lshift_1 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = 0 To dwords - 1
v1 = *mantissa\mantissa[i] % 10000000
v2 = *mantissa\mantissa[i + 1] / 10000000
*mantissa\mantissa[i] = v1 * 10 + v2
*mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 10000000
Next
*mantissa\mantissa[dwords] = 10 * (*mantissa\mantissa[dwords] % 10000000)
EndProcedure
Procedure rshift_2 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = dwords To 1 Step -1
v1 = *mantissa\mantissa[i] / 100
v2 = *mantissa\mantissa[i - 1] % 100
v2 = v2 * 1000000 + v1
*mantissa\mantissa[i] = v2
Next
*mantissa\mantissa[0] = *mantissa\mantissa[0] / 100
EndProcedure
Procedure lshift_2 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = 0 To dwords - 1
v1 = *mantissa\mantissa[i] % 1000000
v2 = *mantissa\mantissa[i + 1] / 1000000
*mantissa\mantissa[i] = v1 * 100 + v2
*mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 1000000
Next
*mantissa\mantissa[dwords] = 100 * (*mantissa\mantissa[dwords] % 1000000)
EndProcedure
Procedure rshift_3 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = dwords To 1 Step -1
v1 = *mantissa\mantissa[i] / 1000
v2 = *mantissa\mantissa[i - 1] % 1000
v2 = v2 * 100000 + v1
*mantissa\mantissa[i] = v2
Next
*mantissa\mantissa[0] = *mantissa\mantissa[0] / 1000
EndProcedure
Procedure lshift_3 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = 0 To dwords - 1
v1 = *mantissa\mantissa[i] % 100000
v2 = *mantissa\mantissa[i + 1] / 100000
*mantissa\mantissa[i] = v1 * 1000 + v2
*mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 100000
Next
*mantissa\mantissa[dwords] = 1000 * (*mantissa\mantissa[dwords] % 100000)
EndProcedure
Procedure rshift_4 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = dwords To 1 Step -1
v1 = *mantissa\mantissa[i] / 10000
v2 = *mantissa\mantissa[i - 1] % 10000
v2 = v2 * 10000 + v1
*mantissa\mantissa[i] = v2
Next
*mantissa\mantissa[0] = *mantissa\mantissa[0] / 10000
EndProcedure
Procedure lshift_4 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = 0 To dwords - 1
v1 = *mantissa\mantissa[i] % 10000
v2 = *mantissa\mantissa[i + 1] / 10000
*mantissa\mantissa[i] = v1 * 10000 + v2
*mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 10000
Next
*mantissa\mantissa[dwords] = 10000 * (*mantissa\mantissa[dwords] % 10000)
EndProcedure
Procedure rshift_5 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = dwords To 1 Step -1
v1 = *mantissa\mantissa[i] / 100000
v2 = *mantissa\mantissa[i - 1] % 100000
v2 = v2 * 1000 + v1
*mantissa\mantissa[i] = v2
Next
*mantissa\mantissa[0] = *mantissa\mantissa[0] / 100000
EndProcedure
Procedure lshift_5 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = 0 To dwords - 1
v1 = *mantissa\mantissa[i] % 1000
v2 = *mantissa\mantissa[i + 1] / 1000
*mantissa\mantissa[i] = v1 * 100000 + v2
*mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 1000
Next
*mantissa\mantissa[dwords] = 100000 * (*mantissa\mantissa[dwords] % 1000)
EndProcedure
Procedure rshift_6 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = dwords To 1 Step -1
v1 = *mantissa\mantissa[i] / 1000000
v2 = *mantissa\mantissa[i - 1] % 1000000
v2 = v2 * 100 + v1
*mantissa\mantissa[i] = v2
Next
*mantissa\mantissa[0] = *mantissa\mantissa[0] / 1000000
EndProcedure
Procedure lshift_6 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = 0 To dwords - 1
v1 = *mantissa\mantissa[i] % 100
v2 = *mantissa\mantissa[i + 1] / 100
*mantissa\mantissa[i] = v1 * 1000000 + v2
*mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 100
Next
*mantissa\mantissa[dwords] = 1000000 * (*mantissa\mantissa[dwords] % 100)
EndProcedure
Procedure rshift_7 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = dwords To 1 Step -1
v1 = *mantissa\mantissa[i] / 10000000
v2 = *mantissa\mantissa[i - 1] % 10000000
v2 = v2 * 10 + v1
*mantissa\mantissa[i] = v2
Next
*mantissa\mantissa[0] = *mantissa\mantissa[0] / 10000000
EndProcedure
Procedure lshift_7 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v1, v2, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = 0 To dwords - 1
v1 = *mantissa\mantissa[i] % 10
v2 = *mantissa\mantissa[i + 1] / 10
*mantissa\mantissa[i] = v1 * 10000000 + v2
*mantissa\mantissa[i + 1] = *mantissa\mantissa[i + 1] % 10
Next
*mantissa\mantissa[dwords] = 10000000 * (*mantissa\mantissa[dwords] % 10)
EndProcedure
Procedure rshift_8 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = dwords To 1 Step -1
*mantissa\mantissa[i] = *mantissa\mantissa[i - 1]
Next
*mantissa\mantissa[0] = 0
EndProcedure
Procedure lshift_8 (*mantissa.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, i
dwords.l = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
For i = 0 To dwords - 1
*mantissa\mantissa[i] = *mantissa\mantissa[i + 1]
Next
*mantissa\mantissa[dwords] = 0
EndProcedure
Procedure.l fpcmp (*x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
Define.l dwords
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
Define.l c, i
If *x\sign < *y\sign
ProcedureReturn -1
EndIf
If *x\sign > *y\sign
ProcedureReturn 1
EndIf
If *x\exponent<*y\exponent
If *x\sign=0
ProcedureReturn -1
Else
ProcedureReturn 1
EndIf
EndIf
If *x\exponent>*y\exponent
If *x\sign=0
ProcedureReturn 1
Else
ProcedureReturn -1
EndIf
EndIf
For i=0 To dwords
c=*x\mantissa[i]-*y\mantissa[i]
If c<>0
Break
EndIf
Next
If c=0
ProcedureReturn 0
EndIf
If c<0
If *x\sign=0
ProcedureReturn -1
Else
ProcedureReturn 1
EndIf
EndIf
If c>0
If *x\sign=0
ProcedureReturn 1
Else
ProcedureReturn -1
EndIf
EndIf
EndProcedure
Procedure norm_fac1 (*fac1.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, i, er, f
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
; normalize the number in fac1
; all routines exit through this one.
;see if the mantissa is all zeros.
;if so, sit the exponent and sign equal to 0.
er = 0: f = 0
For i = 0 To dwords
If *fac1\mantissa[i] > 0
f = 1
EndIf
Next
If f = 0
*fac1\exponent = 0
*fac1\sign = 0
ProcedureReturn
;if the highmost digit in fac1_man is nonzero,
;shift the mantissa right 1 digit and
;increment the exponent
ElseIf *fac1\mantissa[0] > 99999999
rshift_1 (*fac1, dwords)
*fac1\exponent = *fac1\exponent + 1
Else
;now shift fac1_man 1 to the left until a
;nonzero digit appears in the next-to-highest
;digit of fac1_man. decrement exponent for
;each shift.
While *fac1\mantissa[0] = 0
lshift_8 (*fac1, dwords)
*fac1\exponent = *fac1\exponent - 8
If *fac1\exponent = 0
PrintN( "norm_fac1=expu_err")
ProcedureReturn
EndIf
Wend
If *fac1\mantissa[0] < 10
lshift_7 (*fac1, dwords)
*fac1\exponent = *fac1\exponent - 7
ElseIf *fac1\mantissa[0] < 100
lshift_6 (*fac1, dwords)
*fac1\exponent = *fac1\exponent - 6
ElseIf *fac1\mantissa[0] < 1000
lshift_5 (*fac1, dwords)
*fac1\exponent = *fac1\exponent - 5
ElseIf *fac1\mantissa[0] < 10000
lshift_4 (*fac1, dwords)
*fac1\exponent = *fac1\exponent - 4
ElseIf *fac1\mantissa[0] < 100000
lshift_3 (*fac1, dwords)
*fac1\exponent = *fac1\exponent - 3
ElseIf *fac1\mantissa[0] < 1000000
lshift_2 (*fac1, dwords)
*fac1\exponent = *fac1\exponent - 2
ElseIf *fac1\mantissa[0] < 10000000
lshift_1 (*fac1, dwords)
*fac1\exponent = *fac1\exponent - 1
EndIf
EndIf
;check for overflow/underflow
If *fac1\exponent < 0
PrintN( "norm_fac1=#expo_err")
EndIf
EndProcedure
Procedure fpadd_aux (*fac1.decfloat, *fac2.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v, c, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
c = 0
For i = dwords To 1 Step -1
v = *fac2\mantissa[i] + *fac1\mantissa[i] + c
If v > 99999999
v = v - 100000000
c = 1
Else
c = 0
EndIf
*fac1\mantissa[i] = v
Next
v = *fac1\mantissa[0] + *fac2\mantissa[0] + c
*fac1\mantissa[0] = v
norm_fac1 (*fac1, dwords)
EndProcedure
Procedure fpsub_aux (*fac1.decfloat, *fac2.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, v, c, i
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
c = 0
For i = dwords To 1 Step -1
v = *fac1\mantissa[i] - *fac2\mantissa[i] - c
If v < 0
v = v + 100000000
c = 1
Else
c = 0
EndIf
*fac1\mantissa[i] = v
Next
v = *fac1\mantissa[0] - *fac2\mantissa[0] - c
*fac1\mantissa[0] = v
norm_fac1 (*fac1, dwords)
EndProcedure
Procedure fpadd (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, i, t, c, xsign, ysign
Define.decfloat fac1, fac2
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
xsign = *x\sign: *x\sign = 0
ysign = *y\sign: *y\sign = 0
c = fpcmp(*x, *y, dwords_in)
*x\sign = xsign
*y\sign = ysign
If c < 0
copydec(*y, fac1) ;fac1 = y
copydec(*x, fac2) ;fac2 = x
Else
copydec(*x, fac1) ;fac1 = x
copydec(*y, fac2) ;fac2 = y
EndIf
t = fac1\exponent - fac2\exponent
t = ((fac1\exponent & $7fffffff) - #bias - 1) - ((fac2\exponent & $7fffffff) - #bias - 1)
If t < (dwords*8 + 8)
;the difference between the two
;exponents indicate how many times
;we have to multiply the mantissa
;of fac2 by 10 (i.e., shift it right 1 place).
;if we have to shift more times than
;we have digits, the result is already in fac1.
t = fac1\exponent - fac2\exponent
If t > 0 And t < (dwords*8 + 8) ;shift
i = t / 8
While i > 0
rshift_8 (fac2, dwords)
t = t - 8
i = i - 1
Wend
If t = 7
rshift_7 (fac2, dwords)
ElseIf t = 6
rshift_6 (fac2, dwords)
ElseIf t = 5
rshift_5 (fac2, dwords)
ElseIf t = 4
rshift_4 (fac2, dwords)
ElseIf t = 3
rshift_3 (fac2, dwords)
ElseIf t = 2
rshift_2 (fac2, dwords)
ElseIf t = 1
rshift_1 (fac2, dwords)
EndIf
EndIf
;see if the signs of the two numbers
;are the same. if so, add; if not, subtract.
If fac1\sign = fac2\sign ;add
fpadd_aux (fac1, fac2, dwords)
Else
fpsub_aux (fac1, fac2, dwords)
EndIf
EndIf
copydec(fac1, *result) ;result = fac1
EndProcedure
Procedure fpsub (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
Define.l dwords
Define.decfloat fac1, fac2
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
copydec(*x, fac1) ;fac1 = x
copydec(*y, fac2) ;fac2 = y
fac2\sign = fac2\sign ! $8000
fpadd (*result, fac1, fac2, dwords)
EndProcedure
Procedure fpmul (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, i, j, ex, er, ndw, den, num
Define.q digit, carry
Define.decfloat fac1, fac2
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
Dim fac3.q(dwords)
copydec(*x, fac1) ;fac1 = x
copydec(*y, fac2) ;fac2 = y
;check exponents. if either is zero,
;the result is zero
If fac1\exponent = 0 Or fac2\exponent = 0 ;result is zero...clear fac1.
fac1\sign = 0
fac1\exponent = 0
For i = 0 To dwords
fac1\mantissa[i] = 0
Next
copydec(fac1, *result) ;result = fac1
ProcedureReturn
Else
If ex < 0
er = #expo_err
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
;clear fac3 mantissa
For i = 0 To dwords
fac3(i) = 0
Next
den = dwords
While fac2\mantissa[den] = 0
den = den - 1
Wend
num = dwords
While fac1\mantissa[num] = 0
num = num - 1
Wend
If num < den
;Swap fac1, fac2
copydec(*y, fac1) ;;fac1=y
copydec(*x, fac2) ;;fac2=x
Swap den, num
EndIf
For j = den To 0 Step -1
carry = 0
digit = fac2\mantissa[j]
For i = num To 0 Step -1
fac3(i) = fac3(i) + digit * fac1\mantissa[i]
Next
For i = num To 0 Step -1
fac3(i) = fac3(i) + carry
carry = fac3(i) / 100000000
fac3(i) = (fac3(i) % 100000000)
Next
For i = dwords To 1 Step -1
fac3(i) = fac3(i - 1)
Next
fac3(0) = carry
Next
For i = 0 To dwords
ndw = fac3(i)
fac1\mantissa[i] = ndw
Next
EndIf
;now determine exponent of result.
;as you do...watch for overflow.
ex = fac2\exponent - #bias + fac1\exponent
fac1\exponent = ex
;determine the sign of the product
fac1\sign = fac1\sign ! fac2\sign
norm_fac1 (fac1, dwords)
copydec(fac1, *result) ;result = fac1
EndProcedure
Procedure fpmul_si (*result.decfloat, *x.decfloat, y.q, dwords_in.l=#num_dwords)
Define.l dwords, ex, er, i, ndw, count
Define.decfloat fac1, fac2
Define.q carry, digit, prod, value
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
copydec(*x, fac1) ;fac1 = x
digit = Abs(y)
If digit > 99999999
si2fp (fac2, y, dwords)
fpmul (fac1, fac1, fac2, dwords)
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
;check exponents. if either is zero,
;the result is zero
If fac1\exponent = 0 Or y = 0 ;result is zero...clear fac1.
fac1\sign = 0
fac1\exponent = 0
For count = 0 To dwords
fac1\mantissa[count] = 0
Next
norm_fac1 (fac1, dwords)
copydec(fac1, *result) ;result = fac1
ProcedureReturn
Else
If digit = 1
If y < 0
fac1\sign = fac1\sign ! $8000
EndIf
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
;now determine exponent of result.
;as you do...watch for overflow.
If ex < 0
er = #expo_err
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
carry = 0
For i = dwords To 0 Step -1
prod = digit * fac1\mantissa[i] + carry
value = (prod % 100000000)
ndw = value
fac1\mantissa[i] = ndw
carry = prod / 100000000
Next
If carry < 10
rshift_1 (fac1, dwords)
fac1\exponent = fac1\exponent + 1
fac1\mantissa[0] = fac1\mantissa[0] + carry * 10000000
ElseIf carry < 100
rshift_2 (fac1, dwords)
fac1\exponent = fac1\exponent + 2
fac1\mantissa[0] = fac1\mantissa[0] + carry * 1000000
ElseIf carry < 1000
rshift_3 (fac1, dwords)
fac1\exponent = fac1\exponent + 3
fac1\mantissa[0] = fac1\mantissa[0] + carry * 100000
ElseIf carry < 10000
rshift_4 (fac1, dwords)
fac1\exponent = fac1\exponent + 4
fac1\mantissa[0] = fac1\mantissa[0] + carry * 10000
ElseIf carry < 100000
rshift_5 (fac1, dwords)
fac1\exponent = fac1\exponent + 5
fac1\mantissa[0] = fac1\mantissa[0] + carry * 1000
ElseIf carry < 1000000
rshift_6 (fac1, dwords)
fac1\exponent = fac1\exponent + 6
fac1\mantissa[0] = fac1\mantissa[0] + carry * 100
ElseIf carry < 10000000
rshift_7 (fac1, dwords)
fac1\exponent = fac1\exponent + 7
fac1\mantissa[0] = fac1\mantissa[0] + carry * 10
ElseIf carry < 100000000
rshift_8 (fac1, dwords)
fac1\exponent = fac1\exponent + 8
fac1\mantissa[0] = fac1\mantissa[0] + carry
EndIf
EndIf
norm_fac1 (fac1, dwords)
If y < 0
fac1\sign = fac1\sign ! $8000
EndIf
copydec(fac1, *result) ;result = fac1
EndProcedure
Procedure fpmul_ll (*result.decfloat, *x.decfloat, y.q, dwords_in.l=#num_dwords)
Define.l dwords, count, ex, er, i, ndw
Define.decfloat fac1, fac2
Define.q n0, n1, carry, digit, prod, value
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
copydec(*x, fac1) ;fac1 = x
digit = Abs(y)
If digit > 99999999 And digit < 10000000000000000
n0 = digit / 100000000
n1 = digit % 100000000
fpmul_si (fac2, fac1, n1, dwords)
fac1\exponent = fac1\exponent + 8
fpmul_si (fac1, fac1, n0, dwords)
fpadd (fac1, fac1, fac2, dwords)
If y < 0
fac1\sign = fac1\sign ! $8000
EndIf
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
If digit > 9999999999999999
str2fp (fac2, Str(y))
fpmul (fac1, fac1, fac2, dwords)
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
si2fp (fac2, y, dwords)
;check exponents. if either is zero,
;the result is zero
If fac1\exponent = 0 Or y = 0 ;result is zero...clear fac1.
fac1\sign = 0
fac1\exponent = 0
For count = 0 To dwords
fac1\mantissa[count] = 0
Next
norm_fac1 (fac1, dwords)
copydec(fac1, *result) ;result = fac1
ProcedureReturn
Else
If digit = 1
If y < 0
fac1\sign = fac1\sign ! $8000
EndIf
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
;now determine exponent of result.
;as you do...watch for overflow.
If ex < 0
er = #expo_err
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
carry = 0
For i = dwords To 0 Step -1
prod = digit * fac1\mantissa[i] + carry
value = (prod % 100000000)
ndw = value
fac1\mantissa[i] = ndw
carry = prod / 100000000
Next
If carry < 10
rshift_1 (fac1, dwords)
fac1\exponent = fac1\exponent + 1
fac1\mantissa[0] = fac1\mantissa[0] + carry * 10000000
ElseIf carry < 100
rshift_2 (fac1, dwords)
fac1\exponent = fac1\exponent + 2
fac1\mantissa[0] = fac1\mantissa[0] + carry * 1000000
ElseIf carry < 1000
rshift_3 (fac1, dwords)
fac1\exponent = fac1\exponent + 3
fac1\mantissa[0] = fac1\mantissa[0] + carry * 100000
ElseIf carry < 10000
rshift_4 (fac1, dwords)
fac1\exponent = fac1\exponent + 4
fac1\mantissa[0] = fac1\mantissa[0] + carry * 10000
ElseIf carry < 100000
rshift_5 (fac1, dwords)
fac1\exponent = fac1\exponent + 5
fac1\mantissa[0] = fac1\mantissa[0] + carry * 1000
ElseIf carry < 1000000
rshift_6 (fac1, dwords)
fac1\exponent = fac1\exponent + 6
fac1\mantissa[0] = fac1\mantissa[0] + carry * 100
ElseIf carry < 10000000
rshift_7 (fac1, dwords)
fac1\exponent = fac1\exponent + 7
fac1\mantissa[0] = fac1\mantissa[0] + carry * 10
ElseIf carry < 100000000
rshift_8 (fac1, dwords)
fac1\exponent = fac1\exponent + 8
fac1\mantissa[0] = fac1\mantissa[0] + carry
EndIf
EndIf
norm_fac1 (fac1, dwords)
If y < 0
fac1\sign = fac1\sign ! $8000
EndIf
copydec(fac1, *result) ;result = fac1
EndProcedure
Procedure fpreciprocal (*result.decfloat, *m.decfloat, dwords_in.l=#num_dwords)
Define.l dwords
Define.d x
dwords=dwords_in
If dwords>#NUM_DWORDS
dwords=#NUM_DWORDS
EndIf
Define.l k, l, i, ex, s, prec
Define.decfloat r, r2, one, n
copydec(*m, n) ;n=m
l=Log((dwords*8+8)*0.0625)*1.5
Define.s v,f, ts
If n\exponent>0
ex=(n\exponent & $7FFFFFFF)-#BIAS-1
Else
ex=0
EndIf
If n\sign
v="-"
Else
v=" "
EndIf
ts=Str(n\mantissa[0])
If Len(ts)<8
ts=ts+string(8-Len(ts),"0")
EndIf
v=v+Left(ts,1)+"."+Mid(ts,2)
ts=Str(n\mantissa[1])
If Len(ts)<8
ts=string(8-Len(ts),"0")+ts
EndIf
v=v+ts
x=ValD(v)
If x = 0
PrintN( "Div 0")
copydec(r, *result)
ProcedureReturn
EndIf
If x = 1 And ex=0
str2fp(r,"1")
copydec(r, *result)
ProcedureReturn
ElseIf x=1
x=10
ex=ex-1
EndIf
ex=(-1)-ex
x=1/x
str2fp(r,StrD(x,16))
r\exponent=ex+#BIAS+1
one\mantissa[0]=10000000
one\exponent=#BIAS+1
copydec(r, r2) ;r2=r
prec=3
For k=1 To l
prec=2*prec-1
fpmul(r2,n,r, prec)
fpsub(r2,one, r2, prec)
fpmul(r2,r, r2, prec)
fpadd(r,r,r2, prec)
Next
copydec(r, *result)
EndProcedure
Procedure.l min (a.l, b.l)
If a < b
ProcedureReturn a
Else
ProcedureReturn b
EndIf
EndProcedure
Procedure.d realw (Array w.d(1), j.l)
Define.d wx
wx = ((w(j - 1) * 10000 + w(j)) * 10000 + w(j + 1)) * 10000
If ArraySize(w()) >= (j + 2)
wx = wx + w(j + 2)
EndIf
ProcedureReturn wx
EndProcedure
Procedure subtract (Array w.d(1), q.l, Array d.d(1), ka.l, kb.l)
Define.l j
For j = ka To kb
w(j) = w(j) - q * d(j - ka + 2)
Next
EndProcedure
Procedure normalise (Array w.d(1), ka.l, q.l)
w(ka) = w(ka) + w(ka - 1) * 10000
w(ka - 1) = q
EndProcedure
Procedure finalnorm (Array w.d(1), kb.l)
Define.l carry, j
For j = kb To 3 Step -1
If w(j) < 0
carry = ((-Int(w(j)) - 1) / 10000) + 1
Else
If w(j) >= 10000
carry = -(Int(w(j)) / 10000)
Else
carry = 0
EndIf
EndIf
w(j) = w(j) + carry * 10000
w(j - 1) = w(j - 1) - carry
Next
EndProcedure
Procedure fpdiv (*result.decfloat, *x.decfloat, *y.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, i, er, is_power_of_ten, tmp, b, j, stp, last, laststep, q, t
Define.d xd, xn, rund
Define.decfloat fac1, fac2
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
copydec(*x, fac1) ;fac1 = x
copydec(*y, fac2) ;fac2 = y
If fac2\exponent = 0 ; if fac2 = 0, return
; a divide-by-zero error and
; bail out.
For i = 0 To dwords
fac1\mantissa[i] = 99999999
Next
fac1\exponent = 99999 + #bias + 1
er = #divz_err
copydec(fac1, *result) ;result = fac1
ProcedureReturn
ElseIf fac1\exponent = 0 ;fact1=0, just return
er = 0
copydec(fac1, *result) ;result = fac1
ProcedureReturn
Else
;check to see if fac2 is a power of ten
is_power_of_ten = 0
If fac2\mantissa[0] = 10000000
is_power_of_ten = 1
For i = 1 To dwords
If fac2\mantissa[i] <> 0
is_power_of_ten = 0
Break
EndIf
Next
EndIf
;if fac2 is a power of ten then all we need to do is to adjust the sign and exponent and we are finished
If is_power_of_ten = 1
fac1\sign = fac1\sign ! fac2\sign
fac1\exponent = fac1\exponent - fac2\exponent + #bias + 1
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
Dim result.d(2 * dwords + 3)
Dim n.d(2 * dwords + 3)
Dim d.d(2 * dwords + 3)
b= 10000
Dim w.d(ArraySize(n()) + 4)
For j = 0 To dwords
tmp=fac1\mantissa[j] / 10000
n(2 * j + 2) = tmp
tmp=fac1\mantissa[j] % 10000
n(2 * j + 3) = tmp
tmp=fac2\mantissa[j] / 10000
d(2 * j + 2) = tmp
tmp=fac2\mantissa[j] % 10000
d(2 * j + 3) = tmp
Next
tmp=(fac1\exponent & $7fffffff) - #bias - 1
n(1) = tmp
tmp=(fac2\exponent & $7fffffff) - #bias - 1
d(1) = tmp
For j = ArraySize(n()) To ArraySize(w())
w(j) = 0
Next
t = ArraySize(n()) - 1
w(1) = n(1) - d(1) + 1
w(2) = 0
For j = 2 To ArraySize(n())
w(j + 1) = n(j)
Next
xd = (d(2) * b + d(3)) * b + d(4) + d(5) / b
laststep = t + 2
For stp = 1 To laststep
xn = realw(w(), (stp + 2))
q = Int(xn / xd)
last = min(stp + t + 1, ArraySize(w()))
subtract (w(), q, d(), (stp + 2), last)
normalise (w(), (stp + 2), q)
Next
finalnorm (w(), (laststep + 1))
If w(2) <> 0
laststep = laststep - 1
EndIf
rund = w(laststep + 1) / b
If rund >= 0.5
w(laststep) = w(laststep) + 1
EndIf
If w(2) = 0
For j = 1 To t + 1
result(j) = w(j + 1)
Next
Else
For j = 1 To t + 1
result(j) = w(j)
Next
EndIf
If w(2) = 0
result(1) = w(1) - 1
Else
result(1) = w(1)
EndIf
For j = 0 To dwords
fac1\mantissa[j] = result(2 * j + 2) * 10000 + result(2 * j + 3)
Next
norm_fac1 (fac1, dwords)
fac1\exponent = (result(1) + #bias)
EndIf
fac1\sign = fac1\sign ! fac2\sign
copydec(fac1, *result) ;result = fac1
EndProcedure
; sqrt(num)
; sqrt(num)
Procedure fpsqr (*result.decfloat, *num.decfloat, dwords_in.l=#num_dwords)
Define.l dwords, ex, aex, k, l, prec
Define.decfloat r, r2, tmp, n, half
Define.s ts, v
Define.d x
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
l = Log((dwords*8 + 8) * 0.0625) * 1.5
;l=estimated number of iterations needed
;first estimate is accurate to about 16 digits
;l is approximatly = to log2((#num_digits+9)/16)
;#num_digits+9 because decfloat has an extra 9 guard digits
copydec(*num, n) ;n = num
si2fp (tmp, 0, dwords)
If fpcmp(n, tmp, dwords) = 0
si2fp (r, 0, dwords)
copydec(r, *result) ;result = r
ProcedureReturn
EndIf
si2fp (tmp, 1, dwords)
If fpcmp(n, tmp, dwords) = 0
si2fp (r, 1, dwords)
copydec(r, *result) ;result = r
ProcedureReturn
EndIf
si2fp (tmp, 0, dwords)
If fpcmp(n, tmp, dwords) < 0
si2fp (r, 0, dwords)
copydec(r, *result) ;result = r
ProcedureReturn
EndIf
;=====================================================================
;hack to bypass the limitation of double exponent range
;in case the number is larger than what a double can handle
;for example, if the number is 2e500
;we separate the exponent and mantissa in this case 2
;if the exponent is odd then multiply the mantissa by 10
;take the square root and assign it to decfloat
;divide the exponent in half for square root
;in this case 1.414213562373095e250
If n\exponent > 0
ex = (n\exponent & $7fffffff) - #bias - 1
Else
ex = 0
EndIf
ts = Trim(Str(n\mantissa[0]))
If Len(ts) < 8
ts = ts + String(8 - Len(ts), "0")
EndIf
v = v + Left(ts, 1) + "." + Mid(ts, 2)
ts = Trim(Str(n\mantissa[1]))
If Len(ts) < 8
ts = String(8 - Len(ts), "0") + ts
EndIf
v = v + ts
x = ValD(v)
If x = 0
PrintN( "div 0")
copydec(r, *result) ;result = r
ProcedureReturn
EndIf
If x = 1 And ex = 0
si2fp (r, 1, dwords)
copydec(r, *result) ;result = r
ProcedureReturn
EndIf
aex=ex
If aex<0
aex=-aex
EndIf
If aex & 1
x = x * 10
ex = ex - 1
EndIf
x = Sqr(x) ;approximation
v = Trim(StrD(x,16))
k = FindString(v, ".")
str2fp (r, v)
r\exponent = ex / 2 + #bias + 1
If Len(v) > 1 And k = 0
r\exponent = r\exponent + 1
EndIf
For k = 1 To dwords
half\mantissa[k] = 0
Next
half\mantissa[0] = 50000000
half\exponent = #bias
half\sign = 0
;=====================================================================
;newton-raphson method
prec = 3
For k = 1 To l + 1
prec = 2 * prec - 1
fpdiv (tmp, n, r, prec)
fpadd (r2, r, tmp, prec)
fpmul (r, r2, half, prec)
Next
copydec(r, *result) ;result = r
EndProcedure
Procedure fpdiv_si (*result.decfloat, *num.decfloat, den.l, dwords_in.l=#num_dwords)
Define.l dwords, ndw, i
Define.q carry, remder, divisor, quotient
Define.decfloat fac1
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
remder = 0
divisor = Abs(den)
copydec(*num, fac1) ;fac1 = num
If divisor = 0
PrintN( "error: divisor = 0")
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
If divisor > 99999999
PrintN( "error: divisor too large")
copydec(fac1, *result) ;result = fac1
ProcedureReturn
EndIf
For i = 0 To dwords
ndw = i
quotient = fac1\mantissa[ndw] + remder * 100000000
remder = quotient % divisor
fac1\mantissa[ndw] = quotient / divisor
Next
quotient = remder * 100000000
quotient = quotient / divisor
carry = fac1\mantissa[0]
If carry = 0
lshift_8 (fac1, dwords)
fac1\exponent = fac1\exponent - 8
fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient
ElseIf carry < 10
lshift_7 (fac1, dwords)
fac1\exponent = fac1\exponent - 7
fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 10
ElseIf carry < 100
lshift_6 (fac1, dwords)
fac1\exponent = fac1\exponent - 6
fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 100
ElseIf carry < 1000
lshift_5 (fac1, dwords)
fac1\exponent = fac1\exponent - 5
fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 1000
ElseIf carry < 10000
lshift_4 (fac1, dwords)
fac1\exponent = fac1\exponent - 4
fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 10000
ElseIf carry < 100000
lshift_3 (fac1, dwords)
fac1\exponent = fac1\exponent - 3
fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 100000
ElseIf carry < 1000000
lshift_2 (fac1, dwords)
fac1\exponent = fac1\exponent - 2
fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 1000000
ElseIf carry < 10000000
lshift_1 (fac1, dwords)
fac1\exponent = fac1\exponent - 1
fac1\mantissa[dwords] = fac1\mantissa[dwords] + quotient / 10000000
EndIf
;norm_fac1(fac1)
If den < 0
fac1\sign = fac1\sign ! $8000
EndIf
copydec(fac1, *result) ;result = fac1
EndProcedure
;fractional part of num
Procedure fpfrac (*result.decfloat, *num.decfloat, dwords_in.l=#num_dwords)
Define.l dwords
Define.decfloat n
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
fptrunc (n, *num)
fpsub (n, *num, n, dwords)
copydec(n, *result) ;result = n
EndProcedure
;returns the positive of n
Procedure fpabs (*result.decfloat, *n.decfloat)
Define.decfloat x
copydec(*n, x) ;x = n
x\sign = 0
copydec(x, *result) ;result = x
EndProcedure
;changes the sign of n, if n is positive then n will be negative & vice versa
Procedure fpneg (*result.decfloat, *n.decfloat)
Define.decfloat x
copydec(*n, x) ;x = n
x\sign = x\sign ! $8000
copydec(x, *result) ;result = x
EndProcedure
Procedure fpfmod (*quotient.decfloat, *f_mod.decfloat, *num.decfloat, *denom.decfloat, dwords_in.l=#num_dwords)
Define.l dwords
Define.decfloat q, fm
If dwords > #num_dwords
dwords = #num_dwords
EndIf
fpdiv (fm, *num, *denom, dwords)
fptrunc (q, fm)
fpmul(q, q, *denom, dwords)
fpsub (fm, *num, q, dwords)
copydec(q, *quotient) ;quotient = q
copydec(fm, *f_mod) ;f_mod = fm
EndProcedure
Procedure fpeps (*result.decfloat, digits_in.l=#number_of_digits)
Define.l digits, dwords
Define.decfloat ep
digits=digits_in
If digits > #number_of_digits
digits = #number_of_digits
EndIf
dwords=digits/8
si2fp (ep, 5, dwords)
ep\exponent = (-(digits) + #bias + 1)
copydec(ep, *result) ;result = ep
EndProcedure
Procedure fpipow (*result.decfloat, *x.decfloat, e.q, dwords_in.l=#num_dwords)
Define.l dwords, i
Define.decfloat y, z
Define.q n, c
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
;take x to an long power
c = 0
copydec(*x, y) ;y = x
n = Abs(e)
z\sign = 0
z\exponent = (#bias + 1)
z\mantissa[0] = 10000000
For i = 1 To #num_dwords
z\mantissa[i] = 0
Next
While n > 0
While (n & 1) = 0
n = n / 2
fpmul (y, y, y, dwords)
c = c + 1
Wend
n = n - 1
fpmul (z, y, z, dwords)
c = c + 1
Wend
If e < 0
si2fp (y, 1, dwords)
fpdiv (z, y, z, dwords)
EndIf
copydec(z, *result) ;result = z
EndProcedure
Procedure fpfactorial (*result.decfloat, n.l, dwords_in.l=#num_dwords)
Define.l dwords
Define.q i
Define.decfloat f
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
If n < 0
PrintN( "inf")
copydec(f, *result) ;result = f
ProcedureReturn
EndIf
If n = 0 Or n = 1
si2fp (f, 1, dwords)
copydec(f, *result) ;result = f
ProcedureReturn
EndIf
si2fp (f, 2, dwords_in)
If n = 2
copydec(f, *result) ;result = f
ProcedureReturn
EndIf
For i = 3 To n
fpmul_si (f, f, i, dwords)
Next
copydec(f, *result) ;result = f
EndProcedure
Procedure fpnroot (*result.decfloat, *x.decfloat, p.l, dwords_in.l=#num_dwords)
Define.l dwords, i, ex, l, prec
Define.d t, t2
Define.decfloat ry, tmp, tmp2
dwords = dwords_in
If dwords > #num_dwords
dwords = #num_dwords
EndIf
l = Log((dwords*8 + 8) * 0.0625) * 1.5
ex = (*x\exponent & $7fffffff) - #bias - 1
t = *x\mantissa[0] + *x\mantissa[1] / 100000000
t = t / 10000000
t2 = Log(t) / p + Log(10) * ex / p
t2 = Exp(t2)
str2fp (ry, StrD(t2,16))
prec = 3
fpipow (tmp, ry, p - 1, prec)
fpdiv (tmp, *x, tmp, prec)
fpmul_si (tmp2, ry, p - 1, prec)
fpadd (tmp2, tmp2, tmp, prec)
fpdiv_si (ry, tmp2, p, prec)
For i = 1 To l
prec = 2 * prec - 1
fpipow (tmp, ry, p - 1, prec)
fpdiv (tmp, *x, tmp, prec)
fpmul_si (tmp2, ry, p - 1, prec)
fpadd (tmp2, tmp2, tmp, prec)
fpdiv_si (ry, tmp2, p, prec)
Next
copydec(ry, *result) ;result = ry
EndProcedure
Procedure pi_brent_salamin (*pi_bs.decfloat, digits_in.l= #number_of_digits)
Define.l dwords, digits, limit2
Define.decfloat c0, c1, c2, c05, a, b, sum, ak, bk, ck, ab, asq
Define.decfloat pow2, tmp
digits=digits_in
If digits > #number_of_digits
digits = #number_of_digits
EndIf
dwords = digits/8
If (dwords*8)<digits
dwords+1
EndIf
limit2 = -(digits) + #bias + 1
si2fp (c0, 0, dwords)
copydec(c0, ak) ;ak = c0
copydec(c0, bk) ;bk = c0
copydec(c0, ab) ;ab = c0
copydec(c0, asq);asq = c0
si2fp (c1, 1, dwords)
copydec(c1, a) ;a = c1
copydec(c1, ck);ck = c1
copydec(c1, pow2) ;pow2 = c1
si2fp (c2, 2, dwords)
copydec(c2, b) ;b = c2
str2fp (c05, ".5")
copydec(c05, sum) ;sum = c05
si2fp (*pi_bs, 3, dwords)
fpsqr (b, b, dwords)
fpdiv (b, c1, b, dwords)
While fpcmp(ck, c0, dwords) <> 0 And ck\exponent > limit2
fpadd (ak, a, b, dwords)
fpmul (ak, c05, ak, dwords)
fpmul (ab, a, b, dwords)
fpsqr (bk, ab, dwords)
fpmul (asq, ak, ak, dwords)
fpsub (ck, asq, ab, dwords)
fpmul (pow2, pow2, c2, dwords)
fpmul (tmp, pow2, ck, dwords)
fpsub (sum, sum, tmp, dwords)
copydec(ak, a) ;a = ak
copydec(bk, b) ;b = bk
Wend
fpdiv (tmp, asq, sum, dwords)
fpmul (*pi_bs, c2, tmp, dwords)
EndProcedure
Code: Select all
#number_of_digits = 100 ; need to set the constant before including decfloat.pbi
XIncludeFile "decfloat.pbi"
OpenConsole()
Define.decfloat m, x, y, z, v
Define.s s
Define.d t
str2fp (x, "3."+String(#number_of_digits-1, "3"))
str2fp (y, "4."+String(#number_of_digits-1, "4"))
t=ElapsedMilliseconds()
PrintN("x = "+fp2str(x))
PrintN("y = "+fp2str(y))
fpadd(z, x, y)
PrintN("x + y = "+fp2str(z))
; c = fpcmp(x, y)
; PrintN("c="+Str(c))
fpsub(z, x, y)
PrintN("x - y = "+fp2str(z))
fpmul(z, x, y)
PrintN("x * y = "+fp2str(z))
fpdiv(z, x, y)
PrintN("x / y = "+fp2str(z))
si2fp (x, 2)
PrintN("x = "+fp2str(x))
fpsqr (z, x)
PrintN("sqrt(x) = "+fp2str(z))
fpipow (z, x, -3)
PrintN("x^(-3) = "+fp2str(z))
fpnroot (z, x, 3)
PrintN("x^(1/3) = "+fp2str(z))
fpfactorial (z, 10000)
PrintN("10000! = "+fp2str(z))
pi_brent_salamin (z)
PrintN("Pi = "+fp2str(z))
fpfrac(v, z)
PrintN("frac(Pi) = "+fp2str(v))
fptrunc(z, z)
PrintN("trunc(Pi)= "+fp2str(z))
si2fp (x, 1)
;PrintN("x = "+fp2str(x))
str2fp (y, ".18")
PrintN("y = "+fp2str(y))
;fpdiv (z, x, y)
fpreciprocal(z, y)
PrintN("reciprocal(x)="+fp2str(z))
t=ElapsedMilliseconds()-t
PrintN( "elapsed time "+StrD(t/1000)+" seconds")
Input()
CloseConsole()
output
corrected a bug with str2fpx = 3.3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333
y = 4.4444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444
x + y = 7.7777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777
x - y = -1.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
x * y = 14.814814814814814814814814814814814814814814814814814814814814814814814814814814814814814814814814
x / y = 0.75000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
x = 2.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
sqrt(x) = 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415
x^(-3) = 0.12500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
x^(1/3) = 1.2599210498948731647672106072782283505702514647015079800819751121552996765139594837293965624362550
10000! = 2.84625968091705451890641321211986889014805140170279923079417999427441134000376444377299078675778477e+35659
Pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170
frac(Pi) = 0.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706
trunc(Pi)= 3.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
y = 0.18000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
reciprocal(y)= 5.5555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555
elapsed time 0.004 seconds
corrected fpmod as sugested by STARGÅTE, thank you STARGÅTE
exponent increased to 999999999
added the Procedure.s fp2strT (*n.decfloat, places) for trimmed output - no trailing 0's