Code: Select all
' version 2.0 25 November 2013
'=============================
' by Stephan Brunker (stephanbrunker at web punkt de)
'=============================
' The easy way to use - analog to the "hello world" the 1+1:
' #include "big_integer.bas"
' Dim as bigint a,b,c
' a = "1"
' b = "1"
' c = a + b
' print c
' a and b can be so long as you want - up to 1.7E10 bits
' all the operators: + - * \ mod ^ += -= *= \= mod= are overloaded
' also all the comparators: = < > <> >= <=
' for additional functions look in the code, it should be
' well commented to see for what the functions are for
' you can also input and output in various formats, see section
' "Conversion Functions"
'===============================================================
'known bug: Please dont input "integer" or "uinteger".
'Integer<32> is internally handled as long, and Integer<64> as longint.
'But simply Integer can be both, and I couldn't get the constructors
'to work properly - you can try the commented sections below.
'same with undefined assignments:
'(bigint) = (not bigint) {operator} (not bigint) doesn't work too.
'please use a cuint<32> in the syntax, then the conversion is clear.
'maybe someone else can fix that ..!
'===============================================================
' Changelog
'-----------
'-removed comments @shift_right
'-bugfix constructor and operator uinteger, added ulong
'-conversion functions hex_bigint, uhex_bigint,
' bin_bigint, bigint_bin; ubin_bigint, bigint_ubin
'-added mod_power, removed Divide and Remainder
'-integer and uinteger defined for 32-bit/64-bit
'-bug exploitation in negate:
' culngint(not(uinteger)) = cuint(not(uinteger) fixed
'-rewritten divmod(30x faster now)
'-rewritten comparation - 20x faster
'-general speed optimisation: Variables are faster than pointers,
' pointers faster than string operations:
' Square is 10% faster now, Multiply 5% faster
'-using the internal functions is slightly faster than the overloaded
' operators because of copying the input twice
'
' version 1.2 3 July 2011
'================================================================
' Big_Integer_2.bas Arbitrary Precision, Two's Complement Integer.
' by Richard @ freebasic.net/forum
'================================================================
' supported functions are;
' (short constants -1, 0, 1, 2 ) Bigint_m1, Bigint_0, Bigint_1, Bigint_2
' Msbit, Prune, Neg, Absolute, pack, unpack, Addition, Subtract
' Multiply, Square, Mul2, Div2, Power, Factorial
' Shift_Left, Shift_Right, Complement, Bit_value, Bit_set, Bit_reset
' Bit_And, Bit_Or, Bit_Xor, Bit_Eqv, Bit_Imp, Show_Bin, Show_Hex, Show_Oct
' is_ZERO, is_NZ, is_POS, is_NEG, is_EQ, is_NEQ
' is_LT, is_GT, is_LT_EQ, is_GT_EQ, BigInt_Sgn, DivMod, Divide, Remainder
' CInt, CLongInt, CDbl
'----------------------------------------------------------------
' This package only handles integers encoded in a two's complement format.
' The first bit in the two's complement format is always the sign which
' takes the value of 1 = negative, 0 = positive or zero. Byte examples are;
' +127 = 0111 1111 most positive
' +8 = 0000 1000
' +7 = 0000 0111
' +4 = 0000 0100
' +2 = 0000 0010
' +1 = 0000 0001
' 0 = 0000 0000 zero
' -1 = 1111 1111
' -2 = 1111 1110
' -4 = 1111 1100
' -7 = 1111 1001
' -8 = 1111 1000
' -127 = 1000 0001
' -128 = 1000 0000 most negative
'----------------------------------------------------------------
' Each successive 32 bits are packed into a 4 byte block.
' Each block is stored in 4 successive bytes of a string.
'----------------------------------------------------------------
' Strings in FreeBASIC are indexed and printed from left to right. Big
' integers appear to be stored in strings backwards which can be confusing.
' The Left side of the string is the Right side of the number and vice versa.
' Beware: Shift_Left moves a string to the right, Shift_Right moves it left.
'----------------------------------------------------------------
' The first byte in the string is the least significant byte of the number.
' The last block in the string is the most significant block of the number.
' String s indexing has s[0] as the LS byte and s[len(s)-1] as the MS byte.
' These big integer strings are always multiples of 4 bytes long.
' The msb of a number is sign extended to the MS bit of the MS block.
'----------------------------------------------------------------
' A number is always stored in a way that correctly represents that number.
' Where an operation would overflow into the MSB, a sign block is
' appended to the number so as to prevent overflow or sign change.
' Unnecessary leading zero or all ones blocks are removed by prune.
'----------------------------------------------------------------
' String pointers may change if a string is created or any length is changed.
'================================================================
Type bigint ' a little endian, two's complement binary number
s As String ' packed into a string, in blocks four bytes long
' constructor
Declare Constructor ' default constructor
Declare Constructor (Byref As String)
Declare Constructor (Byref As Long)
Declare Constructor (Byref As Ulong)
' Declare Constructor (Byref As Integer)
' Declare Constructor (Byref As Uinteger)
Declare Constructor (Byref As Longint)
Declare Constructor (Byref As Ulongint)
' Let
Declare Operator Let(Byref rhs As String)
Declare Operator Let(Byval rhs As Byte)
Declare Operator Let(Byval rhs As Ubyte)
Declare Operator Let(Byval rhs As Short)
Declare Operator Let(Byval rhs As Ushort)
Declare Operator Let(Byval rhs As Long)
Declare Operator Let(Byval rhs As Ulong)
Declare Operator Let(Byval rhs As Longint)
Declare Operator Let(Byval rhs As Ulongint)
Declare Operator Let(Byval rhs As Double)
' cast
Declare Operator Cast() As String
Declare Operator Cast() As Integer ' CInt
Declare Operator Cast() As LongInt ' CLongInt
Declare Operator Cast() As Double ' Cdbl
' implicit step versions
Declare Operator For ()
Declare Operator Step()
Declare Operator Next(Byref end_cond As bigint) As Integer
' explicit step versions
Declare Operator For (Byref step_var As bigint)
Declare Operator Step(Byref step_var As bigint)
Declare Operator Next(Byref end_cond As bigint, Byref step_var As bigint ) As Integer
' operate and assign
Declare Operator += (Byref rhs As bigint)
Declare Operator -= (Byref rhs As bigint)
Declare Operator *= (Byref rhs As bigint)
Declare Operator \= (Byref rhs As bigint)
Declare Operator mod= (Byref rhs As bigint)
End Type
'================================================================
' only the necessary declarations are here
Declare Function pack(Byref As String) As bigint
Declare Function unpack(Byref As bigint) As String
Declare Function negate(Byref As bigint) As bigint
Declare Function bigint_compare(byref a as bigint, byref b as bigint) as integer
Declare Operator < (Byref As bigint, Byref As bigint) As Integer
Declare Operator > (Byref As bigint, Byref As bigint) As Integer
Declare Operator = (Byref As bigint, Byref As bigint) As Integer
Declare Operator <> (Byref a As bigint, Byref b As bigint) As Integer
Declare Operator <= (Byref a As bigint, Byref b As bigint) As Integer
Declare Operator >= (Byref a As bigint, Byref b As bigint) As Integer
Declare Operator + (Byref x As bigint) As bigint
Declare Operator - (Byref x As bigint) As bigint
Declare Operator + (Byref x As bigint, Byref y As bigint) As bigint
Declare Operator - (Byref x As bigint, Byref y As bigint) As bigint
Declare Operator * (Byref x As bigint, Byref y As bigint) As bigint
Declare Operator \ (Byref x As bigint, Byref y As bigint) As bigint
Declare Operator mod (Byref x As bigint, Byref y As bigint) as bigint
Declare Operator ^ (Byref x As bigint, Byref y As Longint) As bigint
Declare Function is_LT(Byref As bigint, Byref As bigint) As Integer
Declare Function show_hex(Byref s As bigint) As String
Declare Function msbit(Byref a As bigint) As Integer
Declare Function BigInt_2_Double ( Byref a As BigInt ) As Double
Declare Function Bit_Value(Byref v As bigint, Byval b As Ulongint) As Integer
'================================================================
Constructor bigint ' default constructor
this.s = chr(0,0,0,0) ' zero
End Constructor
Constructor bigint (Byref rhs As Long)
this = rhs
End Constructor
Constructor bigint (Byref rhs As ULong)
this = rhs
End Constructor
'Constructor bigint (Byref rhs As Integer)
'this = rhs
'End Constructor
'
'Constructor bigint (Byref rhs As Uinteger)
'this = rhs
'End Constructor
Constructor bigint (Byref rhs As Longint)
this = rhs
End Constructor
Constructor bigint (Byref rhs As ULongint)
this = rhs
End Constructor
Constructor bigint (Byref rhs As String)
this = rhs
End Constructor
'================================================================
' assignment operator
Operator bigint.let (Byref rhs As String)
this = pack(rhs)
End Operator
'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Byte)
If (128 And rhs) Then
this.s = Chr(rhs,-1,-1,-1)
Else
this.s = Chr(rhs,0,0,0)
End If
End Operator
'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Ubyte)
this.s = Chr(rhs,0,0,0)
End Operator
'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Short)
If (32768 And rhs) Then
this.s = Chr(Lobyte(rhs), Hibyte(rhs), -1, -1 )
Else
this.s = Chr(Lobyte(rhs), Hibyte(rhs), 0, 0 )
End If
End Operator
'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Ushort)
this.s = Chr(Lobyte(rhs), Hibyte(rhs), 0, 0 )
End Operator
'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Long)
this.s = Chr(0,0,0,0)
Dim As Integer<32> Ptr bip = Cptr(Integer<32> Ptr, Strptr(this.s))
Dim As Integer<32> Ptr ip = Cptr(Integer<32> Ptr, @rhs)
*bip = *ip
End Operator
'----------------------------------------------------------------
Operator bigint.let (Byval rhs As ULong)
this.s = Chr(0,0,0,0)
Dim As Uinteger<32> Ptr bip = Cptr(Uinteger<32> Ptr, Strptr(this.s))
Dim As Uinteger<32> Ptr uip = Cptr(Uinteger<32> Ptr, @rhs)
*bip = *uip
If (128 And this.s[3]) Then this.s &= Chr(0,0,0,0) ' make it positive
End Operator
''----------------------------------------------------------------
'Operator bigint.let (Byval rhs As Integer)
'this.s = Chr(0,0,0,0)
'if rhs > 2147483467 then print "64-bit Integers not yet supported" : sleep : end
'Dim As Integer<32> Ptr bip = Cptr(Integer<32> Ptr, Strptr(this.s))
'Dim As Integer<32> Ptr ip = Cptr(Integer<32> Ptr, @rhs)
'*bip = *ip
'End Operator
'
''----------------------------------------------------------------
'Operator bigint.let (Byval rhs As UInteger)
'this.s = Chr(0,0,0,0)
'if rhs > 4294967295 then print "64-bit Integers not yet supported" : sleep : end
'Dim As Uinteger<32> Ptr bip = Cptr(Uinteger<32> Ptr, Strptr(this.s))
'Dim As Uinteger<32> Ptr uip = Cptr(Uinteger<32> Ptr, @rhs)
'*bip = *uip
'If (128 And this.s[3]) Then this.s &= Chr(0,0,0,0) ' make it positive
'End Operator
'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Longint)
this.s = String(8, Chr(0) )
Dim As Longint Ptr bip = Cptr(Longint Ptr, Strptr(this.s))
Dim As Longint Ptr lip = Cptr(Longint Ptr, @rhs)
*bip = *lip
End Operator
'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Ulongint)
this.s = String(8, Chr(0) )
Dim As Ulongint Ptr bip = Cptr(Ulongint Ptr, Strptr(this.s))
Dim As Ulongint Ptr ulip = Cptr(Ulongint Ptr, @rhs)
*bip = *ulip
If (128 And this.s[7]) Then this.s &= Chr(0,0,0,0) ' make it positive
End Operator
'----------------------------------------------------------------
Operator bigint.let (Byval d As Double)
Const As Ulongint implicit_bit = 2^52 ' 4503599627370496
Const As Ulongint mant_mask = implicit_bit - 1 ' 4503599627370495
Dim As Ulongint u, mant
Dim As Integer negative, sign, expo
Dim As bigint a
'----------------------------------------------------
If d < 0 Then negative = -1 ' remember sign
d = Int(Abs(d) + 0.5) ' rectify and round to closest integer
'----------------------------------------------------
u = *(Cptr(Ulongint Ptr, @d)) ' copy Double into a Ulongint frame
'sign = (u Shr 63) ' extract the sign bit
expo = (u Shr 52) And 2047 ' the 11 bit exponent
mant = (u And mant_mask ) ' 52 mantissa bits
'----------------------------------------------------
If expo = 0 Then ' the double has zero value or is de-normalized
this.s = Chr(0,0,0,0) ' de-normalized is very very small
Else
mant = mant + implicit_bit ' insert the missing implicit bit
expo = expo - 1023 ' remove the excess exponent
If expo < 53 Then
mant = mant Shr (52 - expo) ' reduce it to integer
a.s = String(8, Chr(0))
*(Cptr(Ulongint Ptr, Strptr(a.s))) = mant ' make bigint
If negative Then a = negate(a)
Else
Print "WARNING: Double argument was unreliable."
Stop
End If
this = a
End If
End Operator
'================================================================
' create a string for printing a bigint
Operator bigint.cast() As String
Dim As String t
t = unpack(this)
Return t
End Operator
Operator bigint.Cast() As Double ' Cdbl
Return BigInt_2_Double(this)
End Operator
'================================================================
' common small constants
'Dim Shared As bigint bigint_m1, bigint_0, bigint_1, bigint_2
'bigint_m1.s = Chr(255,255,255,255)' minus 1
'bigint_0.s = Chr(0,0,0,0) ' zero
'bigint_1.s = Chr(1,0,0,0) ' one
'bigint_2.s = Chr(2,0,0,0) ' two
Dim Shared As bigint bigint_m1=-1
Dim Shared As bigint bigint_0=0
Dim Shared As bigint bigint_1=1
Dim Shared As bigint bigint_2=2
'-----------------------------------------------------------------------
' remove unnecessary leading blocks, prune time easily earns it's keep
Sub prune(Byref a As bigint)
a.s = Left(a.s, (((msbit(a) + 1) Shr 5) + 1 ) Shl 2)
End Sub
' Times and space are improved through the sensible use of prune.
' Addition of opposite signs or subtraction can cancel many blocks.
' A redundant block can appear during multiplication. Square or div2.
' Mul2, Complement, Negate and Absolute do not generate unnecessary blocks.
' Power is pruned internally by the prune in multiply and square.
'================================================================
' ARITHMETRIC OPERATIONS
'================================================================
' Negate the twos complement binary number in a bigint
Function negate(Byref a As bigint) As bigint
Dim As bigint s = a
Dim As Integer blocks = Len(s.s) Shr 2
Dim As Ulongint sum
Dim As Uinteger<32> carry
Dim As Uinteger<32> Ptr ps
ps = Cast( Uinteger<32> Ptr, Strptr(s.s))' the Uinteger data in bigint
carry = 1 ' set carry
Do ' slow ahead until clear of the carry
sum = cuint(Not *ps) + carry
*ps = sum
carry = sum shr 32
ps += 1
blocks -= 1
Loop Until (carry = 0) Or (blocks = 0)
If blocks > 0 Then
Do ' no carry, so full speed ahead
*ps = Not *ps
ps +=1
blocks -= 1
Loop Until blocks = 0
End If
' Negating the most negative integer is a problem because carry propagation
' flips the sign which should have become positive. But negation of zero
' does not change the sign either, so we need to differentiate between zero
' and one by quickly examining the final carry bit from the two's complement.
If carry = 0 Then ' carry was not generated by the most negative number
If (128 And a.s[Len(a.s)-1]) = (128 And s.s[Len(s.s)-1]) Then s.s = s.s + bigint_0.s
End If ' this prevents a negated zero being extended by an extra zero block
Return s
End Function
' digits time on 2 GHz Pentium
' 10 2.2 usec
' 100 2.5 usec
' 1k 4.2 usec
' 10k 16. usec
' 100k 200. usec
'================================================================
' absolute value is positive
Function absolute(Byref a As bigint) As bigint
Dim As bigint b = a
If 128 And b.s[Len(b.s)-1] Then b = negate(b)
Return b
End Function
'================================================================
' addition
Function addition(Byref aa As bigint, Byref bb As bigint) As bigint
Dim As bigint a = aa, b = bb
Dim As Integer blocks, i, j, lena, lenb, sa, sb, delta
'------------------------------------------------------------
' test to see if the two most significant digits differ which
lena = Len(a.s) - 1 ' might change the sign without carry
If a.s[lena] And 128 Then sa = 255 ' sign as a byte
i = a.s[lena] And 192 ' if MSBs differ then extend the sign
If (i = 64) Or (i = 128) Then a.s = a.s + String(4, Chr(sa) )
'------------------------------------------------------------
lenb = Len(b.s) - 1
If b.s[lenb] And 128 Then sb = 255
i = b.s[lenb] And 192
If (i = 64) Or (i = 128) Then b.s = b.s + String(4, Chr(sb) )
'------------------------------------------------------------
' align the two bigints to be added
delta = Len(a.s) - Len(b.s) 'new values
If delta <> 0 Then ' sign extend the shorter
If delta > 0 Then
' a = a
If b.s[Len(b.s)-1] And 128 Then i = 255 Else i = 0
b.s = b.s + String(delta, Chr(i) ) ' extend b
Else
If aa.s[Len(aa.s)-1] And 128 Then i = 255 Else i = 0
a.s = a.s + String(-delta, Chr(i) ) ' extend a
' b = b
End If
End If ' a and b are now the same length
'------------------------------------------------------------
' accumulate b into a
blocks = Len(a.s) Shr 2
Dim As Ulongint sum = 0 ' clear carry
Dim As Uinteger<32> carry
Dim As Uinteger<32> Ptr pa, pb
pa = Cast(Uinteger<32> Ptr, Strptr(a.s) )
pb = Cast(Uinteger<32> Ptr, Strptr(b.s) )
For i = 0 To blocks-1
sum = Clngint(*(pa+i)) + Clngint(*(pb+i)) + carry
*(pa+i) = sum
carry = sum shr 32
Next i
prune(a)
Return a
End Function
' 100 digits in 4.2 usec
'===============================================================
' subtract
Function subtract(Byref aa As bigint, Byref bb As bigint) As bigint
Dim As bigint cc = Negate(bb)
cc = addition(aa, cc)
Return cc
End Function
'=======================================================================
' square a number, approaches twice the speed of multiply for big numbers
Function square(Byref aa As bigint) As bigint
'------------------------------------------------------------
Dim As bigint a = aa, c
If (128 And a.s[Len(a.s)-1]) Then a = negate(a)
'------------------------------------------------------------
' find the dimension of the problem
Dim As Integer i, j, asize = Len(a.s) ' number of bytes in a
c.s = String(2 * asize, Chr(0) ) ' initialise accumulator
asize = (asize Shr 2) - 1 ' highest block number in a
'------------------------------------------------------------
' pointers into all the bigints
Dim As Uinteger<32> Ptr pa, pc
pa = Cast(Uinteger<32> Ptr, Strptr(a.s) ) ' the base addresses of bigints
pc = Cast(Uinteger<32> Ptr, Strptr(c.s) )
Dim As Ulongint product ' bottom half is data, top half will be carry
Dim As Uinteger<32> carry, sum
'------------------------------------------------------------
' multiply one triangular corner only
For i = 1 To asize
' pa starts at 1 not zero because 0,0 is on the diagonal
' the second element in a starts at zero
' pc is the accumulator ic = icz + i + j
carry = 0 ' clear carry
For j = 0 To i - 1
product = Culngint(*(pa+i)) * Culngint(*(pa+j)) + *(pc+j+i) + carry
*(pc+j+i) = product
carry = product shr 32
Next j
*(pc+j+i) = carry ' line of blocks gets one longer each loop
Next i
'------------------------------------------------------------
' double the triangle, cannot do it at same time as squaring diagonal
' because it can cause the product to overflow
carry = 0 ' clear carry
For i = 0 To (2 * asize) + 1
sum = *(pc+i)
product = Culngint(sum) + sum + carry
*(pc+i) = product
carry = product shr 32
Next i
'------------------------------------------------------------
' square and accumulate the diagonal elements
carry = 0 ' clear carry
For i = 0 To asize
' square the diagonal entry, while propagating carry
sum = *(pa+i)
product = Culngint(sum) * Culngint(sum) + *(pc+i+i) + carry
*(pc+i+i) = product
carry = product shr 32
' propagate carry through the following odd block of C
product = Culngint(*(pc+i+i+1)) + carry
*(pc+i+i+1) = product
carry = product shr 32
Next i
'------------------------------------------------------------
If 128 And c.s[Len(c.s)-1] Then c.s = c.s + bigint_0.s ' sign is positive
prune(c)
Return c
End Function
' digits multiply(x,x) square(x)
' 10 3.9 usec 3.9 usec
' 100 11.5 usec 9.0 usec
' 1k 640 usec 340 usec
' 10k 64 msec 33 msec
' 100k 6.35 sec 3.21 sec
' 1M 320 sec
'================================================================
' multiply
Function multiply(Byref aa As bigint, Byref bb As bigint) As bigint
' If @aa = @bb Then ' detect if square, half time if very big
' Return square(aa)
' Else
'------------------------------------------------------------
' sort out the signs and rectify the inputs
Dim As bigint a = aa, b = bb, c
Dim As Integer sign_a, sign_b, sign_c
sign_a = a.s[Len(a.s)-1] And 128
sign_b = b.s[Len(b.s)-1] And 128
sign_c = sign_a Xor sign_b
If sign_a Then a = negate(a)
If sign_b Then b = negate(b)
'------------------------------------------------------------
' find the dimensions of the problem
Dim As Integer i, j, asize, bsize
asize = Len(a.s) ' number of bytes in a
bsize = Len(b.s) ' number of bytes in b
c.s = String(asize + bsize, Chr(0)) ' initialise accumulator
asize = asize Shr 2 - 1 ' number of blocks in a
bsize = bsize Shr 2 - 1 ' Shr 2 is faster than \4
'------------------------------------------------------------
' pointers into all the bigints
Dim As Uinteger<32> Ptr ia, ib, ic
ia = Cast(Uinteger<32> Ptr, Strptr(a.s) )
ib = Cast(Uinteger<32> Ptr, Strptr(b.s) )
ic = Cast(Uinteger<32> Ptr, Strptr(c.s) )
Dim As Ulongint product
Dim As Uinteger<32> carry
'------------------------------------------------------------
For i = 0 To asize
carry = 0 ' clear carry
For j = 0 To bsize
product = Culngint(*(ia+i)) * Culngint(*(ib+j)) + *(ic+i+j) + carry
*(ic+i+j) = product
carry = product shr 32
Next j
*(ic+i+j) = carry
Next i
'------------------------------------------------------------
If sign_c = 128 Then c = negate(c)
prune(c)
Return c
' End If
End Function
' digits time
' 10 3.9 usec
' 100 11.5 usec
' 1k 640. usec
' 10k 64. msec
' 100k 6.35 sec
'=======================================================================
' shift up one bit, low towards high
Function mul2(Byref a As bigint) As bigint
Dim As Integer i, sign, sum, carry = 0
Dim As bigint b = a
sign = b.s[Len(b.s) - 1] And 128 ' snag the msb of highest byte
For i = 0 To Len(b.s) - 1
sum = b.s[i] + b.s[i] + carry
carry = Hibyte(sum)
b.s[i] = Lobyte(sum)
Next i
If ( b.s[Len(b.s) - 1] And 128 ) <> sign Then
carry = carry * 255
b.s = b.s + Chr(carry,carry,carry,carry) ' sign extend four bytes
End If
Return b
End Function
'=======================================================================
' shift down one bit, high towards low
Function div2(Byref a As bigint) As bigint
Dim As Integer i, carry = 0
Dim As bigint b = a
For i = 0 To Len(a.s)-2 ' all except the top byte of four
b.s[i] = ( b.s[i] Shr 1 ) + (128 * (b.s[i+1] And 1))
Next i
i = Len(b.s) - 1
b.s[i] = b.s[i] Shr 1
b.s[i] = b.s[i] Or (2 * ( b.s[i] And 64)) ' sign extend the msb
prune(b)
Return b
End Function
'=======================================================================
' raise a number to a positive power
Function power(Byref x As bigint, Byval n As Longint) As bigint
If n < 0 Then
Print "Cannot raise a big integer to a negative power."
Stop
End If
Dim As Integer i = 64
Do ' find first set bit
i = i - 1
Loop Until Bit(n, i) Or (i = 0)
i = i + 1
Dim As bigint pwr
pwr = bigint_1 ' one
Do
i = i - 1
pwr = square(pwr) ' this was a multiply but square is faster
If Bit(n, i) Then pwr = multiply(pwr, x)
Loop Until i = 0
Return pwr ' pwr was pruned by square and by multiply
End Function
' times with multiply square prune&square
' 2 ^ 3 25 usec 13.5 usec 27.4 usec
' 2 ^ 35 300 usec 147 usec 46.4 usec
'123456789^127 1.43 msec 780 usec 260 usec
' 10 ^ 1000 82 msec 41 msec 186 usec
' 10 ^ 10000 14 sec 7 sec 10.9 msec
' 10 ^ 100000 1.07 sec
'=======================================================================
' integer divide, a \ b, return div and mod
Sub divmod(_
Byref aa As bigint, Byref bb As bigint,_
Byref q As bigint, Byref r As bigint)
Dim As Integer lenaa = Len(aa.s), lenbb = Len(bb.s)
'------------------------------------------------------------
'Test if Division by zero
Dim as uinteger<32> ptr pbb
pbb = Cast(Uinteger<32> Ptr, Strptr(bb.s))
If *pbb = 0 and lenbb = 4 then
Print " Division by zero. "
sleep: End
End If
'------------------------------------------------------------
'Test if Longint Division possible
If (lenaa <= 8) And (lenbb <= 8) Then ' arguments are one or two blocks
q.s = Chr(0,0,0,0,0,0,0,0) ' use 64 bit long integers
r.s = Chr(0,0,0,0,0,0,0,0) ' allocate space for results
If lenaa = 4 Then aa.s = aa.s + String(4, Bit(aa.s[lenaa - 1], 7))
If lenbb = 4 Then bb.s = bb.s + String(4, Bit(bb.s[lenbb - 1], 7))
Dim As Longint Ptr pa1, pb1, pq1, pr1
pa1 = Cast(Longint Ptr, Strptr(aa.s)) '2 block bigint -> Longint
pb1 = Cast(Longint Ptr, Strptr(bb.s))
pq1 = Cast(Longint Ptr, Strptr(q.s))
pr1 = Cast(Longint Ptr, Strptr(r.s))
If *pb1 = 0 then print " Divison by zero. " : sleep : end
*pq1 = *pa1 \ *pb1
*pr1 = *pa1 Mod *pb1
Exit Sub
End If
'------------------------------------------------------------
'Test if divisor is bigger than dividend
If bb > aa then
q = bigint_0
r = aa
exit sub
End if
'------------------------------------------------------------
'Read Signs
Dim As Integer sa, sb, sq
sa = 128 And aa.s[lenaa-1]
sb = 128 And bb.s[lenbb-1]
sq = sa Xor sb ' sign of the result
'---------------------------------------------------------------------
' Setup variables and pointers
' r=dividend and remainder
' b=divisor
' q=quotient
' sum=interim result for subtraction
Dim As bigint b = bb
If sb Then b = negate(b)
r = aa
If sa Then r = negate(r)
Dim As Ushort Ptr pr, pb, pq
Dim As Uinteger<32> sum, high48b, offset=&hFFFF0000
Dim as Uinteger rblocks,bblocks,blockshift,rounds,blocks,i,j
Dim as Uinteger lenr=len(r.s), lenb=len(b.s)
Dim as Ushort qi,carry
Dim as Integer substract
Dim As Ulongint high48r=0
pr = Cast(Ushort Ptr, Strptr(r.s))
pb = Cast(Ushort Ptr, Strptr(b.s))
rblocks=(lenr \ 2)-1
bblocks=(lenb \ 2)-1
'------------------------------------------------------------
'convert to 16bit blocklength for effective testdivison
'append exactly two zero blocks and ignore the second leading zeroblock
if *(pr+rblocks) <> 0 then
r.s = r.s & chr(0,0,0,0)
rblocks += 1
elseif *(pr+rblocks)=0 and *(pr+rblocks-1)=0 then
rblocks -= 1
end if
if *(pb+bblocks) <> 0 then
b.s = b.s & chr(0,0,0,0)
bblocks += 1
elseif *(pb+bblocks)=0 and *(pb+bblocks-1)=0 then
bblocks -= 1
end if
'new pointer after changing string
pr = Cast(Ushort Ptr, Strptr(r.s))
pb = Cast(Ushort Ptr, Strptr(b.s))
blockshift=rblocks-bblocks
rounds=blockshift + 1
'------------------------------------------------------------
' setup quotient
if bit(rounds,0)=0 then
q.s = string(rounds*2,0) 'quotient
pq = Cast(Ushort Ptr, Strptr(q.s))
pq += (len(q.s) \ 2) - 1 'start at msb
else
q.s = string((rounds+1)*2,0)
pq = Cast(Ushort Ptr, Strptr(q.s))
pq += (len(q.s) \ 2) - 2 'start at msb
end if
'------------------------------------------------------------
' start calculation
'most significant 48bits of divisor are constant
if bblocks = 1 then
high48b = Cuint(*(pb+bblocks-1) shl 16) 'second block is empty
else
high48b = Cuint(*(pb+bblocks-1) shl 16) + *(pb+bblocks-2) + 1
end if
for blocks = 1 to rounds
'testdivision: calculate highest three blocks and divide
high48r = Culngint(*(pr+rblocks)) shl 32 + _
Cuint(*(pr+rblocks-1) shl 16) + *(pr+rblocks-2)
qi=(high48r \ high48b ) 'estimate quotient, max error 1
if qi > 0 then
'r= r - q*b for every block, begin at lsblock, with carry
carry=0
for i=blockshift to rblocks
if carry=0 then 'don't apply the offset!
sum = *(pr+i) - (*(pb+i-blockshift) * qi) + carry
else 'offset because the overflow 16->32 bit
sum = *(pr+i) - (*(pb+i-blockshift) * qi) + carry + offset
end if
*(pr+i)=sum
carry=sum shr 16
next i
end if
'------------------------------------------------------------
'if c>b, additional substraction needed
substract=0
for i=rblocks to blockshift step -1
substract= *(pr+i) - *(pb+i-blockshift)
if substract <> 0 then exit for
next i
if substract > 0 then
carry=1
qi += 1
for i=blockshift to rblocks
'substract from remainder from first step
sum = *(pr+i) + Cushort(Not(*(pb+i-blockshift))) + carry
*(pr+i)=sum
carry=sum shr 16
next i
End if
'------------------------------------------------------------
rblocks -= 1
blockshift -= 1
*pq=qi 'write quotient with pointer
pq -= 1 'next block of quotient
next blocks
'------------------------------------------------------------
'finalisation
prune(r) ' trim result
prune(q)
If sq Then q = negate(q) ' Xor of input signs
If sa Then r = negate(r) ' sign of original input A
'------------------------------------------------------------
End Sub
'=======================================================================
Function Factorial Overload(Byref a As bigint) As bigint
Dim As bigint f=bigint_1, n=a
Do Until is_LT(n, bigint_2)
f = multiply(f, n)
n = subtract(n, bigint_1)
Loop
Return f
End Function
'=======================================================================
'exponentiation modulus
Function mod_power(byref bb as bigint, byref ee as bigint, byref m as bigint) as bigint
if (ee.s[Len(ee.s)-1] And 128) then print "Cannot raise a bigint to a negative power"
if ee = bigint_0 then return 1
dim as bigint r= bigint_1, b=bb,x 'x is dump for the quotient
dim as integer bitlen,i,z
Dim as Uinteger<32> ptr pee
pee=cast(Uinteger<32> ptr, strptr(ee.s))
Dim as Uinteger<32> spee
spee = *pee 'load first block of exponent in variable
bitlen=msbit(ee) 'the highest set bit
divmod(b,m,x,b) 'initial reduction
'i counts from the lsb to msb, z counts the 32 bits in a block
for i=0 to bitlen-1
if Bit(spee,z) then 'if bit is set then multiply
r *= b
divmod(r,m,x,r)
end if
b=square(b)
divmod(b,m,x,b)
if z=31 then 'reset z, next block
z=0
pee += 1
spee = *pee
else
z += 1
end if
next i
r *= b 'bitvalue for highest bit=1
divmod(r,m,x,r)
return r
End Function
'=======================================================================
' BIT FUNCTIONS
'=======================================================================
' find the bit position of the first bit that differs from the sign bit
Function msbit(Byref a As bigint) As Integer
Dim As Integer i, j, k = 0
i = Len(a.s)
If 128 And a.s[Len(a.s)-1] Then ' negative
Do ' find the highest non-255 byte in the string
i = i - 1
j = a.s[i]
Loop Until (j < 255) Or (i = 0)
j = 255 - j
Else ' positive
Do ' find the highest non-zero byte in the string
i = i - 1
j = a.s[i]
Loop Until (j > 0) Or (i = 0)
End If
' find the highest non-sign bit in the byte
If j And 1 Then k = 1 ' straight code is faster than a loop
If j And 2 Then k = 2
If j And 4 Then k = 3
If j And 8 Then k = 4
If j And 16 Then k = 5
If j And 32 Then k = 6
If j And 64 Then k = 7
If j And 128 Then k = 8
k = k + (i * 8) - 1 ' if no bits differ (-1 or 0) then return -1
Return k
End Function
'================================================================
' NOT. Invert all the bits in a bigint
Function complement(Byref aa As bigint) As bigint
Dim As bigint a = aa
For i As Integer = 0 To Len(a.s)-1
a.s[i] = 255 - a.s[i]
Next i
Return a
End Function
'=======================================================================
' get the value of a specified bit in a big integer
Function Bit_Value(Byref v As bigint, Byval b As Ulongint) As Integer
Dim As Integer bitval, by = b \ 8
If by < Len(v.s) Then
bitval = Bit(v.s[by], b Mod 8)
Else
If v.s[Len(v.s)-1] And 128 Then bitval = -1 ' the extended sign bit
End If
Return bitval
End Function
'================================================================
' set a specified bit in a big integer
Function Bit_Set(Byref vv As bigint, Byval b As Ulongint) As bigint
Dim As bigint v = vv
Dim As Integer by, bi, delta, sign
by = b \ 8 ' byte number
bi = b Mod 8 ' bit number
delta = by - Len(v.s) + 1
If bi = 7 Then delta = delta + 1 ' protect the sign bit
If delta > 0 Then ' lengthen the number
delta = ((delta + 3)\ 4) * 4
If v.s[Len(v.s)-1] And 128 Then sign = -1 ' the extended sign bit
v.s = v.s + String(delta, sign)
End If
v.s[by] = Bitset(v.s[by], bi)
Return v
End Function
'================================================================
' clear a specified bit in a big integer
Function Bit_Reset(Byref vv As bigint, Byval b As Ulongint) As bigint
Dim As bigint v = vv
Dim As Integer by, bi, delta, sign
by = b \ 8 ' byte number
bi = b Mod 8 ' bit number
delta = by - Len(v.s) + 1
If bi = 7 Then delta = delta + 1 ' protect the sign bit
If delta > 0 Then ' lengthen the number
delta = ((delta + 3)\ 4) * 4
If v.s[Len(v.s)-1] And 128 Then sign = -1 ' the extended sign bit
v.s = v.s + String(delta, sign)
End If
v.s[by] = Bitreset(v.s[by], bi)
Return v
End Function
'================================================================
' shift bigint n bits left
Function shift_left(Byref a As bigint, Byval n As Uinteger) As bigint
If n = 0 Then Return a
Dim As Integer nblocks = n \ 32, nbits = n Mod 32
Dim As bigint s, m
s.s = String(nblocks * 4, Chr(0)) + a.s ' put zeros on the rhs
m = bit_set(bigint_0, nbits)
s = multiply(m, s)
Return s
End Function
'================================================================
' shift bigint n bits right, by shifting left nbits and right nblocks
Function shift_right(Byref a As bigint, Byval n As Uinteger) As bigint
If n = 0 Then Return a
If n > (8*Len(a.s)) Then Return bigint_0
Dim As Integer nblocks = n \ 32, nbits = n Mod 32
Dim As bigint s = a
Dim As bigint m = bit_set(bigint_0, 32 - nbits )
s = multiply(m, s) ' move bits left
s.s = Right(s.s, Len(s.s) - (nblocks+1)*4 )
Return s
End Function
'================================================================
' bitwise AND
Function Bit_And(Byref aa As bigint, Byref bb As bigint) As bigint
Dim As bigint a = aa, b = bb, c
Dim As Integer lena, lenb, i
lena = Len(a.s)
lenb = Len(b.s)
If lena > lenb Then
b.s = b.s + String(lena - lenb, Bit(b.s[lenb - 1], 7))
Else
a.s = a.s + String(lenb - lena, Bit(a.s[lena - 1], 7))
End If
c = b
For i = 0 To Len(c.s) - 1
c.s[i] = c.s[i] And a.s[i]
Next i
Return c
End Function
'================================================================
' bitwise Or
Function Bit_Or(Byref aa As bigint, Byref bb As bigint) As bigint
Dim As bigint a = aa, b = bb, c
Dim As Integer lena, lenb, i
lena = Len(a.s)
lenb = Len(b.s)
If lena > lenb Then
b.s = b.s + String(lena - lenb, Bit(b.s[lenb - 1], 7))
Else
a.s = a.s + String(lenb - lena, Bit(a.s[lena - 1], 7))
End If
c = b
For i = 0 To Len(c.s) - 1
c.s[i] = c.s[i] Or a.s[i]
Next i
Return c
End Function
'================================================================
' bitwise Xor
Function Bit_Xor(Byref aa As bigint, Byref bb As bigint) As bigint
Dim As bigint a = aa, b = bb, c
Dim As Integer lena, lenb, i
lena = Len(a.s)
lenb = Len(b.s)
If lena > lenb Then
b.s = b.s + String(lena - lenb, Bit(b.s[lenb - 1], 7))
Else
a.s = a.s + String(lenb - lena, Bit(a.s[lena - 1], 7))
End If
c = b
For i = 0 To Len(c.s) - 1
c.s[i] = c.s[i] Xor a.s[i]
Next i
Return c
End Function
'================================================================
' bitwise Eqv, equivalence is the complement of Xor
Function Bit_Eqv(Byref aa As bigint, Byref bb As bigint) As bigint
Dim As bigint a = aa, b = bb, c
Dim As Integer lena, lenb, i
lena = Len(a.s)
lenb = Len(b.s)
If lena > lenb Then
b.s = b.s + String(lena - lenb, Bit(b.s[lenb - 1], 7))
Else
a.s = a.s + String(lenb - lena, Bit(a.s[lena - 1], 7))
End If
c = b
For i = 0 To Len(c.s) - 1
c.s[i] = (c.s[i] Eqv a.s[i])
Next i
Return c
End Function
'================================================================
' bitwise Imp, implication
Function Bit_Imp(Byref aa As bigint, Byref bb As bigint) As bigint
Dim As bigint a = aa, b = bb, c
Dim As Integer lena, lenb, i
lena = Len(a.s)
lenb = Len(b.s)
If lena > lenb Then
b.s = b.s + String(lena - lenb, Bit(b.s[lenb - 1], 7))
Else
a.s = a.s + String(lenb - lena, Bit(a.s[lena - 1], 7))
End If
c = b
For i = 0 To Len(c.s) - 1
c.s[i] = (c.s[i] Imp a.s[i])
Next i
Return c
End Function
'================================================================
' Relational functions - Don't know why these are outcommented
'----------------------------------------------------------------
Function is_ZERO(Byref a As bigint) As Integer
If Len(Trim(a.s, Chr(0))) Then Return 0 Else Return -1
'if a=bigint_0 would be faster 942:1764 cycles
End Function
'----------------------------------------------------------------
'Function is_NZ(Byref a As bigint) As Integer
' If Len(Trim(a.s, Chr(0))) Then Return -1 Else Return 0
'End Function
'
''----------------------------------------------------------------
'Function is_POS(Byref a As bigint) As Integer
' ' note that zero is deemed to be positive
' If (128 And a.s[Len(a.s)-1]) Then Return 0 Else Return -1
'End Function
'
''----------------------------------------------------------------
'Function is_NEG(Byref a As bigint) As Integer
' If (128 And a.s[Len(a.s)-1]) Then Return -1 Else Return 0
'End Function
'
''----------------------------------------------------------------
'Function is_EQ(Byref a As bigint, Byref b As bigint) As Integer
' Dim As bigint c
' c = subtract(a, b)
' If Len(Trim(c.s, Chr(0))) Then Return 0 Else Return -1
'End Function
'
''----------------------------------------------------------------
'Function is_NEQ(Byref a As bigint, Byref b As bigint) As Integer
' Dim As bigint c
' c = subtract(a, b)
' If Len(Trim(c.s, Chr(0))) Then Return -1 Else Return 0
'End Function
'
'----------------------------------------------------------------
Function is_LT(Byref a As bigint, Byref b As bigint) As Integer
Dim As bigint c = subtract(a,b)
If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0
End Function
''----------------------------------------------------------------
'Function is_GT(Byref a As bigint, Byref b As bigint) As Integer
' Dim As bigint c = subtract(b, a)
' If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0
'End Function
'
''----------------------------------------------------------------
'Function is_LT_EQ(Byref a As bigint, Byref b As bigint) As Integer
' Dim As bigint c
' c = subtract(b, a)
' If (128 And c.s[Len(c.s)-1]) Then Return 0 Else Return -1
'End Function
''----------------------------------------------------------------
'Function is_GT_EQ(Byref a As bigint, Byref b As bigint) As Integer
' Dim As bigint c
' c = subtract(a, b)
' If (128 And c.s[Len(c.s)-1]) Then Return 0 Else Return -1
'End Function
'
'----------------------------------------------------------------
Function BigInt_Sgn(Byref a As bigint) As Integer
Dim As Integer i = 128 And a.s[Len(a.s)-1]
If i Then Return -1 ' is negative
If a=bigint_0 Then Return 0 ' is zero
Return 1 ' is positive
End Function
'=====================================================================
' Overload arithmetic operators
'=====================================================================
' unary plus +
Operator + (Byref x As bigint) As bigint
Return x
End Operator
' unary minus -
Operator - (Byref x As bigint) As bigint
Return negate(x)
End Operator
' addition
Operator + (Byref x As bigint, Byref y As bigint) As bigint
Return addition(x, y)
End Operator
' subtraction
Operator - (Byref x As bigint, Byref y As bigint) As bigint
Return subtract(x, y)
End Operator
Operator * (Byref x As bigint, Byref y As bigint) As bigint
Return multiply(x, y)
End Operator
Operator \ (Byref x As bigint, Byref y As bigint) As bigint
dim as bigint a,b
divmod(x,y,a,b)
Return a
End Operator
Operator / (Byref x As bigint, Byref y As bigint) As bigint
print "no floating point division in Big_Integer"
sleep: end
return 0
End Operator
Operator mod (Byref x As bigint, Byref y As bigint) as bigint
dim as bigint c,d
divmod(x,y,c,d)
return d
End Operator
'----------------------------------------------------------------
' operate and assign
Operator bigint.+= (Byref rhs As bigint)
this = addition(this, rhs)
End Operator
Operator bigint.-= (Byref rhs As bigint)
this = subtract(this, rhs)
End Operator
Operator bigint.*= (Byref rhs As bigint)
this = multiply(this, rhs)
End Operator
Operator bigint.\= (Byref rhs As bigint)
dim as bigint d
divmod(this,rhs,this,d)
End Operator
Operator bigint.mod= (Byref rhs as bigint)
dim as bigint c
divmod(this,rhs,c,this)
End Operator
'===================================================================
' relational operators, return a boolean as an integer, 0 is false.
'===================================================================
'Main Comparation function - faster than substract ...
Function bigint_compare(byref a as bigint, byref b as bigint) as integer
'return -1 for a>b; 1 for a<b and 0 for equal
dim as integer i,z
dim as ubyte vala, valb
dim as byte signa, signb
signa=128 And a.s[Len(a.s)-1] '-1 (true) negative, 0 (false) positive
signb=128 and b.s[len(b.s)-1]
'-------------------------------------------
'sign is different - easy:
If signa = 0 and signb = -128 Then
return -1
exit function
elseif signa = -128 and signb = 0 then
return 1
exit function
end if
'-------------------------------------------
'len is different - easy:
if len(a.s) > len(b.s) then
if signa = 0 then
return -1
else
return 1
end if
exit function
elseif len(a.s) < len(b.s) then
if signa = 0 then
return 1
else
return -1
end if
exit function
end if
'-------------------------------------------
'compare block for block:
for i=len(a.s)-1 to 0 step -1
z=a.s[i] - b.s[i]
if z > 0 then
return -1
exit function
elseif z < 0 then
return 1
exit function
end if
next i
return 0
end function
'----------------------------------------------------------------------
Operator = (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a, b)
if c=0 then return -1 else return 0
End Operator
'----------------------------------------------------------------
Operator <> (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a, b)
if c=0 then return 0 else return -1
End Operator
'----------------------------------------------------------------
Operator < (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a,b)
if c = 1 then return -1 else return 0
End Operator
'----------------------------------------------------------------
Operator > (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a,b)
if c = -1 then return -1 else return 0
End Operator
'----------------------------------------------------------------
Operator <= (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a,b)
if c = 1 or c = 0 then return -1 else return 0
End Operator
'----------------------------------------------------------------
Operator >= (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a, b)
if c = -1 or c = 0 then return -1 else return 0
End Operator
'----------------------------------------------------------------
' raise a bigint to the power of a long integer, result = x^y
Operator ^ (Byref x As bigint, Byref y As Longint) As bigint
Return power(x, y)
End Operator
'================================================================
' FOR bigint, NEXT bigint and STEP bigint
'=====================================================================
' implicit step versions. Implicit step is 1
Operator bigint.for( )
End Operator
Operator bigint.step( )
this += bigint_1
End Operator
Operator bigint.next( Byref end_cond As bigint ) As Integer
Return this <= end_cond
End Operator
'----------------------------------------------------------------
' explicit step versions
Operator bigint.for( Byref step_var As bigint )
End Operator
Operator bigint.step( Byref step_var As bigint )
this += step_var
End Operator
Operator bigint.next( Byref end_cond As bigint, Byref step_var As bigint ) As Integer
If step_var < bigint_0 Then
Return this >= end_cond
Else
Return this <= end_cond
End If
End Operator
'=================================================================
' more intrinsic numerical functions
'=================================================================
Function Abs_(Byref x As Double) As Double
Return Abs(x)
End Function
#undef Abs
Function Abs Overload (Byref x As Double) As Double
Return abs_(x)
End Function
Function Abs(Byref x As bigint) As bigint
Return Absolute(x)
End Function
'----------------------------------------------------------------
Function Sgn_(Byref x As Double) As Double
Return Sgn(x)
End Function
#undef Sgn
Function Sgn Overload (Byref x As Double) As Double
Return Sgn_(x)
End Function
Function Sgn(Byref x As bigint) As Integer
Return BigInt_Sgn(x)
End Function