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
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
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
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
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
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