Big Number Wrapper (GMP_INT.BI)

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Post Reply
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

Big Number Wrapper (GMP_INT.BI)

Post by frisian »

updated the GMP_INT.BI code (22 may 2017). (typo's, removed some redundant comments).
updated Cipolla.bas to work with the Type code that srvaldez provided.

updated the GMP_INT.BI code and added fdiv_cdiv.bas (20 may 2017).

I extended Yetifoot's Big number wrapper with lots of new stuff making working with GMP integer much easier.
For the moment I have no plans to add/change things. It seems to me that it's fairly complete.

Since GMP 4.0 there are some things have been added so this version will not work. Best to use the newer version.

gmp 6.1.2 from srvaldez (top post)
http://www.freebasic.net/forum/viewtopi ... 55#p230255

GMP_INT.BI

Code: Select all

' version 22/05/2017
' Original by yetifoot http://www.freebasic.net/forum/viewtopic.php?t=7173
' May 2017: srvaldez
' fixed the loop code

' May 2017: frisian
' Use on your own risk
' WARNING: for this file to work you need GMP 6.0 32bit (tested).
'          GMP 4.0 will not work (some functions where later added).
'          GMP 5.0 not tested
'
' changed mpz_mod into mpz_tdiv_r to have mod mimic the default FB behavior
' \ (integer divide) and mod both use tdiv (truncate)
' this behavior does not change when you change the way FB handles float to integer conversion
' you can use mpz_fdiv if you need to floor the result
'  "   "   "  mpz_cdiv if you need to ceil the result
'
' to get the remainder from the floating point division (/) you need some old school basic
' remainder = number - (number / divisor) * divisor
'
' you can create array's and type's
' swap also works
'
' corrected and improved the floating point division (/)
' added operators Shr and Shl
' corrected some errors in comments i made, cleaned the listing
' added function gmp_int_2_base (gmp_int, n = 2) returns a string of a gmp_int number in base 2 to 62 (default 10, decimal)
' added setbit (gmp_int, n) set n'th bit to 1,
' added clrbit (gmp_int, n) set n'th bit to 0
' added combit (gmp_int, n) flips the n'th bit (1 -> 0, 0 -> 1)
' added tstbit (gmp_int, n) tests the n'th bit, returns the value of the bit (0 or 1)
' added root (gmp_int, n = 2) returns the truncated part of the n'th root (default = 2, square root)
'       root (gmp_int, 3) for cube root
' added rootrem (gmp_int, n = 2) returns gmp_int - root(gmp_int, n) ^ n
' added powm (gmp_int1, n, gmp_int2) returns (gmp_int1 ^ n) mod gmp_int2
' added isprime (gmp_int, reps) returns 2 if n is definitely prime
'                               returns 1 if n is probably prime (without being certain)
'                               returns 0 if n is definitely non-prime
'                               reasonable values of reps are between 15 and 50. (default = 15)
' added nextprime (gmp_int) finds the next prime that is greater then gmp_int
' added legendre (gmp_int1, gmp_int2) returns the Legendre symbol (-1, 0, 1)
' added gcd (gmp_int1, gmp_int2) returns the greatest common divisor
' added lcm (gmp_int1, gmp_int2) returns the least common multiple
' added fac (n)  factorial, returns n!
' added mfac (n, m) m-multi-factorial, returns n!^(m) = n*(n-m)*(n-2m)*...
' added primorial (n) primorial, returns the product of all positive prime numbers <= n
' added fib (n) returns the n'th Fibonacci number
'
' added rnd (gmp_int = 0) mersenne twister, returns a random number between 0 and gmp_int -1
'                         if no number is given or the number <= 0 then the number is set to 2 ^ 32
'                         the random generator is seeded when the routine is called for the first time
' if you do a randomize timer the seed will be different every time
' if you don't use randomize or do a randomize "number" the sequence will be the same (just as the FB rnd)
' setbit, clrbit and combit do not return any value
' tstbit, isprime, and legendre return a long
' root, rootrem, powm, nextprime, gcd, lcm, fac, mfac, primorial, fib, gmp_rnd return a gmp_int
'
' overloaded operators +=, -=, *=, \=, /=, ^=, Mod=, And=, Or=, Xor=, Shl= ans Shr= to work with gmp_int's
'
' overloaded abs() and sgn() to work with gmp_int's

#Ifndef __GMP_INT_BI__
#Define __GMP_INT_BI__

#Include Once "gmp.bi"

Type gmp_int
    Declare Constructor ( )
    Declare Constructor ( ByVal i As Long )
    Declare Constructor ( ByRef s As String, ByVal _base As Long = 10 )
    Declare Constructor ( ByRef g As gmp_int )
    Declare Destructor ( )
    Declare Operator Let ( ByRef g As gmp_int )
    Declare Operator Let ( ByVal i As Long)
    Declare Operator Let ( ByRef s As String )
    Declare Operator Cast ( ) As Long
    Declare Operator Cast ( ) As String
    ' For Next Implicit step = +1
    Declare Operator For ( )
    Declare Operator Step( )
    Declare Operator For ( ByRef stp As gmp_int )
    Declare Operator Step ( ByRef stp As gmp_int )
    Declare Operator Next ( ByRef end_cond As gmp_int ) As Integer
    Declare Operator Next ( ByRef end_cond As gmp_int, ByRef step_var As gmp_int ) As Integer
    Declare Function getString ( ByVal _base As Long = 10 ) As String
    num As mpz_ptr = 0

    Declare Operator += (ByRef rhs As gmp_int)
    Declare Operator -= (ByRef rhs As gmp_int)
    Declare Operator *= (ByRef rhs As gmp_int)
    Declare Operator \= (ByRef rhs As gmp_int)
    Declare Operator /= (ByRef rhs As gmp_int)
    Declare Operator ^= (ByVal rhs As Long)
    Declare Operator Mod= (ByRef rhs As gmp_int)
    Declare Operator And= (ByRef rhs As gmp_int)
    Declare Operator Or= (ByRef rhs As gmp_int)
    Declare Operator Xor= (ByRef rhs As gmp_int)
    Declare Operator Shl= (ByVal rhs As Long)
    Declare Operator Shr= (ByVal rhs As Long)

End Type

Declare Operator + ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int
Declare Operator - ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int
Declare Operator * ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int
Declare Operator \ ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int
Declare Operator / ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int
Declare Function root ( ByRef op As gmp_int, ByVal n As Long = 2 ) As gmp_int
Declare Function rootrem ( ByRef op As gmp_int, ByVal n As Long = 2 ) As gmp_int
Declare Operator ^ ( ByRef lhs As gmp_int, ByVal rhs As Long ) As gmp_int
Declare Operator Mod ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int
Declare Function powm (ByRef op As gmp_int, ByRef power As gmp_int, ByRef modulus As gmp_int ) As gmp_int
Declare Operator And ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int
Declare Operator Or ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int
Declare Operator Xor ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int
Declare Operator Shl ( ByRef lhs As gmp_int, ByVal rhs As Long ) As gmp_int
Declare Operator Shr ( ByRef lhs As gmp_int, ByVal rhs As Long ) As gmp_int
Declare Operator Sgn ( ByRef rhs As gmp_int ) As Long
Declare Operator Abs ( ByRef rhs As gmp_int ) As gmp_int
Declare Sub setbit ( ByRef op As gmp_int, ByVal n As Long )
Declare Sub clrbit ( ByRef op As gmp_int, ByVal n As Long )
Declare Sub combit ( ByRef op As gmp_int, ByVal n As Long )
Declare Function tstbit ( ByRef op As gmp_int, ByVal n As Long ) As Long
Declare Operator - ( ByRef rhs As gmp_int ) As gmp_int
Declare Operator Not ( ByRef rhs As gmp_int ) As gmp_int
Declare Operator = ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long
Declare Operator < ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long
Declare Operator > ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long
Declare Operator <= ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long
Declare Operator >= ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long
Declare Operator <> ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long
Declare Function gmp_int_2_base (ByRef bigint As gmp_int, ByVal _base As Long = 10 ) As String
Declare Function isprime ( ByRef op As gmp_int, ByVal reps As Long = 15 ) As Long
Declare Function nextprime ( ByRef op As gmp_int ) As gmp_int
Declare Function legendre ( ByRef a As gmp_int , ByRef p As gmp_int ) As Long
Declare Function gcd ( ByRef a As gmp_int, ByRef b As gmp_int ) As gmp_int
Declare Function lcm ( ByRef a As gmp_int, ByRef b As gmp_int ) As gmp_int
Declare Function fac (ByVal n As Long ) As gmp_int
Declare Function mfac (ByVal n As Long, ByVal m As Long ) As gmp_int ' result = n!^(m) = n*(n-m)*(n-2m)*...
Declare Function primorial (ByVal n As Long ) As gmp_int
Declare Function fib ( ByVal n As Long ) As gmp_int
Declare Function gmp_rnd ( ByRef max As gmp_int = 0) As gmp_int

#EndIf '__GMP_INT_BI__

'::::::::
Constructor gmp_int ( )

num = Callocate( SizeOf( __mpz_struct ) )
mpz_init( num )

End Constructor

'::::::::
Constructor gmp_int ( ByVal i As Long )

num = Callocate( SizeOf( __mpz_struct ) )
mpz_init_set_si( num, i )

End Constructor

'::::::::
Constructor gmp_int ( ByRef s As String, ByVal _base As Long = 10 )

num = Callocate( SizeOf( __mpz_struct ) )
mpz_init_set_str( num, StrPtr( s ), _base )

End Constructor

'::::::::
Constructor gmp_int ( ByRef g As gmp_int )

num = Callocate( SizeOf( __mpz_struct ) )
mpz_init_set( num, g.num )

End Constructor

'::::::::
Destructor gmp_int ( )

mpz_clear( num )
DeAllocate( num )

End Destructor

'::::::::
Operator gmp_int.let ( ByRef g As gmp_int )

    mpz_set( num, g.num )

End Operator

'::::::::
Operator gmp_int.let ( ByVal i As Long )

    mpz_set_si( num, i )

End Operator

'::::::::
Operator gmp_int.let ( ByRef s As String )

    mpz_set_str( num, StrPtr(s), 10 )

End Operator

Operator gmp_int.cast ( ) As Long
Operator = mpz_get_si(num)
End Operator

'::::::::
Operator gmp_int.cast ( ) As String

Operator = getString( 10 )

End Operator

'=====================================================================
' Implicit step   gmp_int.For   gmp_int.Next   gmp_int.Step  is +1

Operator gmp_int.for ( )
End Operator

Operator gmp_int.step ( )
    This += 1 'this = this+1 '
End Operator

Operator gmp_int.next ( ByRef end_cond As gmp_int ) As Integer
    Return This <= end_cond
End Operator

'' explicit step versions
''
Operator gmp_int.for ( ByRef step_var As gmp_int )
End Operator

Operator gmp_int.step ( ByRef step_var As gmp_int )
    This += step_var 'this = this + step_var '
End Operator

Operator gmp_int.next ( ByRef end_cond As gmp_int, ByRef step_var As gmp_int ) As Integer
    If step_var < 0 Then
        Return This >= end_cond
    Else
        Return This <= end_cond
    End If
End Operator

'::::::::
Function gmp_int.getString ( ByVal _base As Long = 10 ) As String

    Dim As ZString Ptr s = mpz_get_str( 0, _base, num )

    If s Then
        Function = *s
        DeAllocate( s )
    End If

End Function

'::::::::
Operator + ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_add( result.num, lhs.num, rhs.num )

    Operator = result

End Operator

'::::::::
Operator - ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_sub( result.num, lhs.num, rhs.num )

    Operator = result

End Operator

'::::::::
Operator * ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_mul( result.num, lhs.num, rhs.num )

    Operator = result

End Operator

'::::::::
Operator \ ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_tdiv_q( result.num, lhs.num, rhs.num )

    Operator = result

End Operator

'::::::::
Operator / ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int

    ' This really is quite nasty, but it's the best way I could come up with
    ' to do an FP divide, while rounding results the way people may expect them.
    ' 0.5 becomes 1, 0.4 becomes 0, -0.5 becomes -1, -0.4 becomes 0
    ' I'm not 100% that temp is needed, but I used it just to be safe, for some
    ' GMP functions it seems a bad idea to pass an mpz as dest and one of the sources.

    Dim As mpf_ptr l, r, q, half, temp
    Dim As gmp_int result

    ' size for gmp floats is fixed, this could lead to loss of precision
    ' the size is set by mpf_set_default_prec or by using mpf_init2( 'var', size)
    ' find the size of lhs and rhs, take the longest size and add 32 bit (9 digits)
    Dim As UInteger l1 = mpz_sizeinbase ( lhs.num, 2 )
    Dim As UInteger r1 = mpz_sizeinbase ( rhs.num, 2 )
    If r1 > l1 Then l1 = r1
    l1 = l1 +32

    l = Callocate( SizeOf( __mpf_struct ) )
    r = Callocate( SizeOf( __mpf_struct ) )
    q = Callocate( SizeOf( __mpf_struct ) )
    half = Callocate( SizeOf( __mpf_struct ) )
    temp = Callocate( SizeOf( __mpf_struct ) )

    mpf_init2( l, l1 )          ' set the size for l to l1
    mpf_init2( r, l1 )          ' set the size for r to l1
    mpf_init2( q, l1 )          ' set the size for q to l1
    mpf_init_set_d( half, .5 )
    mpf_init2( temp, l1 )       ' set the size for temp to l1

    mpf_set_z( l, lhs.num )
    mpf_set_z( r, rhs.num )
    mpf_div( q, l, r )

    Dim As Long cmp_val = mpf_cmp_d( q, 0 )

    If cmp_val > 0 Then
        mpf_add( temp, q, half )
    ElseIf cmp_val < 0 Then
        mpf_sub( temp, q, half )
    Else
        mpf_set( temp, q )
    End If

    mpz_set_f( result.num, temp)

    mpf_clear( l )
    mpf_clear( r )
    mpf_clear( q )
    mpf_clear( half )
    mpf_clear( temp )

    DeAllocate( l )
    DeAllocate( r )
    DeAllocate( q )
    DeAllocate( half )
    DeAllocate( temp )

    Operator = result

End Operator

Function root ( ByRef op As gmp_int, ByVal n As Long = 2 ) As gmp_int

    Dim As gmp_int result

    mpz_root( result.num, op.num, n )

    Function = result

End Function

Function rootrem ( ByRef op As gmp_int, ByVal n As Long = 2 ) As gmp_int

    Dim As gmp_int result, tmp

    mpz_rootrem( tmp.num, result.num, op.num, n )

    Function = result

End Function

'::::::::
Operator ^ ( ByRef lhs As gmp_int, ByVal rhs As Long ) As gmp_int

    Dim As gmp_int result

    mpz_pow_ui( result.num, lhs.num, rhs )

    Operator = result

End Operator


Operator Mod ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_tdiv_r( result.num, lhs.num, rhs.num )

    Operator = result

End Operator

Function powm ( ByRef _base As gmp_int, ByRef power As gmp_int, ByRef modulus As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_powm( result.num, _base.num, power.num, modulus.num )

    Return result

End Function

'::::::::
Operator And ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_and( result.num, lhs.num, rhs.num )

    Operator = result

End Operator

'::::::::
Operator Or ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_ior( result.num, lhs.num, rhs.num )

    Operator = result

End Operator

'::::::::
Operator Xor ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_xor( result.num, lhs.num, rhs.num )

    Operator = result

End Operator

Operator Shl ( ByRef lhs As gmp_int, ByVal rhs As Long ) As gmp_int

    Dim As gmp_int result

    mpz_mul_2exp( result.num, lhs.num, rhs )

    Operator = result

End Operator

Operator Shr ( ByRef lhs As gmp_int, ByVal rhs As Long ) As gmp_int

    Dim As gmp_int result

    mpz_tdiv_q_2exp( result.num, lhs.num, rhs )

    Operator = result

End Operator

Operator Sgn ( ByRef rhs As gmp_int ) As Long

    Return mpz_sgn( rhs.num )

End Operator

Operator Abs ( ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_abs(result.num, rhs.num )

    Return result

End Operator

Sub setbit ( ByRef op As gmp_int, ByVal n As Long)

    mpz_setbit( op.num, n )

End Sub

Sub clrbit ( ByRef op As gmp_int, ByVal n As Long)

    mpz_clrbit( op.num, n )

End Sub

Sub combit ( ByRef op As gmp_int, ByVal n As Long)

    mpz_combit( op.num, n )

End Sub

Function tstbit ( ByRef op As gmp_int, ByVal n As Long) As Long
    ' return 0 if the bit is 0 and 1 if the bit is 1
    Return mpz_tstbit( op.num, n )

End Function

'::::::::
Operator - ( ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_neg( result.num, rhs.num )

    Operator = result

End Operator

'::::::::
Operator Not ( ByRef rhs As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_com( result.num, rhs.num )

    Operator = result

End Operator

'::::::::
Operator = ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long

Operator = (mpz_cmp( lhs.num, rhs.num ) = 0)

End Operator

'::::::::
Operator < ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long

Operator = (mpz_cmp( lhs.num, rhs.num ) < 0)

End Operator

'::::::::
Operator > ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long

Operator = (mpz_cmp( lhs.num, rhs.num ) > 0)

End Operator

'::::::::
Operator <= ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long

Operator = (mpz_cmp( lhs.num, rhs.num ) <= 0)

End Operator

'::::::::
Operator >= ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long

Operator = (mpz_cmp( lhs.num, rhs.num ) >= 0)

End Operator

'::::::::
Operator <> ( ByRef lhs As gmp_int, ByRef rhs As gmp_int ) As Long

Operator = (mpz_cmp( lhs.num, rhs.num ) <> 0)

End Operator

Operator gmp_int.+= (ByRef rhs As gmp_int)
    mpz_add(This.num, this.num, rhs.num)
End Operator

Operator gmp_int.-= (ByRef rhs As gmp_int)
    mpz_sub(This.num, this.num, rhs.num)
End Operator

Operator gmp_int.*= (ByRef rhs As gmp_int)
    mpz_mul(This.num, this.num, rhs.num)
End Operator

Operator gmp_int.\= (ByRef rhs As gmp_int)
    mpz_tdiv_q(This.num, this.num, rhs.num)
End Operator

Operator gmp_int./= (ByRef rhs As gmp_int)
    This = This / rhs
End Operator

Operator gmp_int.^= (ByVal rhs As Long)
    mpz_pow_ui(This.num, this.num, rhs)
End Operator

Operator gmp_int.mod= (ByRef rhs As gmp_int)
    mpz_tdiv_r( this.num, this.num, rhs.num )
End Operator

Operator gmp_int.And= (ByRef rhs As gmp_int)
    mpz_and(This.num, this.num, rhs.num)
End Operator

Operator gmp_int.Or= (ByRef rhs As gmp_int)
    mpz_ior(This.num, this.num, rhs.num)
End Operator

Operator gmp_int.Xor= (ByRef rhs As gmp_int)
    mpz_xor(This.num, this.num, rhs.num)
End Operator

Operator gmp_int.shl= (ByVal rhs As Long)
    mpz_mul_2exp(This.num, this.num, rhs)
End Operator

Operator gmp_int.Shr= (ByVal rhs As Long)
    mpz_tdiv_q_2exp(This.num, this.num, rhs)
End Operator

' return a string of a gmp_int number in base 2 to 62 (default 10, decimal)
Function gmp_int_2_base (ByRef bigint As gmp_int, ByVal _base As Long = 10 ) As String

    Dim As ZString Ptr s = mpz_get_str( 0, _base, bigint.num )

    If s Then
        Function = *s
        DeAllocate( s )
    End If

End Function

Function isprime ( ByRef op As gmp_int, ByVal reps As Long = 15 ) As Long
    ' returns if n is definitely prime
    ' returns 1 if n is probably prime (without being certain)
    ' returns 0 if n is definitely non-prime
    ' reasonable values of reps are between 15 and 50.

    If reps < 15 Then reps = 15

    Function = mpz_probab_prime_p ( op.num, reps )

End Function

Function nextprime ( ByRef op As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_nextprime( result.num , op.num)

    Function = result

End Function

Function legendre ( ByRef a As gmp_int , ByRef p As gmp_int ) As Long
    ' p needs to be a odd positive prime
    Function = mpz_legendre ( a.num, p.num)

End Function

Function gcd ( ByRef a As gmp_int, ByRef b As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_gcd( result.num, a.num, b.num )

    Function = result

End Function

Function lcm ( ByRef a As gmp_int, ByRef b As gmp_int ) As gmp_int

    Dim As gmp_int result

    mpz_lcm( result.num, a.num, b.num )

    Function = result

End Function

Function fac ( ByVal n As Long ) As gmp_int

    Dim As gmp_int result

    mpz_fac_ui( result.num , n )

    Function = result

End Function

Function mfac ( ByVal n As Long, ByVal m As Long ) As gmp_int
    ' result = n!^(m) = n*(n-m)*(n-2m)*...
    Dim As gmp_int result

    mpz_mfac_uiui ( result.num, n, m )

    Function = result

End Function

Function primorial ( ByVal n As Long ) As gmp_int

    Dim As gmp_int result

    mpz_primorial_ui (result.num, n )

    Function = result

End Function

Function fib ( ByVal n As Long ) As gmp_int

    Dim As gmp_int result

    mpz_fib_ui ( result.num, n)

    Function  = result

End Function

' -= Random number generator stuff (integer) =-
Type gmp_random
    Declare Constructor ( )
    Declare Destructor ( )
    g_rnd As Any Ptr
End Type

Constructor gmp_random ( )

g_rnd = Callocate( SizeOf ( __gmp_randstate_struct ) )
gmp_randinit_mt ( g_rnd )   ' Mersenne Twister

End Constructor

Destructor gmp_random ( )

gmp_randclear ( g_rnd )
DeAllocate ( g_rnd )

End Destructor

Function gmp_rnd ( ByRef max As gmp_int = 0 ) As gmp_int

    Static state As gmp_random
    Static flag As Long

    If flag = 0 Then

        Dim As String seed_str
        Dim As gmp_int seed

        For i As Long = 0 To 200 ' create seed
            seed_str += Str( Int( Rnd * 10 ) )
        Next

        seed = seed_str
        Gmp_randseed ( state.g_rnd, seed.num )   ' seed the random generator

        flag = 1

    End If

    If max <= 0 Then
        max = 2
        max ^= 32        ' 2 ^ 32
    End If
    Dim As gmp_int result
    Mpz_urandomm ( result.num , state.g_rnd, max.num )

    Function = result

End Function
Test.bas

Code: Select all

#Include Once "gmp_int.bi"

Dim As gmp_int i, j, two = 2 'if the number is small
' if the number is larger than FB integer then use: two = gmp_int("2")
Dim As Long k

' Powers of two
Print "Some powers of two"
For i = 0 To 128
    Print "2 ^ " & i & " = " & two ^ i
Next i

Print
Print "casting a gmp_int to a FreeBASIC Integer"
Print "this should only be done if you know the gmp_int will fit in a Integer"

k = two
Print "k = two ->";k

Print
i = "1" + String(35, "0")
i = i + 1
Print i; " / 3 = "; i / 3
Print

i = gmp_int(String(40,"1"), 2)
Print "&B" + String(40,"1") + " = ";i

Print
Print i; " shr 39 = "; i Shr 39
Print

i = 1
Print i; " shl 128 = "; i Shl 128

Print
i = "1234567890"
Print i; " in base 62 = "; gmp_int_2_base(i, 62)
Print i; " in base 16 = "; gmp_int_2_base(i, 16)
Print i; " in base 8  = "; gmp_int_2_base(i, 8)
Print i; " in base 2  = "; gmp_int_2_base(i, 2)

Print
i = gmp_int("09azAZ", 62)
Print "09azAZ(62) = "; i

Print
i = 10000
Print "The (truncated) square root of "; i; " = "; root(i)
Print "The remainder = "; rootrem(i)
Print
Print "The (truncated) 5th root of "; i; " = "; root(i, 5)
Print "The remainder = "; rootrem(i, 5)
Print root(i, 5); " ^ 5 + "; root(i, 5) ^ 5; " + "; rootrem(i, 5); " = ";
Print root(i, 5) ^ 5 + rootrem(i ,5); ", i = "; i

Print
i = 40
Print "  starting value "; gmp_int_2_base(i, 2)
setbit (i, 0)
Print "       set bit 0 "; gmp_int_2_base(i, 2)
clrbit (i, 0)
Print "     clear bit 0 "; gmp_int_2_base(i, 2)
combit (i, 1)
Print "compliment bit 1 "; gmp_int_2_base(i, 2)
combit (i, 1)
Print "compliment bit 1 "; gmp_int_2_base(i, 2)
Print

For k = 0 To 5
    Print "bit"; k;
    If tstbit(i, k) = 0 Then Print " is 0"
    If tstbit(i, k) = 1 Then Print " is 1"
Next

Print
i = 10
Print i; " ^ 2 mod 49 = "; powm(i, 2, 49)

Print
Print i;
k = isprime(i)
If k = 2 Then
    Print " is prime"
ElseIf k = 1 Then
    Print " is probably prime (without being certain)"
ElseIf k = 0 Then
    Print " is not a prime"
End If

Print
i = 2 : i = i ^ 61 -1 ' mersenne prime
Print i; " run test 15 times (default)";
k = isprime(i)
If k = 2 Then
    Print " is prime"
ElseIf k = 1 Then
    Print " is probably prime (without being certain)"
ElseIf k = 0 Then
    Print " is not a prime"
End If

Print i; " ,run test 50 times,";
k = isprime(i , 50)
If k = 2 Then
    Print " is prime"
ElseIf k = 1 Then
    Print " is probably prime (without being certain)"
ElseIf k = 0 Then
    Print " is not a prime"
End If

Print
Print i; ", next prime = "; nextprime(i)

Print
i = 20 : j = 100
Print "gcd("; i; ", "; j; ") = "; gcd(i, j)
Print "lcm("; i; ", "; j; ") = "; lcm(i, j)
Print
i = "1234567890123456789" : j = "2345612345623456123456"
Print "gcd("; i; ", "; j; ") = "; gcd(i, j)
Print "lcm("; i; ", "; j; ") = "; lcm(i, j)

Print
Print "fac(10) = "; fac(10)
Print ' mfac(n, m)  = n!^(m) = n*(n-m)*(n-2m)*...
Print "mfac(10, 7) = "; mfac(10, 7) ' 10 * 3

Print
Print "primorial(10) = "; primorial(10) ' 2 * 3 * 5 * 7

Print
Print "the 100'th Fibonacci number = "; fib(100)

Print

Dim As gmp_int array(1 To 2, 1 To 10)

For k = 1 To 10
    array(1, k) = gmp_rnd(10)
    array(2, k) = gmp_rnd ' returns a number between 0 and 2^32-1
Next

For k = 1 To 10
    Print array(1, k), array(2, k)
Next

Print
Print gcd(gmp_int("123456789012345678901234567890123456789012345678901234567890"), gmp_int("123456789012345678901234567890"))

Print
Print  gmp_int("123456789012345678901234567890123456789012345678901234567890")-gmp_int("123456789012345678901234567890")+1

Dim As gmp_int c
Print
Dim As gmp_int a = "123456789123456789123456789123456789123456789123456789"
Dim As gmp_int b = "987654321987654321987654321987654321987654321987654321"
c = a + b +1
Print c
Print gmp_int_2_base(c, 62); " (base 62)"

Print
c = a
c += b +1
Print c

Print
c = 2
c ^= 128
Print c

Print 
'using the GMP print routine for formatted output
Gmp_printf( !"%70Zd\n", c.num )  ' right adjusted 70 char. width, padded with spaces
Gmp_printf( !"%070Zd\n", c.num ) ' right adjusted 70 char. width, padded with zero's

Print : Print "done"
Sleep
End
Using original GMP code for fdiv and cdiv (floor or ciel the result of a integer division).
fdiv_cdiv.bas

Code: Select all

#Include Once "gmp_int.bi"

' the q functions calculate only the quotient
' the r functions calculate only the remainder
' the qr functions calculates the quotient and remainder

Dim As gmp_int q, r, n1 = 9, n2 = -9, d = 3

print
Print "method      n1     n1\d  n1 mod d   n2     n2\d  n2 mod d" 

For i As Long = 1 To 3

Print "  FB  ";
Print Using "     ###"; n1; n1 \ d; n1 Mod d; n2; n2 \ d; n2 Mod d

' floor, round towards - infinity
mpz_fdiv_q(q.num, n1.num, d.num)
mpz_fdiv_r(r.num, n1.num, d.num)
gmp_printf( !"  fdiv     %#3Zd     %#3Zd     %#3Zd", n1.num, q.num, r.num)
mpz_fdiv_qr(q.num, r.num, n2.num, d.num)
gmp_printf( !"     %#3Zd     %#3Zd     %#3Zd\n", n2.num, q.num, r.num)

' ceil, round towards + infinity
mpz_cdiv_q(q.num, n1.num, d.num)
mpz_cdiv_r(r.num, n1.num, d.num)
gmp_printf( !"  cdiv     %#3Zd     %#3Zd     %#3Zd", n1.num, q.num, r.num)
mpz_cdiv_qr(q.num, r.num, n2.num, d.num)
gmp_printf( !"     %#3Zd     %#3Zd     %#3Zd\n\n", n2.num, q.num, r.num)

n1 += 1
n2 -= 1
    
Next

Print "done"
Sleep
End
Some normal FB programs rewritten with GMP_INT.BI to use big number

Cipolla.bas http://rosettacode.org/wiki/Cipolla%27s_algorithm

Code: Select all

#Include Once "gmp_int.bi"

' Cipolla's algorithm
' type fp2 stuff by srvaldez, thanks 

Type fp2
    x As gmp_int
    y As gmp_int
    Declare Constructor ()
    Declare Constructor ( ByRef lhs As gmp_int, ByRef rhs As gmp_int )
End Type

Constructor fp2 ()
    this.x = 0
    this.y = 0
End Constructor

Constructor fp2 ( ByRef lhs As gmp_int, ByRef rhs As gmp_int )
    this.x = lhs
    this.y = rhs
End Constructor

Function fp2mul(a As fp2, b As fp2, p As gmp_int, w2 As gmp_int) As fp2

    Dim As fp2 answer

    answer.x = (a.x * b.x + a.y * b.y * w2) Mod p
    answer.y = (a.x * b.y + a.y * b.x) Mod p

    Return answer

End Function

Function fp2square(a As fp2, p As gmp_int, w2 As gmp_int) As fp2

    Return fp2mul(a, a, p, w2)

End Function

Function fp2pow(a As fp2, n As gmp_int, p As gmp_int, w2 As gmp_int) As fp2

    If n = 0 Then Return Type (1, 0)
    If n = 1 Then Return a
    If n = 2 Then Return fp2square(a, p, w2)
    If (n And 1) = 0 Then
        Return fp2square(fp2pow(a, n \ 2, p, w2), p , w2)
    Else
        Return fp2mul(a, fp2pow(a, n -1, p, w2), p, w2)
    End If

End Function

' ------=< MAIN >=------

Data "10", "13", "56", "101", "8218", "10007", "8219", "10007"
Data "331575", "1000003", "665165880", "1000000007"
Data "881398088036", "1000000000039"
Data "34035243914635549601583369544560650254325084643201"
Data "100000000000000000000000000000000000000000000000151"

Dim As gmp_int n, p, a, w2, x1, x2
Dim As Long i
Dim As fp2 answer
Dim As String n_str, p_str

For i = 1 To 8

    Read n_str, p_str
    n = n_str : p = p_str
    Print
    Print "Find solution for n = "; n; " and p = "; p

    If p = 2 OrElse isprime(p) = 0 Then
        Print "No solution, p is not a odd prime"
        Continue For
    End If

    ' p is checked and is a odd prime
    If legendre(n, p) <> 1 Then
        Print n; " is not a square in F"; Str(p)
        Continue For
    End If

    a = 1
    Do
        Do
            a += 1
            w2 = a * a - n + p
        Loop Until legendre(w2, p) = -1

        answer = Type (a, 1)
        'answer = fp2pow(answer, (p +1) \ 2, p, w2)
        answer = fp2pow(answer, (p +1) Shr 1, p, w2)
        If answer.y <> 0 Then Continue Do

        x1 = answer.x : x2 = p - x1
        If powm(x1, 2, p) = n AndAlso powm(x2, 2, p) = n Then
            Print "Solution found: x1 = "; x1; ", "; "x2 = "; x2
            Exit Do
        End If
    Loop            ' loop until solution is found

Next

Print : Print "done"
Sleep
End
Tonelli-Shanks.bas http://rosettacode.org/wiki/Tonelli-Shanks_algorithm

Code: Select all

#Include Once "gmp_int.bi"

' Tonelli-Shanks algorithm

Data "10", "13", "56", "101", "1030", "10009", "1032", "10009", "44402", "100049"
Data "665820697", "1000000009", "881398088036", "1000000000039"
Data "41660815127637347468140745042827704103445750172002"
Data "100000000000000000000000000000000000000000000000577"         ' 10^50 + 577


Dim As gmp_int b, c, i, m, n, p, q, r, s, t, z
Dim As Long k
Dim As String n_str, p_str

For k = 1 To 8
    Read n_str, p_str
    n = n_str : p = p_str
    Print "Find solution for n = "; n; " and p = ";p

    If p = 2 OrElse Isprime(p, 15) = 0 Then
        Print p;" is not a odd prime"
        Print
        Continue For
    End If

    If legendre(n, p) = -1 Then
        Print n;" is not a quadratic residue"
        Print
        Continue For
    End If

    s = 0 : q = p -1
    Do
        s += 1
        q \= 2
    Loop Until (q And 1) = 1

    If s = 1 And (p Mod 4) = 3 Then
        ' r = powm(n, ((p +1) \ 4), p)
        r = powm(n, ((p +1) Shr 2), p)
        Print "Solution found: "; r; " and "; p - r
        Print
        Continue For
    End If

    z = 1
    Do
        z += 1
    Loop Until legendre(z, p) = -1
    c = powm(z, q, p)
    ' r = powm(n, ((q +1) \ 2), p)
    r = powm(n, ((q +1) shr 1), p)
    t = powm(n, q, p)
    m = s

    Do

        i = 0
        If (t Mod p) = 1 Then
            Print "Solution found: "; r; " and "; p - r
            Print
            Continue For
        End If

        Do
            i += 1
            If i >= m Then Continue For
        Loop Until powm(t, (2 ^ i), p) = 1

        b = powm(c, (2 ^ (m - i -1)), p)
        c = powm(b, 2, p)
        r = (r * b) Mod p
        t = (t * c) Mod p
        m = i

    Loop

Next

Print : Print "done"
Sleep
End
Rank_of_a_permutation.bas http://rosettacode.org/wiki/Permutation ... ermutation

Code: Select all

#Include Once "gmp_int.bi"

' Myrvold and Ruskey

Sub unrank1(n As ULong, r As gmp_int , pi() As UByte)

    If n > 0 Then
        Swap pi(n -1), pi(r Mod n)
        unrank1(n -1, (r \ n), pi())
    End If

End Sub

Function rank1(n As ULong, pi() As UByte, pi_inv() As UByte) As gmp_int

    If n = 1 Then Return 0

    Dim As UByte s = pi(n -1)

    Swap pi(n -1), pi(pi_inv(n -1))
    Swap pi_inv(s), pi_inv(n -1)

    Return (s + n * rank1(n -1, pi(), pi_inv()))

End Function

Sub unrank2(n As ULong, r As gmp_int, pi() As UByte)

    If n > 0 Then
        Dim As gmp_int f1 = fac(n - 1)
        Dim As gmp_int s = r \ f1
        Swap pi(n -1), pi(s)
        unrank2(n -1, r - s * f1, pi())
    End If

End Sub

Function rank2(n As ULong, pi() As UByte, pi_inv() As UByte) As gmp_int

    If n = 1 Then Return 0
    Dim As UByte s = pi(n -1)
    Swap pi(n -1), pi(pi_inv(n -1))
    Swap pi_inv(s), pi_inv(n -1)
    Return (s * fac(n -1) + rank2(n -1, pi(), pi_inv()))

End Function

' ------=< MAIN >=------

Dim As gmp_int i1, j, n1
Dim As ULong i, n
Dim As UByte pi(), pi_inv()
Dim As String frmt1, frmt2
Randomize Timer

n = 3 : n1 = fac(n)
ReDim pi(n -1), pi_inv(n - 1)
frmt1 = " ###"
frmt2 = "##"

Print "Rank:     unrank1       rank1"

For i = 0 To n1 -1
    For j = 0 To n -1
        pi(j) = j
    Next
    Print Using frmt1 & ": --> "; i;
    unrank1(n, i, pi())
    For j = 0 To n -1
        Print Using frmt2; pi(j);
        pi_inv(pi(j))= j
    Next

    Print Using "  -->" & frmt1; rank1(n, pi(), pi_inv())

Next

n = 12 : n1 = fac(n)
ReDim pi(n -1), pi_inv(n - 1)
frmt1 = "###########"
frmt2 = "###"
Print : Print "4 random samples of permutations: 12 objects"
Print "  Rank:                     unrank1                         rank1"

For i = 1 To 4
i1 = gmp_rnd(n1)
For j = 0 To n -1 : pi(j) = j : Next
Print Using frmt1 & ": --> "; i1; : unrank1(n, i1, pi())
For j = 0 To n -1 : Print Using frmt2; pi(j);
    pi_inv(pi(j))= j : Next
    Print Using "  -->" & frmt1; rank1(n, pi(), pi_inv())
Next

Print : Print String(69,"-") : Print
Print "Rank:     unrank2       rank2"

n = 3 : n1 = fac(n)
ReDim pi(n -1), pi_inv(n - 1)
frmt1 = " ###"
frmt2 = "##"

For i = 0 To n1 -1
    For j = 0 To n -1
        pi(j) = j
    Next
    Print Using  frmt1 & ": --> "; i;
    unrank2(n, i, pi())
    For j = 0 To n -1
        Print Using frmt2; pi(j);
        pi_inv(pi(j))= j
    Next

    Print Using "  -->" & frmt1; rank2(n, pi(), pi_inv())

Next

n = 12 : n1 = fac(n)
ReDim pi(n -1), pi_inv(n - 1)
frmt1 = "###########"
frmt2 = "###"
Print : Print "4 random samples of permutations: 12 objects"
Print "  Rank:                     unrank2                         rank2"

For i = 1 To 4
i1 = gmp_rnd(n1)
For j = 0 To n -1 : pi(j) = j : Next
Print Using frmt1 & ": --> "; i1; : unrank2(n, i1, pi())
For j = 0 To n -1 : Print Using frmt2; pi(j);
    pi_inv(pi(j))= j : Next
    Print Using "  -->" & frmt1; rank2(n, pi(), pi_inv())
Next

n = 144 : n1 = fac(n) ' n1 is 251 digits
ReDim pi(n -1), pi_inv(n - 1)
frmt2 = " ###"
Print : Print "10 random samples of permutations: 144 objects"

For i = 1 To 10
    i1 = gmp_rnd(n1)
    For j = 0 To n -1 : pi(j) = j : Next
        unrank2(n, i1, pi())
        For j = 0 To n -1 : Print Using frmt2; pi(j);
        Next
        Print
Next

Print : Print "done"
Sleep
End
Last edited by frisian on May 22, 2017 21:09, edited 2 times in total.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: Big Number Wrapper (GMP_INT.BI)

Post by srvaldez »

hello frisian, you need a constructor for the fp2 type.

Code: Select all

type fp2
	x as gmp_int
	y as gmp_int
	declare constructor()
	declare Constructor ( Byref lhs As gmp_int, Byref rhs As gmp_int )
end type

constructor fp2 ()
	this.x=0
	this.y=0
End Constructor

Constructor fp2 ( Byref lhs As gmp_int, Byref rhs As gmp_int )
	this.x=lhs
	this.y=rhs
End Constructor
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

Re: Big Number Wrapper (GMP_INT.BI)

Post by frisian »

srvaldez wrote:hello frisian, you need a constructor for the fp2 type.
Thanks for the code example, the only thing I knew that it needed a constructor through the error messages I got from my try's. But I had no idea how and where to put it.
Makoto WATANABE
Posts: 231
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

Re: Big Number Wrapper (GMP_INT.BI)

Post by Makoto WATANABE »

Dear frisian

Thank you for extending the "Big Number Wrapper" and adding a lot of useful functions.
And thank you for publishing example programs showing usage.
I will translate them into Japanese and introduce them to Japanese people.
Post Reply