Save this file as “Big_Integer_0.bas”
Code: Select all
' version 1.0 17 March 2010
'================================================================
' Big_Integer_0.bas Arbitrary Precision, Two's Complement Integer.
'================================================================
' 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
'----------------------------------------------------------------
' 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 can change if a string is created or any length is changed.
'================================================================
'Type bigint ' this will be used to overload FreeBASIC operators later
' s As String ' a little endian, two's complement number binary
'End Type ' number packed into a string in blocks of four bytes
'================================================================
' common small constants
Dim Shared As String bigint_m1, bigint_0, bigint_1, bigint_2
bigint_m1 = Chr(255,255,255,255)' minus 1
bigint_0 = Chr(0,0,0,0) ' zero
bigint_1 = Chr(1,0,0,0) ' one
bigint_2 = Chr(2,0,0,0) ' two
'=======================================================================
' find the bit position of the first bit that differs from the sign bit
Function msbit(Byref a As String) As Integer
Dim As Integer i, j, k = 0
i = Len(a)
If 128 And a[Len(a)-1] Then ' negative
Do ' find the highest non-255 byte in the string
i = i - 1
j = a[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[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
'-----------------------------------------------------------------------
' remove unnecessary leading blocks, prune time easily earns it's keep
Sub prune(Byref a As String)
a = Left(a, (((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.
'================================================================
' Negate the twos complement binary number in a string
Function negate(Byref a As String) As String
Dim As String s = a
Dim As Integer blocks = Len(s) Shr 2
Dim As Ulongint sum
Dim As Uinteger Ptr uipd, uipc, uips
uipd = Cast( Uinteger Ptr, Strptr(s))' the Uinteger data in string
uips = Cast( Uinteger Ptr, @sum) ' the high part of sum is carry
uipc = uips + 1 ' the low part of sum is sum mod 2^32
*uipc = 1 ' set carry
Do ' slow ahead until clear of the carry
sum = Clngint(Not *uipd) + *uipc
*uipd = *uips
uipd += 1
blocks -= 1
Loop Until (*uipc = 0) Or (blocks = 0)
If blocks > 0 Then
Do ' no carry, so full speed ahead
*uipd = Not *uipd
uipd +=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 *uipc = 0 Then ' carry was not generated by the most negative number
If (128 And a[Len(a)-1]) = (128 And s[Len(s)-1]) Then s = s + Chr(0,0,0,0)
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 String) As String
Dim As String b = a
If 128 And b[Len(b)-1] Then b = negate(b)
Return b
End Function
'================================================================
' pack ascii numeric string to a straight binary number in a string
' the output string will have a length that is a multiple of 4 bytes
Function pack( Byref ascii As String) As String
Dim As String a, b ' a is the ascii input, b is the binary output
Dim As Integer p, i, j, ch, blocks, sign
a = ascii
If Instr(a, "-") <> 0 Then sign = -1
'------------------------------------------------------------
' extract numerics from input string
j = 0 ' in-place write pointer
For i = 0 To Len(a) - 1
ch = a[i]
If (ch >= Asc("0")) And (ch <= Asc("9")) Then
a[j] = ch
j += 1
End If
Next i
a = Left(a, j) ' j is one ahead from string index = length of a
'------------------------------------------------------------
' extend to next multiple of 9 digits
i = Len(a)
blocks = Int(0.99 + i / 9) ' number of 9 digit blocks needed
p = 9 * blocks
a = String(p - i, "0") + a ' pad to next multiple of 9 digits
'------------------------------------------------------------
' decimal to binary conversion
i = ( 8 + Len(a) * 3.32192809488) \ 8 ' bytes needed for binary
blocks = 1 + (i \ 4) ' adjust to multiple of 4
b = String(blocks * 4, Chr(0) ) ' binary destination string
'------------------------------------------------------------
Dim As Uinteger Ptr bp, bpz, bpcarry, bpdata
bpz = Cast(Uinteger Ptr, Strptr(b)) ' binary output string[0]
Dim As Ulongint product, carry, multiplier = 1e9
bpdata = Cast(Uinteger Ptr, @product) ' points to bottom of product
bpcarry = bpdata + 1 ' points to top half of product
'------------------------------------------------------------
blocks = 1 ' blocks will be advanced as required by carry
For i = 1 To Len(a)-8 Step 9 ' msd to lsd in blocks of 9
bp = bpz ' point back to the start
carry = Val(Mid(a, i, 9)) ' take the next 9 digit block
For j = 1 To blocks
product = Clngint(*bp) * multiplier + carry
*bp = Cuint(*bpdata)
carry = Culngint(*bpcarry)
bp += 1
Next j
' advancing blocks only as needed doubles the speed of conversion
If Carry Then
*bp = carry
blocks += 1
End If
Next i
b = Left(b, blocks * 4) ' keep only used blocks
'-------------------------------------------------------------
If b[Len(b)-1] And 128 Then b = b + Chr(0,0,0,0) ' MSB must be 0
If sign Then b = negate(b)
'-------------------------------------------------------------
Return b
End Function
' on 2 GHz Pentium original time
' digits new time without blocks
' 100 45. usec 46 us
' 1000 520. usec 800 us
' 10k 29. msec 58 ms
' 100k 2.7 sec 5.7 seconds
' 1M 4.5 min 9.42 minutes
'===============================================================
' unpack a straight binary string to a decimal ascii string
Function unpack(Byref b As String) As String
Dim As String bs, ds
bs = b
ds = Chr(0) ' initial decimal output string
Dim As Integer i, j, product, carry, sign
' if negative then negate, append the sign later
If bs[Len(bs)-1] And 128 Then ' negative
sign = -1
bs = negate(bs)
End If
' change from base 2 to base 10
For j = Len(bs)-1 To 0 Step -1 ' each byte in string is base 256
carry = bs[j] ' the next byte to add after multiply
For i = Len(ds)-1 To 0 Step -1
product = 256 * ds[i] + carry
ds[i] = product Mod 10
carry = product \ 10
Next i
Do While carry > 0 ' output string overflows
ds = Chr(carry Mod 10) + ds ' extend output string
carry = carry \ 10 ' as needed
Loop
Next j
' change from Ubyte to ASCII
For i = 0 To Len(ds) - 1
ds[i] = ds[i] + Asc("0")
Next i
If sign Then ds = "-" + ds Else ds = "+" + ds
Return ds
End Function
'===============================================================
Function addition(Byref aa As String, Byref bb As String) As String
Dim As String 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)-1 ' might change the sign without carry
If a[lena] And 128 Then sa = 255 ' sign as a byte
i = a[lena] And 192 ' if MSBs differ then extend the sign
If (i = 64) Or (i = 128) Then a = a + String(4, Chr(sa) )
'------------------------------------------------------------
lenb = Len(b)-1
If b[lenb] And 128 Then sb = 255
i = b[lenb] And 192
If (i = 64) Or (i = 128) Then b = b + String(4, Chr(sb) )
'------------------------------------------------------------
' align the two strings to be added
delta = Len(a) - Len(b) 'new values
If delta <> 0 Then ' sign extend the shorter
If delta > 0 Then
a = a
If b[Len(b)-1] And 128 Then i = 255 Else i = 0
b = b + String(delta, Chr(i) ) ' extend b
Else
If aa[Len(aa)-1] And 128 Then i = 255 Else i = 0
a = a + 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) Shr 2
Dim As Ulongint sum = 0 ' clear carry
Dim As Uinteger Ptr uia, uib, low, carry
uia = Cast(Uinteger Ptr, Strptr(a) )
uib = Cast(Uinteger Ptr, Strptr(b) )
low = Cast(Uinteger Ptr, @sum ) ' data is low half of sum
carry = low + 1 ' carry is the high half of sum
For i = 1 To blocks
sum = Clngint(*uia) + Clngint(*uib) + Clngint(*carry)
*uia = *low
uia += 1
uib += 1
Next i
prune(a)
Return a
End Function
' 100 digits in 4.2 usec
'===============================================================
' subtract
Function subtract(Byref aa As String, Byref bb As String) As String
Dim As String cc = Negate(bb)
cc = addition(aa, cc)
Return cc
End Function
'================================================================
' multiply
Function multiply(Byref aa As String, Byref bb As String) As String
'------------------------------------------------------------
' sort out the signs and rectify the inputs
Dim As String a = aa, b = bb
Dim As Integer sign_a, sign_b, sign_c
sign_a = a[Len(a)-1] And 128
sign_b = b[Len(b)-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) ' number of bytes in a
bsize = Len(b) ' number of bytes in b
Dim As String c = String(asize + bsize, Chr(0)) ' initialise accumulator
asize = asize Shr 2 ' number of blocks in a
bsize = bsize Shr 2 ' Shr 2 is faster than \4
'------------------------------------------------------------
' pointers into all the strings
Dim As Uinteger Ptr ia, ib, ic, iaz, ibz, icz
iaz = Cast(Uinteger Ptr, Strptr(a) )
ibz = Cast(Uinteger Ptr, Strptr(b) )
icz = Cast(Uinteger Ptr, Strptr(c) )
Dim As Ulongint product
Dim As Uinteger Ptr uipcarry, uipdata = Cast(Uinteger Ptr, @product)
uipcarry = uipdata + 1 ' top half of product is carry
'------------------------------------------------------------
ia = iaz
For i = 0 To asize - 1
ib = ibz
product = 0 ' clear carry
For j = 0 To bsize - 1
ic = icz + i + j
product = Culngint(*ia) * Culngint(*ib) + Culngint(*ic) + Culngint(*uipcarry)
*ic = *uipdata
ib += 1
Next j
*(ic+1) = *uipcarry
ia += 1
Next i
'------------------------------------------------------------
If sign_c = 128 Then c = negate(c)
prune(c)
Return c
End Function
' digits time
' 10 3.9 usec
' 100 11.5 usec
' 1k 640. usec
' 10k 64. msec
' 100k 6.35 sec
'=======================================================================
' square a number, approaches twice the speed of multiply for big numbers
Function square(Byref aa As String) As String
'------------------------------------------------------------
Dim As String a = aa
If (128 And a[Len(a)-1]) Then a = negate(a)
'------------------------------------------------------------
' find the dimension of the problem
Dim As Integer i, j, asize = Len(a) ' number of bytes in a
Dim As String c = String(2 * asize, Chr(0) ) ' initialise accumulator
asize = (asize Shr 2) - 1 ' highest block number in a
'------------------------------------------------------------
' pointers into all the strings
Dim As Uinteger Ptr ia, ib, ic, iaz, icz
iaz = Cast(Uinteger Ptr, Strptr(a) ) ' the base addresses of strings
icz = Cast(Uinteger Ptr, Strptr(c) )
Dim As Ulongint product ' bottom half is data, top half will be carry
Dim As Uinteger Ptr uipcarry, uipdata = Cast(Uinteger Ptr, @product)
uipcarry = uipdata + 1 ' top half of product is carry
'------------------------------------------------------------
' multiply one triangular corner only
ia = iaz
For i = 1 To asize
ia += 1 ' ia starts at 1 not zero because 0,0 is on the diagonal
ib = iaz ' ib is second pointer into a, ib starts at zero
ic = icz + i ' ic is the accumulator ic = icz + i + j
product = 0 ' clear carry
For j = 0 To i - 1
product = Culngint(*ia) * Culngint(*ib) + Culngint(*ic) + Culngint(*uipcarry)
*ic = *uipdata
ib += 1
ic += 1
Next j
*ic = *uipcarry ' line of blocks gets one longer each loop
Next i
'------------------------------------------------------------
' double the triangle, cannot do it at same time as squaring diagonal
ic = icz ' because it can cause the product to overflow
product = 0 ' clear carry
For i = 0 To (2 * asize) + 1
product = Culngint(*ic) + Culngint(*ic) + Culngint(*uipcarry)
*ic = *uipdata
ic += 1
Next i
'------------------------------------------------------------
' square and accumulate the diagonal elements
ia = iaz
ic = icz
product = 0 ' clear carry
For i = 0 To asize
' square the diagonal entry, while propagating carry
product = Culngint(*ia) * Culngint(*ia) + Culngint(*ic) + Culngint(*uipcarry)
*ic = *uipdata
ic += 1
' propagate carry through the following odd block of C
product = Culngint(*ic) + Culngint(*uipcarry)
*ic = *uipdata
ic += 1
ia += 1
Next i
'------------------------------------------------------------
If 128 And c[Len(c)-1] Then c = c + Chr(0,0,0,0) ' 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
'=======================================================================
' shift up one bit, low towards high
Function mul2(Byref a As String) As String
Dim As Integer i, sign, sum, carry = 0
Dim As String b = a
sign = b[Len(b) - 1] And 128 ' snag the msb of highest byte
For i = 0 To Len(b) - 1
sum = b[i] + b[i] + carry
carry = Hibyte(sum)
b[i] = Lobyte(sum)
Next i
If ( b[Len(b) - 1] And 128 ) <> sign Then
carry = carry * 255
b = b + 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 String) As String
Dim As Integer i, carry = 0
Dim As String b = a
For i = 0 To Len(a)-2 ' all except the top byte of four
b[i] = ( b[i] Shr 1 ) + (128 * (b[i+1] And 1))
Next i
i = Len(b) - 1
b[i] = b[i] Shr 1
b[i] = b[i] Or (2 * ( b[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 String, Byval n As Ulongint) As String
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 String pwr = Chr(1,0,0,0) ' 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
'=======================================================================
Declare Function is_LT(Byref As String, Byref As String) As Integer
'=======================================================================
Function Factorial(Byref a As String) As String
Dim As String f = Chr(1,0,0,0), n = a
Do Until is_LT(n, Chr(2,0,0,0))
f = multiply(f, n)
n = subtract(n, Chr(1,0,0,0))
Loop
Return f
End Function
'=======================================================================
' bit functions
'=======================================================================
' NOT. Invert all the bits in a string
Function complement(Byref aa As String) As String
Dim As String a = aa
For i As Integer = 0 To Len(a)-1
a[i] = 255 - a[i]
Next i
Return a
End Function
'=======================================================================
' get the value of a specified bit in a big integer
Function Bit_Value(Byref v As String, Byval b As Ulongint) As Integer
Dim As Integer bitval, by = b \ 8
If by < Len(v) Then
bitval = Bit(v[by], b Mod 8)
Else
If v[Len(v)-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 String, Byval b As Ulongint) As String
Dim As String v = vv
Dim As Integer by, bi, delta, sign
by = b \ 8 ' byte number
bi = b Mod 8 ' bit number
delta = by - Len(v) + 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[Len(v)-1] And 128 Then sign = -1 ' the extended sign bit
v = v + String(delta, sign)
End If
v[by] = Bitset(v[by], bi)
Return v
End Function
'================================================================
' clear a specified bit in a big integer
Function Bit_Reset(Byref vv As String, Byval b As Ulongint) As String
Dim As String v = vv
Dim As Integer by, bi, delta, sign
by = b \ 8 ' byte number
bi = b Mod 8 ' bit number
delta = by - Len(v) + 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[Len(v)-1] And 128 Then sign = -1 ' the extended sign bit
v = v + String(delta, sign)
End If
v[by] = Bitreset(v[by], bi)
Return v
End Function
'================================================================
' shift string n bits left
Function shift_left(Byref a As String, Byval n As Uinteger) As String
If n = 0 Then Return a
Dim As Integer nblocks = n \ 32, nbits = n Mod 32
Dim As String s = String(nblocks * 4, Chr(0)) + a ' put zeros on the rhs
Dim As String m = bit_set(Chr(0,0,0,0), nbits)
s = multiply(m, s)
Return s
End Function
'================================================================
' shift string n bits right, by shifting left nbits and right nblocks
Function shift_right(Byref a As String, Byval n As Uinteger) As String
If n = 0 Then Return a
If n > (8*Len(a)) Then Return Chr(0,0,0,0)
Dim As Integer nblocks = n \ 32, nbits = n Mod 32
Dim As String s = a
Dim As String m = bit_set(Chr(0,0,0,0), 32 - nbits )
s = multiply(m, s) ' move bits left
s = Right(s, Len(s) - (nblocks+1)*4 )
Return s
End Function
'================================================================
' bitwise AND
Function Bit_And(Byref aa As String, Byref bb As String) As String
Dim As String a = aa, b = bb
Dim As Integer lena, lenb, i
lena = Len(a)
lenb = Len(b)
If lena > lenb Then
b = b + String(lena - lenb, Bit(b[lenb - 1], 7))
Else
a = a + String(lenb - lena, Bit(a[lena - 1], 7))
End If
Dim As String c = b
For i = 0 To Len(c) - 1
c[i] = c[i] And a[i]
Next i
Return c
End Function
'================================================================
' bitwise Or
Function Bit_Or(Byref aa As String, Byref bb As String) As String
Dim As String a = aa, b = bb
Dim As Integer lena, lenb, i
lena = Len(a)
lenb = Len(b)
If lena > lenb Then
b = b + String(lena - lenb, Bit(b[lenb - 1], 7))
Else
a = a + String(lenb - lena, Bit(a[lena - 1], 7))
End If
Dim As String c = b
For i = 0 To Len(c) - 1
c[i] = c[i] Or a[i]
Next i
Return c
End Function
'================================================================
' bitwise Xor
Function Bit_Xor(Byref aa As String, Byref bb As String) As String
Dim As String a = aa, b = bb
Dim As Integer lena, lenb, i
lena = Len(a)
lenb = Len(b)
If lena > lenb Then
b = b + String(lena - lenb, Bit(b[lenb - 1], 7))
Else
a = a + String(lenb - lena, Bit(a[lena - 1], 7))
End If
Dim As String c = b
For i = 0 To Len(c) - 1
c[i] = c[i] Xor a[i]
Next i
Return c
End Function
'================================================================
' bitwise Eqv, equivalence is the complement of Xor
Function Bit_Eqv(Byref aa As String, Byref bb As String) As String
Dim As String a = aa, b = bb
Dim As Integer lena, lenb, i
lena = Len(a)
lenb = Len(b)
If lena > lenb Then
b = b + String(lena - lenb, Bit(b[lenb - 1], 7))
Else
a = a + String(lenb - lena, Bit(a[lena - 1], 7))
End If
Dim As String c = b
For i = 0 To Len(c) - 1
c[i] = (c[i] Eqv a[i])
Next i
Return c
End Function
'================================================================
' bitwise Imp, implication
Function Bit_Imp(Byref aa As String, Byref bb As String) As String
Dim As String a = aa, b = bb
Dim As Integer lena, lenb, i
lena = Len(a)
lenb = Len(b)
If lena > lenb Then
b = b + String(lena - lenb, Bit(b[lenb - 1], 7))
Else
a = a + String(lenb - lena, Bit(a[lena - 1], 7))
End If
Dim As String c = b
For i = 0 To Len(c) - 1
c[i] = (c[i] Imp a[i])
Next i
Return c
End Function
'================================================================
' convert a string to binary
Function show_bin(Byref s As String) As String
Dim As Integer i
Dim As String h ' lsb is string[0] = little endian
For i = Len(s)-1 To 0 Step -1
h = h + Bin(s[i], 8)
Next i
Return h
End Function
'================================================================
' convert a string to hexadecimal
Function show_hex(Byref s As String) As String
Dim As Integer i
Dim As String h ' lsb is string[0] = little endian
For i = Len(s)-1 To 0 Step -1
h = h + Hex(s[i], 2)
Next i
Return h
End Function
'================================================================
' convert a string to octal
Function show_oct(Byref a As String) As String
Dim As String b, c, s = a
If 128 And a[Len(a)-1] Then ' extend the sign
s = s + Chr(255,255,255,255)
Else
s = s + Chr(0,0,0,0)
End If
Dim As Integer i
Dim As Ulongint u
For i = 1 To Len(a) Step 3
b = Mid(s, i, 3) ' take three bytes = 24 bits
u = b[0] + 256 * (b[1] + 256 * b[2])
c = Oct(u, 8) + c ' eight symbols per three bytes
Next i
Return c
End Function
'================================================================
' Relational functions
'----------------------------------------------------------------
Function is_ZERO(Byref a As String) As Integer
If Len(Trim(a, Chr(0))) Then Return 0 Else Return -1
End Function
'----------------------------------------------------------------
Function is_NZ(Byref a As String) As Integer
If Len(Trim(a, Chr(0))) Then Return -1 Else Return 0
End Function
'----------------------------------------------------------------
Function is_POS(Byref a As String) As Integer
' note that zero is deemed to be positive
If (128 And a[Len(a)-1]) Then Return 0 Else Return -1
End Function
'----------------------------------------------------------------
Function is_NEG(Byref a As String) As Integer
If (128 And a[Len(a)-1]) Then Return -1 Else Return 0
End Function
'----------------------------------------------------------------
Function is_EQ(Byref a As String, Byref b As String) As Integer
If Len(Trim(subtract(a, b), Chr(0))) Then Return 0 Else Return -1
End Function
'----------------------------------------------------------------
Function is_NEQ(Byref a As String, Byref b As String) As Integer
If Len(Trim(subtract(a, b), Chr(0))) Then Return -1 Else Return 0
End Function
'----------------------------------------------------------------
Function is_LT(Byref a As String, Byref b As String) As Integer
Dim As String c = subtract(a, b)
If (128 And c[Len(c)-1]) Then Return -1 Else Return 0
End Function
'----------------------------------------------------------------
Function is_GT(Byref a As String, Byref b As String) As Integer
Dim As String c = subtract(b, a)
If (128 And c[Len(c)-1]) Then Return -1 Else Return 0
End Function
'----------------------------------------------------------------
Function is_LT_EQ(Byref a As String, Byref b As String) As Integer
Dim As String c = subtract(b, a)
If (128 And c[Len(c)-1]) Then Return 0 Else Return -1
End Function
'----------------------------------------------------------------
Function is_GT_EQ(Byref a As String, Byref b As String) As Integer
Dim As String c = subtract(a, b)
If (128 And c[Len(c)-1]) Then Return 0 Else Return -1
End Function
'=======================================================================
Function BigInt_Sgn(Byref a As String) As Integer
Dim As Integer i = 128 And a[Len(a)-1]
If i Then Return -1 ' is negative
If is_zero(a) Then Return 0 ' is zero
Return 1 ' is positive
End Function
'=======================================================================
' integer divide, a \ b, return div and mod
Sub divmod(_
Byref aa As String, Byref bb As String,_
Byref intdiv As String, Byref intmod As String)
Dim As String a = aa, b = bb, c
If is_zero(b) Then
Print " Division by zero. "
Stop
End If
Dim As Integer lena = Len(a), lenb = Len(b)
If (lena <= 8) And (lenb <= 8) Then ' arguments are one or two blocks
intdiv = Chr(0,0,0,0,0,0,0,0) ' use 64 bit long integers
intmod = Chr(0,0,0,0,0,0,0,0) ' allocate space for results
If lena = 4 Then a = a + String(4, Bit(a[lena - 1], 7))
If lenb = 4 Then b = b + String(4, Bit(b[lenb - 1], 7))
Dim As Longint Ptr pa, pb, pc, pd
pa = Cast(Longint Ptr, Strptr(a))
pb = Cast(Longint Ptr, Strptr(b))
pc = Cast(Longint Ptr, Strptr(intdiv))
pd = Cast(Longint Ptr, Strptr(intmod))
*pc = *pa \ *pb
*pd = *pa Mod *pb
'------------------------------------------------------------
Else ' more than two blocks so do a long division
'------------------------------------------------------------
' sign of inputs and output
Dim As Integer sa, sb, sq
sa = 128 And a[Len(a)-1]
sb = 128 And b[Len(b)-1]
sq = sa Xor sb ' sign of the result
' rectify the inputs
If sa Then a = negate(a)
If sb Then b = negate(b)
'------------------------------------------------------------
' prepare to normalise the divisor
Dim As Integer ja = msbit(a), jb = msbit(b), bitshifts = ja - jb
If ja < jb Then ' abandon the division as result is zero
intdiv = Chr(0,0,0,0)
intmod = aa
Else ' proceed with the long slow division
b = shift_Left(b, bitshifts) ' align the first bits
lena = Len(a)
lenb = Len(b)
If lena <> lenb Then
If lena > lenb Then
b = b + String(lena - lenb, Bit(b[lenb - 1], 7) )
Else
a = a + String(lenb - lena, Bit(b[lenb - 1], 7) )
lena = Len(a) ' update length of lena only
End If
End If
'------------------------------------------------------------
Dim As String qu = String(lena, Chr(0)) ' quotient
Dim As Integer i, j, qbit, blocks = lena Shr 2
c = String(lena, Chr(0) ) ' prepare to receive the result
'------------------------------------------------------------
' setup all pointers for 32 bit block arithmetic
Dim As Uinteger Ptr pa, pb, pc, pcarry, pdata
Dim As Ulongint sum
pdata = Cast(Uinteger Ptr, @sum) ' bottom half of sum is data
pcarry = pdata + 1 ' point to carry in top half of sum
For i = 0 To bitshifts
'------------------------------------
' setup and perform subtraction
pa = Cast(Uinteger Ptr, Strptr(a))
pb = Cast(Uinteger Ptr, Strptr(b)) ' strings can move
pc = Cast(Uinteger Ptr, Strptr(c))
j = blocks
*pcarry = 1
Do ' inline subtract, c = a - b
sum = Culngint(*pa) + Culngint(Not(*pb)) + Culngint(*pcarry)
*pc = *pdata ' save this result
pa += 1 ' advance all three pointers
pb += 1
pc += 1
j = j - 1 ' count the 32 bit blocks
Loop Until j = 0
'------------------------------------
' test sign of result
If *pcarry Then
' result of subtraction was positive
pa = Cast(Uinteger Ptr, Strptr(a))
pc = Cast(Uinteger Ptr, Strptr(c)) ' strings can move
For j = 1 To blocks ' copy c to a
*pa = *pc
pa += 1
pc += 1
Next j
qbit = 1
Else ' result of subtraction was negative
qbit = 0
End If
'------------------------------------
' shift quotient one bit left, insert qbit on right
Dim As Uinteger sum16
For j = 0 To Len(b) - 1
sum16 = qu[j] + qu[j] + qbit
qbit = Hibyte(sum16)
qu[j] = Lobyte(sum16)
Next j
'------------------------------------
' shift b one bit right
For j = 0 To Len(a)-2 ' all except the top byte
b[j] = ( b[j] Shr 1 ) + (128 * (b[j+1] And 1))
Next j
b[lena - 1] = b[lena - 1] Shr 1
'------------------------------------
Next i ' do next bit
prune(a) ' finished division
prune(qu)
If sq Then qu = negate(qu) ' Xor of input signs
If sa Then a = negate(a) ' sign of original input A
intdiv = qu ' the quotient
intmod = a ' the remainder
'------------------------------------------------------------
End If
End If
End Sub
' 100 digits divided by 1 takes 775. usec
'===============================================================
'' divide a by b, return the integer result
Function divide(Byref a As String, Byref b As String) As String
Dim As String c, d
divmod(a, b, c, d)
Return c
End Function
'===============================================================
' divide a by b and return only the integer remainder
Function remainder(Byref a As String, Byref b As String) As String
Dim As String c, d
divmod(a, b, c, d)
Return d
End Function
'===============================================================
' E n d o f B i g I n t e g e r C o d e
'===============================================================
Code: Select all
#include "big_integer_0.bas" ' or whatever you call the file
dim as string a, b, c, d
a = pack("123456789012345678901234567890")
print " input ="; unpack(a)
b = square(a)
print " square ="; unpack(b)
c = divide(b, a)
print " divide ="; unpack(c)
d = subtract(c, a)
print "subtract="; unpack(d)
sleep
Overloaded code is now available at
http://www.freebasic.net/forum/viewtopi ... 704#158704