new header file GMP

General FreeBASIC programming questions.
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

Re: new header file GMP

Postby frisian » Dec 16, 2016 15:54

I have corrected two memory leaks in the Greater and Equals routine's. And added detection for floats being used in the integer _div and _mod routine's. Also added output trim for the divide and sqrroot routine's

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 Function

Makoto 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
dodicat
Posts: 6720
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: new header file GMP

Postby dodicat » Dec 16, 2016 20:02

Why not have primefactors as a function in it's own right?
Example (it takes a few seconds for the 111111 ...

Code: Select all

#Include Once "gmp.bi"

 
 Type Pfactors
    As String n
    As String f(Any)
    Declare Operator Cast() As String
End Type

Operator Pfactors.cast() As String   'printout
Print n,
If Ubound(f)=1 Then Return "(prime)"
For z As Integer=1 To Ubound(f)
    Print f(z);" ";
Next z
End Operator

Function primefactors(number As String) As pfactors
    print number,
    Dim As Long i
    Dim As zstring * 10000 outtext
    Dim As __mpz_struct zero,one,two,three,six,Zmod,divisor,Zdiv,tmp,root,num
    mpz_init(@one)
    mpz_init(@tmp)
    mpz_init(@six)
    mpz_init(@Zdiv)
    mpz_init(@Zmod)
    mpz_init(@zero)
    mpz_init(@two)
    mpz_init(@three)
    mpz_init(@divisor)
    mpz_init(@root)
    mpz_init(@num)
    mpz_set_str(@num, number, 10)
    mpz_set_str(@zero,"0", 10)
    mpz_set_str(@two,"2", 10)
    mpz_set_str(@three,"3", 10)
    mpz_set_str(@six,"6", 10)
    mpz_set_str(@one,"1", 10)
    Dim As pfactors p:p.n=number
    #macro fill(x)
    i+=1
    Redim Preserve p.f(1 To i)
    gmp_sprintf( @outtext,"%Zi", @x )
    p.f(i)=Trim(outtext)
    print p.f(i);" ";
    #endmacro
    #macro Eliminate(x)
    mpz_mod(@Zmod,@num,@x)
   
    While Mpz_cmp(@Zmod, @zero) = 0 'number mod x=0
       
        mpz_div(@Zdiv, @num, @x)
        mpz_init_set(@num, @Zdiv)
        fill(x)
        mpz_sqrt( @root,@num)
        mpz_add(@root,@root,@one) 'sqr(num)+1
        mpz_mod(@Zmod,@num,@x)
    Wend
   
    #endmacro
   
    mpz_sqrt( @root,@num)
    mpz_add(@root,@root,@one)
    Eliminate(two)
    Eliminate(three)
    While  Mpz_cmp(@divisor,@root)<0'divisor<root
        mpz_add(@divisor,@divisor,@six)'divisor+=6
        mpz_sub(@tmp,@divisor,@one)
        Eliminate (tmp)'((divisor-1))
        mpz_add(@tmp,@divisor,@one)
        Eliminate(tmp)'((divisor+1))
    Wend
   
    ' If number>1 Then
    If  Mpz_cmp(@num,@one)>0 Then
        fill(num)
    End If
    'clean up
    mpz_clear(@one): mpz_clear(@tmp): mpz_clear(@six): mpz_clear(@Zdiv): mpz_clear(@Zmod)
    mpz_clear(@zero): mpz_clear(@two): mpz_clear(@three): mpz_clear(@divisor): mpz_clear(@root): mpz_clear(@num)
    if ubound(p.f)=1 then print " (prime)";
    print
    Return p
End Function

Function SAR(s0 As String,Find As String,Replacement As String) As String
    Dim s As String=s0
    var position=Instr(s,Find)
    While position>0
        s=Mid(s,1,position-1) & Replacement & Mid(s,position+Len(Find))
        position=Instr(position+Len(Replacement),s,Find)
    Wend
    SAR=s
End Function

dim as pfactors x
Print "number",,"primefactors"
x= primefactors(sar("18,446,744,073,709,551,617",",",""))
x= primefactors(sar("2,305,843,008,139,952,128",",",""))
x= primefactors("12345678909876543212345678909876822")
x= primefactors(sar( "70,000,104,900,007",",",""))
x= primefactors(sar( "1,111,111,111,111,111,111",",",""))
x= primefactors(str(2*3*5*7*11*13*17*19ull))
print
print
'the x contains the last set of prime fcators
print x
print "Done"
sleep
dodicat
Posts: 6720
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: new header file GMP

Postby dodicat » Dec 17, 2016 0:12

Here with a check on the factors (--OK written).
And an option not to print from within the function.
Some prime factors of factorials:

Code: Select all

#Include Once "gmp.bi"

dim shared as zstring * 1000000 outtext
dim shared as long Iprint 'for printing within the function

Type Pfactors
    As String n
    As String f(Any)
    Declare Operator Cast() As String
End Type

Operator Pfactors.cast() As String   'printout
If Ubound(f)=1 Then Return "(prime)"
Dim As __mpz_struct tot,fz,n0
mpz_init(@tot):mpz_init(@fz):mpz_init(@n0):mpz_set_str(@n0,n,10)
mpz_set_str(@tot,"1", 10)
Print n;" = ",
For z As Integer=1 To Ubound(f)
    mpz_set_str(@fz,f(z),10)
    mpz_mul(@tot,@tot,@fz)
   if z<ubound(f) then Print f(z);"*"; else print f(z);
Next z
If Mpz_cmp(@tot,@n0) = 0 Then Print " --OK"
mpz_clear(@tot): mpz_clear(@fz): mpz_clear(@n0)
End Operator

Function primefactors(number As String) As pfactors
   if Iprint then Print number,
    Dim As Long i
    Dim As zstring * 10000 outtext
    Dim As __mpz_struct zero,one,two,three,six,Zmod,divisor,Zdiv,tmp,root,num
   Dim As __mpz_struct tot,fz,n0
    mpz_init(@tot):mpz_init(@fz):mpz_init(@n0):mpz_set_str(@n0,number,10)
    mpz_set_str(@tot,"1", 10)
    mpz_init(@one)
    mpz_init(@tmp)
    mpz_init(@six)
    mpz_init(@Zdiv)
    mpz_init(@Zmod)
    mpz_init(@zero)
    mpz_init(@two)
    mpz_init(@three)
    mpz_init(@divisor)
    mpz_init(@root)
    mpz_init(@num)
    mpz_set_str(@num, number, 10)
    mpz_set_str(@zero,"0", 10)
    mpz_set_str(@two,"2", 10)
    mpz_set_str(@three,"3", 10)
    mpz_set_str(@six,"6", 10)
    mpz_set_str(@one,"1", 10)
    Dim As pfactors p:p.n=number
    #macro fill(x)
    i+=1
    Redim Preserve p.f(1 To i)
    gmp_sprintf( @outtext,"%Zi", @x )
    p.f(i)=Trim(outtext)
    if Iprint then
    Print p.f(i);" ";
    mpz_set_str(@fz,p.f(i),10)
    mpz_mul(@tot,@tot,@fz)
    end if
    #endmacro
    #macro Eliminate(x)
    mpz_mod(@Zmod,@num,@x)
   
    While Mpz_cmp(@Zmod, @zero) = 0 'number mod x=0
       
        mpz_div(@Zdiv, @num, @x)
        mpz_init_set(@num, @Zdiv)
        fill(x)
        mpz_sqrt( @root,@num)
        mpz_add(@root,@root,@one) 'sqr(num)+1
        mpz_mod(@Zmod,@num,@x)
    Wend
   
    #endmacro
   
    mpz_sqrt( @root,@num)
    mpz_add(@root,@root,@one)
    Eliminate(two)
    Eliminate(three)
    While  Mpz_cmp(@divisor,@root)<0'divisor<root
        mpz_add(@divisor,@divisor,@six)'divisor+=6
        mpz_sub(@tmp,@divisor,@one)
        Eliminate (tmp)'((divisor-1))
        mpz_add(@tmp,@divisor,@one)
        Eliminate(tmp)'((divisor+1))
    Wend
   
    ' If number>1 Then
    If  Mpz_cmp(@num,@one)>0 Then
        fill(num)
    End If
   if Iprint then If Mpz_cmp(@tot,@n0) = 0 Then Print " --OK";
    '
    'clean up
    mpz_clear(@one): mpz_clear(@tmp): mpz_clear(@six): mpz_clear(@Zdiv): mpz_clear(@Zmod)
    mpz_clear(@zero): mpz_clear(@two): mpz_clear(@three): mpz_clear(@divisor): mpz_clear(@root): mpz_clear(@num)
   
    mpz_clear(@tot): mpz_clear(@fz): mpz_clear(@n0)
    if Iprint then If Ubound(p.f)=1 Then Print " (prime)";
    if Iprint then Print
    Return p
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

#define range(f,l) int(Rnd*((l+1)-(f))+(f))
Iprint=0
randomize
do
Dim As pfactors x
dim as ulong u=range(10,1000)
dim as string f=factorial(u)
Print "number = ";"factorial(";u;")"
x=primefactors(f)
Print
Print "from cast()"
Print x
print "______________________________"
sleep
loop until inkey=chr(27)

Print "Done"
Sleep
Makoto WATANABE
Posts: 196
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

Re: new header file GMP

Postby Makoto WATANABE » Dec 17, 2016 13:08

Dear frisian
Dear dodicat
Dear srvaldez

Thank you very much for giving big Christmas presents.
You all are Jesus-like.
I am sure all of you are in God's zone.


Dear frisian

Thank you for providing the 16 December 2016 version of the GMP user-defined functions and the revised prime_factors program.
I tried prime factorization of 21 digit number using your programs.
And I got the following result. (It took time a little.)
121,439,531,096,594,251,777 = 3,930,785,153 * 30,894,471,809
I am very grateful to you.


Dear dodicat

Thank you for providing the program of primefactors as a function .
I tried prime factorization of 21 digit number(121,439,531,096,594,251,777) using your program.
It showed the results in just a few minutes (only 2 minutes 38 seconds !).

I would like to translate your program into Japanese and introduce it to Japanese people on my website.
Please consent to this.


Dear srvaldez

Thank you for your posting.
I am amazed at your "gmpf-overload.bi".

I realized the benefits of using "user-defined types" and "operator overloading".
I was able to increase the number of the processing digits without modifying the simple Basic program of integer type.

I would like to translate your "gmpf-overload.bi" into Japanese and introduce it to Japanese people on my website.
Please consent to this.

If you add your "gmpf-overload.bi" to the FreeBASIC document (e.g. following page) ,
I think it will be a Christmas present for people all over the world.
Http://www.freebasic.net/wiki/wikka.php?wakka=ExtLibgmp
http://www.freebasic.net/wiki/wikka.php?wakka=CodeLibrary

Code: Select all

#Include "gmpf-overload.bi"
# include "vbcompat.bi"

Rem Factorization into Prime factors
 
Declare Sub STRreplace(ByRef Expression As String, ByRef Find As String, ByRef Replacement As String, ByVal Start As Integer = 1)
 
 
Dim Suuchi As String
Dim N As gmpf
Dim STARTT As Long
Dim ENDTIME As Long
Dim Minut As Integer
Dim A As gmpf
Dim HANTEI As String
Dim Hairetsu() As gmpf
Dim I As Integer
Dim J As Integer
Dim Kekka As gmpf
Dim Counter As LongInt
 
set_precision(20) 'set precision to at least 20 digits
 
 KURIKAESHI:
 
 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 "The followings are beyond ULongInt, but ..."
 PRINT " 18,446,744,073,709,551,617 = 274,177 * 67,280,421,310,721 "
 Print "121,439,531,096,594,251,777 = 3,930,785,153 * 30,894,471,809  (This takes 2 hours?)"
 Print

 I=0
 Kekka=1
 Counter=0
 
  Line Input Suuchi
 
 STRreplace(Suuchi, ",", "")
 N = Suuchi

 STARTT=VAL(LEFT$(TIME$,2))*3600+VAL(MID$(TIME$,4,2))*60+VAL(RIGHT$(TIME$,2))
 
 Cls
 PRINT Suuchi ;"=";
 
Gusu:
   A=2

Tsugi:
   Counter= Counter+1

   IF A>SQR(N) THEN GOTO Owari
   IF N-INT(N/A)*A=0 THEN GOTO Tsuzuku
   IF A>2 THEN A=A+2 :GOTO Tsugi
   A=A+1 :GOTO Tsugi

Tsuzuku:
 PRINT A; "*" ;
 I=I+1
 ReDim Preserve Hairetsu(I)
 Hairetsu(I)=A
 N=N/A : GOTO Gusu

Owari:
 PRINT N
 
 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

 Print "prove the answer"
 For J=1 To I
    Print Hairetsu(J);"*";
    Kekka=Kekka*Hairetsu(J)
 Next
 Print N;"="
 Print Kekka*N
 Print
 Print "Try again? If Y (or y) again."
 Input "End, other than y. ",HANTEI
 If HANTEI="Y" Or HANTEI="y" Then
    GoTo KURIKAESHI
 EndIf
 
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
dodicat
Posts: 6720
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: new header file GMP

Postby dodicat » Dec 17, 2016 16:59

For a 3 Ghz box (Dell), dual core.
Results for factors of 121439531096594251777

Code: Select all

GMP version 6.1.0

FreeBASIC Compiler - Version 1.05.0 (01-31-2016), built for win32 (32bit)
OS:    Windows NT 6.2 (build 9200) (win 10 in other words)

'______________________________________________________
fbc 32 bit -gen gas
121439531096594251777 =     3930785153*30894471809 --OK

 99.4743767332111 seconds
 '______________________________________________________
 
 fbc 32 bit -gen gcc
 121439531096594251777 =     3930785153*30894471809 --OK

 104.3219764681999 seconds
 '______________________________________________________
 fbc 32 bit -gen gcc -Ofast
 121439531096594251777 =     3930785153*30894471809 --OK

 99.90048057315289 seconds
 
 '======================================================
 
 GMP version 6.1.1
 FreeBASIC Compiler - Version 1.05.0 (01-31-2016), built for win64 (64bit)
 OS:    Windows NT 6.2 (build 9200) (win 10 in other words)
'_______________________________________________________
 fbc 64 bit -gen gcc -Ofast
 121439531096594251777 =     3930785153*30894471809 --OK

 129.5779795260169 seconds
 '______________________________________________________
 
 
 
 

So once again, 32 bit -gen gas beats the rest.
.. and the 64 bit compiler is sluggish (once again).

Makoto WATANABE.
Anything you find usable from my code, please use it, alter it, or whatever, it's up to you!
srvaldez
Posts: 2529
Joined: Sep 25, 2005 21:54

Re: new header file GMP

Postby srvaldez » Dec 17, 2016 20:28

hello Makoto WATANABE
you can use my code anyway you want :-)
Makoto WATANABE
Posts: 196
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

Re: new header file GMP

Postby Makoto WATANABE » Dec 18, 2016 13:17

Dear dodicat
Dear srvaldez

Thanks for your quick reply. Thank you for your consent.
Japan is winter now, the outside is cold and there is ice-skin.
However, my heart is warm.
I am very grateful for your help.

Return to “General”

Who is online

Users browsing this forum: No registered users and 6 guests