Dest Console
Const NUMBER_OF_DIGITS
= 512 Const NUM_DIGITS
= NUMBER_OF_DIGITS
Const NUM_DWORDS
= NUM_DIGITS \
8 Const BIAS
= 1073741824 '2 ^ 30 Const MANTISSA_BYTES
= (NUM_DWORDS
+ 1) * 8 * 4
' Error definitions
Const DIVZ_ERR
= 1 'Divide by zero Const EXPO_ERR
= 2 'Exponent overflow error Const EXPU_ERR
= 3 'Exponent underflow error Call pi_brent_salamin
(pi_dec
, NUMBER_OF_DIGITS
)
Call fpfactorial
(z
, 10000, NUMBER_OF_DIGITS
) Print "factorial(10000) to "; NUMBER_OF_DIGITS;
" digits" Print fp2str
(z
, NUMBER_OF_DIGITS
) Print "in "; t;
" seconds" Call fpdiv
(z
, n
, m
, NUMBER_OF_DIGITS
) Print "3/9 to "; NUMBER_OF_DIGITS;
" digits" Print fp2str
(z
, NUMBER_OF_DIGITS
) Print "in "; t;
" seconds" Call fpmul
(z
, z
, z
, NUMBER_OF_DIGITS
) Print fp2str
(z
, NUMBER_OF_DIGITS
) Call fplog
(z
, n
, NUMBER_OF_DIGITS
) Print "ln(2) to "; NUMBER_OF_DIGITS;
" digits" Print fp2str
(z
, NUMBER_OF_DIGITS
) Print "in "; t;
" seconds" Call fpexp
(z
, z
, NUMBER_OF_DIGITS
) Print fp2str
(z
, NUMBER_OF_DIGITS
) Print "in "; t;
" seconds" Call si2fp
(n
, 3, NUMBER_OF_DIGITS
) Call fpinv
(z
, n
, NUMBER_OF_DIGITS
) Print fp2str
(z
, NUMBER_OF_DIGITS
) Print "in "; t;
" seconds"
ex
= (n.exponent
And &H7FFFFFFF) - BIAS
- 1 ex = 0
ts
= Str$(git
(n.mantissa
, 0)) ts = Trim$(ts)
For i
= 1 To NUM_DWORDS
- 1 ts
= Str$(git
(n.mantissa
, i
)) ts = Trim$(ts)
v = v + ts
v = v + f
fp2str_exp$ = v
'long part of num
Sub fpfix
(result
As decfloat
, num
As decfloat
)
ex
= (num.exponent
And &H7FFFFFFF) - BIAS
ex2 = ex \ 8
k = ex2
ex2 = ex2 - 1
Call set
(ip.mantissa
, ex2
, git
(num.mantissa
, ex2
)) Call set
(ip.mantissa
, k
, 10000000 * (git
(num.mantissa
, k
) \
10000000)) Call set
(ip.mantissa
, k
, 1000000 * (git
(num.mantissa
, k
) \
1000000)) Call set
(ip.mantissa
, k
, 100000 * (git
(num.mantissa
, k
) \
100000)) Call set
(ip.mantissa
, k
, 10000 * (git
(num.mantissa
, k
) \
10000)) Call set
(ip.mantissa
, k
, 1000 * (git
(num.mantissa
, k
) \
1000)) Call set
(ip.mantissa
, k
, 100 * (git
(num.mantissa
, k
) \
100)) Call set
(ip.mantissa
, k
, 10 * (git
(num.mantissa
, k
) \
10)) Call set
(ip.mantissa
, k
, git
(num.mantissa
, k
)) ip.exponent = ex + BIAS
ip.sign = num.sign
result = ip
ex
= (n.exponent
And &H7FFFFFFF) - BIAS
- 1 ex = 0
ts
= Trim$
(Str$(git
(n.mantissa
, 0))) ts
= Trim$
(Str$(git
(n.mantissa
, 1))) v = v + ts
v = v + f
ex
= (n.exponent
And &H7FFFFFFF) - BIAS
- 1 ex = 0
ts
= Trim$
(Str$(git
(n.mantissa
, 0))) v = ts
For i
= 1 To NUM_DWORDS
- 1 ts
= Trim$
(Str$(git
(n.mantissa
, i
))) v = v + ts
fp2str_fix$ = s + v
ex
= (n.exponent
And &H7FFFFFFF) - BIAS
- 1 ex = 0
fp2str = fp2str_fix(n, places)
fp2str = fp2str_exp(n, places)
Dim As Long j
, s
, d
, e
, ep
, ex
, es
, i
, f
, fp
, fln
j = 1
s = 1
d = 0
e = 0
ep = 0
ex = 0
es = 1
i = 0
f = 0
fp = 0
f1 = ""
f2 = ""
f3 = ""
j = j + 1
es = -es
c = ""
j = j + 1
j = j + 1
If (c
> "/") And (c
< ":") Then 'c is digit between 0 and 9 f3 = f3 + c
ex
= 10 * ex
+ (Asc(c
) - 48) j = j + 1
j = j + 1
s = -s
j = j + 1
j = j + 1
j = j + 1
d = 1
If (c
> "/") And (c
< ":") Then 'c is digit between 0 and 9 j = j + 1
e = e - 1
j = j + 1
f1 = f1 + c
i = i + 1
fp = 1
f2 = f2 + c
f = f + 1
ep = 1
j = j + 1
skip_while:
f = 0
f2 = ""
n.sign = s
ex = es * ex - 1 + i + e
f1 = f1 + f2
f1
= Mid$(f1
, 1, (NUM_DIGITS
+ 1 + 8)) f1 = f1 + "0"
j = 1
Call set
(n.mantissa
, i
, ulng
) j = j + 8
If fp
Then n.exponent
= (ex
+ BIAS
+ 1) Else n.exponent
= 0 result = n
Sub si2fp
(result
As decfloat
, m
As Integer64
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Call set
(fac1.mantissa
, i
, 0)
fac1.exponent = 0
fac1.sign = 0
fac1.exponent = BIAS
Call set
(fac1.mantissa
, 0, n
* 10000000) fac1.exponent = fac1.exponent + 1
Call set
(fac1.mantissa
, 0, n
* 1000000) fac1.exponent = fac1.exponent + 2
Call set
(fac1.mantissa
, 0, n
* 100000) fac1.exponent = fac1.exponent + 3
Call set
(fac1.mantissa
, 0, n
* 10000) fac1.exponent = fac1.exponent + 4
Call set
(fac1.mantissa
, 0, n
* 1000) fac1.exponent = fac1.exponent + 5
Call set
(fac1.mantissa
, 0, n
* 100) fac1.exponent = fac1.exponent + 6
Call set
(fac1.mantissa
, 0, n
* 10) fac1.exponent = fac1.exponent + 7
Call set
(fac1.mantissa
, 0, n
) fac1.exponent = fac1.exponent + 8
fac1.exponent = fac1.exponent + 8
Call set
(fac1.mantissa
, 0, n \
10) Call set
(fac1.mantissa
, 1, (n
Mod 10) * 10000000) fac1.exponent = fac1.exponent + 1
Call set
(fac1.mantissa
, 0, n \
100) Call set
(fac1.mantissa
, 1, (n
Mod 100) * 1000000) fac1.exponent = fac1.exponent + 2
Call set
(fac1.mantissa
, 0, n \
1000) Call set
(fac1.mantissa
, 1, (n
Mod 1000) * 100000) fac1.exponent = fac1.exponent + 3
Call set
(fac1.mantissa
, 0, n \
10000) Call set
(fac1.mantissa
, 1, (n
Mod 10000) * 10000) fac1.exponent = fac1.exponent + 4
Call set
(fac1.mantissa
, 0, n \
100000) Call set
(fac1.mantissa
, 1, (n
Mod 100000) * 1000) fac1.exponent = fac1.exponent + 5
Call set
(fac1.mantissa
, 0, n \
1000000) Call set
(fac1.mantissa
, 1, (n
Mod 1000000) * 100) fac1.exponent = fac1.exponent + 6
Call set
(fac1.mantissa
, 0, n \
10000000) Call set
(fac1.mantissa
, 1, (n
Mod 10000000) * 10) fac1.exponent = fac1.exponent + 7
Call set
(fac1.mantissa
, 0, n \
100000000) Call set
(fac1.mantissa
, 1, n
Mod 100000000) fac1.exponent = fac1.exponent + 8
fac1.sign = &H8000
fac1.sign = 0
result = fac1
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1 = git(mantissa.mantissa, i) \ 10
v2
= git
(mantissa.mantissa
, i
- 1) Mod 10 v2 = v2 * 10000000 + v1
Call set
(mantissa.mantissa
, i
, v2
) Call set
(mantissa.mantissa
, 0, git
(mantissa.mantissa
, 0) \
10)
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1
= git
(mantissa.mantissa
, i
) Mod 10000000 v2 = git(mantissa.mantissa, i + 1) \ 10000000
Call set
(mantissa.mantissa
, i
, v1
* 10 + v2
) Call set
(mantissa.mantissa
, i
+ 1, git
(mantissa.mantissa
, i
+ 1) Mod 10000000) Call set
(mantissa.mantissa
, digits
, 10 * (git
(mantissa.mantissa
, digits
) Mod 10000000))
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1 = git(mantissa.mantissa, i) \ 100
v2
= git
(mantissa.mantissa
, i
- 1) Mod 100 v2 = v2 * 1000000 + v1
Call set
(mantissa.mantissa
, i
, v2
) Call set
(mantissa.mantissa
, 0, git
(mantissa.mantissa
, 0) \
100)
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1
= git
(mantissa.mantissa
, i
) Mod 1000000 v2 = git(mantissa.mantissa, i + 1) \ 1000000
Call set
(mantissa.mantissa
, i
, v1
* 100 + v2
) Call set
(mantissa.mantissa
, i
+ 1, git
(mantissa.mantissa
, i
+ 1) Mod 1000000) Call set
(mantissa.mantissa
, digits
, 100 * (git
(mantissa.mantissa
, digits
) Mod 1000000))
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1 = git(mantissa.mantissa, i) \ 1000
v2
= git
(mantissa.mantissa
, i
- 1) Mod 1000 v2 = v2 * 100000 + v1
Call set
(mantissa.mantissa
, i
, v2
) Call set
(mantissa.mantissa
, 0, git
(mantissa.mantissa
, 0) \
1000)
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1
= git
(mantissa.mantissa
, i
) Mod 100000 v2 = git(mantissa.mantissa, i + 1) \ 100000
Call set
(mantissa.mantissa
, i
, v1
* 1000 + v2
) Call set
(mantissa.mantissa
, i
+ 1, git
(mantissa.mantissa
, i
+ 1) Mod 100000) Call set
(mantissa.mantissa
, digits
, 1000 * (git
(mantissa.mantissa
, digits
) Mod 100000))
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1 = git(mantissa.mantissa, i) \ 10000
v2
= git
(mantissa.mantissa
, i
- 1) Mod 10000 v2 = v2 * 10000 + v1
Call set
(mantissa.mantissa
, i
, v2
) Call set
(mantissa.mantissa
, 0, git
(mantissa.mantissa
, 0) \
10000)
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1
= git
(mantissa.mantissa
, i
) Mod 10000 v2 = git(mantissa.mantissa, i + 1) \ 10000
Call set
(mantissa.mantissa
, i
, v1
* 10000 + v2
) Call set
(mantissa.mantissa
, i
+ 1, git
(mantissa.mantissa
, i
+ 1) Mod 10000) Call set
(mantissa.mantissa
, digits
, 10000 * (git
(mantissa.mantissa
, digits
) Mod 10000))
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1 = git(mantissa.mantissa, i) \ 100000
v2
= git
(mantissa.mantissa
, i
- 1) Mod 100000 v2 = v2 * 1000 + v1
Call set
(mantissa.mantissa
, i
, v2
) Call set
(mantissa.mantissa
, 0, git
(mantissa.mantissa
, 0) \
100000)
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1
= git
(mantissa.mantissa
, i
) Mod 1000 v2 = git(mantissa.mantissa, i + 1) \ 1000
Call set
(mantissa.mantissa
, i
, v1
* 100000 + v2
) Call set
(mantissa.mantissa
, i
+ 1, git
(mantissa.mantissa
, i
+ 1) Mod 1000) Call set
(mantissa.mantissa
, digits
, 100000 * (git
(mantissa.mantissa
, digits
) Mod 1000))
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1 = git(mantissa.mantissa, i) \ 1000000
v2
= git
(mantissa.mantissa
, i
- 1) Mod 1000000 v2 = v2 * 100 + v1
Call set
(mantissa.mantissa
, i
, v2
) Call set
(mantissa.mantissa
, 0, git
(mantissa.mantissa
, 0) \
1000000)
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1
= git
(mantissa.mantissa
, i
) Mod 100 v2 = git(mantissa.mantissa, i + 1) \ 100
Call set
(mantissa.mantissa
, i
, v1
* 1000000 + v2
) Call set
(mantissa.mantissa
, i
+ 1, git
(mantissa.mantissa
, i
+ 1) Mod 100) Call set
(mantissa.mantissa
, digits
, 1000000 * (git
(mantissa.mantissa
, digits
) Mod 100))
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1 = git(mantissa.mantissa, i) \ 10000000
v2
= git
(mantissa.mantissa
, i
- 1) Mod 10000000 v2 = v2 * 10 + v1
Call set
(mantissa.mantissa
, i
, v2
) Call set
(mantissa.mantissa
, 0, git
(mantissa.mantissa
, 0) \
10000000)
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
v1
= git
(mantissa.mantissa
, i
) Mod 10 v2 = git(mantissa.mantissa, i + 1) \ 10
Call set
(mantissa.mantissa
, i
, v1
* 10000000 + v2
) Call set
(mantissa.mantissa
, i
+ 1, git
(mantissa.mantissa
, i
+ 1) Mod 10) Call set
(mantissa.mantissa
, digits
, 10000000 * (git
(mantissa.mantissa
, digits
) Mod 10))
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Call set
(mantissa.mantissa
, i
, git
(mantissa.mantissa
, i
- 1)) Call set
(mantissa.mantissa
, 0, 0)
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Call set
(mantissa.mantissa
, i
, git
(mantissa.mantissa
, i
+ 1)) Call set
(mantissa.mantissa
, digits
, 0)
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
If x.exponent
= y.exponent
Then c = git(x.mantissa, i) - git(y.mantissa, i)
' if we reach this point it means that the signs are different
' and if the sign of x is set meaning that x is negative then x < y
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
' normalize the number in fac1
' all routines exit through this one.
'see if the mantissa is all zeros.
'if so, set the exponent and sign equal to 0.
er = 0: f = 0
If git
(fac1.mantissa
, i
) > 0 Then f
= 1 fac1.exponent = 0
fac1.sign = 0
'if the highmost Digit in fac1_man is nonzero,
'shift the mantissa right 1 Digit and
'increment the exponent
Call RSHIFT_1
(fac1
, digits
) fac1.exponent = fac1.exponent + 1
'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 git
(fac1.mantissa
, 0) = 0 Call LSHIFT_8
(fac1
, digits
) fac1.exponent = fac1.exponent - 8
Print "NORM_FAC1=EXPU_ERR" If git
(fac1.mantissa
, 0) < 10 Then Call LSHIFT_7
(fac1
, digits
) fac1.exponent = fac1.exponent - 7
Call LSHIFT_6
(fac1
, digits
) fac1.exponent = fac1.exponent - 6
Call LSHIFT_5
(fac1
, digits
) fac1.exponent = fac1.exponent - 5
Call LSHIFT_4
(fac1
, digits
) fac1.exponent = fac1.exponent - 4
Call LSHIFT_3
(fac1
, digits
) fac1.exponent = fac1.exponent - 3
Call LSHIFT_2
(fac1
, digits
) fac1.exponent = fac1.exponent - 2
Call LSHIFT_1
(fac1
, digits
) fac1.exponent = fac1.exponent - 1
'check for overflow/underflow
Print "NORM_FAC1=EXPO_ERR"
Sub fpadd_aux
(fac1
As decfloat
, fac2
As decfloat
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
c = 0
v = git(fac2.mantissa, i) + git(fac1.mantissa, i) + c
v = v - 100000000
c = 1
c = 0
Call set
(fac1.mantissa
, i
, v
) v = git(fac1.mantissa, 0) + git(fac2.mantissa, 0) + c
Call set
(fac1.mantissa
, 0, v
)
Call NORM_FAC1
(fac1
, digits
)
Sub fpsub_aux
(fac1
As decfloat
, fac2
As decfloat
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
c = 0
v = git(fac1.mantissa, i) - git(fac2.mantissa, i) - c
v = v + 100000000
c = 1
c = 0
Call set
(fac1.mantissa
, i
, v
) v = git(fac1.mantissa, 0) - git(fac2.mantissa, 0) - c
Call set
(fac1.mantissa
, 0, v
)
Call NORM_FAC1
(fac1
, digits
)
Sub fpadd
(result
As decfloat
, x
As decfloat
, y
As decfloat
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As decfloat fac1
, fac2
xsign = x.sign: x.sign = 0
ysign = y.sign: y.sign = 0
c = fpcmp(x, y, digits)
x.sign = xsign
y.sign = ysign
fac1 = y
fac2 = x
fac1 = x
fac2 = y
t = fac1.exponent - fac2.exponent
t
= ((fac1.exponent
And &H7FFFFFFF) - BIAS
- 1) - ((fac2.exponent
And &H7FFFFFFF) - BIAS
- 1)
'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
< (NUM_DIGITS
+ 8) Then 'shift
i = t \ 8
Call RSHIFT_8
(fac2
, digits
) t = t - 8
i = i - 1
Call RSHIFT_7
(fac2
, digits
) Call RSHIFT_6
(fac2
, digits
) Call RSHIFT_5
(fac2
, digits
) Call RSHIFT_4
(fac2
, digits
) Call RSHIFT_3
(fac2
, digits
) Call RSHIFT_2
(fac2
, digits
) Call RSHIFT_1
(fac2
, digits
) 'See if the signs of the two numbers
'are the same. If so, add; if not, subtract.
If fac1.sign
= fac2.sign
Then 'add Call fpadd_aux
(fac1
, fac2
, digits
) Call fpsub_aux
(fac1
, fac2
, digits
) result = fac1
Sub fpsub
(result
As decfloat
, x
As decfloat
, y
As decfloat
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As decfloat fac1
, fac2
fac1 = x
fac2 = y
fac2.sign
= fac2.sign
Xor &H8000 Call fpadd
(fac1
, fac1
, fac2
, digits
) result = fac1
Sub fpmul
(result
As decfloat
, x
As decfloat
, y
As decfloat
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As decfloat fac1
, fac2
Dim As Integer64 digit
, carry
Dim As Unsigned Integer64 fac3
(0 To digits
)
fac1 = x
fac2 = y
'check exponents. if either is zero,
'the result is zero
If fac1.exponent
= 0 Or fac2.exponent
= 0 Then 'result is zero...clear fac1. fac1.sign = 0
fac1.exponent = 0
Call set
(fac1.mantissa
, i
, 0) 'NORM_FAC1(fac1)
er = EXPO_ERR
'clear fac3 mantissa
fac3(i) = 0
den = digits
While git
(fac2.mantissa
, den
) = 0 den = den - 1
num = digits
While git
(fac1.mantissa
, num
) = 0 num = num - 1
'fac1=y
'fac2=x
carry = 0
digit = git(fac2.mantissa, j)
fac3(i) = fac3(i) + digit * git(fac1.mantissa, i)
fac3(i) = fac3(i) + carry
carry = fac3(i) \ 100000000
fac3
(i
) = (fac3
(i
) Mod 100000000)
fac3(i) = fac3(i - 1)
fac3(0) = carry
Call set
(fac1.mantissa
, i
, fac3
(i
)) '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
Xor fac2.sign
Call NORM_FAC1
(fac1
, digits
) result = fac1
Sub fpmul_si
(result
As decfloat
, x
As decfloat
, y
As Integer64
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As decfloat fac1
, fac2
Dim As Integer64 carry
, digit
, prod
, value
fac1 = x
Call si2fp
(fac2
, y
, digits
) Call fpmul
(fac1
, fac1
, fac2
, digits
) 'check exponents. if either is zero,
'the result is zero
If fac1.exponent
= 0 Or y
= 0 Then 'result is zero...clear fac1. fac1.sign = 0
fac1.exponent = 0
Call set
(fac1.mantissa
, count
, 0) Call NORM_FAC1
(fac1
, digits
) fac1.sign
= fac1.sign
Xor &H8000 'now determine exponent of result.
'as you do...watch for overflow.
er = EXPO_ERR
carry = 0
prod = digit * git(fac1.mantissa, i) + carry
value
= (prod
Mod 100000000) Call set
(fac1.mantissa
, i
, value
) carry = prod \ 100000000
Call RSHIFT_1
(fac1
, digits
) fac1.exponent = fac1.exponent + 1
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 10000000) Call RSHIFT_2
(fac1
, digits
) fac1.exponent = fac1.exponent + 2
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 1000000) Call RSHIFT_3
(fac1
, digits
) fac1.exponent = fac1.exponent + 3
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 100000) Call RSHIFT_4
(fac1
, digits
) fac1.exponent = fac1.exponent + 4
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 10000) Call RSHIFT_5
(fac1
, digits
) fac1.exponent = fac1.exponent + 5
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 1000) Call RSHIFT_6
(fac1
, digits
) fac1.exponent = fac1.exponent + 6
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 100) Call RSHIFT_7
(fac1
, digits
) fac1.exponent = fac1.exponent + 7
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 10) Call RSHIFT_8
(fac1
, digits
) fac1.exponent = fac1.exponent + 8
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
)
Call NORM_FAC1
(fac1
, digits
)
fac1.sign
= fac1.sign
Xor &H8000 result = fac1
Sub fpmul_ll
(result
As decfloat
, x
As decfloat
, y
As Integer64
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As decfloat fac1
, fac2
Dim As Integer64 carry
, digit
, prod
, value
, n0
, n1
fac1 = x
If digit
> 99999999 And digit
< 10000000000000000 Then n0 = digit \ 100000000
Call fpmul_si
(fac2
, fac1
, n1
, digits
) fac1.exponent = fac1.exponent + 8
Call fpmul_si
(fac1
, fac1
, n0
, digits
) Call fpadd
(fac1
, fac1
, fac2
, digits
) If y
< 0 Then fac1.sign
= fac1.sign
Xor &H8000 If digit
> 9999999999999999 Then Call fpmul
(fac1
, fac1
, fac2
, digits
) Call si2fp
(fac2
, y
, digits
)
'check exponents. if either is zero,
'the result is zero
If fac1.exponent
= 0 Or y
= 0 Then 'result is zero...clear fac1. fac1.sign = 0
fac1.exponent = 0
Call set
(fac1.mantissa
, count
, 0) Call NORM_FAC1
(fac1
, digits
) fac1.sign
= fac1.sign
Xor &H8000 'now determine exponent of result.
'as you do...watch for overflow.
er = EXPO_ERR
carry = 0
prod = digit * git(fac1.mantissa, i) + carry
value
= (prod
Mod 100000000) Call set
(fac1.mantissa
, i
, value
) carry = prod \ 100000000
Call RSHIFT_1
(fac1
, digits
) fac1.exponent = fac1.exponent + 1
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 10000000) Call RSHIFT_2
(fac1
, digits
) fac1.exponent = fac1.exponent + 2
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 1000000) Call RSHIFT_3
(fac1
, digits
) fac1.exponent = fac1.exponent + 3
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 100000) Call RSHIFT_4
(fac1
, digits
) fac1.exponent = fac1.exponent + 4
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 10000) Call RSHIFT_5
(fac1
, digits
) fac1.exponent = fac1.exponent + 5
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 1000) Call RSHIFT_6
(fac1
, digits
) fac1.exponent = fac1.exponent + 6
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 100) Call RSHIFT_7
(fac1
, digits
) fac1.exponent = fac1.exponent + 7
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
* 10) Call RSHIFT_8
(fac1
, digits
) fac1.exponent = fac1.exponent + 8
Call set
(fac1.mantissa
, 0, git
(fac1.mantissa
, 0) + carry
)
Call NORM_FAC1
(fac1
, digits
)
fac1.sign
= fac1.sign
Xor &H8000 result = fac1
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As decfloat r
, r2
, one
, n
n = m
l
= Log((NUM_DIGITS
+ 8) * 0.0625) * 1.5 ex
= (n.exponent
And &H7FFFFFFF) - BIAS
- 1 ex = 0
ts
= Str$(git
(n.mantissa
, 0)) ts
= Str$(git
(n.mantissa
, 1)) v = v + ts
x = 10
ex = ex - 1
ex = (-1) - ex
x = 1 / x
r.exponent = ex + BIAS + 1
Call set
(one.mantissa
, 0, 10000000) one.exponent = BIAS + 1
r2 = r
prec = 3
prec = 2 * prec - 1
Call fpmul
(r2
, n
, r
, prec
) Call fpsub
(r2
, one
, r2
, prec
) Call fpmul
(r2
, r
, r2
, prec
) Call fpadd
(r
, r
, r2
, prec
) result = r
wx = ((w(j - 1) * 10000 + w(j)) * 10000 + w(j + 1)) * 10000
RealW# = wx
w(j) = w(j) - q * d(j - ka + 2)
w(ka) = w(ka) + w(ka - 1) * 10000
w(ka - 1) = q
carry = ((-w(j) - 1) \ 10000) + 1
carry = -(w(j) \ 10000)
carry = 0
w(j) = w(j) + carry * 10000
w(j - 1) = w(j - 1) - carry
Sub fpdiv
(result
As decfloat
, x
As decfloat
, y
As decfloat
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As decfloat fac1
, fac2
fac1 = x
fac2 = y
If fac2.exponent
= 0 Then ' if fac2 = 0, return ' a divide-by-zero error and
' bail out.
Call set
(fac1.mantissa
, i
, 99999999) fac1.exponent = 99999 + BIAS + 1
er = DIVZ_ERR
result = fac1
ElseIf fac1.exponent
= 0 Then 'fact1=0, just return er = 0
result = fac1
'check to see if fac2 is a power of ten
is_power_of_ten = 0
If git
(fac2.mantissa
, 0) = 10000000 Then is_power_of_ten = 1
If git
(fac2.mantissa
, i
) <> 0 Then is_power_of_ten = 0
'if fac2 is a power of ten then all we need to do is to adjust the sign and exponent and we are finished
fac1.sign
= fac1.sign
Xor fac2.sign
fac1.exponent = fac1.exponent - fac2.exponent + BIAS + 1
result = fac1
Dim As Double result
(1 To 2 * digits
+ 3), n
(1 To 2 * digits
+ 3), d
(1 To 2 * digits
+ 3)
n(2 * j + 2) = git(fac1.mantissa, j) \ 10000
n
(2 * j
+ 3) = git
(fac1.mantissa
, j
) Mod 10000 d(2 * j + 2) = git(fac2.mantissa, j) \ 10000
d
(2 * j
+ 3) = git
(fac2.mantissa
, j
) Mod 10000 n
(1) = (fac1.exponent
And &H7FFFFFFF) - BIAS
- 1 d
(1) = (fac2.exponent
And &H7FFFFFFF) - BIAS
- 1 w(j) = 0
w(1) = n(1) - d(1) + 1
w(2) = 0
w(j + 1) = n(j)
xd = (d(2) * b + d(3)) * b + d(4) + d(5) / b
laststep = t + 2
xn = RealW(w(), (stp + 2))
last
= min
(stp
+ t
+ 1, UBound(w
)) Call subtract
(w
(), q
, d
(), (stp
+ 2), last
) Call normalize
(w
(), (stp
+ 2), q
) Call finalnorm
(w
(), (laststep
+ 1)) If w
(2) <> 0 Then laststep
= laststep
- 1 rund = w(laststep + 1) / b
If rund
>= 0.5 Then w
(laststep
) = w
(laststep
) + 1 result(j) = w(j + 1)
result(j) = w(j)
If w
(2) = 0 Then result
(1) = w
(1) - 1 Else result
(1) = w
(1)
Call set
(fac1.mantissa
, j
, result
(2 * j
+ 2) * 10000 + result
(2 * j
+ 3)) Call NORM_FAC1
(fac1
, digits
) fac1.exponent = (result(1) + BIAS)
fac1.sign
= fac1.sign
Xor fac2.sign
result = fac1
' sqrt(num)
Sub fpsqr
(result
As decfloat
, num
As decfloat
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As decfloat r
, r2
, tmp
, n
, half
l
= Log((NUM_DIGITS
+ 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
n = num
Call si2fp
(tmp
, 0, digits
) If fpcmp
(n
, tmp
, digits
) = 0 Then result = r
Call si2fp
(tmp
, 1, digits
) If fpcmp
(n
, tmp
, digits
) = 0 Then result = r
Call si2fp
(tmp
, 0, digits
) If fpcmp
(n
, tmp
, digits
) < 0 Then result = r
'=====================================================================
'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
ex
= (n.exponent
And &H7FFFFFFF) - BIAS
- 1 ex = 0
ts
= Trim$
(Str$(git
(n.mantissa
, 0))) ts
= Trim$
(Str$(git
(n.mantissa
, 1))) v = v + ts
result = r
x = x * 10
ex = ex - 1
x
= Sqr(x
) 'approximation r.exponent = ex \ 2 + BIAS + 1
Call set
(half.mantissa
, k
, 0) Call set
(half.mantissa
, 0, 50000000) half.exponent = BIAS
half.sign = 0
'=====================================================================
'Newton-Raphson method
prec = 3
prec = 2 * prec - 1
Call fpdiv
(tmp
, n
, r
, prec
) Call fpadd
(r2
, r
, tmp
, prec
) Call fpmul
(r
, r2
, half
, prec
) result = r
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As Unsigned Integer64 carry
, remder
Dim As Integer64 i
, divisor
Dim As Integer64 quotient
remder = 0
fac1 = num
Print "error: divisor = 0" Print "error: divisor too large"
quotient = git(fac1.mantissa, i) + remder * 100000000
remder
= quotient
Mod divisor
Call set
(fac1.mantissa
, i
, quotient \ divisor
) quotient = remder * 100000000
quotient = quotient \ divisor
carry = git(fac1.mantissa, 0)
Call LSHIFT_8
(fac1
, digits
) fac1.exponent = fac1.exponent - 8
Call set
(fac1.mantissa
, digits
, git
(fac1.mantissa
, digits
) + quotient
) Call LSHIFT_7
(fac1
, digits
) fac1.exponent = fac1.exponent - 7
Call set
(fac1.mantissa
, digits
, git
(fac1.mantissa
, digits
) + quotient \
10) Call LSHIFT_6
(fac1
, digits
) fac1.exponent = fac1.exponent - 6
Call set
(fac1.mantissa
, digits
, git
(fac1.mantissa
, digits
) + quotient \
100) Call LSHIFT_5
(fac1
, digits
) fac1.exponent = fac1.exponent - 5
Call set
(fac1.mantissa
, digits
, git
(fac1.mantissa
, digits
) + quotient \
1000) Call LSHIFT_4
(fac1
, digits
) fac1.exponent = fac1.exponent - 4
Call set
(fac1.mantissa
, digits
, git
(fac1.mantissa
, digits
) + quotient \
10000) Call LSHIFT_3
(fac1
, digits
) fac1.exponent = fac1.exponent - 3
Call set
(fac1.mantissa
, digits
, git
(fac1.mantissa
, digits
) + quotient \
100000) Call LSHIFT_2
(fac1
, digits
) fac1.exponent = fac1.exponent - 2
Call set
(fac1.mantissa
, digits
, git
(fac1.mantissa
, digits
) + quotient \
1000000) Call LSHIFT_1
(fac1
, digits
) fac1.exponent = fac1.exponent - 1
Call set
(fac1.mantissa
, digits
, git
(fac1.mantissa
, digits
) + quotient \
10000000)
'NORM_FAC1(fac1)
fac1.sign
= fac1.sign
Xor &H8000 result = fac1
'fractional part of num
Sub fpfrac
(result
As decfloat
, num
As decfloat
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Call fpsub
(n
, num
, n
, digits
) result = n
'returns the positive of n
Sub fpabs
(result
As decfloat
, n
As decfloat
) x = n
x.sign = 0
result = x
'changes the sign of n, if n is positive then n will be negative & vice versa
Sub fpneg
(result
As decfloat
, n
As decfloat
) x = n
x.sign
= x.sign
Xor &H8000 result = x
'returns the negative of n regardless of the sign of n
Sub fpnegative
(result
As decfloat
, n
As decfloat
) x = n
x.sign = &H8000
result = x
Sub fpfmod
(quotient
As decfloat
, f_mod
As decfloat
, num
As decfloat
, denom
As decfloat
, digits
As Long) If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As decfloat q
, fm
'make copy in case the destination and source are the same Call fpdiv
(fm
, num
, denom
, digits
):
Call fpfix
(q
, fm
) Call fpsub
(fm
, fm
, q
, digits
):
Call fpmul
(fm
, fm
, denom
, digits
) quotient = q
f_mod = fm
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Call si2fp
(ep
, 1, digits
) ep.exponent = (-(NUM_DIGITS) + BIAS + 1)
result = ep
Sub fpipow
(result
As decfloat
, x
As decfloat
, e
As Integer64
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
'take x to an Long power
c = 0
y = x
z.sign = 0
z.exponent = (BIAS + 1)
Call set
(z.mantissa
, 0, 10000000) Call set
(z.mantissa
, i
, 0) n = n \ 2
Call fpmul
(y
, y
, y
, digits
) c = c + 1
n = n - 1
Call fpmul
(z
, y
, z
, digits
) c = c + 1
result = z
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Call fpmul_si
(f
, f
, i
, digits
) result = f
Sub fplogTaylor
(result
As decfloat
, x
As decfloat
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
'taylor series
'====================Log Guard==================
zero.sign = 0
zero.exponent = 0
Call set
(zero.mantissa
, i
, 0) '=============================================
Dim As decfloat XX
, Term
, Accum
, x9
, tmp
, tmp2
Dim As decfloat T
, B
, one
, Q
, two
one.sign = 0
one.exponent = (BIAS + 1)
Call set
(one.mantissa
, 0, 10000000) two.sign = 0
two.exponent = (BIAS + 1)
Call set
(two.mantissa
, 0, 20000000) Call set
(one.mantissa
, i
, 0) Call set
(two.mantissa
, i
, 0) x9 = x
If fpcmp
(x
, one
, digits
) < 0 Then invflag = 1
Call fpdiv
(x9
, one
, x9
, digits
) Call fpsub
(T
, x9
, one
, digits
) Call fpadd
(B
, x9
, one
, digits
) Call fpdiv
(Accum
, T
, B
, digits
) Call fpdiv
(Q
, T
, B
, digits
) tmp = Q
Call fpmul
(XX
, Q
, Q
, digits
) c = 1
c = c + 2
tmp2 = tmp
Call fpmul
(Q
, Q
, XX
, digits
) Call fpdiv_si
(Term
, Q
, c
, digits
) Call fpadd
(Accum
, tmp
, Term
, digits
) Call fpmul_si
(Accum
, Accum
, 2, digits
) result = Accum
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
'====================Log Guard==================
Dim As decfloat g
, one
, zero
zero.sign = 0
zero.exponent = 0
Call set
(zero.mantissa
, i
, 0) Call si2fp
(one
, 1, digits
) '=============================================
Dim As decfloat approx
, ans
, logx
approx = x
factor = 4096
Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fpsqr
(approx
, approx
, digits
) Call fplogTaylor
(logx
, approx
, digits
) Call fpmul_si
(ans
, logx
, factor
, digits
) result = ans
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
'taylor series
Dim As decfloat fac
, x9
, temp
, accum
, p
, term
Call si2fp
(temp
, 0, digits
)
fac.sign = 0
fac.exponent = (BIAS + 1)
Call set
(fac.mantissa
, 0, 10000000) Call set
(fac.mantissa
, i
, 0) c = 1
Call fpdiv_si
(x9
, x
, 8192, digits
) 'fpdiv_si(x, 67108864) ' p = x9
Call fpadd
(accum
, fac
, x9
, digits
) '1 + x
c = c + 1
temp = accum
Call fpdiv_si
(fac
, fac
, c
, digits
) Call fpmul
(p
, p
, x9
, digits
) Call fpmul
(term
, p
, fac
, digits
) Call fpadd
(accum
, temp
, term
, digits
) Call fpmul
(accum
, accum
, accum
, digits
) result = accum
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
' x
' cf = ------------------------------
' x^2
' 2 + -------------------------
' x^2
' 6 + ----------------------
' x^2
' 10 + ------------------
' x^2
' 14 + --------------
' x^2
' 18 + ---------
' x^2
' 22 + -----
' 26 + ...
' (1 + cf)
' exp(x) = ----------
' (1 - cf)
Dim As decfloat s
, x1
, x2
, tmp
, tmp2
Call fpdiv_si
(x1
, x
, 8192, digits
) ' 2^13 Call fpmul
(x2
, x1
, x1
, digits
) Call si2fp
(s
, 4 * n
+ 6, digits
) Call si2fp
(tmp
, 4 * i
+ 2, digits
) Call fpdiv
(s
, x2
, s
, digits
) Call fpadd
(s
, s
, tmp
, digits
) Call fpdiv
(s
, x1
, s
, digits
) Call si2fp
(tmp
, 1, digits
) tmp2 = tmp
Call fpadd
(tmp
, tmp
, s
, digits
) Call fpsub
(tmp2
, tmp2
, s
, digits
) Call fpdiv
(s
, tmp
, tmp2
, digits
) Call fpmul
(s
, s
, s
, digits
) result = s
Sub fplog_r
(result
As decfloat
, x
As decfloat
, digits_in
As Long) digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
'====================Log Guard==================
Dim As decfloat x2
, y
, tmp
tmp.sign = 0
tmp.exponent = (BIAS + 1)
Call set
(tmp.mantissa
, 0, 10000000) Call set
(tmp.mantissa
, i
, 0) ex
= (x.exponent
And &H7FFFFFFF) - BIAS
- 1 t = git(x.mantissa, 0) + git(x.mantissa, 1) / 100000000
t = t / 10000000
n
= Log((digits
+ 1) * 0.5) * 1.5
prec = 3
terms = 2
prec = 2 * prec - 1
Call fpexp_cf
(y
, x2
, terms
, prec
) Call fpsub
(tmp
, y
, x
, prec
) Call fpdiv
(tmp
, tmp
, y
, prec
) Call fpsub
(x2
, x2
, tmp
, prec
) prec = 2 * prec - 1
terms = 2 * terms + 1
Call fpexp_cf
(y
, x2
, terms
, prec
) Call fpsub
(tmp
, y
, x
, prec
) Call fpdiv
(tmp
, tmp
, y
, prec
) Call fpsub
(x2
, x2
, tmp
, prec
) prec = 2 * prec - 1
terms = 2 * terms + 1
Call fpexp_cf
(y
, x2
, terms
, prec
) Call fpsub
(tmp
, y
, x
, prec
) Call fpdiv
(tmp
, tmp
, y
, prec
) Call fpsub
(x2
, x2
, tmp
, prec
) prec = 2 * prec - 1
terms = 2 * terms + 1
Call fpexp_cf
(y
, x2
, terms
, prec
) Call fpsub
(tmp
, y
, x
, prec
) Call fpdiv
(tmp
, tmp
, y
, prec
) Call fpsub
(x2
, x2
, tmp
, prec
) prec = 2 * prec - 1
terms = 2 * terms + 1
Call fpexp_cf
(y
, x2
, terms
, prec
) Call fpsub
(tmp
, y
, x
, prec
) Call fpdiv
(tmp
, tmp
, y
, prec
) Call fpsub
(x2
, x2
, tmp
, prec
) prec = 2 * prec - 1
terms = 2 * terms + 1
Call fpexp_cf
(y
, x2
, terms
, prec
) Call fpsub
(tmp
, y
, x
, prec
) Call fpdiv
(tmp
, tmp
, y
, prec
) Call fpsub
(x2
, x2
, tmp
, prec
) result = x2
Sub fppow
(result
As decfloat
, lhs
As decfloat
, rhs
As decfloat
) Call fplog
(lhs2
, lhs
, NUM_DWORDS
) Call fpmul
(lhs2
, lhs2
, rhs
, NUM_DWORDS
) Call fpexp
(lhs2
, lhs2
, NUM_DWORDS
) result = lhs2
digits = digits_in
If digits
> NUM_DWORDS
Then digits
= NUM_DWORDS
Dim As decfloat ry
, tmp
, tmp2
l
= Log((NUM_DIGITS
+ 8) * 0.0625) * 1.5 ex
= (x.exponent
And &H7FFFFFFF) - BIAS
- 1 t = git(x.mantissa, 0) + git(x.mantissa, 1) / 100000000
t = t / 10000000
t2
= Log(t
) / p
+ Log(10) * ex
/ p
prec = 3
Call fpipow
(tmp
, ry
, p
- 1, prec
) Call fpdiv
(tmp
, x
, tmp
, prec
) Call fpmul_si
(tmp2
, ry
, p
- 1, prec
) Call fpadd
(tmp2
, tmp2
, tmp
, prec
) Call fpdiv_si
(ry
, tmp2
, p
, prec
) prec = 2 * prec - 1
Call fpipow
(tmp
, ry
, p
- 1, prec
) Call fpdiv
(tmp
, x
, tmp
, prec
) Call fpmul_si
(tmp2
, ry
, p
- 1, prec
) Call fpadd
(tmp2
, tmp2
, tmp
, prec
) Call fpdiv_si
(ry
, tmp2
, p
, prec
) result = ry
Sub pi_brent_salamin
(pi_bs
As decfloat
, digits_in
As Unsigned
Long) digits = digits_in
If digits
> NUMBER_OF_DIGITS
Then digits
= NUMBER_OF_DIGITS
Dim As decfloat c0
, c1
, c2
, c05
Dim As decfloat a
, b
, sum
Dim As decfloat ak
, bk
, ck
Dim As decfloat pow2
, tmp
limit2 = -digits + BIAS + 1
Call si2fp
(c0
, 0, NUMBER_OF_DIGITS
): ak
= c0: bk
= c0: ab
= c0: asq
= c0
Call si2fp
(c1
, 1, NUMBER_OF_DIGITS
): a
= c1: ck
= c1: pow2
= c1
Call si2fp
(c2
, 2, NUMBER_OF_DIGITS
): b
= c2
Call str2fp
(c05
, ".5"): sum
= c05
Call si2fp
(pi_bs
, 3, NUMBER_OF_DIGITS
)
Call fpsqr
(b
, b
, NUMBER_OF_DIGITS
) Call fpdiv
(b
, c1
, b
, NUMBER_OF_DIGITS
) While fpcmp
(ck
, c0
, NUMBER_OF_DIGITS
) <> 0 And ck.exponent
> limit2
Call fpadd
(ak
, a
, b
, NUMBER_OF_DIGITS
) Call fpmul
(ak
, c05
, ak
, NUMBER_OF_DIGITS
) Call fpmul
(ab
, a
, b
, NUMBER_OF_DIGITS
) Call fpsqr
(bk
, ab
, NUMBER_OF_DIGITS
) Call fpmul
(asq
, ak
, ak
, NUMBER_OF_DIGITS
) Call fpsub
(ck
, asq
, ab
, NUMBER_OF_DIGITS
) Call fpmul
(pow2
, pow2
, c2
, NUMBER_OF_DIGITS
) Call fpmul
(tmp
, pow2
, ck
, NUMBER_OF_DIGITS
) Call fpsub
(sum
, sum
, tmp
, NUMBER_OF_DIGITS
) a = ak: b = bk
Call fpdiv
(tmp
, asq
, sum
, NUMBER_OF_DIGITS
) Call fpmul
(pi_bs
, c2
, tmp
, NUMBER_OF_DIGITS
)
Sub fpsin
(result
As decfloat
, x
As decfloat
, digits_in
As Unsigned
Long) Dim As decfloat XX
, Term
, Accum
, p
, temp2
, fac
, x_2
Dim As decfloat pi2
, circ
, Ab
precision = digits_in
x_2 = x
pi2 = pi_dec
Call fpmul_si
(circ
, pi2
, 2, precision
) If fpcmp
(Ab
, circ
, precision
) > 0 Then '======== CENTRALIZE ==============
'floor/ceil to centralize
Dim As decfloat tmp
, tmp2
pi2 = pi_dec
Call fpmul_si
(pi2
, pi2
, 2, precision
) 'got 2*pi Call fpdiv
(tmp
, x_2
, pi2
, precision
) tmp2 = tmp
Call fpfix
(tmp
, tmp
) 'int part Call fpsub
(tmp
, tmp2
, tmp
, precision
) 'frac part Call fpmul
(tmp
, tmp
, pi2
, precision
) x_2 = tmp
lm = precision
limit2
= Int(-0.45344993886092585968 + 0.022333002852398072433 * lm
+ 5.0461814408333079844E-7 * lm
* lm
- 4.2338453039804235772E-11 * lm
* lm
* lm
)
Call si2fp
(factor
, 5, precision
) Call fpipow
(factor
, factor
, limit2
, precision
) Call fpdiv
(x_2
, x_2
, factor
, precision
) 'x_=x_/5^limit2 '==================================
Accum = x_2
Call si2fp
(fac
, 1, precision
) p = x_2
Call fpmul
(XX
, x_2
, x_2
, precision
) c = c + 2
temp2 = Accum
Call fpmul_si
(fac
, fac
, c
* (c
- 1), precision
) Call fpmul
(p
, p
, XX
, precision
) Call fpdiv
(Term
, p
, fac
, precision
) Call fpsub
(Accum
, temp2
, Term
, precision
) Call fpadd
(Accum
, temp2
, Term
, precision
) Loop Until fpcmp
(Accum
, temp2
, precision
) = 0 'multiply the result by 5^limit2
Call fpmul
(p
, Accum
, Accum
, precision
) Call fpmul
(temp2
, Accum
, p
, precision
) '*** sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
Call fpmul_si
(Accum
, Accum
, 5, precision
) Call fpmul_si
(Term
, temp2
, 20, precision
) Call fpmul_si
(XX
, temp2
, 16, precision
) Call fpmul
(XX
, XX
, p
, precision
) Call fpsub
(Accum
, Accum
, Term
, precision
) Call fpadd
(Accum
, Accum
, XX
, precision
) result = Accum
Sub fpcos
(result
As decfloat
, z
As decfloat
, digits_in
As Unsigned
Long) precision = digits_in
Call fpdiv_si
(pi2
, pi_dec
, 2, precision
) Call fpsub
(x_2
, pi2
, z
, precision
)
Call fpsin
(result
, x_2
, precision
)
Sub fptan
(result
As decfloat
, z
As decfloat
, digits_in
As Unsigned
Long) Dim As decfloat x_2
, s
, c
precision = digits_in
x_2 = z
Call fpsin
(s
, x_2
, precision
) x_2 = z
Call fpcos
(c
, x_2
, precision
) Call fpdiv
(result
, s
, c
, precision
)
Sub fpatn
(result
As decfloat
, x
As decfloat
, digits_in
As Unsigned
Long) precision = digits_in
Dim As decfloat XX
, Term
, Accum
, strC
, x_2
, mt
, mt2
, p
Dim As decfloat decnum
, one
, decnum2
, factor
decnum2 = x
decnum2.sign = 0
Call si2fp
(one
, 1, precision
) If fpcmp
(decnum2
, one
, precision
) = 0 Then Call fpdiv_si
(result
, pi_dec
, 4, precision
) result.sign = x.sign
decnum2.sign = x.sign
Call si2fp
(factor
, SHL
(2, limit2
- 1), precision
) Call fpmul
(decnum
, decnum2
, decnum2
, precision
) Call fpadd
(decnum
, decnum
, one
, precision
) Call fpsqr
(decnum
, decnum
, precision
) Call fpadd
(decnum
, decnum
, one
, precision
) Call fpdiv
(decnum
, decnum2
, decnum
, precision
) decnum2 = decnum
mt = decnum
x_2 = decnum
p = decnum
Call fpmul
(XX
, x_2
, x_2
, precision
) c = c + 2
mt2 = mt
Call si2fp
(strC
, c
, precision
) Call fpmul
(p
, p
, XX
, precision
) Call fpdiv
(Term
, p
, strC
, precision
) Call fpsub
(Accum
, mt
, Term
, precision
) Call fpadd
(Accum
, mt
, Term
, precision
) Call fpmul
(result
, factor
, mt
, precision
)
Sub fpasin
(result
As decfloat
, x
As decfloat
, digits_in
As Unsigned
Long) precision = digits_in
' ARCSIN = ATN(x / SQR(-x * x + 1))
'============= ARCSIN GUARD =========
num = fp2dbl(x)
'========================
Dim As decfloat one
, T
, B
, term1
, minusone
Call si2fp
(one
, 1, precision
) Call si2fp
(minusone
, -1, precision
) T = x
Call fpmul
(B
, x
, x
, precision
) 'x*x 'for 1 and -1
If fpcmp
(B
, one
, precision
) = 0 Then Dim As decfloat two
, atn1
Call si2fp
(two
, 2, precision
) Call fpatn
(atn1
, one
, precision
) If fpcmp
(x
, minusone
, precision
) = 0 Then Call fpmul
(two
, two
, atn1
, precision
) Call fpmul
(result
, two
, minusone
, precision
) Call fpmul
(result
, two
, atn1
, precision
) Call fpsub
(B
, one
, B
, precision
) '1-x*x Call fpsqr
(B
, B
, precision
) 'sqr(1-x*x) Call fpdiv
(term1
, T
, B
, precision
) Call fpatn
(result
, term1
, precision
)
Sub fpacos
(result
As decfloat
, x
As decfloat
, digits_in
As Unsigned
Long) precision = digits_in
Dim As decfloat one
, minusone
, two
, atn1
, tail
, T
, B
, term1
, atnterm1
',_x,temp 'ARCCOS = ATN(-x / SQR(-x * x + 1)) + 2 * ATN(1)
'============= ARCCOS GUARD =========
num = fp2dbl(x)
'========================
Call si2fp
(one
, 1, precision
) Call si2fp
(minusone
, -1, precision
) Call si2fp
(two
, 2, precision
) Call fpatn
(atn1
, one
, precision
) Call fpmul
(tail
, two
, atn1
, precision
) '2*atn(1) Call fpmul
(T
, minusone
, x
, precision
) '-x Call fpmul
(B
, x
, x
, precision
) 'x*x If fpcmp
(B
, one
, precision
) = 0 Then 'for 1 and -1
If fpcmp
(x
, minusone
, precision
) = 0 Then Call si2fp
(four
, 4, precision
) Call fpmul
(result
, four
, atn1
, precision
) Call si2fp
(result
, 0, precision
) Call fpsub
(B
, one
, B
, precision
) '1-x*x Call fpsqr
(B
, B
, precision
) 'sqr(1-x*x) Call fpdiv
(term1
, T
, B
, precision
) Call fpatn
(atnterm1
, term1
, precision
) Call fpadd
(result
, atnterm1
, tail
, precision
)
Show Posts