Code: Select all
' 1 December 2016 - original by dodicat
' 11 December 2016 - frisian, added mpf_clear/mpz_clear statements
' changed some code by using GMP code, cleaned up some init/init_set stuff
' reduced the number of big integer/float variables needed
' 16 December 2016 - frisian added error routine
' added trim for divide and sqrroot
' added cleanup for greater and equals
#Include Once "gmp.bi"
Type mpf_t As __mpf_struct
' floating point routine's
'approx list
'equals,greater,less,absolute
'sin,cos,tan,logtaylor,log,exp,power,atn,acos,asin,greater,equals,less
'absolute,Pi_ui
'
' Integer routine's : _mod, _div
Dim Shared As ULongInt PRECISION
'========= Just in case you forget to set_precision ==========
precision=60
mpf_set_default_prec(PRECISION * 4)
'========================================================
Sub big_num_error(n As UInteger)
Select Case n
Case 1
Print "Number is Floating Point, expected Integer"
Case Else
Print "Unknow Error"
End Select
Print
Print "Program stops"
Sleep
End
End Sub
Sub set_precision(n As UInteger)
PRECISION = n
mpf_set_default_prec(PRECISION * 4)
End Sub
Function equals Overload(a As mpf_t, b As mpf_t) As Integer
If Mpf_cmp(@a, @b) = 0 Then Return -1
Return 0
End Function
Function greater Overload(a As mpf_t, b As mpf_t) As Integer 'a>b
If mpf_cmp(@a, @b) > 0 Then Return -1
Return 0
End Function
Function less Overload(a As mpf_t, b As mpf_t) As Integer 'a<b
If mpf_cmp(@a, @b) < 0 Then Return -1
Return 0
End Function
Function Absolute Overload(a As mpf_t) As mpf_t
Dim As mpf_t Ab
mpf_init(@Ab)
mpf_abs(@Ab, @a)
Return Ab
End Function
Function Pi_ui Overload(places As UInteger) As mpf_t
' Dim As __mpf_struct a,b,t,p,aa,bb,tt,pp,pi
Dim As __mpf_struct a, b, t, aa, bb, tt, pi
mpf_init2(@a, 4*places)
mpf_init2(@b, 4*places)
mpf_init2(@t, 4*places)
' mpf_init2( @p, 4*places)
Dim As UInteger p
mpf_init2(@aa,4*places)
mpf_init2(@bb,4*places)
mpf_init2(@tt,4*places)
' mpf_init2( @pp,4*places)
mpf_init2(@pi,4*places)
mpf_set_ui(@a, 1)
mpf_set_ui(@b, 2) : mpf_sqrt(@b, @b)
mpf_set_str(@t,".25",10)
' mpf_set_str(@p,"1",10)
mpf_ui_div(@b,1,@b)
Do
mpf_add(@aa, @a, @b)
' mpf_div_ui(@aa, @aa, 2)
mpf_div_2exp(@aa, @aa, 1)
mpf_mul(@bb, @a, @b)
mpf_sqrt(@bb, @bb)
mpf_sub(@tt, @a, @aa)
mpf_mul(@tt,@tt,@tt)
' mpf_mul(@tt, @tt,@p)
mpf_mul_2exp(@tt, @tt, p)
p += 1
mpf_sub(@tt, @t, @tt)
' mpf_mul_ui(@pp, @p, 2)
mpf_swap(@a, @aa)
mpf_swap(@b, @bb)
mpf_swap(@t, @tt)
' mpf_swap(@p, @pp)
Loop Until Mpf_cmp(@a, @aa) = 0
mpf_add(@pi, @a, @b)
mpf_mul(@pi, @pi, @pi)
' mpf_div_ui(@pi, @pi, 4)
mpf_div_2exp(@pi, @pi, 2)
mpf_div(@pi, @pi, @t)
' remove big int's from memory
mpf_clear(@a) : mpf_clear(@aa)
mpf_clear(@b) : mpf_clear(@bb)
mpf_clear(@t) : mpf_clear(@tt)
Return pi
End Function
Function _sin Overload(x As mpf_t) As mpf_t
Dim As mpf_t XX, Term, Accum, _x, temp2, fac, pi2
mpf_init_set(@_x, @x)
mpf_init(@pi2)
mpf_init(@temp2)
If mpf_cmp_d(@x, 6.283185) >= 0 OrElse mpf_cmp_ui(@x, 0) < 0 Then
'======== CENTRALIZE ==============
'floor/ceil to centralize
pi2 = Pi_ui(precision)
mpf_mul_2exp(@pi2, @pi2, 1) ' pi2 = pi * 2
mpf_div(@temp2, @_x, @pi2) ' temp2 = _x / pi2
mpf_floor(@temp2, @temp2) ' temp2 = floor(temp2) 'rounds temp2 down towards minus infinity
mpf_mul(@temp2, @temp2, @pi2) ' temp2 = temp2 * p2
mpf_sub(@_x, @_x, @temp2) ' _x = _x - temp2 (_x = _x mod 2*pi)
End If
'==================================
' based on Richard methode for sin/cos in squares
Dim As Integer c = 1
mpf_init(@XX)
mpf_init(@Term)
mpf_init_set(@accum, @_x)
mpf_init_set_ui(@fac, 1)
mpf_mul(@XX, @_x, @_x)
mpf_neg(@XX, @XX) ' make XX negative
Do
c = c + 2
Mpf_swap(@temp2, @Accum)
mpf_mul_ui(@fac, @fac, c) ' avoid an overflow when c > 65536
mpf_mul_ui(@fac, @fac, (c -1)) ' split in two separate multiply's
mpf_mul(@_x , @_x, @XX) ' _x alternates between negative an positive
mpf_div(@Term, @_x, @fac)
mpf_add(@Accum, @temp2, @term) ' temp2 and accum are swapped, temp2 holds the previous content of accum
Loop Until Mpf_cmp(@Accum, @temp2) = 0
' clean up
mpf_clear(@XX) : mpf_clear(@Term) : mpf_clear(@fac)
mpf_clear(@_x) : mpf_clear(@temp2) : mpf_clear(@pi2)
Return accum
End Function
Function _cos Overload(x As mpf_t) As mpf_t
Dim As mpf_t XX, Term, Accum, _x, p, temp2, fac, pi2
mpf_init_set(@_x,@x)
mpf_init(@pi2)
mpf_init(@temp2)
If mpf_cmp_d(@x, 6.283185) >= 0 OrElse mpf_cmp_ui(@x, 0) < 0 Then
'======== CENTRALIZE ==============
'floor/ceil to centralize
pi2 = Pi_ui(precision)
mpf_mul_2exp(@pi2, @pi2, 1) ' pi2 = pi * 2
mpf_div(@temp2, @_x, @pi2) ' temp2 = _x / pi2
mpf_floor(@temp2, @temp2) ' temp2 = floor(temp2) 'rounds temp2 down towards minus infinity
mpf_mul(@temp2, @temp2, @pi2) ' temp2 = temp2 * p2
mpf_sub(@_x, @_x, @temp2) ' _x = _x - temp2 (_x = _x mod 2*pi)
End If
'==================================
' based on Richard methode for sin/cos in squares
Dim As Integer c
mpf_init(@XX)
mpf_init(@Term)
mpf_init_set_ui(@Accum, 1)
mpf_init_set_ui(@fac, 1)
mpf_init_set_ui(@p, 1)
mpf_mul(@XX, @_x, @_x)
mpf_neg(@XX, @XX)
Do
c += 2
Mpf_swap(@temp2, @accum)
mpf_mul_ui(@fac, @fac, c)
mpf_mul_ui(@fac, @fac, (c-1))
mpf_mul(@p, @p, @XX)
mpf_div(@term, @p, @fac)
mpf_add(@accum, @temp2, @term)
Loop Until Mpf_cmp(@accum, @temp2) = 0
' clean up
mpf_clear(@XX) : mpf_clear(@Term) : mpf_clear(@fac)
mpf_clear(@_x) : mpf_clear(@temp2) : mpf_clear(@pi2)
mpf_clear(@p)
Return accum
End Function
Function _tan Overload(x As mpf_t) As mpf_t
Dim As mpf_t s, c, ans
mpf_init(@ans)
s = _sin(x)
c = _cos(x)
mpf_div(@ans, @s, @c)
mpf_clear(@c) : mpf_clear(@s)
Return ans
End Function
Function _logTaylor(x As mpf_t) As mpf_t
'taylor series
'====================Log Guard==================
If mpf_cmp_ui(@x, 0) <= 0 Then Exit Function ' exit if x = 0 or x = negative
'===============================================
Dim As Integer invflag
Dim As Mpf_t Q, tmp, _x, accum, term, XX
mpf_init(@XX)
mpf_init(@Q)
mpf_init(@tmp)
mpf_init(@accum)
mpf_init(@term)
mpf_init_set(@_x, @x)
If mpf_cmp_ui(@_x, 1) < 0 Then
invflag = 1
mpf_ui_div(@_x, 1, @_x)
End If
mpf_sub_ui(@tmp, @_x, 1)
mpf_add_ui(@Q, @_x, 1) ' q = b
mpf_div(@accum,@tmp, @Q)
Mpf_set(@Q , @accum)
Mpf_mul(@XX, @Q, @Q)
Dim As Integer c=1
Do
c += 2
Mpf_swap(@tmp,@accum)
mpf_mul(@Q, @Q, @XX)
mpf_div_ui(@term, @Q, c)
mpf_add(@Accum, @tmp, @Term)
Loop Until Mpf_cmp(@tmp, @accum) = 0
mpf_mul_2exp(@accum, @accum, 1)
If invflag Then
mpf_neg(@accum, @accum)
End If
mpf_clear(@_x) : mpf_clear(@tmp) : mpf_clear(@Q)
mpf_clear(@XX) : mpf_clear(@term)
Return accum
End Function
Function _log Overload(x As mpf_t) As mpf_t
/'
'====================Log Guard==================
If Mpf_cmp_ui(@x, 0) <= 0 Then Exit Function ' exit if x = 0 or x is negative
'===============================================
Dim As mpf_t approx, ans, logx ',factor
Mpf_init_set(@approx, @x)
mpf_init(@ans)
mpf_sqrt(@approx, @approx) ' 1
mpf_sqrt(@approx, @approx) ' 2
mpf_sqrt(@approx, @approx) ' 3
logx = _logTaylor(approx)
mpf_mul_2exp(@ans, @logx, 3)
' clean up
mpf_clear(@approx) : mpf_clear(@logx)
'/
'====================Log Guard==================
If Mpf_cmp_ui(@x, 0) <= 0 Then Exit Function ' exit if x = 0 or x is negative
'===============================================
Dim As mpf_t approx, ans, logx ',factor
Mpf_init_set(@approx, @x)
mpf_init(@ans)
Dim As Integer c
While mpf_cmp_d(@approx, 1.25) > 0
mpf_sqrt(@approx, @approx)
c += 1
Wend
logx = _logTaylor(approx)
mpf_mul_2exp(@ans, @logx, c)
' clean up
mpf_clear(@approx) : mpf_clear(@logx)
Return ans
End Function
Function _exp Overload(x As mpf_t) As mpf_t
'taylor series
Dim As mpf_t fac, temp2, accum, p, term
mpf_init(@temp2)
mpf_init(@term)
mpf_init_set_ui(@fac, 1)
mpf_init_set_ui(@p, 1)
mpf_init_set_ui(@accum, 1)
Dim As Integer c
Do
c += 1
Mpf_swap(@temp2, @accum)
mpf_mul_ui(@fac, @fac, c)
mpf_mul(@p, @p, @x)
mpf_div(@term, @p, @fac)
mpf_add(@Accum, @temp2, @Term)
Loop Until Mpf_cmp(@accum, @temp2) = 0
' clean up
mpf_clear(@temp2) : mpf_clear(@fac)
mpf_clear(@term) : mpf_clear(@p)
Return accum
End Function
Function power Overload(a As mpf_t,p As mpf_t) As mpf_t
'a^p= exp(p*log(a))
'====================Power Guard==================
If Mpf_cmp_ui(@a, 0) <= 0 Then Exit Function ' exit if x = 0 or x is negative
'=================================================
Dim As mpf_t loga, product, ans
mpf_init(@product)
loga = _log(a)
mpf_mul(@product, @p, @loga)
ans = _exp(product)
'clean up
mpf_clear(@loga) : mpf_clear(@product)
Return ans
End Function
Function _Atn Overload(x As mpf_t) As mpf_t
#Macro ArctanTaylor(decnum)
mpf_init(@XX)
mpf_init(@Term)
mpf_init(@Accum)
mpf_init_set(@mt, @decnum)
mpf_init_set(@p, @decnum)
mpf_mul(@XX, @mt, @mt)
mpf_neg(@XX, @XX)
Do
c += 2
mpf_mul(@p, @p, @XX)
mpf_div_ui(@Term, @p, c)
mpf_add(@Accum, @mt, @Term)
Mpf_swap(@mt, @Accum)
Loop Until Mpf_cmp(@mt, @accum) = 0
#EndMacro
Dim As UInteger c = 1
Dim As mpf_t XX, Term, Accum, mt, p, _temp, _temp2
mpf_init(@_temp)
mpf_init_set(@_temp2, @x)
Dim As Integer limit = 16
For z As Integer = 1 To limit
mpf_mul( @_temp, @_temp2, @_temp2)
mpf_add_ui(@_temp ,@_temp, 1)
mpf_sqrt( @_temp, @_temp)
mpf_add_ui(@_temp, @_temp, 1)
mpf_div( @_temp, @_temp2, @_temp)
Mpf_swap( @_temp, @_temp2)
Next z
ArctanTaylor(_temp)
mpf_mul_2exp(@mt, @mt, (limit -1))
mpf_clear(@_temp) : mpf_clear(@_temp2) : mpf_clear(@XX)
mpf_clear(@accum) : mpf_clear(@term) : mpf_clear(@p)
Return mt
End Function
Function _Acos Overload(x As mpf_t) As mpf_t
'ARCCOS = ATN(-x / SQR(-x * x + 1)) + 2 * ATN(1)
'============= ARCCOS GUARD =========
Dim As Mpf_t B : Mpf_init(@B)
Mpf_mul(@B, @x, @x) 'x*x
If Mpf_cmp_ui(@B, 1) > 0 Then
mpf_clear(@B) ' need to clean up B
Exit Function
End If
'====================================
Dim As Mpf_t atn1, term1, ans
Mpf_init_set_ui(@term1, 1)
Mpf_init(@ans)
atn1=_Atn(term1)
If Mpf_cmp_ui(@B, 1) = 0 Then
'for 1 and -1
If Mpf_cmp_si(@x, -1) = 0 Then
Mpf_mul_2exp(@ans, @atn1, 2)
End If
mpf_clear(@B) : mpf_clear(@atn1) : mpf_clear(@term1)
' ans = pi or 0
Return ans
End If
Dim As Mpf_t tail, T, atnterm1
mpf_init(@tail)
mpf_init(@T)
mpf_mul_2exp(@tail, @atn1, 1) ' 2*atn(1)
mpf_neg(@T, @x) ' -x
'mpf_mul(@B,@x,@x) ' x*x ' done at the begin
mpf_ui_sub(@B, 1, @B) ' 1 - x*x
mpf_sqrt(@B, @B) ' sqr(1 - x*x)
mpf_div(@term1, @T, @B)
atnterm1 = _Atn(term1)
mpf_add(@ans, @atnterm1, @tail)
'clean up
mpf_clear(@B) : mpf_clear(@atn1) : mpf_clear(@term1)
mpf_clear(@T) : mpf_clear(@tail) : mpf_clear(@atnterm1)
Return ans
End Function
Function _Asin Overload(x As mpf_t) As mpf_t
' ARCSIN = ATN(x / SQR(-x * x + 1))
'============= ARCSIN GUARD =========
Dim As Mpf_t B : Mpf_init(@B)
Mpf_mul(@B, @x, @x) 'x*x
If Mpf_cmp_ui(@B, 1) > 0 Then
mpf_clear(@B)
Exit Function
End If
'====================================
Dim As mpf_t term1
mpf_init_set_ui(@term1, 1)
'for 1 and -1
If Mpf_cmp_ui(@B, 1) = 0 Then
Dim As mpf_t atn1
atn1 = _Atn(term1)
mpf_mul_2exp(@atn1, @atn1, 1)
If mpf_cmp_si(@x, -1) = 0 Then
mpf_neg(@atn1, @atn1)
End If
' clean up
mpf_clear(@B) : mpf_clear(@term1)
Return atn1
End If
Dim As Mpf_t T, atnterm1
Mpf_init_set(@T, @x)
mpf_ui_sub(@B, 1, @B) '1 - x*x
mpf_sqrt(@B, @B) 'sqr(1 - x*x)
mpf_div(@term1, @T, @B)
atnterm1 = _Atn(term1)
' clean up
mpf_clear(@B) : mpf_clear(@T) : mpf_clear(@term1)
Return atnterm1
End Function
'===========================================================================
'======================= OVERLOADS FOR STRINGS =============================
Dim Shared As ZString * 100000000 outtext
Function Pi_ui Overload(places As Integer) As String
Dim As Mpf_t ans
Var pl=CUInt(places)
ans=Pi_ui(pl)
gmp_sprintf(@outtext, "%." & pl & "Ff", @ans )
mpf_clear(@ans)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function _sin Overload(x As String) As String
Dim As Mpf_t _x, ans
mpf_init_set_str(@_x, x, 10)
ans = _sin(_x)
gmp_sprintf(@outtext, "%." & precision & "Ff", @ans )
mpf_clear(@_x) : mpf_clear(@ans)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function _cos Overload(x As String) As String
Dim As Mpf_t _x, ans
mpf_init_set_str(@_x,x,10)
ans = _cos(_x)
gmp_sprintf(@outtext, "%." & precision & "Ff", @ans )
mpf_clear(@_x) : mpf_clear(@ans)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function _tan Overload(x As String) As String
Dim As Mpf_t _x, ans
mpf_init_set_str(@_x,x,10)
ans = _tan(_x)
gmp_sprintf( @outtext, "%." & precision & "Ff", @ans )
mpf_clear(@_x) : mpf_clear(@ans)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function _log Overload(x As String) As String
Dim As Mpf_t _x, ans
mpf_init_set_str(@_x,x,10)
ans = _log(_x)
gmp_sprintf( @outtext,"%." & precision & "Ff",@ans )
mpf_clear(@_x) : mpf_clear(@ans)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function _exp Overload(x As String) As String
Dim As Mpf_t _x, ans
mpf_init_set_str(@_x, x, 10)
ans = _exp(_x)
gmp_sprintf( @outtext,"%." & precision & "Ff",@ans )
mpf_clear(@_x) : mpf_clear(@ans)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function power Overload(a As String,p As String) As String
Dim As Mpf_t _x, ans, pow
mpf_init_set_str(@_x, a, 10)
mpf_init_set_str(@pow, p, 10)
ans = power(_x, pow)
gmp_sprintf( @outtext,"%." & precision & "Ff", @ans )
mpf_clear(@_x) : mpf_clear(@ans) : mpf_clear(@pow)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function _Atn Overload(x As String) As String
Dim As Mpf_t _x, ans
mpf_init_set_str(@_x, x, 10)
ans = _Atn(_x)
gmp_sprintf( @outtext,"%." & precision & "Ff", @ans )
mpf_clear(@_x) : mpf_clear(@ans)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function _Acos Overload(x As String) As String
Dim As Mpf_t _x, ans
mpf_init_set_str(@_x, x, 10)
ans = _Acos(_x)
gmp_sprintf( @outtext,"%." & precision & "Ff", @ans )
mpf_clear(@_x) : mpf_clear(@ans)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function _Asin Overload(x As String) As String
Dim As Mpf_t _x, ans
mpf_init_set_str(@_x, x, 10)
ans = _Asin(_x)
gmp_sprintf( @outtext,"%." & precision & "Ff", @ans )
mpf_clear(@_x) : mpf_clear(@ans)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function factorial(n As UInteger) As String 'Automatic precision
Dim As __mpz_struct Intanswer
mpz_init( @Intanswer)
mpz_fac_ui(@Intanswer,n)
gmp_sprintf( @outtext,"%Zi", @Intanswer )
mpz_clear(@Intanswer)
Return Trim(outtext)
End Function
Function _mod(n1 As String,n2 As String) As String
If InStr(n1, ".") <> 0 OrElse InStr(n2, ".") <> 0 Then
big_num_error(1)
End If
Dim As __mpz_struct answer, mn1, mn2
mpz_init(@answer)
mpz_init_set_str(@mn1, n1, 10)
mpz_init_set_str(@mn2, n2, 10)
mpz_mod(@answer, @mn1, @mn2)
gmp_sprintf( @outtext,"%Zi", @answer )
mpz_clear(@answer) : mpz_clear(@mn1) : mpz_clear(@mn2)
Return Trim(outtext)
End Function
Function _div(n1 As String,n2 As String) As String
If InStr(n1, ".") <> 0 OrElse InStr(n2, ".") <> 0 Then
big_num_error(1)
End If
Dim As __mpz_struct answer, mn1, mn2
mpz_init(@answer)
mpz_init_set_str(@mn1, n1, 10)
mpz_init_set_str(@mn2, n2, 10)
mpz_div(@answer, @mn1, @mn2)
gmp_sprintf( @outtext,"%Zi", @answer )
mpz_clear(@answer) : mpz_clear(@mn1) : mpz_clear(@mn2)
Return Trim(outtext)
End Function
Function Sqrroot(number As String,decimals As UInteger=PRECISION) As String'precision parameter
If InStr(number,"-") Then Exit Function
Dim As __mpf_struct num, FloatAnswer
Dim As Integer LN = Len(number)
mpf_init2(@num, 4*Ln) : mpf_init2(@FloatAnswer, 4*Ln+4*decimals)
mpf_set_str(@num, number, 10)
mpf_sqrt( @FloatAnswer, @num)
gmp_sprintf( @outtext,"%." & Str(decimals) & "Ff",@FloatAnswer )
mpf_clear(@num) : mpf_clear(@FloatAnswer)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function mult(number1 As String,number2 As String) As String'Automatic precision
Dim As Integer Ln1 = Len(number1), Ln2 = Len(number2)
Dim As __mpf_struct num1,num2,FloatAnswer
mpf_init2(@num1, 4*(Ln1+Ln2+1) )
mpf_init2(@num2, 4*(Ln1+Ln2+1) )
mpf_init2(@Floatanswer, 4*(Ln1+Ln2+1))
Ln1=InStr(1, number1, ".") : Ln2 = InStr(1, number2, ".")
Var decimals = Len(Mid(number1, Ln1+1))+Len(Mid(number2, Ln2+1))+1
mpf_set_str(@num1, number1, 10)
mpf_set_str(@num2, number2, 10)
mpf_mul(@Floatanswer, @num1, @num2)
gmp_sprintf( @outtext,"%." & Str(decimals) & "Ff",@FloatAnswer )
mpf_clear(@num1) : mpf_clear(@num2) : mpf_clear(@FloatAnswer)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
'precision parameter
Function divide(number1 As String,number2 As String,decimals As UInteger=PRECISION) As String
Dim As Integer Ln1=Len(number1),Ln2=Len(number2),Ln
If Ln1>=Ln2 Then Ln=Ln1 Else Ln=Ln2
Dim As __mpf_struct num1,num2,FloatAnswer
mpf_init2(@num1, 4*(Ln+1) )
mpf_init2(@num2, 4*(Ln+1) )
mpf_init2(@Floatanswer, 4*(Ln+1)+4*decimals)
mpf_set_str(@num1, number1, 10)
mpf_set_str(@num2, number2, 10)
mpf_div(@Floatanswer, @num1, @num2)
gmp_sprintf( @outtext,"%." & Str(decimals) & "Ff",@FloatAnswer)
mpf_clear(@num1) : mpf_clear(@num2) : mpf_clear(@FloatAnswer)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function power Overload(number As String, n As UInteger) As String'automatic precision
#Define dp 3321921
Dim As __mpf_struct _number, FloatAnswer
Dim As ULongInt ln=Len(number)*(n)*4
If ln>dp Then ln=dp
mpf_init2(@FloatAnswer, ln)
mpf_init2(@_number, ln) 'or 4*len(number)
mpf_set_str(@_number, number, 10)
mpf_pow_ui(@Floatanswer, @_number, n)
gmp_sprintf( @outtext,"%." & Str(n) & "Ff",@FloatAnswer )
mpf_clear(@_number) : mpf_clear(@FloatAnswer)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function plus(number1 As String,number2 As String) As String'automatic precision
Dim As Integer Ln1=Len(number1),Ln2=Len(number2),decimals,Ln
If Ln1>=Ln2 Then Ln=Ln1 Else Ln=Ln2
Ln=ln+1
Dim As __mpf_struct num1,num2,FloatAnswer
mpf_init2(@num1, 4*(Ln1+1) )
mpf_init2(@num2, 4*(Ln2+1) )
mpf_init2(@Floatanswer, 4*(Ln))
mpf_set_str(@num1,number1,10)
mpf_set_str(@num2,number2,10)
Ln1=InStr(1,number1,"."):Ln2=InStr(1,number2,".")
If Ln1 Or Ln2 Then
decimals=Len(Mid(number1,Ln1+1))+Len(Mid(number2,Ln2+1))+1
End If
mpf_add(@Floatanswer, @num1, @num2)
gmp_sprintf( @outtext,"%." & Str(decimals) & "Ff",@FloatAnswer )
mpf_clear(@num1) : mpf_clear(@num2) : mpf_clear(@FloatAnswer)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function minus(number1 As String,number2 As String) As String'automatic precision
Dim As Integer Ln1=Len(number1),Ln2=Len(number2),decimals,Ln
If Ln1>=Ln2 Then Ln=Ln1 Else Ln=Ln2
Ln=ln+1
Dim As __mpf_struct num1, num2, FloatAnswer
mpf_init2(@num1, 4*(Ln1+1) )
mpf_init2(@num2, 4*(Ln2+1) )
mpf_init2(@Floatanswer, 4*(Ln))
mpf_set_str(@num1, number1, 10)
mpf_set_str(@num2, number2, 10)
Ln1=InStr(1,number1,"."):Ln2=InStr(1,number2,".")
If Ln1 Or Ln2 Then
decimals=Len(Mid(number1,Ln1+1))+Len(Mid(number2,Ln2+1))+1
End If
mpf_sub(@Floatanswer, @num1, @num2)
gmp_sprintf( @outtext,"%." & Str(decimals) & "Ff",@FloatAnswer )
mpf_clear(@num1) : mpf_clear(@num2) : mpf_clear(@FloatAnswer)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function Absolute Overload(a As String) As String
Dim As __mpf_struct Ab, Floatanswer
mpf_init2(@FloatAnswer, 4*precision )
mpf_init2(@Ab, 4*precision )
mpf_set_str(@Ab, a, 10)
mpf_abs(@FloatAnswer, @Ab)
gmp_sprintf( @outtext,"%." & Str(precision) & "Ff",@FloatAnswer )
mpf_clear(@Ab) : mpf_clear(@FloatAnswer)
Var outtxt=Trim(outtext)
If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
Return Trim(outtxt)
End Function
Function greater Overload(a As String,b As String) As Integer
Dim As mpf_t ma, mb
mpf_init2(@ma, 4*precision)
mpf_init2(@mb, 4*precision)
mpf_set_str(@ma, a, 10)
mpf_set_str(@mb, b, 10)
greater = greater(ma, mb)
mpf_clear(@ma) : mpf_clear(@mb)
End Function
Function equals Overload(a As String,b As String) As Integer
Dim As mpf_t ma, mb
mpf_init2(@ma, 4*precision)
mpf_init2(@mb, 4*precision)
mpf_set_str(@ma, a, 10)
mpf_set_str(@mb, b, 10)
equals = equals(ma, mb)
mpf_clear(@ma) : mpf_clear(@mb)
End Function
Function less Overload(a As String,b As String) As Integer
If equals(a,b) Then Return 0
If greater(a,b) Then Return 0
Return -1
End FunctionMakoto WATANABE
Please update your gmpFunctionsFrisian+DivIntJP.bas file with the new one, the old has still some memory leaks.
I have added a deterministic prime detection routine to your program and some other small changes to speed things up.
Code: Select all
#Include "gmpFunctionsFrisian+DivIntJP.bas"
#include "vbcompat.bi"
Declare Sub STRreplace(ByRef Expression As String, ByRef Find As String, ByRef Replacement As String, ByVal Start As Integer = 1)
Rem Factorization into Prime factors
Dim As String Suuchi, sqrroot_suuchi
Dim STARTT As Long
Dim ENDTIME As Long
Dim Minut As Integer
Dim A As String
Dim HANTEI As String
Dim Hairetsu() As String
Dim I As Integer
Dim J As Integer
Dim Kekka As String
' Dim Counter As Double
Dim As ULongInt Counter ' FreeBASIC 64bit integer
Dim f As Integer 'File number
' Find the first free file number
f = FreeFile
Open "Counter.txt" For Output As #f
If Err>0 Then Print "Error opening the file":End
'Function used
'greater,_mod,plus,divide,Sqrroot
KURIKAESHI:
set_precision(30)
Cls
Print "Factorization into Prime factors"
Print "Please input any integer."
Print "For example, input 12,319 or 4,294,967,297 and Enter key."
Print "For example, 99,400,891 or 6,041,375,340 or 123,456,789,013."
Print "For example, 2,305,843,008,139,952,128 ."
Print "For example, 70,000,104,900,007 or 70,000,005,000,000,007."
Print
Print "exceeds ULongInt: 18,446,744,073,709,551,617 = 274,177 * 67,280,421,310,721 ." '
Print "Try 1,111,111,111,111,111,111 (19 one's) = prime"
I=0
Kekka="1"
Counter=0
Line Input Suuchi
STRreplace(Suuchi, ",", "")
STARTT=Val(Left$(Time$,2))*3600+Val(Mid$(Time$,4,2))*60+Val(Right$(Time$,2))
' ------ prime test ------
' Determine whether n is prime. Return 2 if n is definitely prime, return 1 if n is probably
' prime (without being certain), or return 0 if n is definitely non-prime.
' Reasonable values of reps are between 15 and 50
If Len(suuchi) > 10 Then
Dim As __mpz_struct num2test
mpz_set_str(@num2test, suuchi, 10)
j = mpz_probab_prime_p (@num2test , 50 )
mpz_clear(@num2test)
If j = 2 Then Print suuchi; " = prime"
If j = 1 Then Print suuchi; " = probaly prime"
' j = 0, suuchi has factors
If j <> 0 Then GoTo label1
End If
Cls
Print Suuchi;"=";
'Gusu:
A = "2"
gusu: ' no need to start from 2 when we find a factor, start with the factor we have found
sqrroot_suuchi = sqrroot(suuchi) ' calculate at the start, no need to calculate it every loop
Tsugi:
Counter= Counter+1
If Counter Mod 1000000 =0 Then Print #f, Counter 'Write strings to a file using the number "f"
If greater(A,Sqrroot_Suuchi)=-1 Then GoTo Owari
If _mod (Suuchi,A)="0" Then GoTo Tsuzuku
If greater(A,"2")=-1 Then A=plus(A,"2") :GoTo Tsugi
A=plus(A,"1") :GoTo Tsugi
Tsuzuku:
Print A; "*" ;
I=I+1
ReDim Preserve Hairetsu(I)
Hairetsu(I)=A
Suuchi= _div(Suuchi,A)
' If InStr(Suuchi,".") Then Suuchi= RTrim(Suuchi,"0"):Suuchi=RTrim(Suuchi,".")
GoTo Gusu
Owari:
Print Suuchi
ENDTIME = Val(Left$(Time$,2))*3600+Val(Mid$(Time$,4,2))*60+Val(Right$(Time$,2))
Minut=(ENDTIME-STARTT)\60
Print Using "Processing time is #### minutes ## seconds"; Minut; (ENDTIME-STARTT)-Minut*60
Print
Print "Number of loops:"; Format(Counter, "#,##0")
Print
'Sleep
Print "prove the answer"
For J=1 To I
Print Hairetsu(J);"*";
Kekka = mult(Kekka,Hairetsu(J))
Next
Print Suuchi;"="
Print mult(Kekka,Suuchi)
label1:
Print
Print "Try again? If Y (or y) again."
Input "End, other than y . ",HANTEI
Print
If HANTEI="Y" Or HANTEI="y" Then
GoTo KURIKAESHI
EndIf
Close #f
Sub STRreplace(ByRef Expression As String, ByRef Find As String, ByRef Replacement As String, ByVal Start As Integer = 1)
Var p = InStr(Start, Expression, Find), li = Len(Find), ls = Len(Replacement) : If li = ls Then li = 0
While p
If li Then Expression = Left(Expression, p - 1) & Replacement & Mid(Expression, p + li) Else Mid(Expression, p) = Replacement
p = InStr(p + ls, Expression, Find)
Wend
End Sub