about bignum
about bignum
I am still working on my big-float package and while examining the multiplication routine by Stephan Brunker https://github.com/stephanbrunker/big_integer I noticed that he was assigning a ulongint to a ulong, it turns out that the least significant portion of the ulongint is stored in the ulong so he saves from having to do a mod operation, he uses base 2^32 I think
he's multiplication routine is over 2 times faster than mine partially because of the larger base and the trick I mentioned, I plan on rewriting my bigfloat package using base 2^32 to take advantage of this trick, I also plan on implementing the toom 3 algorithm
he's multiplication routine is over 2 times faster than mine partially because of the larger base and the trick I mentioned, I plan on rewriting my bigfloat package using base 2^32 to take advantage of this trick, I also plan on implementing the toom 3 algorithm
Re: about bignum
Hi srvaldez.
Many years ago in topic circles, Richard, Albert and myself (to a lesser degree) messed about with bignums.
Richard finally used base 1e9 for his method.
Albert used various bases (as those mults I posted a while back).
I stuck to high school methods using string indexing base 10.
My method could never get up to speed, but I stuck to it anyway, I still post bits and pieces now and then.
Maybe you can test your right/left shifts (> fb ulongint) with this, It looks and tests out OK.
I have kept the test numbers smallish to suit 32 bit gas.
Many years ago in topic circles, Richard, Albert and myself (to a lesser degree) messed about with bignums.
Richard finally used base 1e9 for his method.
Albert used various bases (as those mults I posted a while back).
I stuck to high school methods using string indexing base 10.
My method could never get up to speed, but I stuck to it anyway, I still post bits and pieces now and then.
Maybe you can test your right/left shifts (> fb ulongint) with this, It looks and tests out OK.
I have kept the test numbers smallish to suit 32 bit gas.
Code: Select all
Function plus(_num1 As String,_num2 As String) As String
Static As integer ADDQmod(19)={48,49,50,51,52,53,54,55,56,57,48,49,50,51,52,53,54,55,56,57}
Static As integer ADDbool(19)={0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1}
Var _flag=0,n_=0
Dim As integer addup=Any,addcarry=Any
#macro finish()
answer=Ltrim(answer,"0")
If _flag=1 Then Swap _num2,_num1
Return answer
#endmacro
If Len(_num2)>Len(_num1) Then
Swap _num2,_num1
_flag=1
Endif
Var diff=Len(_num1)-Len(_num2)
Var answer="0"+_num1
addcarry=0
For n_=Len(_num1)-1 To diff Step -1
addup=_num2[n_-diff]+_num1[n_]-96
answer[n_+1]=ADDQmod(addup+addcarry)
addcarry=ADDbool(addup+addcarry)
Next n_
If addcarry=0 Then
finish()
End If
If n_=-1 Then
answer[0]=addcarry+48
finish()
Endif
For n_=n_ To 0 Step -1
addup=_num1[n_]-48
answer[n_+1]=ADDQmod(addup+addcarry)
addcarry=ADDbool(addup+addcarry)
If addcarry=0 Then Exit For
Next n_
answer[0]=addcarry+48
finish()
End Function
' ============== BINARY TO DECIMAL =================================
Function Base10(BinaryNumber As String) As String
Dim As String sum
sum=Left(BinaryNumber,1)
For n As Integer=2 To Len(BinaryNumber)
sum=plus(plus(sum,sum),Chr(BinaryNumber[n-1]))
Next n
Return sum
End Function
' ==================== DECIMAL TO BINARY ==================================
Function base2(DecimalNumber As String) As String
Dim As String starter=DecimalNumber,ans,m,b
Dim As Ubyte main,carry,temp',c
Dim As Integer c
#macro reverse(s)
Scope
Var lens=Len(s)
For n As Integer=0 To Int((lens-1)/2):Swap s[n],s[lens-1-n]:Next
End Scope
#endmacro
#macro div2(s,m,c)
carry=0:ans=s
For z As Integer=0 To Len(s)-1
temp=(s[z]-48+carry)
main=temp Shr 1
carry=(temp And 1) Shl 3 +(temp And 1) Shl 1
ans[z]=main+48
Next z
c= carry\10
m= Ltrim(ans,"0")
#endmacro
Do
c=0
div2(starter,m,c)
b=b+Str(c)
starter=m
Loop Until m=""
reverse(b)
b=Str(m)+b
Return b
End Function
Function shiftr(sd As String,k As Long ) As String
Var b2=base2(sd)
Var l=Len(b2)
Dim As String s=String(k,"0")+b2
var ans=base10(Left(s,l))
Return iif(ans="","0",ans)
End Function
Function shiftl(sd As String,k As Long ) As String
Var b2=base2(sd)
Dim As String s=b2+String(k,"0")
var ans=base10(s)
Return iif(ans="","0",ans)
End Function
dim as integer a,b
print "loop test . . ."
for n as long=1 to 10000
a=20+rnd*10000
b=1+rnd*15
if a shl b<>val(shiftl(str(a),b)) then print n
if a shr b<>val(shiftr(str(a),b)) then print n
next
Print
Print shiftl("2021",2021)
print "done"
sleep
Re: about bignum
Hi dodicat :-)
thank you very much, it's just what I needed
I got started on the 32-bit binary math routines and found out that the bit shifter with the asm code is way too slow so I used fxm suggestion and it's working ok so far, your conversion routines are just what I needed :-)
thank you very much, it's just what I needed
I got started on the 32-bit binary math routines and found out that the bit shifter with the asm code is way too slow so I used fxm suggestion and it's working ok so far, your conversion routines are just what I needed :-)
Re: about bignum
so far I have implemented add, subtract, multiply and divide_by_long, a big hurdle to overcome is the conversion from string to bigfloat and back
I tested the routines by calculating sin(0.5) to 1030 decimal digits verified with the pari/gp calculator
the binary routines were about 38% faster than my decimal routines, I am somewhat disappointed
I tested the routines by calculating sin(0.5) to 1030 decimal digits verified with the pari/gp calculator
the binary routines were about 38% faster than my decimal routines, I am somewhat disappointed
Re: about bignum
just for fun, here's my try at multi-precision binary floating-point math
the input/output is crude but it's better than nothing, was having trouble figuring what was wrong with the conversion to string until it dawned on me that since I use double to extract the digits the value could easily be out of range
the input/output is crude but it's better than nothing, was having trouble figuring what was wrong with the conversion to string until it dawned on me that since I use double to extract the digits the value could easily be out of range
Code: Select all
Extern "c"
Declare Function __builtin_clz (Byval x As Ulong) As Long
Declare Function __builtin_smull_overflow (Byval As Long, Byval As Long, Byref As Long) As Boolean
Declare Function __builtin_smulll_overflow (Byval As Longint, Byval As Longint, Byref As Longint) As Boolean
Declare Function __builtin_umull_overflow (Byval As Ulong, Byval As Ulong, Byref As Ulong) As Boolean
Declare Function __builtin_umulll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As Boolean
Declare Function __builtin_uaddll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As Boolean
End Extern
Const NUMBER_OF_DIGITS = 96
Const NUMBER_OF_BITS = NUMBER_OF_DIGITS*3.321928094887362
#If (NUMBER_OF_BITS Mod 32)<>0
Const NUM_BITS = (NUMBER_OF_BITS\32+1)*32
#Else
Const NUM_BITS = NUMBER_OF_BITS
#Endif
Const NUM_DWORDS = NUM_BITS\32
Const BIAS = 1073741824 '2 ^ 30
' Error definitions
#Define DIVZ_ERR 1 'Divide by zero
#Define EXPO_ERR 2 'Exponent overflow error
#Define EXPU_ERR 3 'Exponent underflow error
Type bigfloat_struct
Declare Constructor ( )
Declare Constructor ( Byref rhs As bigfloat_struct )
Declare Destructor ( )
Declare Operator Let ( Byref rhs As bigfloat_struct )
As Short sign
As Ulong exponent
As Ulong Ptr mantissa
End Type
Constructor bigfloat_struct ( )
If this.mantissa<>0 Then
Print "pointer is not 0"
End(1)
End If
this.mantissa=Callocate((NUM_DWORDS +1), 4)
If this.mantissa=0 Then
Print "unable to allocate memory"
End(4)
End If
this.sign=0
this.exponent=0
End Constructor
Constructor bigfloat_struct ( Byref rhs As bigfloat_struct )
If this.mantissa<>0 Then
Print "pointer is not 0"
End(1)
End If
this.mantissa=Callocate((NUM_DWORDS +1), 4)
If this.mantissa=0 Then
Print "unable to allocate memory"
End(4)
End If
this.sign=rhs.sign
this.exponent=rhs.exponent
For i As Long=0 To NUM_DWORDS
this.mantissa[i]=rhs.mantissa[i]
Next
End Constructor
Operator bigfloat_struct.let ( Byref rhs As bigfloat_struct )
this.sign=rhs.sign
this.exponent=rhs.exponent
For i As Long=0 To NUM_DWORDS
this.mantissa[i]=rhs.mantissa[i]
Next
End Operator
Destructor bigfloat_struct ( )
If this.mantissa=0 Then
Print "unable to de-allocate memory"
End(-4)
Else
Deallocate(this.mantissa)
End If
this.mantissa=0
End Destructor
Dim Shared As bigfloat_struct Log102
Log102.sign=0
Log102.exponent=(-2)+BIAS+1
If NUM_DWORDS>=0 Then Log102.mantissa[0]=&h13441350
If NUM_DWORDS>=1 Then Log102.mantissa[1]=&h9F79FEF3
If NUM_DWORDS>=2 Then Log102.mantissa[2]=&h11F12B35
If NUM_DWORDS>=3 Then Log102.mantissa[3]=&h816F922F
Private Function log2_32(Byval value As Ulong) As Long
'https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers
Static tab32(0 To 31) As Const Long = {0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30, 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31}
value Or= Culng(value Shr 1)
value Or= Culng(value Shr 2)
value Or= Culng(value Shr 4)
value Or= Culng(value Shr 8)
value Or= Culng(value Shr 16)
Return tab32(Culng(Culng(value * &h07C4ACDD) Shr 27))
End Function
Private Function log2_64(Byval value As Ulongint) As Long
'https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers
Static tab64(0 To 63) As Const Long ={63, 0, 58, 1, 59, 47, 53, 2, 60, 39, 48, 27, 54, 33, 42, 3, 61, 51, 37, 40, 49, 18, 28, 20, 55, 30, 34, 11, 43, 14, 22, 4, 62, 57, 46, 52, 38, 26, 32, 41, 50, 36, 17, 19, 29, 10, 13, 21, 56, 45, 25, 31, 35, 16, 9, 12, 44, 24, 15, 8, 23, 7, 6, 5}
value Or= Culngint(value Shr 1)
value Or= Culngint(value Shr 2)
value Or= Culngint(value Shr 4)
value Or= Culngint(value Shr 8)
value Or= Culngint(value Shr 16)
value Or= Culngint(value Shr 32)
Return tab64((Culngint((value - (value Shr 1)) * &h07EDD5E59A4E28C2) Shr 58))
End Function
Function ipower (Byval x As Ulongint, Byval e As Ulongint) As Ulongint
'take x to an integer power
Dim As Ulongint z, y, n
y = x
n = e
z = 1
While n > 0
While (n And 1) = 0
n = n \ 2
y = y * y
Wend
n = n - 1
z = y * z
Wend
Return z
End Function
Function shl32(Byval n As Ulong, Byval k As Ubyte, Byref c As Ulong) As Ulong
If k>0 And k<32 Then
Dim As Ulong carry=n
Dim As Ubyte k32=32-k
Asm
mov cl, [k32]
Shr dword Ptr [carry], cl
mov cl, [k]
Shl dword Ptr [n], cl
End Asm
c=carry
End If
Return n
End Function
Function shr32(Byval n As Ulong, Byval k As Ubyte, Byref c As Ulong) As Ulong
If k>0 And k<32 Then
Dim As Ulong carry=n
Dim As Ubyte k32=32-k
Asm
mov cl, [k32]
Shl dword Ptr [carry], cl
mov cl, [k]
Shr dword Ptr [n], cl
End Asm
c=carry
End If
Return n
End Function
Sub shiftl(Byref fac1 As bigfloat_struct, Byval k As Long, Byval dwords As Long=NUM_DWORDS)
If k>0 And k<32 Then
Dim As Long i
Dim As Ulong n, carry, c=0
Dim As Long k32=32-k
For i=dwords To 0 Step -1
n=fac1.mantissa[i]
carry=n
carry=Cast(Ulong, (carry Shr k32))
n=Cast(Ulong, (n Shl k))
fac1.mantissa[i]=n+c
c=carry
Next
Elseif k=32 Then
Dim As Long i
For i=0 To dwords-1
fac1.mantissa[i]=fac1.mantissa[i+1]
Next
fac1.mantissa[dwords]=0
End If
End Sub
Sub shiftr(Byref fac1 As bigfloat_struct, Byval k As Long, Byval dwords As Long=NUM_DWORDS)
If k>0 And k<32 Then
Dim As Long i
Dim As Ulong n, carry, c=0
Dim As Long k32=32-k
For i=0 To dwords
n=fac1.mantissa[i]
carry=n
carry=Cast(Ulong, (carry Shl k32))
n=Cast(Ulong, (n Shr k))
fac1.mantissa[i]=c+n
c=carry
Next
Elseif k=32 Then
Dim As Long i
For i=dwords To 1 Step -1
fac1.mantissa[i]=fac1.mantissa[i-1]
Next
fac1.mantissa[0]=0
End If
End Sub
Function fpcmp(Byref x As bigfloat_struct, Byref y As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As Long c, i
If x.sign < y.sign Then fpcmp = -1: Exit Function
If x.sign > y.sign Then fpcmp = 1: Exit Function
If x.sign=y.sign Then
If x.exponent=y.exponent Then
For i=0 To dwords
c=x.mantissa[i]-y.mantissa[i]
If c<>0 Then Exit For
Next
If c<0 Then Return -1
If c=0 Then Return 0
If c>0 Then Return 1
End If
If x.exponent<y.exponent Then Return -1
If x.exponent>y.exponent Then Return 1
End If
End Function
Function NORM_FAC1(Byref fac1 As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As Integer
If dwords>NUM_DWORDS Then dwords=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.
Dim As Integer i,er=0,f=0
For i=0 To dwords
If fac1.mantissa[i]>0 Then f=1
Next
If f=0 Then
fac1.exponent=0
fac1.sign=0
Exit Function
End If
While fac1.mantissa[0]=0
shiftl(fac1, 32, dwords)
fac1.exponent-=32
If fac1.exponent=0 Then
NORM_FAC1=EXPU_ERR
Exit Function
End If
Wend
i=__builtin_clz(fac1.mantissa[0])-3
If i>0 Then
shiftl(fac1, i, dwords)
fac1.exponent-=i
End If
'if the highmost Bit in fac1_man is nonzero,
'shift the mantissa right 1 Bit and
'increment the exponent
If fac1.mantissa[0]>(&h1ffffffful)Then
While fac1.mantissa[0]>&h1ffffffful
shiftr(fac1, 1, dwords)
fac1.exponent+=1
Wend
Elseif fac1.mantissa[0]<&h10000000ul Then
/' the following will probably never be executed '/
'now shift fac1_man 1 to the left until a
'nonzero Bit appears in the next-to-highest
'Bit of fac1_man. decrement exponent for
'each shift.
While fac1.mantissa[0]<&h10000000ul
shiftl(fac1, 1, dwords)
fac1.exponent-=1
Wend
End If
'check for overflow/underflow
If fac1.exponent<0 Then
NORM_FAC1=EXPO_ERR
End If
End Function
Function si2fp(Byval m As Long, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct fac1
Dim As Long n=Abs(m)
If m=0 Then
Return fac1
End If
fac1.mantissa[0]=n
NORM_FAC1(fac1, NUM_DWORDS)
fac1.exponent=log2_32(n)+BIAS+1
If m<0 Then
fac1.sign=&H8000
Else
fac1.sign=0
End If
Return fac1
End Function
Function ui2fp(Byval m As Ulong, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct fac1
If m=0 Then
Return fac1
End If
fac1.mantissa[0]=m
NORM_FAC1(fac1, NUM_DWORDS)
fac1.exponent=log2_32(m)+BIAS+1
Return fac1
End Function
Sub fpadd_aux(Byref fac1 As bigfloat_struct, Byref fac2 As bigfloat_struct, Byval dwords As Long=NUM_DWORDS)
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As Long c, i
Dim As Longint v
c=0
For i=dwords To 1 Step -1
v=Clngint(fac2.mantissa[i])+Clngint(fac1.mantissa[i])+c
If v>(&h100000000ull-1) Then
v=v-&h100000000ull
c=1
Else
c=0
End If
fac1.mantissa[i]=v
Next
v=Clngint(fac1.mantissa[0])+Clngint(fac2.mantissa[0])+c
fac1.mantissa[0]=v
NORM_FAC1(fac1, dwords)
End Sub
Sub fpsub_aux(Byref fac1 As bigfloat_struct, Byref fac2 As bigfloat_struct, Byval dwords As Long=NUM_DWORDS)
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As Long c, i
Dim As Longint v
c=0
For i=dwords To 1 Step -1
v=Clngint(fac1.mantissa[i])-clngint(fac2.mantissa[i])-c
If v<0 Then
v=v+&h100000000ull
c=1
Else
c=0
End If
fac1.mantissa[i]=v
Next
v=Clngint(fac1.mantissa[0])-clngint(fac2.mantissa[0])-c
fac1.mantissa[0]=v
NORM_FAC1(fac1, dwords)
End Sub
Function fpadd(Byref x As bigfloat_struct, Byref y As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct fac1,fac2
Dim As Long i, t, c, xsign, ysign
xsign=x.sign:x.sign=0
ysign=y.sign:y.sign=0
c=fpcmp(x, y, dwords)
x.sign=xsign
y.sign=ysign
If c<0 Then
fac1=y
fac2=x
Else
fac1=x
fac2=y
End If
t=fac1.exponent-fac2.exponent
If t<(NUM_BITS+32) Then
'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 dwords, the result is already in FAC1.
If t>0 And t<(NUM_BITS+32) Then 'shift
i=t\32
While i>0
shiftr(fac2, 32, dwords)
t-=32
i-=1
Wend
'While t>0
If t>0 Then shiftr(fac2, t, dwords)
' t-=1
'Wend
End If
'See if the signs of the two numbers
'are the same. If so, add; if not, subtract.
If fac1.sign=fac2.sign Then 'add
fpadd_aux(fac1,fac2, dwords)
Else
fpsub_aux(fac1,fac2, dwords)
End If
Endif
Return fac1
End Function
Function fpsub(Byref x As bigfloat_struct, Byref y As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct fac1,fac2
Dim As Long s
fac1=x
fac2=y
fac2.sign=fac2.sign Xor &h8000
fac1=fpadd(fac1,fac2, dwords)
Return fac1
End Function
Function fpmul(Byref x As bigfloat_struct, Byref y As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct fac1,fac2
Dim As Integer i, j, ex, er, den, num
Dim As Ulongint digit, carry, tmp, prod
Dim As Ulong fac3(0 To 2*dwords+1)
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
For i=0 To dwords
fac1.mantissa[i]=0
Next
'NORM_FAC1(fac1)
Return fac1
Else
If ex<0 Then
er=EXPO_ERR
Return fac1 'Exit Function
End If
'clear fac3 mantissa
For i=0 To dwords
fac3(i)=0
Next
den=dwords
While fac2.mantissa[den]=0
den-=1
Wend
num=dwords
While fac1.mantissa[num]=0
num-=1
Wend
If num<den Then
Swap fac1, fac2
'fac1=y
'fac2=x
Swap den, num
End If
For j=den To 0 Step -1
carry=0
digit=fac2.mantissa[j]
For i=num To 0 Step -1
/' prod=fac3(i+j+1)+digit*fac1.mantissa[i]+carry '/
tmp=fac1.mantissa[i]
If __builtin_umulll_overflow (digit, tmp, prod) Then Print "mul overflow"
tmp=fac3(i+j+1)
If __builtin_uaddll_overflow (tmp, prod, prod) Then Print "add overflow"
If __builtin_uaddll_overflow (carry, prod, prod) Then Print "add overflow"
carry=prod Shr 32'\&H100000000
fac3(i+j+1)=prod '(prod mod &H100000000)
Next
fac3(j)=carry
Next
For i=0 To dwords
fac1.mantissa[i]=fac3(i)
Next
End If
'now determine exponent of result.
'as you do...watch for overflow.
'ex=fac2.exponent-BIAS+fac1.exponent
'fac1.exponent=ex
ex=(fac2.exponent And &h7FFFFFFFul)-BIAS+1
ex=ex+(fac1.exponent And &h7FFFFFFFul)-BIAS+1
fac1.exponent=ex+BIAS+1
'determine the sign of the product
fac1.sign=fac1.sign Xor fac2.sign
NORM_FAC1(fac1, dwords)
Return fac1
End Function
Function fpmul_si(Byref x As bigfloat_struct, Byval y As Long, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct fac1,fac2
fac1=x
fac2=ui2fp(y, dwords)
fac1=fpmul(fac1, fac2, dwords)
Return fac1
End Function
Function fpdiv_si(Byref num As bigfloat_struct, Byval den As Long, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct fac1
Dim As Ulongint carry, remder=0
Dim As Longint i, divisor
Dim As Longint quotient
Dim As Long j
divisor=Abs(den)
fac1=num
If divisor=1 Then
Return fac1
End If
If divisor = 0 Then
Print "error: divisor = 0"
Return fac1 'Exit function
End If
If divisor>2147483647 Then
Print "error: divisor too large"
Return fac1 'Exit function
End If
For i = 0 To dwords
quotient = fac1.mantissa[i] + remder * &h100000000ull
remder = quotient Mod divisor
fac1.mantissa[i]=quotient \ divisor
Next
quotient = remder * &h100000000ull
quotient=quotient \ divisor
carry=fac1.mantissa[0]
j=__builtin_clz(fac1.mantissa[0])-3
If j>0 Then
shiftl(fac1, j, dwords)
fac1.exponent-=j
End If
While fac1.mantissa[0]>(&h1ffffffful)
shiftr(fac1, 1, dwords)
fac1.exponent+=1
Wend
NORM_FAC1(fac1)
If den<0 Then
fac1.sign=fac1.sign Xor &h8000
End If
Return fac1
End Function
Function fpdiv(Byref x As bigfloat_struct, Byref y As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct fac1, fac2, one
Dim As Integer i, er, is_power_of_two
Dim As Short sign
fac1=x
fac2=y
one.exponent=(BIAS+1)
one.mantissa[0]=&h10000000ul
one=ui2fp(1)
sign=fac2.sign
fac2.sign=0
If fpcmp(fac2, one)=0 Then
Return fac1
Else
fac2.sign=sign
End If
If fac2.exponent=0 Then ' if fac2 = 0, return
' a divide-by-zero error and
' bail out.
fac1.mantissa[i]=&H1FFFFFFFul
For i=1 To dwords
fac1.mantissa[i]=&HFFFFFFFFul
Next
fac1.exponent=&HFFFFF+BIAS+1
er=DIVZ_ERR
Return fac1
Elseif fac1.exponent=0 Then 'fact1=0, just return
er=0
Return fac1
Else
'check to see if fac2 is a power of ten
is_power_of_two=0
If fac2.mantissa[0]=&H10000000ul Then
is_power_of_two=1
For i=1 To dwords
If fac2.mantissa[i]<>0 Then
is_power_of_two=0
Exit For
End If
Next
End If
'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_two=1 Then
fac1.sign=fac1.sign Xor fac2.sign
fac1.exponent=fac1.exponent-fac2.exponent+BIAS+1
Return fac1
End If
#Define min(a,b) Iif(((a)<(b)),(a),(b))
#Macro realw(w, j)
((w(j - 1)*b + w(j))*b + w(j + 1))*b +Iif(Ubound(w)>=j+2,w(j+2),0)
#Endmacro
#Macro subtract(w, q, d, ka, kb)
For j=ka To kb
w(j) = w(j) - q*d(j - ka + 2)
Next
#Endmacro
#Macro normalize(w, ka, q)
w(ka) = w(ka) + w(ka - 1)*b
w(ka - 1) = q
#Endmacro
#Macro finalnorm(w, kb)
For j=kb To 3 Step -1
carry=Iif(w(j)<0, ((-w(j) - 1)\b) + 1, Iif(w(j) >= b, -(w(j)\b), 0))
w(j) = w(j) + carry*b
w(j - 1) = w(j - 1) - carry
Next
#Endmacro
Dim As Double result(1 To 2*dwords+3), n(1 To 2*dwords+3), d(1 To 2*dwords+3)
Const b=&H10000
Dim As Long j, last, laststep, q, t
Dim As Long stp, carry
Dim As Double xd, xn, rund
Dim As Double w(1 To Ubound(n)+4)
For j=0 To dwords
n(2*j+2)=fac1.mantissa[j]\&H10000
n(2*j+3)=fac1.mantissa[j] Mod &H10000
d(2*j+2)=fac2.mantissa[j]\&H10000
d(2*j+3)=fac2.mantissa[j] Mod &H10000
Next
n(1)=(fac1.exponent And &h7FFFFFFF)-BIAS-1
d(1)=(fac2.exponent And &h7FFFFFFF)-BIAS-1
For j=Ubound(n) To Ubound(w)
w(j)=0
Next
t=Ubound(n)-1
w(1)=n(1)-d(1)+1
w(2)=0
For j=2 To Ubound(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, Ubound(W))
subtract(w, q, d, (stp + 2), last)
normalize(w, (stp + 2), q)
Next
FinalNorm(w, (laststep + 1))
laststep = Iif(w(2) = 0, laststep, laststep - 1)
rund = w(laststep + 1)/b
w(laststep) = w(laststep) + Iif(rund >= 0.5, 1, 0)
If w(2)=0 Then
For j=1 To t+1
result(j)=w(j+1)
Next
Else
For j=1 To t+1
result(j)=w(j)
Next
End If
result(1) = Iif(w(2) = 0, w(1) - 1, w(1))
For j=0 To dwords
fac1.mantissa[j]=result(2*j+2)*&H10000+result(2*j+3)
Next
NORM_FAC1(fac1, dwords)
fac1.exponent=(result(1)+BIAS)
End If
fac1.sign=fac1.sign Xor fac2.sign
Return fac1
End Function
Function fpipow(Byref x As bigfloat_struct, Byval e As Longint, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
'take x to an Long power
Dim As bigfloat_struct y=x
Dim As bigfloat_struct z, one
Dim As Longint n, c=0
Dim As Integer i
n = Abs(e)
one.exponent=(BIAS+1)
one.mantissa[0]=&h10000000ul
z=one
While n>0
While (n And 1)=0
n\=2
y=fpmul(y, y, dwords)
c+=1
Wend
n-=1
z=fpmul(y, z, dwords)
c+=1
Wend
If e<0 Then
z=fpdiv(one, z, dwords)
End If
Return z
End Function
Function dbl2fp(Byval x As Double) As bigfloat_struct
Dim As bigfloat_struct bf
Dim As Ulongint n, e
n=Peek ( Ulongint, @x )
e=n
n=Cast(Ulongint, (n Shl 12))
n=Cast(Ulongint, (n Shr 4))
n+=&h1000000000000000ull
bf.mantissa[0]=Cast(Ulongint, (n Shr 32))
n=e
n=Cast(Ulongint, (n Shl 40))
n=Cast(Ulongint, (n Shr 32))
bf.mantissa[1]=n
n=n+&h1000000000000000ull
e=Cast(Ulongint, (e Shr 52))
e=e-&h3FFul
bf.exponent=e+BIAS+1
Return bf
End Function
Function fp2dbl(Byref x As bigfloat_struct) As Double
Dim As Double dbl
Dim As Ulongint n, e
Dim As Ushort ex
Dim As Long sign
If x.sign<>0 Then
sign=-1
Else
sign=1
End If
If x.exponent>0 Then
ex=(x.exponent And &h7FFFFFFF)-BIAS-1
Else
ex=0
End If
n=x.mantissa[0]
n-=&h10000000ul
n=Cast(Ulongint, (n Shl 24))
e=&h3FF+ex
e=e*&h10000000000000ull
n=n+e
e=x.mantissa[1]\&h100
n=n+e
Poke Ulongint,@dbl,n
Return dbl*sign
End Function
Function fp2dbl2(Byref x As bigfloat_struct) As Double
Dim As Double dbl
Dim As Ulongint n, e
n=x.mantissa[0]
n-=&h10000000ul
n=Cast(Ulongint, (n Shl 24))
e=&h3FF+1
e=e*&h10000000000000ull
n=n+e
e=x.mantissa[1]\&h100
n=n+e
Poke Ulongint,@dbl,n
Return dbl
End Function
Function factorial(Byval n As Ulong) As bigfloat_struct
Dim As bigfloat_struct one, f, indx
one.mantissa[0]=&h10000000ul
one.exponent=0+BIAS+1
f=one
indx=one
For i As Long=1 To n
f=fpmul(f, indx)
indx=fpadd(one, indx)
Next
Return f
End Function
Function factorial_inv(Byval n As Ulong) As bigfloat_struct
Dim As bigfloat_struct f
f.mantissa[0]=&h10000000ul
f.exponent=0+BIAS+1
For i As Long=1 To n
f=fpdiv_si(f, i)
Next
Return f
End Function
Function fp2ui(Byref x As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As Ulong
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct fac1
Dim As Double n
fac1.exponent=0
fac1.mantissa[0]=x.mantissa[0]
n=fp2dbl(x)
Return Fix(n)
End Function
Function fp2si(Byref x As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct fac1
Dim As Long sign
Dim As Double n
fac1.exponent=0
fac1.mantissa[0]=x.mantissa[0]
n=fp2dbl(x)
Return Fix(n)
End Function
Function str2fp_aux(Byval value As String, Byref s As Long, Byref ex As Long ) As String
Dim As Integer j,d,e,ep,es,i,f,fp,fln
Dim As String c,f1,f2,f3, ts
Dim As Ulong ulng
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 Then
If c=" " Then
j=j+1
Continue While
Endif
If c="-" Then
es=-es
c=""
Endif
If c="+" Then
j=j+1
Continue While
Endif
If (c="0") And (f3="") Then
j=j+1
Continue While
Endif
If (c>"/") And (c<":") Then 'c is digit between 0 and 9
f3=f3+c
ex=10*ex+(Asc(c)-48)
j=j+1
Continue While
Endif
Endif
If c=" " Then
j=j+1
Continue While
Endif
If c="-" Then
s=-s
j=j+1
Continue While
Endif
If c="+" Then
j=j+1
Continue While
Endif
If c="." Then
If d=1 Then
j=j+1
Continue While
Endif
d=1
Endif
If (c>"/") And (c<":") Then 'c is digit between 0 and 9
If ((c="0") And (i=0)) Then
If d=0 Then
j=j+1
Continue While
Endif
If (d=1) And (f=0) Then
e=e-1
j=j+1
Continue While
Endif
Endif
If d=0 Then
f1=f1+c
i=i+1
Else
If (c>"0") Then
fp=1
Endif
f2=f2+c
f=f+1
Endif
Endif
If c="E" Or c="D" Then
ep=1
Endif
j=j+1
Wend
If fp=0 Then
f=0
f2=""
Endif
If s=-1 Then s=&h8000 Else s=0
f1=f1+f2
f1=Mid(f1,1,1)+Right(f1,Len(f1)-1)
fln=Len(f1)
If fp=0 Andalso ex=0 Then
ex=es*fln-1
Else
ex=es*ex-1+i+e
End If
If fln>((NUM_BITS+32)*0.3010299956639811) Then
f1=Mid(f1,1,((NUM_BITS+32)*0.3010299956639811)-1)
Endif
Return f1
End Function
Function str2fp(Byref x As String) As bigfloat_struct
Dim As Long strlen, i, n, ex, sign
Dim As String s2, strin=x
Dim As bigfloat_struct y, z, pw, ten
strin=str2fp_aux(strin, sign, ex)
strlen=Len(strin)
n=1
ten=ui2fp(Valuint("1"+String(n,"0")))
If strlen>n Then
s2=Left(strin, n)
strin=Mid(strin, n+1)
Else
z=ui2fp(Valuint(strin))
If ex<>0 Then
pw=ui2fp(10)
pw=fpipow(pw, strlen-ex-1)
z=fpdiv(z, pw)
End If
Return z
End If
z=ui2fp(Valuint(s2))
While Len(strin)>=n
s2=Left(strin, n)
strin=Mid(strin, n+1)
y=ui2fp(Valuint(s2))
z=fpmul(z, ten)
z=fpadd(z,y)
Wend
pw=ui2fp(10)
pw=fpipow(pw,strlen-ex-1)
z=fpdiv(z, pw)
Return z
End Function
Function fp2str(Byref zz As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As String
Dim As Long ex, i, n, powten, sign, tn
Dim As Long m
Dim As bigfloat_struct y, ten, z, ep
Dim As String s, s2
If zz.exponent=0 Then
Return "0"
End If
z=zz
z.sign=0
sign=zz.sign
If z.exponent<>0 Then
ex=(z.exponent And &h7FFFFFFFul)-BIAS-1
Else
ex=0
End If
n=NUM_BITS-ex
ep.exponent=-(n+15)+BIAS+1
ep.mantissa[0]=&h10000000ul
z=fpadd(z, ep)
If ex>0 Then
While ex>=3
z=fpdiv_si(z, 10)
powten+=1
ex=(z.exponent And &h7FFFFFFFul)-BIAS-1
Wend
Else
While ex<0
z=fpmul_si(z, 10)
powten-=1
ex=(z.exponent And &h7FFFFFFFul)-BIAS-1
Wend
End If
n=1
m=fp2ui(z)
If m=0 Then
z=fpmul_si(z, 10)
powten-=1
ex=(z.exponent And &h7FFFFFFFul)-BIAS-1
m=fp2ui(z)
End If
s=Trim(Str(m))+"."
If sign<>0 Then
s="-"+s
Else
s=" "+s
End If
ten=ui2fp(Clng("1"+String(n,"0")))
tn=Clng("1"+String(n,"0"))
For i=1 To (dwords*9.63)\n+1
y=ui2fp(m)
z=fpsub(z,y)
z=fpmul(z,ten)
m=fp2ui(z)
If m<0 Or m>tn Then
m=0
End If
s2=Trim(Str(m))
s2=String(n-len(s2),"0")+s2
s=s+s2
Next
s2=Trim(Str(Abs(powten)))
s2=String(5-len(s2),"0")+s2
If powten<0 Then
s2="-"+s2
Else
s2="+"+s2
End If
s=s+"e"+s2
Return s
End Function
'integer part of num
Function fpfix( Byref num As bigfloat_struct ) As bigfloat_struct
Dim As bigfloat_struct ip
Dim As Long ex, ex2, j, k, c
ex=(num.exponent And &h7FFFFFFFul)-BIAS
If ex<1 Then
Return ip
End If
If ex>=(NUM_BITS) Then
Return num
End If
ex2=ex\32
k=ex2
j=ex Mod 32
While ex2>0
ex2-=1
ip.mantissa[ex2]=num.mantissa[ex2]
Wend
ip.mantissa[k]=num.mantissa[k]
ip.mantissa[k]=shr32(ip.mantissa[k], (32-j)-3, c)
ip.mantissa[k]=shl32(ip.mantissa[k], (32-j)-3, c)
ip.exponent=ex+BIAS
ip.sign=num.sign
Return ip
End Function
'fractional part of num
Function fpfrac( Byref num As bigfloat_struct ) As bigfloat_struct
Dim As bigfloat_struct ip, fp
ip=fpfix(num)
fp=fpsub(num, ip)
Return fp
End Function
Dim As bigfloat_struct x, x2, x3, y, z, term, fac, pi, eps
Dim As Long i, n, sign
Dim As String s, spi
Dim As Double t
eps.exponent=-(NUM_BITS)+BIAS+1
eps.mantissa[0]=&h10000000ul
spi="3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821"
pi=str2fp(spi)
Print "Pi+Pi"
y=fpadd(pi, pi)
Print fp2str(y)
Print "Pi/2"
y=fpdiv_si(pi, 2)
Print fp2str(y)
Print "integer part of Pi/2"
z=fpfix(y)
Print fp2str(z)
Print "fractional part of Pi/2"
z=fpfrac(y)
Print fp2str(z)
Print "Pi*Pi"
z=fpmul(pi, pi)
Print fp2str(z)
Print "factorial(10000)"
z=factorial(10000)
Print fp2str(z)
Print "Pi^100"
z=fpipow(pi, 100)
Print fp2str(z)
t=Timer
'' sin(.5)
x=str2fp(".5")
y=x
x3=x
x2=fpmul(x, x)
term=ui2fp(1)
fac=term
sign=-1
i=1
While fpcmp(term, eps)>0 ' fpcmp returns -1 if term<eps, 0 if equal and 1 if greater
x3=fpmul(x3, x2)
fac=fpdiv_si(fac, i+1)
fac=fpdiv_si(fac, i+2)
term=fpmul(x3, fac)
If sign<0 Then
y=fpsub(y, term)
Else
y=fpadd(y, term)
End If
i+=2
sign=-sign
Wend
t=timer-t
Print "sin(.5)"
Print fp2str(y)
Print "time to comput sine function ";t;" seconds"
Print
Print "you can specify the number of dwords to the fp2str function"
Print fp2str(pi,2)
Sleep
Re: about bignum
with overloading
Code: Select all
Extern "c"
Declare Function __builtin_clz (Byval x As Ulong) As Long
Declare Function __builtin_sadd_overflow (Byval As Long, Byval As Long, Byref As Long) As Boolean
Declare Function __builtin_saddll_overflow (Byval As Longint, Byval As Longint, Byref As Longint) As Boolean
Declare Function __builtin_uadd_overflow (Byval As Ulong, Byval As Ulong, Byref As Ulong) As Boolean
Declare Function __builtin_uaddll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As Boolean
Declare Function __builtin_ssub_overflow (Byval As Long, Byval As Long, Byref As Long) As Boolean
Declare Function __builtin_ssubll_overflow (Byval As Longint, Byval As Longint, Byref As Longint) As Boolean
Declare Function __builtin_usub_overflow (Byval As Ulong, Byval As Ulong, Byref As Ulong) As Boolean
Declare Function __builtin_usubll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As Boolean
Declare Function __builtin_smul_overflow (Byval As Long, Byval As Long, Byref As Long) As Boolean
Declare Function __builtin_smulll_overflow (Byval As Longint, Byval As Longint, Byref As Longint) As Boolean
Declare Function __builtin_umul_overflow (Byval As Ulong, Byval As Ulong, Byref As Ulong) As Boolean
Declare Function __builtin_umulll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As Boolean
End Extern
Const NUMBER_OF_DIGITS = 96
Const NUMBER_OF_BITS = NUMBER_OF_DIGITS*3.321928094887362
#If (NUMBER_OF_BITS Mod 32)<>0
Const NUM_BITS = (NUMBER_OF_BITS\32+1)*32
#Else
Const NUM_BITS = NUMBER_OF_BITS
#Endif
Const NUM_DWORDS = NUM_BITS\32
Const BIAS = 1073741824 '2 ^ 30
' Error definitions
#Define DIVZ_ERR 1 'Divide by zero
#Define EXPO_ERR 2 'Exponent overflow error
#Define EXPU_ERR 3 'Exponent underflow error
Type BigFloat_struct
Declare Constructor ( )
Declare Constructor ( Byref rhs As BigFloat_struct )
Declare Destructor ( )
Declare Operator Let ( Byref rhs As BigFloat_struct )
As Short sign
As Ulong exponent
As Ulong Ptr mantissa
End Type
Constructor BigFloat_struct ( )
If this.mantissa<>0 Then
Print "pointer is not 0"
End(1)
End If
this.mantissa=Callocate((NUM_DWORDS +1), 4)
If this.mantissa=0 Then
Print "unable to allocate memory"
End(4)
End If
this.sign=0
this.exponent=0
End Constructor
Constructor BigFloat_struct ( Byref rhs As BigFloat_struct )
If this.mantissa<>0 Then
Print "pointer is not 0"
End(1)
End If
this.mantissa=Callocate((NUM_DWORDS +1), 4)
If this.mantissa=0 Then
Print "unable to allocate memory"
End(4)
End If
this.sign=rhs.sign
this.exponent=rhs.exponent
For i As Long=0 To NUM_DWORDS
this.mantissa[i]=rhs.mantissa[i]
Next
End Constructor
Operator BigFloat_struct.let ( Byref rhs As BigFloat_struct )
this.sign=rhs.sign
this.exponent=rhs.exponent
For i As Long=0 To NUM_DWORDS
this.mantissa[i]=rhs.mantissa[i]
Next
End Operator
Destructor BigFloat_struct ( )
If this.mantissa=0 Then
Print "unable to de-allocate memory"
End(-4)
Else
Deallocate(this.mantissa)
End If
this.mantissa=0
End Destructor
Type BigFloat
Declare Constructor ( )
Declare Constructor ( Byval rhs As Long )
Declare Constructor ( Byval rhs As Double )
Declare Constructor ( Byref rhs As String )
Declare Constructor ( Byref rhs As BigFloat )
Declare Destructor ( )
Declare Operator Let ( Byval rhs As Long )
Declare Operator Let ( Byval rhs As Double )
Declare Operator Let ( Byref rhs As String )
Declare Operator Let ( Byref rhs As BigFloat )
Declare Operator Cast ( ) As String
Declare Operator Cast ( ) As Long
Declare Operator Cast ( ) As Double
'----------------------------------------------
Declare Operator += (Byref rhs As BigFloat)
Declare Operator += (Byval rhs As Long)
Declare Operator += (Byval rhs As Double)
Declare Operator += (Byref rhs As String)
Declare Operator -= (Byref rhs As BigFloat)
Declare Operator -= (Byval rhs As Long)
Declare Operator -= (Byval rhs As Double)
Declare Operator -= (Byref rhs As String)
Declare Operator *= (Byref rhs As BigFloat)
Declare Operator *= (Byval rhs As Long)
Declare Operator *= (Byval rhs As Double)
Declare Operator *= (Byref rhs As String)
Declare Operator /= (Byref rhs As BigFloat)
Declare Operator /= (Byval rhs As Long)
Declare Operator /= (Byval rhs As Double)
Declare Operator /= (Byval rhs As Single)
Declare Operator /= (Byref rhs As String)
' For Next Implicit step = +1
Declare Operator For ( )
Declare Operator Step( )
Declare Operator Next( Byref end_cond As BigFloat ) As Integer
' For Next Exlicit step
Declare Operator For ( Byref stp As BigFloat )
Declare Operator Step( Byref stp As BigFloat )
Declare Operator Next( Byref end_cond As BigFloat, Byref step_var As BigFloat ) As Integer
Declare Function toString( Byval places As Long=NUM_DWORDS*9.63 ) As String
Declare Function toLong ( ) As Long
Declare Function toDouble ( ) As Double
BigNum As BigFloat_struct
End Type
Declare Function fpcmp(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
Declare Function fpadd(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpsub(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpmul(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpmul_si(Byref x As BigFloat_struct, Byval y As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpdiv_si(Byref num As BigFloat_struct, Byval den As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpdiv(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpipow(Byref x As BigFloat_struct, Byval e As Longint, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function dbl2fp(Byval x As Double) As BigFloat_struct
Declare Function fp2dbl(Byref x As BigFloat_struct) As Double
Declare Function ui2fp(Byval m As Ulong, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fp2ui(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Ulong
Declare Function si2fp(Byval m As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fp2si(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
Declare Function str2fp(Byref x As String) As BigFloat_struct
Declare Function fp2str(Byref zz As BigFloat_struct, Byval digits As Long=NUM_DWORDS*9.63) As String
Declare Function fpfix( Byref num As BigFloat_struct ) As BigFloat_struct
Declare Function fpfrac( Byref num As BigFloat_struct ) As BigFloat_struct
Declare Function fpsqr (Byref num As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
Function Bigfloat.toString( Byval places As Long=NUM_DWORDS*9.63 ) As String
Function = fp2str( this.BigNum, places )
End Function
Function Bigfloat.toLong ( ) As Long
Dim As Double x
x=fp2dbl(this.BigNum)
Return Clng(x)
End Function
Function Bigfloat.toDouble ( ) As Double
Function = fp2dbl(this.BigNum)
End Function
Constructor Bigfloat ( )
this.BigNum=si2fp(0, NUM_DWORDS)
End Constructor
Constructor Bigfloat ( Byval rhs As Long )
this.BigNum=si2fp( rhs, NUM_DWORDS )
End Constructor
Constructor Bigfloat ( Byref rhs As String )
this.BigNum=str2fp( rhs )
End Constructor
Constructor Bigfloat ( Byref rhs As Bigfloat)
this.BigNum.sign=rhs.BigNum.sign
this.BigNum.exponent=rhs.BigNum.exponent
For i As Long=0 To NUM_DWORDS
this.BigNum.mantissa[i]=rhs.BigNum.mantissa[i]
Next
End Constructor
Constructor Bigfloat ( Byval rhs As Double )
this.BigNum=str2fp( Str( rhs ) )
End Constructor
Operator Bigfloat.let ( Byref rhs As Bigfloat )
this.BigNum.sign=rhs.BigNum.sign
this.BigNum.exponent=rhs.BigNum.exponent
For i As Long=0 To NUM_DWORDS
this.BigNum.mantissa[i]=rhs.BigNum.mantissa[i]
Next
End Operator
Operator Bigfloat.let ( Byval rhs As Long )
this.BigNum=si2fp( rhs, NUM_DWORDS )
End Operator
Operator Bigfloat.let ( Byref rhs As String )
this.BigNum=str2fp( rhs )
End Operator
Operator Bigfloat.Let ( Byval rhs As Double )
this.BigNum=str2fp( Str(rhs) )
End Operator
Operator Bigfloat.cast ( ) As String
Operator = fp2str(this.BigNum)
End Operator
Operator Bigfloat.cast ( ) As Long
Dim As Double x
x=fp2dbl(this.BigNum)
Operator = Clng(x)
End Operator
Operator Bigfloat.cast ( ) As Double
Operator = fp2dbl(this.BigNum)
End Operator
Destructor Bigfloat ( )
End Destructor
'============================================================================
Operator BigFloat.for ( )
End Operator
Operator BigFloat.step ( )
this.BigNum = fpadd(this.BigNum,si2fp(1, NUM_DWORDS), NUM_DWORDS)
End Operator
Operator BigFloat.next ( Byref end_cond As BigFloat ) As Integer
Return fpcmp(This.BigNum, end_cond.BigNum, NUM_DWORDS)<=0
End Operator
'' explicit step versions
''
Operator BigFloat.for ( Byref step_var As BigFloat )
End Operator
Operator BigFloat.step ( Byref step_var As BigFloat )
this.BigNum = fpadd(this.BigNum, step_var.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.next ( Byref end_cond As BigFloat, Byref step_var As BigFloat ) As Integer
If fpcmp(step_var.BigNum, si2fp(0, NUM_DWORDS), NUM_DWORDS) < 0 Then
Return fpcmp(This.BigNum, end_cond.BigNum, NUM_DWORDS) >= 0
Else
Return fpcmp(This.BigNum, end_cond.BigNum, NUM_DWORDS) <= 0
End If
End Operator
'============================================================================
Operator + ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum= fpadd(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator + ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp( Str(rhs) )
result.BigNum= fpadd(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator + ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp( Str(lhs) )
result.BigNum= fpadd( result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator + (Byref lhs As Bigfloat, Byref rhs As String) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp(rhs)
result.BigNum= fpadd(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - ( Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result = rhs
result.BigNum.sign = result.BigNum.sign Xor &h8000
Operator = result
End Operator
Operator - ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum= fpsub(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp( Str(rhs) )
result.BigNum= fpsub(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp( Str(lhs) )
result.BigNum= fpsub( result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - ( Byref lhs As Bigfloat, Byval rhs As Long ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
result.BigNum= fpsub(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - ( Byval lhs As Long, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(lhs, NUM_DWORDS)
result.BigNum= fpsub(result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - (Byref lhs As Bigfloat, Byref rhs As String) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp(rhs)
result.BigNum= fpsub(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum= fpmul(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp( Str(rhs) )
result.BigNum= fpmul(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp( Str(lhs) )
result.BigNum= fpmul( result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * ( Byref lhs As Bigfloat, Byval rhs As Long ) As Bigfloat
Dim As Bigfloat result
If Abs(rhs)<&h100000000ull Then
result.BigNum=fpmul_si(lhs.BigNum, rhs)
Else
result.BigNum=si2fp(rhs, NUM_DWORDS)
result.BigNum= fpmul(lhs.BigNum, result.BigNum, NUM_DWORDS)
End If
Operator = result
End Operator
Operator * ( Byval lhs As Long, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(lhs, NUM_DWORDS)
result.BigNum= fpmul(result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * (Byref lhs As Bigfloat, Byref rhs As String) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp(rhs)
result.BigNum= fpmul(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum= fpdiv(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp( Str(rhs) )
result.BigNum= fpdiv(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp( Str(lhs) )
result.BigNum= fpdiv( result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / ( Byref lhs As Bigfloat, Byval rhs As Long ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
result.BigNum= fpdiv(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / ( Byval lhs As Long, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(lhs, NUM_DWORDS)
result.BigNum= fpdiv(result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / (Byref lhs As Bigfloat, Byref rhs As String) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp(rhs)
result.BigNum= fpdiv(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator ^ ( Byref lhs As Bigfloat, Byval rhs As Longint ) As Bigfloat
Dim As Bigfloat lhs2
lhs2.BigNum=fpipow(lhs.BigNum, rhs, NUM_DWORDS)
Operator = lhs2
End Operator
Operator BigFloat.+= (Byref rhs As BigFloat)
this.BigNum = fpadd(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.+= (Byval rhs As Long)
Dim As BigFloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.+= (Byval rhs As Double)
Dim As BigFloat result
result.BigNum=str2fp( Str(rhs) )
this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.+= (Byref rhs As String)
Dim As BigFloat result
result.BigNum=str2fp(rhs)
this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.-= (Byref rhs As BigFloat)
this.BigNum = fpsub(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.-= (Byval rhs As Long)
Dim As BigFloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.-= (Byval rhs As Double)
Dim As BigFloat result
result.BigNum=str2fp( Str(rhs) )
this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.-= (Byref rhs As String)
Dim As BigFloat result
result.BigNum=str2fp(rhs)
this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.*= (Byref rhs As BigFloat)
this.BigNum = fpmul(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.*= (Byval rhs As Long)
Dim As BigFloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.*= (Byval rhs As Double)
Dim As BigFloat result
result.BigNum=str2fp( Str(rhs) )
this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.*= (Byref rhs As String)
Dim As BigFloat result
result.BigNum=str2fp(rhs)
this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat./= (Byref rhs As BigFloat)
this.BigNum = fpdiv(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat./= (Byval rhs As Long)
Dim As BigFloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat./= (Byval rhs As Double)
Dim As BigFloat result
result.BigNum=str2fp( Str(rhs) )
this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat./= (Byref rhs As String)
Dim As BigFloat result
result.BigNum=str2fp(rhs)
this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator = ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)=0
End Operator
Operator < ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<0
End Operator
Operator > ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)>0
End Operator
Operator <= ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<=0
End Operator
Operator >= ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)>=0
End Operator
Operator <> ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<>0
End Operator
Operator Abs(Byref rhs As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=rhs.BigNum
result.BigNum.sign = 0
Operator = result
End Operator
Operator Fix (Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum =fpfix(x.BigNum)
Operator = result
End Operator
Operator Frac(Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=fpfrac(x.BigNum)
Operator = result
End Operator
Operator Sqr(Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=fpsqr(x.BigNum)
Operator = result
End Operator
Dim Shared As BigFloat_struct Log102
Log102.sign=0
Log102.exponent=(-2)+BIAS+1
If NUM_DWORDS>=0 Then Log102.mantissa[0]=&h13441350
If NUM_DWORDS>=1 Then Log102.mantissa[1]=&h9F79FEF3
If NUM_DWORDS>=2 Then Log102.mantissa[2]=&h11F12B35
If NUM_DWORDS>=3 Then Log102.mantissa[3]=&h816F922F
Private Function log2_32(Byval value As Ulong) As Long
'https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers
Static tab32(0 To 31) As Const Long = {0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30, 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31}
value Or= Culng(value Shr 1)
value Or= Culng(value Shr 2)
value Or= Culng(value Shr 4)
value Or= Culng(value Shr 8)
value Or= Culng(value Shr 16)
Return tab32(Culng(Culng(value * &h07C4ACDD) Shr 27))
End Function
Private Function log2_64(Byval value As Ulongint) As Long
'https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers
Static tab64(0 To 63) As Const Long ={63, 0, 58, 1, 59, 47, 53, 2, 60, 39, 48, 27, 54, 33, 42, 3, 61, 51, 37, 40, 49, 18, 28, 20, 55, 30, 34, 11, 43, 14, 22, 4, 62, 57, 46, 52, 38, 26, 32, 41, 50, 36, 17, 19, 29, 10, 13, 21, 56, 45, 25, 31, 35, 16, 9, 12, 44, 24, 15, 8, 23, 7, 6, 5}
value Or= Culngint(value Shr 1)
value Or= Culngint(value Shr 2)
value Or= Culngint(value Shr 4)
value Or= Culngint(value Shr 8)
value Or= Culngint(value Shr 16)
value Or= Culngint(value Shr 32)
Return tab64((Culngint((value - (value Shr 1)) * &h07EDD5E59A4E28C2) Shr 58))
End Function
Function ipower (Byval x As Ulongint, Byval e As Ulongint) As Ulongint
'take x to an integer power
Dim As Ulongint z, y, n
y = x
n = e
z = 1
While n > 0
While (n And 1) = 0
n = n \ 2
y = y * y
Wend
n = n - 1
z = y * z
Wend
Return z
End Function
Function shl32(Byval n As Ulong, Byval k As Ubyte, Byref c As Ulong) As Ulong
If k>0 And k<32 Then
Dim As Ulong carry=n
Dim As Ubyte k32=32-k
Asm
mov cl, [k32]
Shr dword Ptr [carry], cl
mov cl, [k]
Shl dword Ptr [n], cl
End Asm
c=carry
End If
Return n
End Function
Function shr32(Byval n As Ulong, Byval k As Ubyte, Byref c As Ulong) As Ulong
If k>0 And k<32 Then
Dim As Ulong carry=n
Dim As Ubyte k32=32-k
Asm
mov cl, [k32]
Shl dword Ptr [carry], cl
mov cl, [k]
Shr dword Ptr [n], cl
End Asm
c=carry
End If
Return n
End Function
Sub shiftl(Byref fac1 As BigFloat_struct, Byval k As Long, Byval dwords As Long=NUM_DWORDS)
If k>0 And k<32 Then
Dim As Long i
Dim As Ulong n, carry, c=0
Dim As Long k32=32-k
For i=dwords To 0 Step -1
n=fac1.mantissa[i]
carry=n
carry=Cast(Ulong, (carry Shr k32))
n=Cast(Ulong, (n Shl k))
fac1.mantissa[i]=n+c
c=carry
Next
Elseif k=32 Then
Dim As Long i
For i=0 To dwords-1
fac1.mantissa[i]=fac1.mantissa[i+1]
Next
fac1.mantissa[dwords]=0
End If
End Sub
Sub shiftr(Byref fac1 As BigFloat_struct, Byval k As Long, Byval dwords As Long=NUM_DWORDS)
If k>0 And k<32 Then
Dim As Long i
Dim As Ulong n, carry, c=0
Dim As Long k32=32-k
For i=0 To dwords
n=fac1.mantissa[i]
carry=n
carry=Cast(Ulong, (carry Shl k32))
n=Cast(Ulong, (n Shr k))
fac1.mantissa[i]=c+n
c=carry
Next
Elseif k=32 Then
Dim As Long i
For i=dwords To 1 Step -1
fac1.mantissa[i]=fac1.mantissa[i-1]
Next
fac1.mantissa[0]=0
End If
End Sub
Function fpcmp(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As Long i
Dim As Longint c
If x.sign < y.sign Then fpcmp = -1: Exit Function
If x.sign > y.sign Then fpcmp = 1: Exit Function
If x.exponent<y.exponent Then Return -1
If x.exponent>y.exponent Then Return 1
For i=0 To dwords
c=x.mantissa[i]-y.mantissa[i]
If c<>0 Then Exit For
Next
If c<0 Then Return -1
If c=0 Then Return 0
If c>0 Then Return 1
End Function
Function NORM_FAC1(Byref fac1 As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Integer
If dwords>NUM_DWORDS Then dwords=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.
Dim As Integer i,er=0,f=0
For i=0 To dwords
If fac1.mantissa[i]>0 Then f=1
Next
If f=0 Then
fac1.exponent=0
fac1.sign=0
Exit Function
End If
While fac1.mantissa[0]=0
shiftl(fac1, 32, dwords)
fac1.exponent-=32
If fac1.exponent=0 Then
NORM_FAC1=EXPU_ERR
Exit Function
End If
Wend
i=__builtin_clz(fac1.mantissa[0])-3
If i>0 Then
shiftl(fac1, i, dwords)
fac1.exponent-=i
End If
'if the highmost Bit in fac1_man is nonzero,
'shift the mantissa right 1 Bit and
'increment the exponent
If fac1.mantissa[0]>(&h1ffffffful)Then
While fac1.mantissa[0]>&h1ffffffful
shiftr(fac1, 1, dwords)
fac1.exponent+=1
Wend
Elseif fac1.mantissa[0]<&h10000000ul Then
/' the following will probably never be executed '/
'now shift fac1_man 1 to the left until a
'nonzero Bit appears in the next-to-highest
'Bit of fac1_man. decrement exponent for
'each shift.
While fac1.mantissa[0]<&h10000000ul
shiftl(fac1, 1, dwords)
fac1.exponent-=1
Wend
End If
'check for overflow/underflow
If fac1.exponent<0 Then
NORM_FAC1=EXPO_ERR
End If
End Function
Function si2fp(Byval m As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1
Dim As Long n=Abs(m)
If m=0 Then
Return fac1
End If
fac1.mantissa[0]=n
NORM_FAC1(fac1, NUM_DWORDS)
fac1.exponent=log2_32(n)+BIAS+1
If m<0 Then
fac1.sign=&H8000
Else
fac1.sign=0
End If
Return fac1
End Function
Function ui2fp(Byval m As Ulong, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1
If m=0 Then
Return fac1
End If
fac1.mantissa[0]=m
NORM_FAC1(fac1, NUM_DWORDS)
fac1.exponent=log2_32(m)+BIAS+1
Return fac1
End Function
Sub fpadd_aux(Byref fac1 As BigFloat_struct, Byref fac2 As BigFloat_struct, Byval dwords As Long=NUM_DWORDS)
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As Long c, i
Dim As Longint v
c=0
For i=dwords To 1 Step -1
v=Clngint(fac2.mantissa[i])+Clngint(fac1.mantissa[i])+c
If v>(&h100000000ull-1) Then
v=v-&h100000000ull
c=1
Else
c=0
End If
fac1.mantissa[i]=v
Next
v=Clngint(fac1.mantissa[0])+Clngint(fac2.mantissa[0])+c
fac1.mantissa[0]=v
NORM_FAC1(fac1, dwords)
End Sub
Sub fpsub_aux(Byref fac1 As BigFloat_struct, Byref fac2 As BigFloat_struct, Byval dwords As Long=NUM_DWORDS)
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As Long c, i
Dim As Longint v
c=0
For i=dwords To 1 Step -1
v=Clngint(fac1.mantissa[i])-clngint(fac2.mantissa[i])-c
If v<0 Then
v=v+&h100000000ull
c=1
Else
c=0
End If
fac1.mantissa[i]=v
Next
v=Clngint(fac1.mantissa[0])-clngint(fac2.mantissa[0])-c
fac1.mantissa[0]=v
NORM_FAC1(fac1, dwords)
End Sub
Function fpadd(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1,fac2
Dim As Long i, t, c, xsign, ysign
xsign=x.sign:x.sign=0
ysign=y.sign:y.sign=0
c=fpcmp(x, y, dwords)
x.sign=xsign
y.sign=ysign
If c<0 Then
fac1=y
fac2=x
Else
fac1=x
fac2=y
End If
t=fac1.exponent-fac2.exponent
If t<(NUM_BITS+32) Then
'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 dwords, the result is already in FAC1.
If t>0 And t<(NUM_BITS+32) Then 'shift
i=t\32
While i>0
shiftr(fac2, 32, dwords)
t-=32
i-=1
Wend
'While t>0
If t>0 Then shiftr(fac2, t, dwords)
't-=1
'Wend
End If
'See if the signs of the two numbers
'are the same. If so, add; if not, subtract.
If fac1.sign=fac2.sign Then 'add
fpadd_aux(fac1,fac2, dwords)
Else
fpsub_aux(fac1,fac2, dwords)
End If
Endif
Return fac1
End Function
Function fpsub(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1,fac2
Dim As Long s
fac1=x
fac2=y
fac2.sign=fac2.sign Xor &h8000
fac1=fpadd(fac1,fac2, dwords)
Return fac1
End Function
Function fpmul(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1,fac2
Dim As Integer i, j, ex, er, den, num
Dim As Ulongint digit, carry, prod ', tmp
Dim As Ulong fac3(0 To 2*dwords+1)
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
For i=0 To dwords
fac1.mantissa[i]=0
Next
'NORM_FAC1(fac1)
Return fac1
Else
If ex<0 Then
er=EXPO_ERR
Return fac1 'Exit Function
End If
'clear fac3 mantissa
For i=0 To dwords
fac3(i)=0
Next
den=dwords
While fac2.mantissa[den]=0
den-=1
Wend
num=dwords
While fac1.mantissa[num]=0
num-=1
Wend
If num<den Then
Swap fac1, fac2
'fac1=y
'fac2=x
Swap den, num
End If
For j=den To 0 Step -1
carry=0
digit=fac2.mantissa[j]
For i=num To 0 Step -1
prod=fac3(i+j+1)+digit*fac1.mantissa[i]+carry
/'tmp=fac1.mantissa[i]
If __builtin_umulll_overflow (digit, tmp, prod) Then Print "mul overflow"
tmp=fac3(i+j+1)
If __builtin_uaddll_overflow (tmp, prod, prod) Then Print "add overflow"
If __builtin_uaddll_overflow (carry, prod, prod) Then Print "add overflow"'/
carry=prod Shr 32'\&H100000000
fac3(i+j+1)=prod '(prod mod &H100000000)
Next
fac3(j)=carry
Next
For i=0 To dwords
fac1.mantissa[i]=fac3(i)
Next
End If
'now determine exponent of result.
'as you do...watch for overflow.
'ex=fac2.exponent-BIAS+fac1.exponent
'fac1.exponent=ex
ex=(fac2.exponent And &h7FFFFFFFul)-BIAS+1
ex=ex+(fac1.exponent And &h7FFFFFFFul)-BIAS+1
fac1.exponent=ex+BIAS+1
'determine the sign of the product
fac1.sign=fac1.sign Xor fac2.sign
NORM_FAC1(fac1, dwords)
Return fac1
End Function
Function fpmul_si(Byref x As BigFloat_struct, Byval y As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1,fac2
fac1=x
fac2=ui2fp(y, dwords)
fac1=fpmul(fac1, fac2, dwords)
Return fac1
End Function
Function fpdiv_si(Byref num As BigFloat_struct, Byval den As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1
Dim As Ulongint carry, remder=0
Dim As Longint i, divisor
Dim As Longint quotient
Dim As Long j
divisor=Abs(den)
fac1=num
If divisor=1 Then
Return fac1
End If
If divisor = 0 Then
Print "error: divisor = 0"
Return fac1 'Exit function
End If
If divisor>2147483647 Then
Print "error: divisor too large"
Return fac1 'Exit function
End If
For i = 0 To dwords
quotient = fac1.mantissa[i] + remder * &h100000000ull
remder = quotient Mod divisor
fac1.mantissa[i]=quotient \ divisor
Next
quotient = remder * &h100000000ull
quotient=quotient \ divisor
carry=fac1.mantissa[0]
j=__builtin_clz(fac1.mantissa[0])-3
If j>0 Then
shiftl(fac1, j, dwords)
fac1.exponent-=j
End If
While fac1.mantissa[0]>(&h1ffffffful)
shiftr(fac1, 1, dwords)
fac1.exponent+=1
Wend
NORM_FAC1(fac1)
If den<0 Then
fac1.sign=fac1.sign Xor &h8000
End If
Return fac1
End Function
Function fpdiv(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1, fac2, one
Dim As Integer i, er, is_power_of_two
Dim As Short sign
fac1=x
fac2=y
one.exponent=(BIAS+1)
one.mantissa[0]=&h10000000ul
one=ui2fp(1)
sign=fac2.sign
fac2.sign=0
If fpcmp(fac2, one)=0 Then
Return fac1
Else
fac2.sign=sign
End If
If fac2.exponent=0 Then ' if fac2 = 0, return
' a divide-by-zero error and
' bail out.
fac1.mantissa[i]=&H1FFFFFFFul
For i=1 To dwords
fac1.mantissa[i]=&HFFFFFFFFul
Next
fac1.exponent=&HFFFFF+BIAS+1
er=DIVZ_ERR
Return fac1
Elseif fac1.exponent=0 Then 'fact1=0, just return
er=0
Return fac1
Else
'check to see if fac2 is a power of ten
is_power_of_two=0
If fac2.mantissa[0]=&H10000000ul Then
is_power_of_two=1
For i=1 To dwords
If fac2.mantissa[i]<>0 Then
is_power_of_two=0
Exit For
End If
Next
End If
'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_two=1 Then
fac1.sign=fac1.sign Xor fac2.sign
fac1.exponent=fac1.exponent-fac2.exponent+BIAS+1
Return fac1
End If
#Define min(a,b) Iif(((a)<(b)),(a),(b))
#Macro realw(w, j)
((w(j - 1)*b + w(j))*b + w(j + 1))*b +Iif(Ubound(w)>=j+2,w(j+2),0)
#Endmacro
#Macro subtract(w, q, d, ka, kb)
For j=ka To kb
w(j) = w(j) - q*d(j - ka + 2)
Next
#Endmacro
#Macro normalize(w, ka, q)
w(ka) = w(ka) + w(ka - 1)*b
w(ka - 1) = q
#Endmacro
#Macro finalnorm(w, kb)
For j=kb To 3 Step -1
carry=Iif(w(j)<0, ((-w(j) - 1)\b) + 1, Iif(w(j) >= b, -(w(j)\b), 0))
w(j) = w(j) + carry*b
w(j - 1) = w(j - 1) - carry
Next
#Endmacro
Dim As Double result(1 To 2*dwords+3), n(1 To 2*dwords+3), d(1 To 2*dwords+3)
Const b=&H10000
Dim As Long j, last, laststep, q, t
Dim As Long stp, carry
Dim As Double xd, xn, rund
Dim As Double w(1 To Ubound(n)+4)
For j=0 To dwords
n(2*j+2)=fac1.mantissa[j]\&H10000
n(2*j+3)=fac1.mantissa[j] Mod &H10000
d(2*j+2)=fac2.mantissa[j]\&H10000
d(2*j+3)=fac2.mantissa[j] Mod &H10000
Next
n(1)=(fac1.exponent And &h7FFFFFFF)-BIAS-1
d(1)=(fac2.exponent And &h7FFFFFFF)-BIAS-1
For j=Ubound(n) To Ubound(w)
w(j)=0
Next
t=Ubound(n)-1
w(1)=n(1)-d(1)+1
w(2)=0
For j=2 To Ubound(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, Ubound(W))
subtract(w, q, d, (stp + 2), last)
normalize(w, (stp + 2), q)
Next
FinalNorm(w, (laststep + 1))
laststep = Iif(w(2) = 0, laststep, laststep - 1)
rund = w(laststep + 1)/b
w(laststep) = w(laststep) + Iif(rund >= 0.5, 1, 0)
If w(2)=0 Then
For j=1 To t+1
result(j)=w(j+1)
Next
Else
For j=1 To t+1
result(j)=w(j)
Next
End If
result(1) = Iif(w(2) = 0, w(1) - 1, w(1))
For j=0 To dwords
fac1.mantissa[j]=result(2*j+2)*&H10000+result(2*j+3)
Next
NORM_FAC1(fac1, dwords)
fac1.exponent=(result(1)+BIAS)
End If
fac1.sign=fac1.sign Xor fac2.sign
Return fac1
End Function
Function fpipow(Byref x As BigFloat_struct, Byval e As Longint, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
'take x to an Long power
Dim As BigFloat_struct y=x
Dim As BigFloat_struct z, one
Dim As Longint n, c=0
Dim As Integer i
n = Abs(e)
one.exponent=(BIAS+1)
one.mantissa[0]=&h10000000ul
z=one
While n>0
While (n And 1)=0
n\=2
y=fpmul(y, y, dwords)
c+=1
Wend
n-=1
z=fpmul(y, z, dwords)
c+=1
Wend
If e<0 Then
z=fpdiv(one, z, dwords)
End If
Return z
End Function
' sqrt(num)
Function fpsqr (Byref num As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct r, r2, tmp, n, half, one, eps, zero
Dim As Long ex, ex2, ex3, k, l, prec
Dim As String ts, v
Dim As Double x
one.exponent=(BIAS+1)
one.mantissa[0]=&h10000000ul
eps.exponent=-(NUM_BITS)+BIAS+1
eps.mantissa[0]=&h10000000ul
l=Log((NUM_DWORDS*9.63+9.63)*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 bigfloat has an extra 9 guard digits
n=num
If fpcmp(n, zero, dwords)=0 Then
Return r ' r is stil 0
End If
If fpcmp(n, one, dwords)=0 Then
Return one
End If
If fpcmp(n, si2fp(0, dwords), dwords)<0 Then
Return r ' r is stil 0
End If
'=====================================================================
'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 bigfloat
'divide the exponent in half for square root
'in this case 1.414213562373095e250
If n.exponent<>0 Then
ex=(n.exponent And &h7FFFFFFFul)-BIAS-1
Else
ex=0
End If
ts=fp2str(n, 16)
k=Instr(ts, "e")
v=Mid(ts,k+1)
ts=Left(ts, k-1)
ex2=Clng(v)
x=Val(ts)
If x = 0 Then Print "Div 0": Return r 'Exit Function
If x = 1 Andalso ex=0 Then
Return one
End If
If Abs(ex2) And 1 Then
x=x*10
ex2-=1
End If
x = Sqr(x) 'approximation
v=Trim(Str(x))
k=Instr(v,".")
If Len(v)>1 And k=0 Then ex2+=1
v=v+"e"+Trim(Str(ex2\2))
r=str2fp(v)
half.mantissa[0]=&h10000000
half.exponent=BIAS
half.sign=0
'=====================================================================
'Newton-Raphson method
prec=3
For k=1 To l+2
prec=2*prec-1
tmp = fpdiv(n, r, prec)
r2 = fpadd(r, tmp, prec)
r = fpmul(r2, half, prec)
Next
Return r
End Function
Function dbl2fp(Byval x As Double) As BigFloat_struct
Dim As BigFloat_struct bf
Dim As Ulongint n, e
n=Peek ( Ulongint, @x )
e=n
n=Cast(Ulongint, (n Shl 12))
n=Cast(Ulongint, (n Shr 4))
n+=&h1000000000000000ull
bf.mantissa[0]=Cast(Ulongint, (n Shr 32))
n=e
n=Cast(Ulongint, (n Shl 40))
n=Cast(Ulongint, (n Shr 32))
bf.mantissa[1]=n
n=n+&h1000000000000000ull
e=Cast(Ulongint, (e Shr 52))
e=e-&h3FFul
bf.exponent=e+BIAS+1
Return bf
End Function
Function fp2dbl(Byref x As BigFloat_struct) As Double
Dim As Double dbl
Dim As Ulongint n, e
Dim As Ushort ex
Dim As Long sign
If x.sign<>0 Then
sign=-1
Else
sign=1
End If
If x.exponent>0 Then
ex=(x.exponent And &h7FFFFFFF)-BIAS-1
Else
ex=0
End If
n=x.mantissa[0]
n-=&h10000000ul
n=Cast(Ulongint, (n Shl 24))
e=&h3FF+ex
e=e*&h10000000000000ull
n=n+e
e=x.mantissa[1]\&h100
n=n+e
Poke Ulongint,@dbl,n
Return dbl*sign
End Function
Function factorial(Byval n As Ulong) As BigFloat_struct
Dim As BigFloat_struct one, f, indx
one.mantissa[0]=&h10000000ul
one.exponent=0+BIAS+1
f=one
indx=one
For i As Long=1 To n
f=fpmul(f, indx)
indx=fpadd(one, indx)
Next
Return f
End Function
Function factorial_inv(Byval n As Ulong) As BigFloat_struct
Dim As BigFloat_struct f
f.mantissa[0]=&h10000000ul
f.exponent=0+BIAS+1
For i As Long=1 To n
f=fpdiv_si(f, i)
Next
Return f
End Function
Function fp2ui(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Ulong
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1
Dim As Double n
fac1.exponent=0
fac1.mantissa[0]=x.mantissa[0]
n=fp2dbl(x)
Return Fix(n)
End Function
Function fp2si(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1
Dim As Long sign
Dim As Double n
fac1.exponent=0
fac1.mantissa[0]=x.mantissa[0]
n=fp2dbl(x)
Return Fix(n)
End Function
Function str2fp_aux(Byval value As String, Byref s As Long, Byref ex As Long ) As String
Dim As Integer j,d,e,ep,es,i,f,fp,fln
Dim As String c,f1,f2,f3, ts
Dim As Ulong ulng
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 Then
If c=" " Then
j=j+1
Continue While
Endif
If c="-" Then
es=-es
c=""
Endif
If c="+" Then
j=j+1
Continue While
Endif
If (c="0") And (f3="") Then
j=j+1
Continue While
Endif
If (c>"/") And (c<":") Then 'c is digit between 0 and 9
f3=f3+c
ex=10*ex+(Asc(c)-48)
j=j+1
Continue While
Endif
Endif
If c=" " Then
j=j+1
Continue While
Endif
If c="-" Then
s=-s
j=j+1
Continue While
Endif
If c="+" Then
j=j+1
Continue While
Endif
If c="." Then
If d=1 Then
j=j+1
Continue While
Endif
d=1
Endif
If (c>"/") And (c<":") Then 'c is digit between 0 and 9
If ((c="0") And (i=0)) Then
If d=0 Then
j=j+1
Continue While
Endif
If (d=1) And (f=0) Then
e=e-1
j=j+1
Continue While
Endif
Endif
If d=0 Then
f1=f1+c
i=i+1
Else
If (c>"0") Then
fp=1
Endif
f2=f2+c
f=f+1
Endif
Endif
If c="E" Or c="D" Then
ep=1
Endif
j=j+1
Wend
If fp=0 Then
f=0
f2=""
Endif
If s=-1 Then s=&h8000 Else s=0
f1=f1+f2
f1=Mid(f1,1,1)+Right(f1,Len(f1)-1)
fln=Len(f1)
If fp=0 Andalso ex=0 Then
ex=es*fln-1
Else
ex=es*ex-1+i+e
End If
If fln>((NUM_BITS+32)*0.3010299956639811) Then
f1=Mid(f1,1,((NUM_BITS+32)*0.3010299956639811)-1)
Endif
Return f1
End Function
Function str2fp(Byref x As String) As BigFloat_struct
Dim As Long strlen, i, n, ex, sign
Dim As String s2, strin=x
Dim As BigFloat_struct y, z, pw, ten
strin=str2fp_aux(strin, sign, ex)
strlen=Len(strin)
n=1
ten=ui2fp(Valuint("1"+String(n,"0")))
If strlen>n Then
s2=Left(strin, n)
strin=Mid(strin, n+1)
Else
z=ui2fp(Valuint(strin))
If ex<>0 Then
pw=ui2fp(10)
pw=fpipow(pw, strlen-ex-1)
z=fpdiv(z, pw)
End If
Return z
End If
z=ui2fp(Valuint(s2))
While Len(strin)>=n
s2=Left(strin, n)
strin=Mid(strin, n+1)
y=ui2fp(Valuint(s2))
z=fpmul(z, ten)
z=fpadd(z,y)
Wend
pw=ui2fp(10)
pw=fpipow(pw,strlen-ex-1)
z=fpdiv(z, pw)
Return z
End Function
Function fp2str(Byref zz As BigFloat_struct, Byval digits As Long=NUM_DWORDS*9.63) As String
if digits>NUM_DWORDS*9.63 or digits<0 then digits=NUM_DWORDS*9.63
Dim As Long ex, i, m, n, powten, sign, tn
Dim As BigFloat_struct y, ten, tenp, z, ep, one, nine, zero
Dim As String s, s2
If zz.exponent=0 Then
Return "0"
End If
digits+=2
n=8
one=ui2fp(1)
nine=ui2fp(9)
ten=ui2fp(10)
tenp=ui2fp(Clng("1"+String(n,"0")))
tn=Clng("1"+String(n,"0"))
z=zz
sign=zz.sign
z.sign=0
If z.exponent<>0 Then
ex=(z.exponent And &h7FFFFFFFul)-BIAS-1
else
ex=0
end if
m=NUM_BITS-ex
ep.exponent=-(m)+BIAS+1
ep.mantissa[0]=&h10000000ul
z=fpadd(z,ep)
If fpcmp(z,nine)>0 Then
While fpcmp(z,nine)>0
z=fpdiv_si(z, 10)
powten+=1
Wend
Elseif fpcmp(z,one)<0 then
While fpcmp(z,one)<0
z=fpmul(z, ten)
powten-=1
Wend
End If
m=fp2ui(z)
z=fpfrac(z)
z=fpmul(z, tenp)
s=Trim(Str(m))+"."
For i=1 To digits
if fpcmp(z, one)<0 then
while fpcmp(z, one)<0
if (len(s))>=digits then exit while
z=fpmul(z, tenp)
s=s+string(n,"0")
wend
else
if (len(s))>=digits then exit for
y=fpfix(z)
m=fp2ui(y)
z=fpfrac(z)
z=fpmul(z, tenp)
s2=Trim(Str(m))
s2=String(n-len(s2),"0")+s2
s=s+s2
end if
Next
s=ltrim(s, "0")
if len(s)>(digits-1) then
s=left(s, digits-1)
end if
i=instr(trim(s),".")
if i>2 then
s=left(s,i-1)+mid(s, i+1)
s=left(s,1)+"."+mid(s, 2)
powten+=1
elseif i=1 then
s=mid(s,2)
s=left(s,1)+"."+mid(s,2)
powten-=1
end if
s2=Trim(Str(Abs(powten)))
if (5-len(s2))>0 then s2=String(5-len(s2),"0")+s2
If powten<0 Then
s2="-"+s2
Else
s2="+"+s2
End If
s=s+"e"+s2
If sign<>0 Then
s="-"+s
Else
s=" "+s
End If
Return s
End Function
'integer part of num
Function fpfix( Byref num As BigFloat_struct ) As BigFloat_struct
Dim As BigFloat_struct ip
Dim As Long ex, ex2, j, k, c
ex=(num.exponent And &h7FFFFFFFul)-BIAS
If ex<1 Then
Return ip
End If
If ex>=(NUM_BITS) Then
Return num
End If
ex2=ex\32
k=ex2
j=ex Mod 32
While ex2>0
ex2-=1
ip.mantissa[ex2]=num.mantissa[ex2]
Wend
ip.mantissa[k]=num.mantissa[k]
ip.mantissa[k]=shr32(ip.mantissa[k], (32-j)-3, c)
ip.mantissa[k]=shl32(ip.mantissa[k], (32-j)-3, c)
ip.exponent=ex+BIAS
ip.sign=num.sign
NORM_FAC1(ip)
Return ip
End Function
'fractional part of num
Function fpfrac( Byref num As BigFloat_struct ) As BigFloat_struct
Dim As BigFloat_struct ip, fp
ip=fpfix(num)
fp=fpsub(num, ip)
Return fp
End Function
Function pi_brent_salamin (Byval digits As Ulong = NUM_DWORDS*9.63) As bigfloat_struct
If digits > NUM_DWORDS*9.63 Then digits = NUM_DWORDS*9.63
Dim As Long limit
Dim As bigfloat_struct c0, c1, c2, c05
Dim As bigfloat_struct a, b, sum
Dim As bigfloat_struct ak, bk, ck
Dim As bigfloat_struct ab, asq
Dim As bigfloat_struct pow2, tmp, eps
eps.exponent=-(NUM_BITS-10)+BIAS+1
eps.mantissa[0]=&h10000000ul
limit = (-digits) + BIAS + 1
c0=si2fp(0, digits): ak = c0: bk = c0: ab = c0: asq = c0
c1=si2fp(1, digits): a = c1: ck = c1: pow2 = c1
c2=si2fp(2, digits): b = c2
c05=ui2fp(1) : c05=fpdiv_si(c05, 2)
sum = c05
b=fpsqr(b, digits)
b=fpdiv(c1, b, digits)
While fpcmp(ck, eps) > 0
ak=fpadd(a, b)
ak=fpmul(c05, ak)
ab=fpmul(a, b)
bk=fpsqr(ab)
asq=fpmul(ak, ak)
ck=fpsub(asq, ab)
pow2=fpmul(pow2, c2)
tmp=fpmul(pow2, ck)
sum=fpsub(sum, tmp)
a = ak: b = bk
Wend
tmp=fpdiv(asq, sum)
tmp=fpmul(c2, tmp)
Return tmp
End Function
Dim As BigFloat x, x2, x3, y, z, term, fac, pi, eps
Dim As Long i, n, sign
Dim As String s, spi
Dim As Double t
eps.BigNum.exponent=-(NUM_BITS-10)+BIAS+1
eps.BigNum.mantissa[0]=&h10000000ul
spi="3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865"
pi=spi
Print "Pi+Pi"
y=pi + pi
Print y
Print "Pi/2"
y=pi / 2
Print y
Print "integer part of Pi/2"
z=Fix(y)
Print z
Print "fractional part of Pi/2"
z=Frac(y)
Print z
Print "Pi*Pi"
z=pi * pi
Print z
Print "factorial(10000)"
z.BigNum=factorial(10000)
Print z
Print "Pi^100"
z=pi ^ 100
Print z
t=Timer
'' sin(.5)
x=".5"
y=x
x3=x
x2=x * x
term=1
fac=term
sign=-1
i=1
While term>eps
x3*= x2
fac/= (i+1)
fac/= (i+2)
term=x3 * fac
If sign<0 Then
y=y - term
Else
y+= term
End If
i+=2
sign=-sign
Wend
t=timer-t
Print "sin(.5)"
Print y
Print "time to comput sine function ";t;" seconds"
Print
Print "you can specify the number of digits to the toString function"
Print pi.toString(16)
Print
Print "Square root"
For j As bigfloat=1 To 20
y=Sqr(j)
? j.toLong, y.toString(60)
Next
Print "press return to exit"
Sleep
Last edited by srvaldez on Dec 11, 2021 15:12, edited 4 times in total.
Re: about bignum
I don't think it's worth pursuing
other than the primitive/inefficient float to string conversion I think it's ok
Last edited by srvaldez on Dec 11, 2021 14:03, edited 1 time in total.
Re: about bignum
Thanks srvaldez, nice routines.
I got a strange warning while testing , I got your factorial(10000) to check against my own, I kept your answer inside a hide macro to compare via fbide quick runs, the compiler finds my comment inside the macro and wants to truncate it?
I got a strange warning while testing , I got your factorial(10000) to check against my own, I kept your answer inside a hide macro to compare via fbide quick runs, the compiler finds my comment inside the macro and wants to truncate it?
Code: Select all
#macro hide
'2.846259680917054518906413212119868890148051401702799230794179994274411340
#endmacro
'get warning
'warning 8(2): Literal number too big, truncated
sleep
Re: about bignum
hi dodicat
thanks for testing, about your macro warning it's strange to say the lest
about the math routines being kosher, I just printed out the value of sqr(1e100+1) in hex and it's ok, it's my bigfloat to string that's at fault
I did a torture test of sorts by calculating Pi to 715000 places using the Brent_Salamin algorithm and they agree with what Maple produces but the same calculation is at least 60% faster using my decfloat routines
thanks for testing, about your macro warning it's strange to say the lest
about the math routines being kosher, I just printed out the value of sqr(1e100+1) in hex and it's ok, it's my bigfloat to string that's at fault
I did a torture test of sorts by calculating Pi to 715000 places using the Brent_Salamin algorithm and they agree with what Maple produces but the same calculation is at least 60% faster using my decfloat routines
Re: about bignum
corrected a couple of bugs
Updated
Updated
Code: Select all
Extern "c"
Declare Function __builtin_clz (Byval x As Ulong) As Long
Declare Function __builtin_sadd_overflow (Byval As Long, Byval As Long, Byref As Long) As Boolean
Declare Function __builtin_saddll_overflow (Byval As Longint, Byval As Longint, Byref As Longint) As Boolean
Declare Function __builtin_uadd_overflow (Byval As Ulong, Byval As Ulong, Byref As Ulong) As Boolean
Declare Function __builtin_uaddll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As Boolean
Declare Function __builtin_ssub_overflow (Byval As Long, Byval As Long, Byref As Long) As Boolean
Declare Function __builtin_ssubll_overflow (Byval As Longint, Byval As Longint, Byref As Longint) As Boolean
Declare Function __builtin_usub_overflow (Byval As Ulong, Byval As Ulong, Byref As Ulong) As Boolean
Declare Function __builtin_usubll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As Boolean
Declare Function __builtin_smul_overflow (Byval As Long, Byval As Long, Byref As Long) As Boolean
Declare Function __builtin_smulll_overflow (Byval As Longint, Byval As Longint, Byref As Longint) As Boolean
Declare Function __builtin_umul_overflow (Byval As Ulong, Byval As Ulong, Byref As Ulong) As Boolean
Declare Function __builtin_umulll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As Boolean
End Extern
Const NUMBER_OF_DIGITS = 1024
Const NUMBER_OF_BITS = NUMBER_OF_DIGITS*3.321928094887362
#If (NUMBER_OF_BITS Mod 32)<>0
Const NUM_BITS = (NUMBER_OF_BITS\32+1)*32
#Else
Const NUM_BITS = NUMBER_OF_BITS
#Endif
Const NUM_DWORDS = NUM_BITS\32
Const BIAS = 1073741824 '2 ^ 30
' Error definitions
#Define DIVZ_ERR 1 'Divide by zero
#Define EXPO_ERR 2 'Exponent overflow error
#Define EXPU_ERR 3 'Exponent underflow error
Type BigFloat_struct
Declare Constructor ( )
Declare Constructor ( Byref rhs As BigFloat_struct )
Declare Destructor ( )
Declare Operator Let ( Byref rhs As BigFloat_struct )
As Short sign
As Ulong exponent
As Ulong Ptr mantissa
End Type
Constructor BigFloat_struct ( )
If this.mantissa<>0 Then
Print "pointer is not 0"
End(1)
End If
this.mantissa=Callocate((NUM_DWORDS +1), 4)
If this.mantissa=0 Then
Print "unable to allocate memory"
End(4)
End If
this.sign=0
this.exponent=0
End Constructor
Constructor BigFloat_struct ( Byref rhs As BigFloat_struct )
If this.mantissa<>0 Then
Print "pointer is not 0"
End(1)
End If
this.mantissa=Callocate((NUM_DWORDS +1), 4)
If this.mantissa=0 Then
Print "unable to allocate memory"
End(4)
End If
this.sign=rhs.sign
this.exponent=rhs.exponent
For i As Long=0 To NUM_DWORDS
this.mantissa[i]=rhs.mantissa[i]
Next
End Constructor
Operator BigFloat_struct.let ( Byref rhs As BigFloat_struct )
this.sign=rhs.sign
this.exponent=rhs.exponent
For i As Long=0 To NUM_DWORDS
this.mantissa[i]=rhs.mantissa[i]
Next
End Operator
Destructor BigFloat_struct ( )
If this.mantissa=0 Then
Print "unable to de-allocate memory"
End(-4)
Else
Deallocate(this.mantissa)
End If
this.mantissa=0
End Destructor
Type BigFloat
Declare Constructor ( )
Declare Constructor ( Byval rhs As Long )
Declare Constructor ( Byval rhs As Double )
Declare Constructor ( Byref rhs As String )
Declare Constructor ( Byref rhs As BigFloat )
Declare Destructor ( )
Declare Operator Let ( Byval rhs As Long )
Declare Operator Let ( Byval rhs As Double )
Declare Operator Let ( Byref rhs As String )
Declare Operator Let ( Byref rhs As BigFloat )
Declare Operator Cast ( ) As String
Declare Operator Cast ( ) As Long
Declare Operator Cast ( ) As Double
'----------------------------------------------
Declare Operator += (Byref rhs As BigFloat)
Declare Operator += (Byval rhs As Long)
Declare Operator += (Byval rhs As Double)
Declare Operator += (Byref rhs As String)
Declare Operator -= (Byref rhs As BigFloat)
Declare Operator -= (Byval rhs As Long)
Declare Operator -= (Byval rhs As Double)
Declare Operator -= (Byref rhs As String)
Declare Operator *= (Byref rhs As BigFloat)
Declare Operator *= (Byval rhs As Long)
Declare Operator *= (Byval rhs As Double)
Declare Operator *= (Byref rhs As String)
Declare Operator /= (Byref rhs As BigFloat)
Declare Operator /= (Byval rhs As Long)
Declare Operator /= (Byval rhs As Double)
Declare Operator /= (Byval rhs As Single)
Declare Operator /= (Byref rhs As String)
' For Next Implicit step = +1
Declare Operator For ( )
Declare Operator Step( )
Declare Operator Next( Byref end_cond As BigFloat ) As Integer
' For Next Exlicit step
Declare Operator For ( Byref stp As BigFloat )
Declare Operator Step( Byref stp As BigFloat )
Declare Operator Next( Byref end_cond As BigFloat, Byref step_var As BigFloat ) As Integer
Declare Function toString( Byval places As Long=NUM_DWORDS*9.63 ) As String
Declare Function toLong ( ) As Long
Declare Function toDouble ( ) As Double
BigNum As BigFloat_struct
End Type
Declare Function fpcmp(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
Declare Function fpadd(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpsub(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpmul(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpmul_si(Byref x As BigFloat_struct, Byval y As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpdiv_si(Byref num As BigFloat_struct, Byval den As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpdiv(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpipow(Byref x As BigFloat_struct, Byval e As Longint, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function dbl2fp(Byval x As Double) As BigFloat_struct
Declare Function fp2dbl(Byref x As BigFloat_struct) As Double
Declare Function ui2fp(Byval m As Ulong, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fp2ui(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Ulong
Declare Function si2fp(Byval m As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fp2si(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
Declare Function str2fp(Byref x As String) As BigFloat_struct
Declare Function fp2str(Byref zz As BigFloat_struct, Byval digits As Long=NUM_DWORDS*9.63) As String
Declare Function fpfix( Byref num As BigFloat_struct ) As BigFloat_struct
Declare Function fpfrac( Byref num As BigFloat_struct ) As BigFloat_struct
Declare Function fpsqr (Byref num As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
Declare Function mp_sin(Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
Declare Function mp_cos (Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
Declare Function mp_tan(Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
Declare Function fpatn (Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
Declare Function fpasin (Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
Declare Function fpacos (Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
Function Bigfloat.toString( Byval places As Long=NUM_DWORDS*9.63 ) As String
Function = fp2str( this.BigNum, places )
End Function
Function Bigfloat.toLong ( ) As Long
Dim As Double x
x=fp2dbl(this.BigNum)
Return Clng(x)
End Function
Function Bigfloat.toDouble ( ) As Double
Function = fp2dbl(this.BigNum)
End Function
Constructor Bigfloat ( )
this.BigNum=si2fp(0, NUM_DWORDS)
End Constructor
Constructor Bigfloat ( Byval rhs As Long )
this.BigNum=si2fp( rhs, NUM_DWORDS )
End Constructor
Constructor Bigfloat ( Byref rhs As String )
this.BigNum=str2fp( rhs )
End Constructor
Constructor Bigfloat ( Byref rhs As Bigfloat)
this.BigNum.sign=rhs.BigNum.sign
this.BigNum.exponent=rhs.BigNum.exponent
For i As Long=0 To NUM_DWORDS
this.BigNum.mantissa[i]=rhs.BigNum.mantissa[i]
Next
End Constructor
Constructor Bigfloat ( Byval rhs As Double )
this.BigNum=str2fp( Str( rhs ) )
End Constructor
Operator Bigfloat.let ( Byref rhs As Bigfloat )
this.BigNum.sign=rhs.BigNum.sign
this.BigNum.exponent=rhs.BigNum.exponent
For i As Long=0 To NUM_DWORDS
this.BigNum.mantissa[i]=rhs.BigNum.mantissa[i]
Next
End Operator
Operator Bigfloat.let ( Byval rhs As Long )
this.BigNum=si2fp( rhs, NUM_DWORDS )
End Operator
Operator Bigfloat.let ( Byref rhs As String )
this.BigNum=str2fp( rhs )
End Operator
Operator Bigfloat.Let ( Byval rhs As Double )
this.BigNum=str2fp( Str(rhs) )
End Operator
Operator Bigfloat.cast ( ) As String
Operator = fp2str(this.BigNum)
End Operator
Operator Bigfloat.cast ( ) As Long
Dim As Double x
x=fp2dbl(this.BigNum)
Operator = Clng(x)
End Operator
Operator Bigfloat.cast ( ) As Double
Operator = fp2dbl(this.BigNum)
End Operator
Destructor Bigfloat ( )
End Destructor
'============================================================================
Operator BigFloat.for ( )
End Operator
Operator BigFloat.step ( )
this.BigNum = fpadd(this.BigNum,si2fp(1, NUM_DWORDS), NUM_DWORDS)
End Operator
Operator BigFloat.next ( Byref end_cond As BigFloat ) As Integer
Return fpcmp(This.BigNum, end_cond.BigNum, NUM_DWORDS)<=0
End Operator
'' explicit step versions
''
Operator BigFloat.for ( Byref step_var As BigFloat )
End Operator
Operator BigFloat.step ( Byref step_var As BigFloat )
this.BigNum = fpadd(this.BigNum, step_var.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.next ( Byref end_cond As BigFloat, Byref step_var As BigFloat ) As Integer
If fpcmp(step_var.BigNum, si2fp(0, NUM_DWORDS), NUM_DWORDS) < 0 Then
Return fpcmp(This.BigNum, end_cond.BigNum, NUM_DWORDS) >= 0
Else
Return fpcmp(This.BigNum, end_cond.BigNum, NUM_DWORDS) <= 0
End If
End Operator
'============================================================================
Operator + ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum= fpadd(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator + ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat
Dim As Bigfloat result
result.BigNum=dbl2fp( rhs )
result.BigNum= fpadd(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator + ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=dbl2fp( lhs )
result.BigNum= fpadd( result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator + ( Byref lhs As Bigfloat, Byval rhs As Long ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp( rhs )
result.BigNum= fpadd(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator + ( Byval lhs As Long, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp( lhs )
result.BigNum= fpadd( result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator + (Byref lhs As Bigfloat, Byref rhs As String) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp(rhs)
result.BigNum= fpadd(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - ( Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result = rhs
result.BigNum.sign = result.BigNum.sign Xor &h8000
Operator = result
End Operator
Operator - ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum= fpsub(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat
Dim As Bigfloat result
result.BigNum=dbl2fp( rhs )
result.BigNum= fpsub(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=dbl2fp( lhs )
result.BigNum= fpsub( result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - ( Byref lhs As Bigfloat, Byval rhs As Long ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
result.BigNum= fpsub(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - ( Byval lhs As Long, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(lhs, NUM_DWORDS)
result.BigNum= fpsub(result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator - (Byref lhs As Bigfloat, Byref rhs As String) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp(rhs)
result.BigNum= fpsub(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum= fpmul(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat
Dim As Bigfloat result
result.BigNum=dbl2fp( rhs )
result.BigNum= fpmul(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=dbl2fp( lhs )
result.BigNum= fpmul( result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * ( Byref lhs As Bigfloat, Byval rhs As Long ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
result.BigNum= fpmul(result.BigNum, lhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * ( Byval lhs As Long, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(lhs, NUM_DWORDS)
result.BigNum= fpmul(result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator * (Byref lhs As Bigfloat, Byref rhs As String) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp(rhs)
result.BigNum= fpmul(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum= fpdiv(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat
Dim As Bigfloat result
result.BigNum=dbl2fp( rhs )
result.BigNum= fpdiv(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=dbl2fp( lhs )
result.BigNum= fpdiv( result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / ( Byref lhs As Bigfloat, Byval rhs As Long ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
result.BigNum= fpdiv(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / ( Byval lhs As Long, Byref rhs As Bigfloat ) As Bigfloat
Dim As Bigfloat result
result.BigNum=si2fp(lhs, NUM_DWORDS)
result.BigNum= fpdiv(result.BigNum, rhs.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator / (Byref lhs As Bigfloat, Byref rhs As String) As Bigfloat
Dim As Bigfloat result
result.BigNum=str2fp(rhs)
result.BigNum= fpdiv(lhs.BigNum, result.BigNum, NUM_DWORDS)
Operator = result
End Operator
Operator ^ ( Byref lhs As Bigfloat, Byval rhs As Longint ) As Bigfloat
Dim As Bigfloat lhs2
lhs2.BigNum=fpipow(lhs.BigNum, rhs, NUM_DWORDS)
Operator = lhs2
End Operator
Operator BigFloat.+= (Byref rhs As BigFloat)
this.BigNum = fpadd(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.+= (Byval rhs As Long)
Dim As BigFloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.+= (Byval rhs As Double)
Dim As BigFloat result
result.BigNum=dbl2fp( rhs )
this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.+= (Byref rhs As String)
Dim As BigFloat result
result.BigNum=str2fp(rhs)
this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.-= (Byref rhs As BigFloat)
this.BigNum = fpsub(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.-= (Byval rhs As Long)
Dim As BigFloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.-= (Byval rhs As Double)
Dim As BigFloat result
result.BigNum=dbl2fp( rhs )
this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.-= (Byref rhs As String)
Dim As BigFloat result
result.BigNum=str2fp(rhs)
this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.*= (Byref rhs As BigFloat)
this.BigNum = fpmul(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.*= (Byval rhs As Long)
Dim As BigFloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.*= (Byval rhs As Double)
Dim As BigFloat result
result.BigNum=dbl2fp( rhs )
this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat.*= (Byref rhs As String)
Dim As BigFloat result
result.BigNum=str2fp(rhs)
this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat./= (Byref rhs As BigFloat)
this.BigNum = fpdiv(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat./= (Byval rhs As Long)
Dim As BigFloat result
result.BigNum=si2fp(rhs, NUM_DWORDS)
this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat./= (Byval rhs As Double)
Dim As BigFloat result
result.BigNum=dbl2fp( rhs )
this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator BigFloat./= (Byref rhs As String)
Dim As BigFloat result
result.BigNum=str2fp(rhs)
this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator
Operator = ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)=0
End Operator
Operator < ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<0
End Operator
Operator > ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)>0
End Operator
Operator <= ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<=0
End Operator
Operator >= ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)>=0
End Operator
Operator <> ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<>0
End Operator
Operator Abs(Byref rhs As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=rhs.BigNum
result.BigNum.sign = 0
Operator = result
End Operator
Operator Fix (Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum =fpfix(x.BigNum)
Operator = result
End Operator
Operator Frac(Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=fpfrac(x.BigNum)
Operator = result
End Operator
Operator Sqr(Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=fpsqr(x.BigNum)
Operator = result
End Operator
Operator Sin(Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=mp_sin(x.BigNum)
Operator = result
End Operator
Operator Cos(Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=mp_cos(x.BigNum)
Operator = result
End Operator
Operator Tan(Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=mp_tan(x.BigNum)
Operator = result
End Operator
Operator Atn(Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=fpatn(x.BigNum)
Operator = result
End Operator
Operator Asin(Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=fpasin(x.BigNum)
Operator = result
End Operator
Operator Acos(Byref x As BigFloat) As BigFloat
Dim As BigFloat result
result.BigNum=fpatn(x.BigNum)
Operator = result
End Operator
Dim Shared As BigFloat_struct Log102
Log102.sign=0
Log102.exponent=(-2)+BIAS+1
If NUM_DWORDS>=0 Then Log102.mantissa[0]=&h13441350
If NUM_DWORDS>=1 Then Log102.mantissa[1]=&h9F79FEF3
If NUM_DWORDS>=2 Then Log102.mantissa[2]=&h11F12B35
If NUM_DWORDS>=3 Then Log102.mantissa[3]=&h816F922F
Private Function log2_32(Byval value As Ulong) As Long
'https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers
Static tab32(0 To 31) As Const Long = {0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30, 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31}
value Or= Culng(value Shr 1)
value Or= Culng(value Shr 2)
value Or= Culng(value Shr 4)
value Or= Culng(value Shr 8)
value Or= Culng(value Shr 16)
Return tab32(Culng(Culng(value * &h07C4ACDD) Shr 27))
End Function
Private Function log2_64(Byval value As Ulongint) As Long
'https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers
Static tab64(0 To 63) As Const Long ={63, 0, 58, 1, 59, 47, 53, 2, 60, 39, 48, 27, 54, 33, 42, 3, 61, 51, 37, 40, 49, 18, 28, 20, 55, 30, 34, 11, 43, 14, 22, 4, 62, 57, 46, 52, 38, 26, 32, 41, 50, 36, 17, 19, 29, 10, 13, 21, 56, 45, 25, 31, 35, 16, 9, 12, 44, 24, 15, 8, 23, 7, 6, 5}
value Or= Culngint(value Shr 1)
value Or= Culngint(value Shr 2)
value Or= Culngint(value Shr 4)
value Or= Culngint(value Shr 8)
value Or= Culngint(value Shr 16)
value Or= Culngint(value Shr 32)
Return tab64((Culngint((value - (value Shr 1)) * &h07EDD5E59A4E28C2) Shr 58))
End Function
Function ipower (Byval x As Ulongint, Byval e As Ulongint) As Ulongint
'take x to an integer power
Dim As Ulongint z, y, n
y = x
n = e
z = 1
While n > 0
While (n And 1) = 0
n = n \ 2
y = y * y
Wend
n = n - 1
z = y * z
Wend
Return z
End Function
Function shl32(Byval n As Ulong, Byval k As Ubyte, Byref c As Ulong) As Ulong
If k>0 And k<32 Then
Dim As Ulong carry=n
Dim As Ubyte k32=32-k
Asm
mov cl, [k32]
Shr dword Ptr [carry], cl
mov cl, [k]
Shl dword Ptr [n], cl
End Asm
c=carry
End If
Return n
End Function
Function shr32(Byval n As Ulong, Byval k As Ubyte, Byref c As Ulong) As Ulong
If k>0 And k<32 Then
Dim As Ulong carry=n
Dim As Ubyte k32=32-k
Asm
mov cl, [k32]
Shl dword Ptr [carry], cl
mov cl, [k]
Shr dword Ptr [n], cl
End Asm
c=carry
End If
Return n
End Function
Sub shiftl(Byref fac1 As BigFloat_struct, Byval k As Long, Byval dwords As Long=NUM_DWORDS)
If k>0 And k<32 Then
Dim As Long i
Dim As Ulong n, carry, c=0
Dim As Long k32=32-k
For i=dwords To 0 Step -1
n=fac1.mantissa[i]
carry=n
carry=Cast(Ulong, (carry Shr k32))
n=Cast(Ulong, (n Shl k))
fac1.mantissa[i]=n+c
c=carry
Next
Elseif k=32 Then
Dim As Long i
For i=0 To dwords-1
fac1.mantissa[i]=fac1.mantissa[i+1]
Next
fac1.mantissa[dwords]=0
End If
End Sub
Sub shiftr(Byref fac1 As BigFloat_struct, Byval k As Long, Byval dwords As Long=NUM_DWORDS)
If k>0 And k<32 Then
Dim As Long i
Dim As Ulong n, carry, c=0
Dim As Long k32=32-k
For i=0 To dwords
n=fac1.mantissa[i]
carry=n
carry=Cast(Ulong, (carry Shl k32))
n=Cast(Ulong, (n Shr k))
fac1.mantissa[i]=c+n
c=carry
Next
Elseif k=32 Then
Dim As Long i
For i=dwords To 1 Step -1
fac1.mantissa[i]=fac1.mantissa[i-1]
Next
fac1.mantissa[0]=0
End If
End Sub
Function fpcmp(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As Long i
Dim As Longint c
If x.sign < y.sign Then
Return -1
End If
If x.sign > y.sign Then
Return 1
End If
If x.exponent<y.exponent Then
If x.sign=0 Then
Return -1
Else
Return 1
End If
End If
If x.exponent>y.exponent Then
If x.sign=0 Then
Return 1
Else
Return -1
End If
End If
For i=0 To dwords
c=Clngint(x.mantissa[i])-clngint(y.mantissa[i])
If c<>0 Then Exit For
Next
If c=0 Then Return 0
If c<0 Then
If x.sign=0 Then
Return -1
Else
Return 1
End If
End If
If c>0 Then
If x.sign=0 Then
Return 1
Else
Return -1
End If
End If
End Function
Function NORM_FAC1(Byref fac1 As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Integer
If dwords>NUM_DWORDS Then dwords=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.
Dim As Integer i,er=0,f=0
For i=0 To dwords
If fac1.mantissa[i]>0 Then f=1
Next
If f=0 Then
fac1.exponent=0
fac1.sign=0
Exit Function
End If
While fac1.mantissa[0]=0
shiftl(fac1, 32, dwords)
fac1.exponent-=32
If fac1.exponent=0 Then
NORM_FAC1=EXPU_ERR
Exit Function
End If
Wend
i=__builtin_clz(fac1.mantissa[0])-3
If i>0 Then
shiftl(fac1, i, dwords)
fac1.exponent-=i
End If
'if the highmost Bit in fac1_man is nonzero,
'shift the mantissa right 1 Bit and
'increment the exponent
If fac1.mantissa[0]>(&h1ffffffful)Then
While fac1.mantissa[0]>&h1ffffffful
shiftr(fac1, 1, dwords)
fac1.exponent+=1
Wend
Elseif fac1.mantissa[0]<&h10000000ul Then
/' the following will probably never be executed '/
'now shift fac1_man 1 to the left until a
'nonzero Bit appears in the next-to-highest
'Bit of fac1_man. decrement exponent for
'each shift.
While fac1.mantissa[0]<&h10000000ul
shiftl(fac1, 1, dwords)
fac1.exponent-=1
Wend
End If
'check for overflow/underflow
If fac1.exponent<0 Then
Print "NORM_FAC1 fac1.exponent<0"
NORM_FAC1=EXPO_ERR
End If
End Function
Function si2fp(Byval m As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1
Dim As Long n=Abs(m)
If m=0 Then
Return fac1
End If
fac1.mantissa[0]=n
NORM_FAC1(fac1, NUM_DWORDS)
fac1.exponent=log2_32(n)+BIAS+1
If m<0 Then
fac1.sign=&H8000
Else
fac1.sign=0
End If
Return fac1
End Function
Function ui2fp(Byval m As Ulong, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1
If m=0 Then
Return fac1
End If
fac1.mantissa[0]=m
NORM_FAC1(fac1, NUM_DWORDS)
fac1.exponent=log2_32(m)+BIAS+1
Return fac1
End Function
Sub fpadd_aux(Byref fac1 As BigFloat_struct, Byref fac2 As BigFloat_struct, Byval dwords As Long=NUM_DWORDS)
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As Long c, i
Dim As Longint v
c=0
For i=dwords To 1 Step -1
v=Clngint(fac2.mantissa[i])+Clngint(fac1.mantissa[i])+c
If v>(&h100000000ull-1) Then
v=v-&h100000000ull
c=1
Else
c=0
End If
fac1.mantissa[i]=v
Next
v=Clngint(fac1.mantissa[0])+Clngint(fac2.mantissa[0])+c
fac1.mantissa[0]=v
NORM_FAC1(fac1, dwords)
End Sub
Sub fpsub_aux(Byref fac1 As BigFloat_struct, Byref fac2 As BigFloat_struct, Byval dwords As Long=NUM_DWORDS)
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As Long c, i
Dim As Longint v
c=0
For i=dwords To 1 Step -1
v=Clngint(fac1.mantissa[i])-clngint(fac2.mantissa[i])-c
If v<0 Then
v=v+&h100000000ull
c=1
Else
c=0
End If
fac1.mantissa[i]=v
Next
v=Clngint(fac1.mantissa[0])-clngint(fac2.mantissa[0])-c
fac1.mantissa[0]=v
NORM_FAC1(fac1, dwords)
End Sub
Function fpadd(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1,fac2
Dim As Long i, t, c, xsign, ysign
xsign=x.sign:x.sign=0
ysign=y.sign:y.sign=0
c=fpcmp(x, y, dwords)
x.sign=xsign
y.sign=ysign
If c<0 Then
fac1=y
fac2=x
Else
fac1=x
fac2=y
End If
t=fac1.exponent-fac2.exponent
If t<(NUM_BITS+32) Then
'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 dwords, the result is already in FAC1.
If t>0 And t<(NUM_BITS+32) Then 'shift
i=t\32
While i>0
shiftr(fac2, 32, dwords)
t-=32
i-=1
Wend
'While t>0
If t>0 Then shiftr(fac2, t, dwords)
't-=1
'Wend
End If
'See if the signs of the two numbers
'are the same. If so, add; if not, subtract.
If fac1.sign=fac2.sign Then 'add
fpadd_aux(fac1,fac2, dwords)
Else
fpsub_aux(fac1,fac2, dwords)
End If
Endif
Return fac1
End Function
Function fpsub(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1,fac2
Dim As Long s
fac1=x
fac2=y
fac2.sign=fac2.sign Xor &h8000
fac1=fpadd(fac1,fac2, dwords)
Return fac1
End Function
Function fpmul(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1,fac2
Dim As Integer i, j, ex, er, den, num
Dim As Ulongint digit, carry, prod ', tmp
Dim As Ulong fac3(0 To 2*dwords+1)
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
For i=0 To dwords
fac1.mantissa[i]=0
Next
'NORM_FAC1(fac1)
Return fac1
Else
If ex<0 Then
er=EXPO_ERR
Return fac1 'Exit Function
End If
'clear fac3 mantissa
For i=0 To dwords
fac3(i)=0
Next
den=dwords
While fac2.mantissa[den]=0
den-=1
Wend
num=dwords
While fac1.mantissa[num]=0
num-=1
Wend
If num<den Then
Swap fac1, fac2
'fac1=y
'fac2=x
Swap den, num
End If
For j=den To 0 Step -1
carry=0
digit=fac2.mantissa[j]
For i=num To 0 Step -1
prod=fac3(i+j+1)+digit*fac1.mantissa[i]+carry
/'tmp=fac1.mantissa[i]
If __builtin_umulll_overflow (digit, tmp, prod) Then Print "mul overflow"
tmp=fac3(i+j+1)
If __builtin_uaddll_overflow (tmp, prod, prod) Then Print "add overflow"
If __builtin_uaddll_overflow (carry, prod, prod) Then Print "add overflow"'/
carry=prod Shr 32'\&H100000000
fac3(i+j+1)=prod '(prod mod &H100000000)
Next
fac3(j)=carry
Next
For i=0 To dwords
fac1.mantissa[i]=fac3(i)
Next
End If
'now determine exponent of result.
'as you do...watch for overflow.
'ex=fac2.exponent-BIAS+fac1.exponent
'fac1.exponent=ex
ex=(fac2.exponent And &h7FFFFFFFul)-BIAS+1
ex=ex+(fac1.exponent And &h7FFFFFFFul)-BIAS+1
fac1.exponent=ex+BIAS+1
'determine the sign of the product
fac1.sign=fac1.sign Xor fac2.sign
NORM_FAC1(fac1, dwords)
Return fac1
End Function
Function fpmul_si(Byref x As BigFloat_struct, Byval y As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1,fac2
fac1=x
fac2=ui2fp(y, dwords)
fac1=fpmul(fac1, fac2, dwords)
Return fac1
End Function
Function fpdiv_si(Byref num As BigFloat_struct, Byval den As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1
Dim As Ulongint carry, remder=0
Dim As Longint i, divisor
Dim As Longint quotient
Dim As Long j
divisor=Abs(den)
fac1=num
If divisor=1 Then
Return fac1
End If
If divisor = 0 Then
Print "error: divisor = 0"
Return fac1 'Exit function
End If
If divisor>2147483647 Then
Print "error: divisor too large"
Return fac1 'Exit function
End If
For i = 0 To dwords
quotient = fac1.mantissa[i] + remder * &h100000000ull
remder = quotient Mod divisor
fac1.mantissa[i]=quotient \ divisor
Next
quotient = remder * &h100000000ull
quotient=quotient \ divisor
carry=fac1.mantissa[0]
j=__builtin_clz(fac1.mantissa[0])-3
If j>0 Then
shiftl(fac1, j, dwords)
fac1.exponent-=j
End If
While fac1.mantissa[0]>(&h1ffffffful)
shiftr(fac1, 1, dwords)
fac1.exponent+=1
Wend
NORM_FAC1(fac1)
If den<0 Then
fac1.sign=fac1.sign Xor &h8000
End If
Return fac1
End Function
Function fpdiv(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1, fac2, one
Dim As Integer i, er, is_power_of_two
Dim As Short sign
fac1=x
fac2=y
one.exponent=(BIAS+1)
one.mantissa[0]=&h10000000ul
one=ui2fp(1)
sign=fac2.sign
fac2.sign=0
If fpcmp(fac2, one)=0 Then
Return fac1
Else
fac2.sign=sign
End If
If fac2.exponent=0 Then ' if fac2 = 0, return
' a divide-by-zero error and
' bail out.
fac1.mantissa[i]=&H1FFFFFFFul
For i=1 To dwords
fac1.mantissa[i]=&HFFFFFFFFul
Next
fac1.exponent=&HFFFFF+BIAS+1
er=DIVZ_ERR
Return fac1
Elseif fac1.exponent=0 Then 'fact1=0, just return
er=0
Return fac1
Else
'check to see if fac2 is a power of ten
is_power_of_two=0
If fac2.mantissa[0]=&H10000000ul Then
is_power_of_two=1
For i=1 To dwords
If fac2.mantissa[i]<>0 Then
is_power_of_two=0
Exit For
End If
Next
End If
'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_two=1 Then
fac1.sign=fac1.sign Xor fac2.sign
fac1.exponent=fac1.exponent-fac2.exponent+BIAS+1
Return fac1
End If
#Define min(a,b) Iif(((a)<(b)),(a),(b))
#Macro realw(w, j)
((w(j - 1)*b + w(j))*b + w(j + 1))*b +Iif(Ubound(w)>=j+2,w(j+2),0)
#Endmacro
#Macro subtract(w, q, d, ka, kb)
For j=ka To kb
w(j) = w(j) - q*d(j - ka + 2)
Next
#Endmacro
#Macro normalize(w, ka, q)
w(ka) = w(ka) + w(ka - 1)*b
w(ka - 1) = q
#Endmacro
#Macro finalnorm(w, kb)
For j=kb To 3 Step -1
carry=Iif(w(j)<0, ((-w(j) - 1)\b) + 1, Iif(w(j) >= b, -(w(j)\b), 0))
w(j) = w(j) + carry*b
w(j - 1) = w(j - 1) - carry
Next
#Endmacro
Dim As Double result(1 To 2*dwords+3), n(1 To 2*dwords+3), d(1 To 2*dwords+3)
Const b=&H10000
Dim As Long j, last, laststep, q, t
Dim As Long stp, carry
Dim As Double xd, xn, rund
Dim As Double w(1 To Ubound(n)+4)
For j=0 To dwords
n(2*j+2)=fac1.mantissa[j]\&H10000
n(2*j+3)=fac1.mantissa[j] Mod &H10000
d(2*j+2)=fac2.mantissa[j]\&H10000
d(2*j+3)=fac2.mantissa[j] Mod &H10000
Next
n(1)=(fac1.exponent And &h7FFFFFFF)-BIAS-1
d(1)=(fac2.exponent And &h7FFFFFFF)-BIAS-1
For j=Ubound(n) To Ubound(w)
w(j)=0
Next
t=Ubound(n)-1
w(1)=n(1)-d(1)+1
w(2)=0
For j=2 To Ubound(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, Ubound(W))
subtract(w, q, d, (stp + 2), last)
normalize(w, (stp + 2), q)
Next
FinalNorm(w, (laststep + 1))
laststep = Iif(w(2) = 0, laststep, laststep - 1)
rund = w(laststep + 1)/b
w(laststep) = w(laststep) + Iif(rund >= 0.5, 1, 0)
If w(2)=0 Then
For j=1 To t+1
result(j)=w(j+1)
Next
Else
For j=1 To t+1
result(j)=w(j)
Next
End If
result(1) = Iif(w(2) = 0, w(1) - 1, w(1))
For j=0 To dwords
fac1.mantissa[j]=result(2*j+2)*&H10000+result(2*j+3)
Next
NORM_FAC1(fac1, dwords)
fac1.exponent=(result(1)+BIAS)
End If
fac1.sign=fac1.sign Xor fac2.sign
Return fac1
End Function
Function fpipow(Byref x As BigFloat_struct, Byval e As Longint, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
'take x to an Long power
Dim As BigFloat_struct y=x
Dim As BigFloat_struct z, one
Dim As Longint n, c=0
Dim As Integer i
n = Abs(e)
one.exponent=(BIAS+1)
one.mantissa[0]=&h10000000ul
z=one
While n>0
While (n And 1)=0
n\=2
y=fpmul(y, y, dwords)
c+=1
Wend
n-=1
z=fpmul(y, z, dwords)
c+=1
Wend
If e<0 Then
z=fpdiv(one, z, dwords)
End If
Return z
End Function
' sqrt(num)
Function fpsqr (Byref num As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As bigfloat_struct r, r2, tmp, n, half, one, eps, zero
Dim As Long ex, ex2, ex3, k, l, prec
Dim As String ts, v
Dim As Double x
one.exponent=(BIAS+1)
one.mantissa[0]=&h10000000ul
eps.exponent=-(NUM_BITS)+BIAS+1
eps.mantissa[0]=&h10000000ul
l=Log((NUM_DWORDS*9.63+9.63)*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 bigfloat has an extra 9 guard digits
n=num
If fpcmp(n, zero, dwords)=0 Then
Return r ' r is stil 0
End If
If fpcmp(n, one, dwords)=0 Then
Return one
End If
If fpcmp(n, si2fp(0, dwords), dwords)<0 Then
Return r ' r is stil 0
End If
'=====================================================================
'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 bigfloat
'divide the exponent in half for square root
'in this case 1.414213562373095e250
If n.exponent<>0 Then
ex=(n.exponent And &h7FFFFFFFul)-BIAS-1
Else
ex=0
End If
ts=fp2str(n, 16)
k=Instr(ts, "e")
v=Mid(ts,k+1)
ts=Left(ts, k-1)
ex2=Clng(v)
x=Val(ts)
If x = 0 Then Print "Div 0": Return r 'Exit Function
If x = 1 Andalso ex=0 Then
Return one
End If
If Abs(ex2) And 1 Then
x=x*10
ex2-=1
End If
x = Sqr(x) 'approximation
v=Trim(Str(x))
k=Instr(v,".")
If Len(v)>1 And k=0 Then ex2+=1
v=v+"e"+Trim(Str(ex2\2))
r=str2fp(v)
half.mantissa[0]=&h10000000
half.exponent=BIAS
half.sign=0
'=====================================================================
'Newton-Raphson method
prec=3
For k=1 To l+2
prec=2*prec-1
tmp = fpdiv(n, r, prec)
r2 = fpadd(r, tmp, prec)
r = fpmul(r2, half, prec)
Next
Return r
End Function
Function dbl2fp(Byval x As Double) As BigFloat_struct
Dim As BigFloat_struct bf
Dim As Ulongint n, e
dim as long sign
if x<0 then sign=&h8000
x=abs(x)
n=Peek ( Ulongint, @x )
e=n
n=Cast(Ulongint, (n Shl 12))
n=Cast(Ulongint, (n Shr 4))
n+=&h1000000000000000ull
bf.mantissa[0]=Cast(Ulongint, (n Shr 32))
n=e
n=Cast(Ulongint, (n Shl 40))
n=Cast(Ulongint, (n Shr 32))
bf.mantissa[1]=n
n=n+&h1000000000000000ull
e=Cast(Ulongint, (e Shr 52))
e=e-&h3FFul
bf.exponent=e+BIAS+1
bf.sign=sign
Return bf
End Function
Function fp2dbl(Byref x As BigFloat_struct) As Double
Dim As Double dbl
Dim As Ulongint n, e
Dim As Ushort ex
Dim As Long sign
If x.sign<>0 Then
sign=-1
Else
sign=1
End If
If x.exponent>0 Then
ex=(x.exponent And &h7FFFFFFF)-BIAS-1
Else
ex=0
End If
n=x.mantissa[0]
n-=&h10000000ul
n=Cast(Ulongint, (n Shl 24))
e=&h3FF+ex
e=e*&h10000000000000ull
n=n+e
e=x.mantissa[1]\&h100
n=n+e
Poke Ulongint,@dbl,n
Return dbl*sign
End Function
Function factorial(Byval n As Ulong) As BigFloat_struct
Dim As BigFloat_struct one, f, indx
one.mantissa[0]=&h10000000ul
one.exponent=0+BIAS+1
f=one
indx=one
For i As Long=1 To n
f=fpmul(f, indx)
indx=fpadd(one, indx)
Next
Return f
End Function
Function factorial_inv(Byval n As Ulong) As BigFloat_struct
Dim As BigFloat_struct f
f.mantissa[0]=&h10000000ul
f.exponent=0+BIAS+1
For i As Long=1 To n
f=fpdiv_si(f, i)
Next
Return f
End Function
Function fp2ui(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Ulong
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1
Dim As Double n
fac1.exponent=0
fac1.mantissa[0]=x.mantissa[0]
n=fp2dbl(x)
Return Fix(n)
End Function
Function fp2si(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
If dwords>NUM_DWORDS Then dwords=NUM_DWORDS
Dim As BigFloat_struct fac1
Dim As Long sign
Dim As Double n
fac1.exponent=0
fac1.mantissa[0]=x.mantissa[0]
n=fp2dbl(x)
Return Fix(n)
End Function
Function str2fp_aux(Byval value As String, Byref s As Long, Byref ex As Long ) As String
Dim As Integer j,d,e,ep,es,i,f,fp,fln
Dim As String c,f1,f2,f3, ts
Dim As Ulong ulng
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 Then
If c=" " Then
j=j+1
Continue While
Endif
If c="-" Then
es=-es
c=""
Endif
If c="+" Then
j=j+1
Continue While
Endif
If (c="0") And (f3="") Then
j=j+1
Continue While
Endif
If (c>"/") And (c<":") Then 'c is digit between 0 and 9
f3=f3+c
ex=10*ex+(Asc(c)-48)
j=j+1
Continue While
Endif
Endif
If c=" " Then
j=j+1
Continue While
Endif
If c="-" Then
s=-s
j=j+1
Continue While
Endif
If c="+" Then
j=j+1
Continue While
Endif
If c="." Then
If d=1 Then
j=j+1
Continue While
Endif
d=1
Endif
If (c>"/") And (c<":") Then 'c is digit between 0 and 9
If ((c="0") And (i=0)) Then
If d=0 Then
j=j+1
Continue While
Endif
If (d=1) And (f=0) Then
e=e-1
j=j+1
Continue While
Endif
Endif
If d=0 Then
f1=f1+c
i=i+1
Else
If (c>"0") Then
fp=1
Endif
f2=f2+c
f=f+1
Endif
Endif
If c="E" Or c="D" Then
ep=1
Endif
j=j+1
Wend
If fp=0 Then
f=0
f2=""
Endif
If s=-1 Then s=&h8000 Else s=0
f1=f1+f2
f1=Mid(f1,1,1)+Right(f1,Len(f1)-1)
fln=Len(f1)
If fp=0 Andalso ex=0 Then
ex=es*fln-1
Else
ex=es*ex-1+i+e
End If
If fln>((NUM_BITS+32)*0.3010299956639811) Then
f1=Mid(f1,1,((NUM_BITS+32)*0.3010299956639811)-1)
Endif
Return f1
End Function
Function str2fp(Byref x As String) As BigFloat_struct
Dim As Long strlen, i, n, ex, sign
Dim As String s2, strin=x
Dim As BigFloat_struct y, z, pw, ten
strin=str2fp_aux(strin, sign, ex)
strlen=Len(strin)
n=1
ten=ui2fp(Valuint("1"+String(n,"0")))
If strlen>n Then
s2=Left(strin, n)
strin=Mid(strin, n+1)
Else
z=ui2fp(Valuint(strin))
If ex<>0 Then
pw=ui2fp(10)
pw=fpipow(pw, strlen-ex-1)
z=fpdiv(z, pw)
End If
Return z
End If
z=ui2fp(Valuint(s2))
While Len(strin)>=n
s2=Left(strin, n)
strin=Mid(strin, n+1)
y=ui2fp(Valuint(s2))
z=fpmul(z, ten)
z=fpadd(z,y)
Wend
pw=ui2fp(10)
pw=fpipow(pw,strlen-ex-1)
z=fpdiv(z, pw)
z.sign=sign
Return z
End Function
Function fp2str(Byref zz As BigFloat_struct, Byval digits As Long=NUM_DWORDS*9.63) As String
If digits>NUM_DWORDS*9.63 Or digits<0 Then digits=NUM_DWORDS*9.63
Dim As Long ex, i, m, n, powten, sign, tn
Dim As BigFloat_struct y, ten, tenp, z, ep, one, nine, zero
Dim As String s, s2
If zz.exponent=0 Then
Return "0"
End If
digits+=2
n=8
one=ui2fp(1)
nine=ui2fp(9)
ten=ui2fp(10)
tenp=ui2fp(Clng("1"+String(n,"0")))
tn=Clng("1"+String(n,"0"))
z=zz
sign=zz.sign
z.sign=0
If z.exponent<>0 Then
ex=(z.exponent And &h7FFFFFFFul)-BIAS-1
Else
ex=0
End If
m=NUM_BITS-ex
ep.exponent=-(m)+BIAS+1
ep.mantissa[0]=&h10000000ul
z=fpadd(z,ep)
If fpcmp(z,nine)>0 Then
While fpcmp(z,nine)>0
z=fpdiv_si(z, 10)
powten+=1
Wend
Elseif fpcmp(z,one)<0 Then
While fpcmp(z,one)<0
z=fpmul(z, ten)
powten-=1
Wend
End If
m=fp2ui(z)
z=fpfrac(z)
z=fpmul(z, tenp)
s=Trim(Str(m))+"."
For i=1 To digits
If fpcmp(z, one)<0 Then
While fpcmp(z, one)<0
If (Len(s))>=digits Then Exit While
z=fpmul(z, tenp)
s=s+String(n,"0")
Wend
Else
If (Len(s))>=digits Then Exit For
y=fpfix(z)
m=fp2ui(y)
z=fpfrac(z)
z=fpmul(z, tenp)
s2=Trim(Str(m))
s2=String(n-len(s2),"0")+s2
s=s+s2
End If
Next
s=Ltrim(s, "0")
If Len(s)>(digits-1) Then
s=Left(s, digits-1)
End If
i=Instr(Trim(s),".")
If i>2 Then
s=Left(s,i-1)+Mid(s, i+1)
s=Left(s,1)+"."+Mid(s, 2)
powten+=1
Elseif i=1 Then
s=Mid(s,2)
s=Left(s,1)+"."+Mid(s,2)
powten-=1
End If
s2=Trim(Str(Abs(powten)))
If (5-len(s2))>0 Then s2=String(5-len(s2),"0")+s2
If powten<0 Then
s2="-"+s2
Else
s2="+"+s2
End If
s=s+"e"+s2
If sign<>0 Then
s="-"+s
Else
s=" "+s
End If
Return s
End Function
'integer part of num
Function fpfix( Byref num As BigFloat_struct ) As BigFloat_struct
Dim As BigFloat_struct ip
Dim As Long ex, ex2, j, k, c
ex=(num.exponent And &h7FFFFFFFul)-BIAS
If ex<1 Then
Return ip
End If
If ex>=(NUM_BITS) Then
Return num
End If
ex2=ex\32
k=ex2
j=ex Mod 32
While ex2>0
ex2-=1
ip.mantissa[ex2]=num.mantissa[ex2]
Wend
ip.mantissa[k]=num.mantissa[k]
ip.mantissa[k]=shr32(ip.mantissa[k], (32-j)-3, c)
ip.mantissa[k]=shl32(ip.mantissa[k], (32-j)-3, c)
ip.exponent=ex+BIAS
ip.sign=num.sign
NORM_FAC1(ip)
Return ip
End Function
'fractional part of num
Function fpfrac( Byref num As BigFloat_struct ) As BigFloat_struct
Dim As BigFloat_struct ip, fp
ip=fpfix(num)
fp=fpsub(num, ip)
Return fp
End Function
'returns 1 if integer part is odd
Function fpfix_is_odd( Byref num As bigfloat_struct ) As Long
Dim As Long ex, j, k, c
Dim As Ulong m
ex=(num.exponent And &h7FFFFFFF)-BIAS
If ex<1 Then
Return 0
End If
If ex>=(NUM_BITS) Then
Print "error in function fpfix_is_odd"
Return -1
End If
k=ex\32
j=ex Mod 32
m=num.mantissa[k]
m=shr32(m, (32-j)-3, c)
Return m And 1
End Function
Function pi_brent_salamin (Byval digits As Ulong = NUM_DWORDS*9.63) As bigfloat_struct
If digits > NUM_DWORDS*9.63 Then digits = NUM_DWORDS*9.63
Dim As Long limit
Dim As bigfloat_struct c0, c1, c2, c05
Dim As bigfloat_struct a, b, sum
Dim As bigfloat_struct ak, bk, ck
Dim As bigfloat_struct ab, asq
Dim As bigfloat_struct pow2, tmp, eps
eps.exponent=-(NUM_BITS-10)+BIAS+1
eps.mantissa[0]=&h10000000ul
limit = (-digits) + BIAS + 1
c0=si2fp(0, digits): ak = c0: bk = c0: ab = c0: asq = c0
c1=si2fp(1, digits): a = c1: ck = c1: pow2 = c1
c2=si2fp(2, digits): b = c2
c05=ui2fp(1) : c05=fpdiv_si(c05, 2)
sum = c05
b=fpsqr(b, digits)
b=fpdiv(c1, b, digits)
While fpcmp(ck, eps) > 0
ak=fpadd(a, b)
ak=fpmul(c05, ak)
ab=fpmul(a, b)
bk=fpsqr(ab)
asq=fpmul(ak, ak)
ck=fpsub(asq, ab)
pow2=fpmul(pow2, c2)
tmp=fpmul(pow2, ck)
sum=fpsub(sum, tmp)
a = ak: b = bk
Wend
tmp=fpdiv(asq, sum)
tmp=fpmul(c2, tmp)
Return tmp
End Function
Last edited by srvaldez on Dec 16, 2021 10:11, edited 7 times in total.
Re: about bignum
paste this after the above post
Updated
Updated
Code: Select all
Dim Shared As BigFloat pi_bf, pi2_bf, pi_bf_half, pi_bf_quarter
pi_bf.BigNum = pi_brent_salamin()
pi2_bf.BigNum = fpadd(pi_bf.BigNum, pi_bf.BigNum)
pi_bf_half.BigNum = fpdiv_si( pi_bf.BigNum, 2)
pi_bf_quarter.BigNum = fpdiv_si( pi_bf.BigNum, 4)
Dim Shared As BigFloat tan_half_num(15), tan_half_den(14)
tan_half_num(0).BigNum=si2fp(992)
tan_half_num(1).BigNum=str2fp("161388480")
tan_half_num(2).BigNum=str2fp("7686610525440")
tan_half_num(3).BigNum=str2fp("167256984742848000")
tan_half_num(4).BigNum=str2fp("2000393537524462080000")
tan_half_num(5).BigNum=str2fp("14467646220791919547392000")
tan_half_num(6).BigNum=str2fp("66677300813465118447390720000")
tan_half_num(7).BigNum=str2fp("201117789910786072985458237440000")
tan_half_num(8).BigNum=str2fp("400342706504764747636935691468800000")
tan_half_num(9).BigNum=str2fp("521967288977995909272835160309760000000")
tan_half_num(10).BigNum=str2fp("435052278687602761865918494187323392000000")
tan_half_num(11).BigNum=str2fp("221463964316902607512694240578598338560000000")
tan_half_num(12).BigNum=str2fp("63663507608965602906315837691661069058048000000")
tan_half_num(13).BigNum=str2fp("8994510946805140308046160658488525397688320000000")
tan_half_num(14).BigNum=str2fp("470550277118574335327341015729793791741132800000000")
tan_half_num(15).BigNum=str2fp("3827142253897737927329040261268989506161213440000000")
tan_half_den(0).BigNum=si2fp(491040)
tan_half_den(1).BigNum=str2fp("39540177600")
tan_half_den(2).BigNum=str2fp("1232419887578880")
tan_half_den(3).BigNum=str2fp("19569067214913216000")
tan_half_den(4).BigNum=str2fp("180435497084706479616000")
tan_half_den(5).BigNum=str2fp("1036847979156754234229760000")
tan_half_den(6).BigNum=str2fp("3857758118493338995884748800000")
tan_half_den(7).BigNum=str2fp("9452536125806945430316537159680000")
tan_half_den(8).BigNum=str2fp("15257505370126034271052104685977600000")
tan_half_den(9).BigNum=str2fp("15972199042726674823748755905478656000000")
tan_half_den(10).BigNum=str2fp("10480804895655884717678945541785518080000000")
tan_half_den(11).BigNum=str2fp("4060172679143214471066061077274302873600000000")
tan_half_den(12).BigNum=str2fp("837419984702547545921539095790310985302016000000")
tan_half_den(13).BigNum=str2fp("75810877980214754024960496978688999780515840000000")
tan_half_den(14).BigNum=str2fp("1913571126948868963664520130634494753080606720000000")
Function fptan_half(Byref x As bigfloat_struct, Byval dwords As Ulong = NUM_DWORDS) As bigfloat_struct
If dwords > NUM_DWORDS Then dwords = NUM_DWORDS
Dim As bigfloat_struct accum, x_, xx, thn, thd
Dim As Long i, sign
xx=fpmul(x, x, dwords)
'the folowing is a rational polynomial for tan(x/2) derived from the continued fraction
' x^2
'tan(x/2) = 2 - -------------------------
' x^2
' 6 - ----------------------
' x^2
' 10 - ------------------
' x^2
' 14 - --------------
' x^2
' 18 - ---------
' x^2
' 22 - -----
' 26 - ...
' and so on
'
'the nice quality of this expansion is that you can calculate sin, cos and tan very easily
'
'if y = tan(x/2) then
'
' 2*x*y
'sin(x) = -------
' y^2+x^2
'
' y^2-x^2
'cos(x) = ---------
' y^2+x^2
'
' 2*x*y
'tan(x) = -------
' y^2-x^2
sign=-1
thn=tan_half_num(0).BigNum
For i=1 To 15
thn=fpmul(xx, thn, dwords)
If sign=1 Then
thn=fpadd(thn, tan_half_num(i).BigNum, dwords)
Else
thn=fpsub(thn, tan_half_num(i).BigNum, dwords)
End If
sign=-sign
Next
thd=fpsub(xx, tan_half_den(0).BigNum, dwords)
sign=1
For i=1 To 14
thd=fpmul(xx, thd, dwords)
If sign=1 Then
thd=fpadd(thd, tan_half_den(i).BigNum, dwords)
Else
thd=fpsub(thd, tan_half_den(i).BigNum, dwords)
End If
sign=-sign
Next
accum=fpdiv(thn,thd, dwords) 'tan(x/2)
Return accum
End Function
Function mp_sin(Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
If dwords > NUM_DWORDS Then dwords = NUM_DWORDS
Dim As BigFloat_struct ab, accum, factor, x_2, xx, p, tmp, tmp2
Dim As BigFloat_struct thn, thd
Dim As Long sign(0 To 3), c, limit, i
x_2=x
ab=x
ab.sign=0
If fpcmp(ab, pi2_bf.BigNum, dwords)=1 Then
'======== centralize ==============
'floor/ceil to centralize
tmp=fpdiv(x_2, pi2_bf.BigNum, dwords)
tmp2=tmp
tmp=fpfix(tmp) 'int part
tmp=fpsub(tmp2, tmp, dwords) 'frac part
tmp=fpmul(tmp, pi2_bf.BigNum, dwords)
x_2=tmp
End If
'==================================
limit=9.63*dwords
'limit = digits of precision, here's a polynomial fit to get an exstimate for limit
'works well for precision from 60 to 10000 digits
limit=1+Int(-0.45344993886092585968+0.022333002852398072433*limit+5.0461814408333079844e-7*limit*limit-4.2338453039804235772e-11*limit*limit*limit)
If limit<0 Then limit=0
factor=si2fp(5, dwords)
For i=1 To limit
factor=fpmul_si(factor, 5, dwords)
Next
'factor=factor^limit
x_2=fpdiv(x_2,factor, dwords) 'x_2=x_2/5^limit
accum=fptan_half(x_2, dwords)
xx=fpmul(x_2, x_2, dwords)
'now convert to sin(x)
tmp=fpmul_si(x_2, 2, dwords)
tmp=fpmul(tmp, accum, dwords)
tmp2=fpmul(accum, accum, dwords)
tmp2=fpadd(tmp2, xx, dwords)
accum=fpdiv(tmp, tmp2, dwords)
'multiply the result by 5^limit
For i=1 To limit+1
tmp=fpmul(accum, accum, dwords)
tmp2=fpmul(accum, tmp, dwords)
'sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
accum=fpmul_si(accum, 5, dwords)
accum=fpsub(accum, fpmul_si(tmp2, 20, dwords), dwords)
tmp=fpmul(tmp2, tmp, dwords)
accum=fpadd(accum, fpmul_si(tmp, 16, dwords), dwords)
Next i
Return accum
End Function
Function mp_cos (Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
If dwords > NUM_DWORDS Then dwords = NUM_DWORDS
Dim As BigFloat_struct ab, accum, factor, x_2, xx, p, tmp, tmp2
Dim As BigFloat_struct thn, thd
Dim As Long sign(0 To 3), c, limit, i
x_2=fpadd(pi_bf_half.BigNum, x, dwords)
ab=x_2 :ab.sign=0
If fpcmp(ab, pi2_bf.BigNum, dwords)=1 Then
'======== centralize ==============
'floor/ceil to centralize
tmp=fpdiv(x_2, pi2_bf.BigNum, dwords)
tmp2=tmp
tmp=fpfix(tmp) 'int part
tmp=fpsub(tmp2, tmp, dwords) 'frac part
tmp=fpmul(tmp, pi2_bf.BigNum, dwords)
x_2=tmp
End If
'==================================
limit=9.63*dwords
'limit = digits of precision, here's a polynomial fit to get an exstimate for limit
'works well for precision from 60 to 10000 digits
limit=1+Int(-0.45344993886092585968+0.022333002852398072433*limit+5.0461814408333079844e-7*limit*limit-4.2338453039804235772e-11*limit*limit*limit)
If limit<0 Then limit=0
factor=si2fp(5, dwords)
For i=1 To limit
factor=fpmul_si(factor, 5, dwords)
Next
'factor=factor^limit
x_2=fpdiv(x_2,factor, dwords) 'x_2=x_2/5^limit
accum=fptan_half(x_2, dwords)
xx=fpmul(x_2, x_2, dwords)
'now convert to sin(x)
tmp=fpmul_si(x_2, 2, dwords)
tmp=fpmul(tmp, accum, dwords)
tmp2=fpmul(accum, accum, dwords)
tmp2=fpadd(tmp2, xx, dwords)
accum=fpdiv(tmp, tmp2, dwords)
'multiply the result by 5^limit
For i=1 To limit+1
tmp=fpmul(accum, accum, dwords)
tmp2=fpmul(accum, tmp, dwords)
'sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
accum=fpmul_si(accum, 5, dwords)
accum=fpsub(accum, fpmul_si(tmp2, 20, dwords), dwords)
tmp=fpmul(tmp2, tmp, dwords)
accum=fpadd(accum, fpmul_si(tmp, 16, dwords), dwords)
Next i
Return accum
End Function
Function mp_tan(Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
If dwords > NUM_DWORDS Then dwords = NUM_DWORDS
Dim As BigFloat_struct ab, accum, factor, x_2, xx, p, pi, pi2, circ, tmp, tmp2
Dim As BigFloat_struct thn, thd
Dim As Long c, limit, i, sign
pi=pi_bf.BigNum
pi2=fpdiv_si(pi_bf.BigNum, 2, dwords)
x_2=x
'calculate the sign of tan(x)
'if integer part of x/(pi/2) is odd then the sign is negative
'unless x is negative then the sign is positive
sign=fpfix_is_odd(fpdiv(x,pi2, dwords))
If sign=1 Then
sign=&h8000
End If
circ=fpmul_si(pi, 2, dwords)
ab=x :ab.sign=0
If fpcmp(ab, circ, dwords)=1 Then
'======== centralize ==============
'floor/ceil to centralize
pi=fpadd(pi, pi, dwords) 'got 2*pi
tmp=fpdiv(x_2, pi, dwords)
tmp2=tmp
tmp=fpfix(tmp) 'int part
tmp=fpsub(tmp2, tmp, dwords) 'frac part
tmp=fpmul(tmp, pi, dwords)
x_2=tmp
End If
'==================================
limit=dwords*9.63
'lm = dwords of precision, here's a polynomial fit to get an exstimate for limit
'works well for precision from 60 to 10000 dwords
limit=1+Int(-0.45344993886092585968+0.022333002852398072433*limit+5.0461814408333079844e-7*limit*limit-4.2338453039804235772e-11*limit*limit*limit)
If limit<0 Then limit=0
factor=si2fp(5, dwords)
For i=1 To limit
factor=fpmul_si(factor, 5, dwords)
Next
'factor=factor^limit
x_2=fpdiv(x_2,factor, dwords) 'x_2=x_2/5^limit
accum=fptan_half(x_2, dwords)
xx=fpmul(x_2, x_2, dwords)
'now convert to sin(x)
tmp=fpmul_si(x_2, 2, dwords)
tmp=fpmul(tmp, accum, dwords)
tmp2=fpmul(accum, accum, dwords)
tmp2=fpadd(tmp2, xx, dwords)
accum=fpdiv(tmp, tmp2, dwords)
'multiply the result by 5^limit
For i=1 To limit+1
tmp=fpmul(accum, accum, dwords)
tmp2=fpmul(accum, tmp, dwords)
'sin(5*x) = 5 * sin(x) - 20 * sin(x)^3 + 16 * sin(x)^5
accum=fpmul_si(accum, 5, dwords)
accum=fpsub(accum, fpmul_si(tmp2, 20, dwords), dwords)
tmp=fpmul(tmp2, tmp, dwords)
accum=fpadd(accum, fpmul_si(tmp, 16, dwords), dwords)
Next i
tmp=fpmul(accum, accum, dwords)
tmp=fpsub(si2fp(1, dwords), tmp, dwords)
tmp=fpsqr(tmp, dwords)
tmp=fpdiv(accum, tmp, dwords)
tmp.sign=sign Xor x.sign
Return tmp
End Function
Function fpatn (Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
If dwords > NUM_DWORDS Then dwords = NUM_DWORDS
Dim As Long sign(3): sign(3) = 1
Dim As Long z, c = 1
Dim As BigFloat_struct XX, Term, Accum, strC, x_2, mt, mt2, p
Dim As BigFloat_struct Bignum, one, Bignum2, factor
Bignum2 = x
Bignum2.sign = 0
one = si2fp(1, dwords)
If fpcmp(Bignum2, one, dwords) = 0 Then
Bignum = pi_bf_quarter.BigNum
Bignum.sign = x.sign
Return Bignum
End If
Bignum2.sign = x.sign
Dim As Long limit = 16
factor = si2fp(2 Shl (limit - 1), dwords)
For z = 1 To limit
Bignum = fpmul(Bignum2, Bignum2, dwords)
Bignum = fpadd(Bignum, one, dwords)
Bignum = fpsqr(Bignum, dwords)
Bignum = fpadd(Bignum, one, dwords)
Bignum = fpdiv(Bignum2, Bignum, dwords)
Bignum2 = Bignum
Next z
mt = Bignum
x_2 = Bignum
p = Bignum
XX = fpmul(x_2, x_2, dwords)
Do
c = c + 2
mt2 = mt
p = fpmul(p, XX, dwords)
Term = fpdiv_si(p, c, dwords)
If sign(c And 3) Then
Accum = fpsub(mt, Term, dwords)
Else
Accum = fpadd(mt, Term, dwords)
End If
Swap mt, Accum
Loop Until fpcmp(mt, mt2, dwords) = 0
Return fpmul(factor, mt, dwords)
End Function
Function fpasin (Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
If dwords > NUM_DWORDS Then dwords = NUM_DWORDS
Dim As Double num
Dim As BigFloat_struct one, T, B, term1, minusone
' ARCSIN = ATN(x / SQR(-x * x + 1))
'============= ARCSIN GUARD =========
T = x
one =si2fp(1, dwords)
minusone = si2fp(-1, dwords)
If fpcmp(T, one)>0 Then
Return one
End If
If fpcmp(T, minusone)<0 Then
Return minusone
End If
B = fpmul(x, x, dwords) 'x*x
'for 1 and -1
If fpcmp(B, one, dwords) = 0 Then
Dim As BigFloat_struct two, Atn1
two = si2fp(2, dwords)
'atn1 = fpatn(one, dwords)
If fpcmp(x, minusone, dwords) = 0 Then
two = pi_bf_half.BigNum 'fpmul(two, atn1, dwords)
two.sign =&h8000
Return two 'fpmul(two, minusone, dwords)
Else
Return pi_bf_half.BigNum 'fpmul(two, atn1, dwords)
End If
End If
B = fpsub(one, B, dwords) '1-x*x
B = fpsqr(B, dwords) 'sqr(1-x*x)
term1 = fpdiv(T, B, dwords)
Return fpatn(term1, dwords)
End Function
Function fpacos (Byref x As BigFloat_struct, Byval dwords As Ulong = NUM_DWORDS) As BigFloat_struct
If dwords > NUM_DWORDS Then dwords = NUM_DWORDS
Dim As BigFloat_struct one, minusone, two, Atn1, tail, T, B, term1, atnterm1, zero, ep
ep.exponent=-(NUM_BITS-20)+BIAS+1
ep.mantissa[0]=&h10000000ul
'ARCCOS = ATN(-x / SQR(-x * x + 1)) + 2 * ATN(1)
'========================
one = si2fp(1, dwords)
minusone = si2fp(-1, dwords)
T=fpsub(x,one)
T.sign=0
If fpcmp(T, ep)<=0 Then
Return zero
End If
If fpcmp(x, one)>=0 Then
Return zero
End If
If fpcmp(x, minusone)<0 Then
Return pi_bf.BigNum
End If
two = si2fp(2, dwords)
Atn1 = pi_bf_quarter.BigNum
tail = pi_bf_half.BigNum '2*atn(1)
T = fpmul(minusone, x, dwords) '-x
B = fpmul(x, x, dwords) 'x*x
If fpcmp(B, one, dwords) = 0 Then
'for 1 and -1
If fpcmp(x, minusone, dwords) = 0 Then
Return pi_bf.BigNum
Else
Return si2fp(0, dwords)
End If
End If
B = fpsub(one, B, dwords) '1-x*x
B = fpsqr(B, dwords) 'sqr(1-x*x)
term1 = fpdiv(T, B, dwords)
atnterm1 = fpatn(term1, dwords)
Return fpadd(atnterm1, tail, dwords)
End Function
Dim As BigFloat x, y, z, t, eps, one, minusone
Dim As Double tm
Dim As Longint i, factor, limit=NUM_DWORDS*9.63
Dim As Long j, ex
Dim As String s
eps.BigNum.exponent=-(NUM_BITS-16)+BIAS+1
eps.BigNum.mantissa[0]=&h10000000ul
Print "current precision is ";NUMBER_OF_DIGITS;" digits"
tm=Timer
x=1
Print "Sin(";x.toLong;") = ";(Sin(x)).toString(20)
Print "Cos(";x.toLong;") = ";(Cos(x)).toString(20)
Print "Tan(";x.toLong;") = ";(Tan(x)).toString(20)
tm=timer-tm
Print "elapsed time ";tm;" seconds"
Sleep
Last edited by srvaldez on Dec 13, 2021 21:43, edited 4 times in total.
Re: about bignum
Thanks srvaldez.
no problems here.
I am inspired finish my own bignum stuff, I have +, -, *, and \ and loads of bits and pieces.
I use string[index], It'll take a bit of time to put it all together, I'll start off with sin/cos/tan like yourself.
no problems here.
I am inspired finish my own bignum stuff, I have +, -, *, and \ and loads of bits and pieces.
I use string[index], It'll take a bit of time to put it all together, I'll start off with sin/cos/tan like yourself.
Re: about bignum
Hi dodicat, I am looking forward to it :-)
Re: about bignum
I did some timings between bigfloat and my decfloat math, bigfloat is over 2 times faster
NUMBER_OF_DIGITS=1040
bigfloat sin 2.80 sec., decfloat 6.26, factor of 2.23
bigfloat cos 2.82 sec., decfloat 6.34, factor of 2.24
bigfloat tan 3.28 sec., decfloat 6.67, factor of 2.03
the inverse
same for acos and atn
bigfloat asin(sin) 17.064 sec., decfloat 23.452, factor of 1.374
bigfloat acos(cos) 16.616 sec., decfloat 24.82, factor of 1.674
bigfloat atan(tan) 16.702 sec., decfloat 23.565, factor of 1.41
NUMBER_OF_DIGITS=1040
Code: Select all
tm=timer
for x=-1000 to 1000
y.BigNum=mp_sin(x.BigNum)
z+=y
next
tm=timer-tm
Print z
Print tm
bigfloat cos 2.82 sec., decfloat 6.34, factor of 2.24
bigfloat tan 3.28 sec., decfloat 6.67, factor of 2.03
the inverse
Code: Select all
tm=timer
for x=-1000 to 1000
t.BigNum=mp_sin(x.BigNum)
y.BigNum=fpasin(t.BigNum)
z+=y
next
tm=timer-tm
bigfloat asin(sin) 17.064 sec., decfloat 23.452, factor of 1.374
bigfloat acos(cos) 16.616 sec., decfloat 24.82, factor of 1.674
bigfloat atan(tan) 16.702 sec., decfloat 23.565, factor of 1.41
Re: about bignum
dodicat, I look forward to your implementation, meanwhile looking at my trig implementation I see that I use 5*x reduction, it may be faster to do 4*x reduction because in binary floating-point a division of a power of 2 is simply decrementing the exponent by the appropriate amount
[edit]
corrected a couple of bugs, please update your copy
[edit]
corrected a couple of bugs, please update your copy