## Arbitrarily Big Integer Routines

User projects written in or related to FreeBASIC.
Richard
Posts: 3029
Joined: Jan 15, 2007 20:44
Location: Australia

### Arbitrarily Big Integer Routines

Here is a full set of fast big integer routines that automatically adjust their size as needed. Hopefully they all work correctly. Please post any bugs or omissions found.

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)
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
'===============================================================

Here is an example program.

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

The next step will be to write them as overload functions for FB.
Overloaded code is now available at
viewtopic.php?p=158704#158704
Last edited by Richard on Jun 06, 2011 8:43, edited 1 time in total.
Ophelius
Posts: 428
Joined: Feb 26, 2006 1:57
Good job! Looks very useful. I'll test it later
PyRoe
Posts: 61
Joined: Sep 07, 2005 16:10
Location: not quite sure, i just woke up here

### This is freaking sweet

With your help I just calculated 2^1,000,000 in no time flat. BRAVO!!!!
duke4e
Posts: 717
Joined: Dec 04, 2005 0:16
Location: Varazdin, Croatia, Europe
Contact:
Great! Can't wait for overloaded version
integer
Posts: 391
Joined: Feb 01, 2007 16:54
Location: usa

### Cannot grasp what is wrong inside the function

I am trying to use a power-mod function:

Code: Select all

#include "big_integer_0.bas"

'----------------------------------------------------------------------
' raise a number to a positive power; take the number mod m
'   result = b^e mod m         [base^exponent : modulus]
'
Function powermod(Byval Number As String, Byval Exponent As string, byval Modulus as string) As String

Dim As String pwr, t, b, e, m

b = Number
e = Exponent
m = Modulus

t = multiply(m,b)

print "inside function powermod()"
print "    b= "; unpack(b)
print "    e= "; unpack(e)
print "    m= "; unpack(m)
print "t=b*m= "; unpack(t)
t = multiply(m,e)
print " b=m*e ";unpack(t)
sleep

pwr = pack("1")

Do while is_NZ( e )
if Bit_Value(e,0) then
pwr = multiply(pwr,b)
pwr = remainder(pwr,m)
endif
e = div2(e)
b = square(b)
b = remainder(b,m)
Loop
Return pwr

End Function

'==========================================
'sample trial program

dim as string b,e,m,r,s

b = pack("53")
e = pack("1023")
m = pack("2047")
s = multiply(b,m)

print "-------------------------"
print "before the function call"
print "    b= ";unpack(b)
print "    e= ";unpack(e)
print "    m= ";unpack(m)
print "s=b*m= ";unpack(s)
print "-------------------------"

r = powermod( b, e, m )

print "b^e mod m= ";unpack(r)
sleep

The multiply operation work fine before the call to the powermod function.
It appears as if the values were passed correctly, as the same values are printed inside the function match the values printed before the call to the powermod function.

HOWEVER, the multiply function inside the powermod() returns a zero (+0 )

This is bugging me.

The values created before the call are the same values printed inside the called function [ powermod (b,e,m) ]

The multipy operation (before the call) gives the desired answer,
The same values inside the function yield zero as the result of the multiply.

Code: Select all

-------------------------
before the function call
b= +53
e= +1023
m= +2047
s=b*m= +108491
-------------------------
inside function powermod()
b= +53
e= +1023
m= +2047
t=b*m= +0
t = multiply(m,e)
b=m*e  +0

Anyone have an idea as to what is happening?
Richard
Posts: 3029
Joined: Jan 15, 2007 20:44
Location: Australia
I don't have an answer to your question now, but I do have the latest Overload version of Bigint.

Code: Select all

' version 1.1   6 June 2011
'================================================================
' Big_Integer_1.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 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 Integer)
Declare Constructor (Byref As Longint)
Declare Constructor (Byref As String)
' 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 Integer)
Declare Operator Let(Byval rhs As Uinteger)
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
' 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)
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 Operator < (Byref As bigint, Byref As bigint) As Integer
Declare Function is_LT(Byref As bigint, Byref As bigint) As Integer

'================================================================
Constructor bigint  ' default constructor
this.s = Chr(0,0,0,0)   ' zero
End Constructor

Constructor bigint (Byref rhs As Integer)
this = rhs
End Constructor

Constructor bigint (Byref rhs As Longint)
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 Integer)
this.s = Chr(0,0,0,0)
Dim As Integer Ptr bip = Cptr(Integer Ptr, Strptr(this.s))
Dim As Integer Ptr  ip = Cptr(Integer Ptr, @rhs)
*bip = *ip
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Uinteger)
this.s = Chr(0,0,0,0)
Dim As Integer Ptr bip = Cptr(Integer Ptr, Strptr(this.s))
Dim As Integer Ptr uip = Cptr(Integer Ptr, @rhs)
*bip = *uip
If (128 And rhs) 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

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

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

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

'================================================================
' 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 Ptr uipd, uipc, uips
uipd = Cast( Uinteger Ptr, Strptr(s.s))' the Uinteger data in bigint
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.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

'================================================================
' pack ascii numeric string to a straight binary number in a bigint
' the output bigint will have a length that is a multiple of 4 bytes
Function pack( Byref ascii As String) As bigint
Dim As String a  ' a is the ascii input
Dim As bigint b  ' 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.s = String(blocks * 4, Chr(0) ) ' binary destination string
'------------------------------------------------------------
Dim As Uinteger Ptr bp, bpz, bpcarry, bpdata
bpz = Cast(Uinteger Ptr, Strptr(b.s)) ' binary output string[0]
Dim As Ulongint product, carry, multiplier = 1e9
bpdata = Cast(Uinteger Ptr, @product) ' bottom half of product
bpcarry = bpdata + 1                ' 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 = Valulng(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 ' an exact count of the blocks used
End If
Next i
b.s = Left(b.s, blocks * 4) ' keep only used blocks
'-------------------------------------------------------------
If b.s[Len(b.s)-1] And 128 Then b.s = b.s + bigint_0.s ' 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 bb As bigint) As String
Dim As bigint b
Dim As String d
b = bb
d = Chr(0)   ' initial decimal output string
Dim As Integer i, j, product, carry, sign
' if negative then negate, append the sign later
If b.s[Len(b.s)-1] And 128 Then   ' negative
sign = -1
b = negate(b)
End If
' change from base 2 to base 10
For j = Len(b.s)-1 To 0 Step -1 ' each byte in string is base 256
carry = b.s[j]   ' the next byte to add after multiply
For i = Len(d)-1 To 0 Step -1
product = 256 * d[i] + carry
d[i] = product Mod 10
carry = product \ 10
Next i
Do While carry > 0  ' output string overflows
d = Chr(carry Mod 10) + d   ' extend output string
carry = carry \ 10            '  as needed
Loop
Next j
' change from Ubyte to ASCII
For i = 0 To Len(d) - 1
d[i] = d[i] + Asc("0")
Next i
If sign Then d = "-" + d Else d = "+" + d
Return d
End Function

'===============================================================
Function addition(Byref aa As bigint, Byref bb As bigint) As bigint
Dim As bigint a, b
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 Ptr uia, uib, low, carry
uia = Cast(Uinteger Ptr, Strptr(a.s) )
uib = Cast(Uinteger Ptr, Strptr(b.s) )
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 bigint, Byref bb As bigint) As bigint
Dim As bigint cc = Negate(bb)
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 Ptr ia, ib, ic, iaz, icz
iaz = Cast(Uinteger Ptr, Strptr(a.s) )  ' the base addresses of bigints
icz = Cast(Uinteger Ptr, Strptr(c.s) )
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.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 ' number of blocks in a
bsize = bsize Shr 2 ' Shr 2 is faster than \4
'------------------------------------------------------------
' pointers into all the bigints
Dim As Uinteger Ptr ia, ib, ic, iaz, ibz, icz
iaz = Cast(Uinteger Ptr, Strptr(a.s) )
ibz = Cast(Uinteger Ptr, Strptr(b.s) )
icz = Cast(Uinteger Ptr, Strptr(c.s) )
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 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

'=======================================================================
Function Factorial Overload(Byref a As bigint) As bigint
Dim As bigint f, n
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

'=======================================================================
' bit functions
'=======================================================================
' 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

'================================================================
' convert a bigint to binary
Function show_bin(Byref s As bigint) As String
Dim As Integer i
Dim As String h     ' lsb is string[0] = little endian
For i = Len(s.s)-1 To 0 Step -1
h = h + Bin(s.s[i], 8)
Next i
Return h
End Function

'================================================================
' convert a bigint to hexadecimal
Function show_hex(Byref s As bigint) As String
Dim As Integer i
Dim As String h     ' lsb is string[0] = little endian
For i = Len(s.s)-1 To 0 Step -1
h = h + Hex(s.s[i], 2)
Next i
Return h
End Function

'================================================================
' convert a bigint to octal
Function show_oct(Byref a As bigint) As String
Dim As String b, c
Dim As bigint s
s = a
If 128 And a.s[Len(a.s)-1] Then ' extend the sign
s.s = s.s + bigint_m1.s
Else
s.s = s.s + bigint_0.s
End If
Dim As Integer i
Dim As Ulongint u
For i = 1 To Len(a.s) Step 3
b = Mid(s.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 bigint) As Integer
If Len(Trim(a.s, Chr(0))) Then Return 0 Else Return -1
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
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 = "0" Then Return 0 ' is zero
If is_zero(a) Then Return 0 ' is zero
Return 1 ' is positive
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 = "0" then return 0    ' 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 bigint, Byref bb As bigint,_
Byref intdiv As bigint, Byref intmod As bigint)
Dim As bigint a = aa, b = bb, c
If b = "0" Then 'If is_zero(b) Then
Print " Division by zero. "
Stop
End If
Dim As Integer lena = Len(a.s), lenb = Len(b.s)
If (lena <= 8) And (lenb <= 8) Then ' arguments are one or two blocks ++++++++++++++++++++++
intdiv.s = Chr(0,0,0,0,0,0,0,0)   ' use 64 bit long integers      signed or unsigned?
intmod.s = Chr(0,0,0,0,0,0,0,0)   ' allocate space for results    is this buggy because of Longint ?
If lena = 4 Then a.s = a.s + String(4, Bit(a.s[lena - 1], 7))
If lenb = 4 Then b.s = b.s + String(4, Bit(b.s[lenb - 1], 7))
Dim As Longint Ptr pa, pb, pc, pd
pa = Cast(Longint Ptr, Strptr(a.s))
pb = Cast(Longint Ptr, Strptr(b.s))
pc = Cast(Longint Ptr, Strptr(intdiv.s))
pd = Cast(Longint Ptr, Strptr(intmod.s))
*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.s[Len(a.s)-1]
sb = 128 And b.s[Len(b.s)-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 = bigint_0
intmod = aa
Else    ' proceed with the long slow division
b = shift_Left(b, bitshifts)    ' align the first bits
lena = Len(a.s)
lenb = Len(b.s)
If lena <> lenb Then
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(b.s[lenb - 1], 7) )
lena = Len(a.s)   ' update length of lena only
End If
End If
'------------------------------------------------------------
Dim As bigint qu
qu.s = String(lena, Chr(0))  ' quotient
Dim As Integer i, j, qbit, blocks = lena Shr 2
c.s = 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.s))
pb = Cast(Uinteger Ptr, Strptr(b.s))  ' strings can move
pc = Cast(Uinteger Ptr, Strptr(c.s))
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.s))
pc = Cast(Uinteger Ptr, Strptr(c.s))  ' 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.s) - 1
sum16 = qu.s[j] + qu.s[j] + qbit
qbit = Hibyte(sum16)
qu.s[j] = Lobyte(sum16)
Next j
'------------------------------------
' shift b one bit right
For j = 0 To Len(a.s)-2   ' all except the top byte
b.s[j] = ( b.s[j] Shr 1 ) + (128 * (b.s[j+1] And 1))
Next j
b.s[lena - 1] = b.s[lena - 1] Shr 1 ' top byte
'------------------------------------
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 bigint, Byref b As bigint) As bigint
Dim As bigint 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 bigint, Byref b As bigint) As bigint
Dim As bigint c, d
divmod(a, b, c, d)
Return d
End Function

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

Operator + (Byref x As bigint, Byref y As bigint) As bigint
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
Return divide(x, y)
End Operator

Operator / (Byref x As bigint, Byref y As bigint) As bigint
Return divide(x, y) ' should this be prohibited ?
End Operator

'----------------------------------------------------------------
' operate and assign
Operator bigint.+= (Byref rhs As bigint)
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)
this = divide(this, rhs)
End Operator

'=================================================================
' relational operators, return a boolean as an integer, 0 is false.
'=================================================================
Operator = (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 Operator

'----------------------------------------------------------------
Operator <> (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 Operator

'----------------------------------------------------------------
Operator < (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 -1 Else Return 0
End Operator

'----------------------------------------------------------------
Operator > (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 Operator

'----------------------------------------------------------------
Operator <= (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 Operator

'----------------------------------------------------------------
Operator >= (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 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

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

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

'===============================================================
'       E n d    o f    B i g    I n t e g e r    C o d e
'===============================================================
Here is an old test program that demonstrates UDT and Overloaded operators. It names above as "big_integer_1.bas"

Code: Select all

#include "big_integer_1.bas"
Dim As Double startime, stoptime
Dim As bigint a, b = 2, c
Dim As String s1, s2
Print

dim as double bits, blocks, bytes, dig = 10000 ' decimal digits
print dig; " decimal digits"
bits = dig * 3.3219281
blocks = int(bits / 32)
bytes = blocks * 4

a.s = string(bytes, chr(0)) + chr(1,0,0,0)
print "length "; len(a.s); " bytes"

startime = timer
b = a * a
stoptime = timer
print " squareing ", 1e3*(stoptime - startime); " millisec"

b = a
startime = timer
b = b * a
stoptime = timer
print "multiplying", 1e3*(stoptime - startime); " millisec"

print "length "; len(b.s); " bytes"
print "done"

Sleep
Richard
Posts: 3029
Joined: Jan 15, 2007 20:44
Location: Australia
I think due to Byval using Zstrings and finding a zero byte ?
Sorry about the delay in posting the latest Bigint, I am still not satisfied.
Bigint Overload seems to handle the Byval mode without the Zstring problem.
@ Integer. Here is your new code, does it now work as you would expect ?

Code: Select all

#include "big_integer_1.bas"
'----------------------------------------------------------------------
' raise a number to a positive power; take the number mod m
'   result = b^e mod m         [base^exponent : modulus]
'----------------------------------------------------------------------
' convert Byval to Byref
function powermod(Byref Number As Bigint, Byref Exponent As Bigint, Byref Modulus As Bigint) As Bigint

Dim As Bigint pwr, t, b, e, m

b = Number
e = Exponent
m = Modulus

t = m * b

Print "inside function powermod()"
Print "    b= "; b
Print "    e= "; e
Print "    m= "; m
Print "t=b*m= "; t
t = m * e
Print " b=m*e "; t
' Sleep
Print
pwr = pack("1")

Do While e <> 0             'is_NZ( e )
If Bit_Value(e,0) Then
pwr = pwr * b
pwr = remainder(pwr, m)
End If
e = div2(e)
b = square(b)
b = remainder(b, m)
Loop
Return pwr

End Function

'==========================================
'sample trial program

Dim As Bigint b,e,m,r,s

b = "53"
e = "1023"
m = "2047"
s = b * m

Print "-------------------------"
Print "before the function call"
Print "    b= "; b
Print "    e= "; e
Print "    m= "; m
Print "s=b*m= "; s
Print "-------------------------"

r = powermod( b, e, m )

Print "b^e mod m= "; r
Sleep
integer
Posts: 391
Joined: Feb 01, 2007 16:54
Location: usa
That is SUPER!
It works as expected.

Thank you, Richard.

This BigInt program is very beneficial.
I appreciate the enormous amount of time you have invested in this.

P.S.

WOW!

Thanks.
integer
Posts: 391
Joined: Feb 01, 2007 16:54
Location: usa
@Richard

The overloaded version is simply awesome!
I know you tire of the adulation from someone who does not grasp the entire intricate complexity of your creation. I am nonplussed -- awestruck.

And now a question:
Is there a function to convert from BigInt back to longint (or integer)?

This is my routine for that:

Code: Select all

'--------------------------------------
'convert a BigInteger to a Long Integer
'--------------------------------------
function BigInt2LongInt(byref BigNumber as BigInt) as longint
dim as string s
dim as longint i
s = mid(unpack(BigNumber),2)    'skip the sign byte
i = vallng(s)
if BigNumber < 0 then i=-i : endif
return (i)
end function

The problem is the time required.
For a few thousand calls the time used is not significant.
However, when dealing with a few (hundred) million conversations, it begins to be important.

Is there a more efficient way to convert a BigInter to a LongInteger?

Thanks.
Richard
Posts: 3029
Joined: Jan 15, 2007 20:44
Location: Australia
Try this. It avoids conversion to decimal and then back again to binary.
Here are two functions with test code.

Code: Select all

#include "Big_Integer_1.bas"

Function BigInt_2_LongInt ( Byref a As BigInt ) As Longint
Dim As String s = a.s
If Len(s) > 8 Then
Print " Overflow in BigInt to LongInteger conversion. "
Print a
Sleep : Stop
End If
If Len(s) = 4 Then  ' sign extend 4 byte integer to 8 byte LongInt
If Bit(s[3], 7) Then
s = s + String(4, Chr(255))
Else
s = s + String(4, Chr(0))
End If
End If
Return *Cptr(Longint Ptr, Strptr(s))
End Function

Function BigInt_2_Integer ( Byref a As BigInt ) As Integer
Dim As String s = a.s
If Len(s) <> 4 Then
Print " Overflow in BigInt to Integer conversion. "
Print a
Sleep : Stop
End If
Return *Cptr(Integer Ptr, Strptr(s))
End Function

Dim As BigInt x
x = 3
For i As Integer = 1 To 20
Print BigInt_2_Integer( x ),
Print x
x = x * 10
Next i

Sleep
Adulation is always welcome from an Integer.
It still surprises me that I wrote BigInt and that it actually works.
Richard
Posts: 3029
Joined: Jan 15, 2007 20:44
Location: Australia
Here is a fast BigInt to Double converter with some test code.

This function may not be much use because it can only accurately represent 15 digit integers.
It never generates overflow since it will return a double precision coded as IEEE signed infinity for BigInts over 308 digits. When used to pass numbers to other FB code it can replace the BigInt_2_Integer code. It is slower but avoids overflow.

It does however give a result that can represent a magnitude of 308 digits which means it may have some application passing numbers to other FB code for human interpretation or graphical display.

Code: Select all

'===================================================================
' Convert a BigInt to a double with 1.#INF overflow if too big
'===================================================================
#include "Big_Integer_1.bas"

Function BigInt_2_Double ( Byref a As BigInt ) As Double
Dim As Bigint b = a
Dim As Ulongint uli, Sign_Bit
Dim As Integer i
If Bit( b.s[Len(b.s) - 1], 7) Then  ' extract the sign bit
Sign_Bit = Culngint(1) Shl 63   ' remember the sign
b = -b  ' rectify
End If  ' b is now a positive BigInt
' overflows double? if mag > 1.797693134862310e308 = signed infinity
If Len(b.s) > 128 Then ' 128 bytes * 8 bits per byte = 1024 bits
If Len(b.s) = 132 Then  ' special case of sign bit = entire block
If ( b.s[131] Or b.s[130] Or b.s[129] Or b.s[128] ) <> 0 Then
uli = Sign_Bit Or &h7FF0000000000000 ' all ones exponent
Return *Cast(Double Ptr, @uli)  ' return the bit pattern
End If
End If
End If
If Len(b.s) = 4 Then    ' test for simple zero
If ( b.s[3] Or b.s[2] Or b.s[1] Or b.s[0] ) = 0 Then Return 0
End If
Dim As Longint expo = 8 * Len(b.s) + 1022 ' = bits - 1 + bias of expo
' if needed for the conversion, extend tail with two zero LS blocks
If Len(b.s) < 12 Then b.s = Chr(0,0,0,0, 0,0,0,0) + b.s
' the ms block still contains the data, so no change to expo
Dim As Ubyte Ptr ubp = Strptr(b.s) + Len(b.s) - 1 ' point to the MSbyte
For i = 0 To 4  ' find the leading non-zero byte, MS block may be zero
If *ubp > 0 Then Exit For
ubp = ubp - 1
expo = expo - 8 ' expo reduction of 8 bits per zero byte skipped
Next i  ' ubp now points to the MS non-zero byte
uli = *Cast(Ulongint Ptr, ubp - 7)  ' normalize bytes, left justify
For i = 63 To 56 Step -1    ' find the MS set bit
If Bit(uli, i) Then Exit For
expo = expo - 1
Next i  ' i now points to MSbit
uli = uli Shr (i - 52)  ' move MSbit into bit 52
uli = Bitreset(uli, 52) ' kill the MSbit implied bit in 52
uli = Sign_Bit Or (expo Shl 52) Or uli  ' build the double
Return *Cast(Double Ptr, @uli)
End Function

'===================================================================
' test the BigInt_2_Double function
Dim As BigInt x
x = -8
Dim As Double d
Dim As Integer i
For i = 0 To 17
Print x; "  ";
Print Tab(8)
Print BigInt_2_Double( x )
x = x +1
Next i
Print

For i = 303 To 308
Dim As String s = "1797" + String(i, "0")
x = s
Print "  "; Left(s, 20); "  "; Len(s); ".  "; Len(x.s); "  ";
Print BigInt_2_Double( s )
Next i
Print

For i = 303 To 308
Dim As String s = "1798" + String(i, "0")
x = s
Print "  "; Left(s, 20); "  "; Len(s); ".  "; Len(x.s); "  ";
Print BigInt_2_Double( s )
Next i
Print

x.s = String(128, Chr(255)) + Chr(0,0,0,0) ' most positive double
Print x
Print Len(x.s),
Print BigInt_2_Double(x)
Print

x = x + 1
Print x
Print Len(x.s),
Print BigInt_2_Double(x)
Print

x.s = Chr(1) + String(127, Chr(0)) + Chr(255,255,255,255) ' most negative
Print x
Print Len(x.s),
Print BigInt_2_Double(x)
Print

x = x - 1
Print x
Print Len(x.s),
Print BigInt_2_Double(x)

'===================================================================
Sleep
'===================================================================
' double       Sign         Exponent         Significand
' Bit count      1             11                 52
' Bit range     63           62..52           51......0
' notes       1 = neg     bias = 1023     plus implied = 53 bits
'
' Signed infinities are represented with Exp being all 1s (2047).
' If the fraction part is zero then it is Infinity, else it is NaN.
'===================================================================
Edited to improve infinity detection.
Richard
Posts: 3029
Joined: Jan 15, 2007 20:44
Location: Australia
CInt(), CLongInt() and CDbl() are now supported in this latest version called Big_Integer_2.bas

Code: Select all

' version 1.2   3 July 2011
'================================================================
' Big_Integer_2.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
' 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 Integer)
Declare Constructor (Byref As Longint)
Declare Constructor (Byref As String)
' 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 Integer)
Declare Operator Let(Byval rhs As Uinteger)
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)
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 Operator < (Byref As bigint, Byref As bigint) As Integer
Declare Function is_LT(Byref As bigint, Byref As bigint) As Integer

'================================================================
Constructor bigint  ' default constructor
this.s = Chr(0,0,0,0)   ' zero
End Constructor

Constructor bigint (Byref rhs As Integer)
this = rhs
End Constructor

Constructor bigint (Byref rhs As Longint)
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 Integer)
this.s = Chr(0,0,0,0)
Dim As Integer Ptr bip = Cptr(Integer Ptr, Strptr(this.s))
Dim As Integer Ptr  ip = Cptr(Integer Ptr, @rhs)
*bip = *ip
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Uinteger)
this.s = Chr(0,0,0,0)
Dim As Integer Ptr bip = Cptr(Integer Ptr, Strptr(this.s))
Dim As Integer Ptr uip = Cptr(Integer Ptr, @rhs)
*bip = *uip
If (128 And rhs) 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

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

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

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

'================================================================
' 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 Ptr uipd, uipc, uips
uipd = Cast( Uinteger Ptr, Strptr(s.s))' the Uinteger data in bigint
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.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

'================================================================
' pack ascii numeric string to a straight binary number in a bigint
' the output bigint will have a length that is a multiple of 4 bytes
Function pack( Byref ascii As String) As bigint
Dim As String a  ' a is the ascii input
Dim As bigint b  ' 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.s = String(blocks * 4, Chr(0) ) ' binary destination string
'------------------------------------------------------------
Dim As Uinteger Ptr bp, bpz, bpcarry, bpdata
bpz = Cast(Uinteger Ptr, Strptr(b.s)) ' binary output string[0]
Dim As Ulongint product, carry, multiplier = 1e9
bpdata = Cast(Uinteger Ptr, @product) ' bottom half of product
bpcarry = bpdata + 1                ' 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 = Valulng(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 ' an exact count of the blocks used
End If
Next i
b.s = Left(b.s, blocks * 4) ' keep only used blocks
'-------------------------------------------------------------
If b.s[Len(b.s)-1] And 128 Then b.s = b.s + bigint_0.s ' 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 bb As bigint) As String
Dim As bigint b
Dim As String d
b = bb
d = Chr(0)   ' initial decimal output string
Dim As Integer i, j, product, carry, sign
' if negative then negate, append the sign later
If b.s[Len(b.s)-1] And 128 Then   ' negative
sign = -1
b = negate(b)
End If
' change from base 2 to base 10
For j = Len(b.s)-1 To 0 Step -1 ' each byte in string is base 256
carry = b.s[j]   ' the next byte to add after multiply
For i = Len(d)-1 To 0 Step -1
product = 256 * d[i] + carry
d[i] = product Mod 10
carry = product \ 10
Next i
Do While carry > 0  ' output string overflows
d = Chr(carry Mod 10) + d   ' extend output string
carry = carry \ 10            '  as needed
Loop
Next j
' change from Ubyte to ASCII
For i = 0 To Len(d) - 1
d[i] = d[i] + Asc("0")
Next i
If sign Then d = "-" + d Else d = "+" + d
Return d
End Function

'===============================================================
Function addition(Byref aa As bigint, Byref bb As bigint) As bigint
Dim As bigint a, b
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 Ptr uia, uib, low, carry
uia = Cast(Uinteger Ptr, Strptr(a.s) )
uib = Cast(Uinteger Ptr, Strptr(b.s) )
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 bigint, Byref bb As bigint) As bigint
Dim As bigint cc = Negate(bb)
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 Ptr ia, ib, ic, iaz, icz
iaz = Cast(Uinteger Ptr, Strptr(a.s) )  ' the base addresses of bigints
icz = Cast(Uinteger Ptr, Strptr(c.s) )
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.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 ' number of blocks in a
bsize = bsize Shr 2 ' Shr 2 is faster than \4
'------------------------------------------------------------
' pointers into all the bigints
Dim As Uinteger Ptr ia, ib, ic, iaz, ibz, icz
iaz = Cast(Uinteger Ptr, Strptr(a.s) )
ibz = Cast(Uinteger Ptr, Strptr(b.s) )
icz = Cast(Uinteger Ptr, Strptr(c.s) )
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 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

'=======================================================================
Function Factorial Overload(Byref a As bigint) As bigint
Dim As bigint f, n
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

'=======================================================================
' bit functions
'=======================================================================
' 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

'================================================================
' convert a bigint to binary
Function show_bin(Byref s As bigint) As String
Dim As Integer i
Dim As String h     ' lsb is string[0] = little endian
For i = Len(s.s)-1 To 0 Step -1
h = h + Bin(s.s[i], 8)
Next i
Return h
End Function

'================================================================
' convert a bigint to hexadecimal
Function show_hex(Byref s As bigint) As String
Dim As Integer i
Dim As String h     ' lsb is string[0] = little endian
For i = Len(s.s)-1 To 0 Step -1
h = h + Hex(s.s[i], 2)
Next i
Return h
End Function

'================================================================
' convert a bigint to octal
Function show_oct(Byref a As bigint) As String
Dim As String b, c
Dim As bigint s
s = a
If 128 And a.s[Len(a.s)-1] Then ' extend the sign
s.s = s.s + bigint_m1.s
Else
s.s = s.s + bigint_0.s
End If
Dim As Integer i
Dim As Ulongint u
For i = 1 To Len(a.s) Step 3
b = Mid(s.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 bigint) As Integer
If Len(Trim(a.s, Chr(0))) Then Return 0 Else Return -1
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
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 = "0" Then Return 0 ' is zero
If is_zero(a) Then Return 0 ' is zero
Return 1 ' is positive
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 = "0" then return 0    ' 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 bigint, Byref bb As bigint,_
Byref intdiv As bigint, Byref intmod As bigint)
Dim As bigint a = aa, b = bb, c
If b = "0" Then 'If is_zero(b) Then
Print " Division by zero. "
Stop
End If
Dim As Integer lena = Len(a.s), lenb = Len(b.s)
If (lena <= 8) And (lenb <= 8) Then ' arguments are one or two blocks ++++++++++++++++++++++
intdiv.s = Chr(0,0,0,0,0,0,0,0)   ' use 64 bit long integers      signed or unsigned?
intmod.s = Chr(0,0,0,0,0,0,0,0)   ' allocate space for results    is this buggy because of Longint ?
If lena = 4 Then a.s = a.s + String(4, Bit(a.s[lena - 1], 7))
If lenb = 4 Then b.s = b.s + String(4, Bit(b.s[lenb - 1], 7))
Dim As Longint Ptr pa, pb, pc, pd
pa = Cast(Longint Ptr, Strptr(a.s))
pb = Cast(Longint Ptr, Strptr(b.s))
pc = Cast(Longint Ptr, Strptr(intdiv.s))
pd = Cast(Longint Ptr, Strptr(intmod.s))
*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.s[Len(a.s)-1]
sb = 128 And b.s[Len(b.s)-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 = bigint_0
intmod = aa
Else    ' proceed with the long slow division
b = shift_Left(b, bitshifts)    ' align the first bits
lena = Len(a.s)
lenb = Len(b.s)
If lena <> lenb Then
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(b.s[lenb - 1], 7) )
lena = Len(a.s)   ' update length of lena only
End If
End If
'------------------------------------------------------------
Dim As bigint qu
qu.s = String(lena, Chr(0))  ' quotient
Dim As Integer i, j, qbit, blocks = lena Shr 2
c.s = 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.s))
pb = Cast(Uinteger Ptr, Strptr(b.s))  ' strings can move
pc = Cast(Uinteger Ptr, Strptr(c.s))
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.s))
pc = Cast(Uinteger Ptr, Strptr(c.s))  ' 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.s) - 1
sum16 = qu.s[j] + qu.s[j] + qbit
qbit = Hibyte(sum16)
qu.s[j] = Lobyte(sum16)
Next j
'------------------------------------
' shift b one bit right
For j = 0 To Len(a.s)-2   ' all except the top byte
b.s[j] = ( b.s[j] Shr 1 ) + (128 * (b.s[j+1] And 1))
Next j
b.s[lena - 1] = b.s[lena - 1] Shr 1 ' top byte
'------------------------------------
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 bigint, Byref b As bigint) As bigint
Dim As bigint 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 bigint, Byref b As bigint) As bigint
Dim As bigint c, d
divmod(a, b, c, d)
Return d
End Function

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

Operator + (Byref x As bigint, Byref y As bigint) As bigint
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
Return divide(x, y)
End Operator

Operator / (Byref x As bigint, Byref y As bigint) As bigint
Return divide(x, y) ' should this WARN or be prohibited ?
End Operator

'----------------------------------------------------------------
' operate and assign
Operator bigint.+= (Byref rhs As bigint)
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)
this = divide(this, rhs)
End Operator

'=================================================================
' relational operators, return a boolean as an integer, 0 is false.
'=================================================================
Operator = (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 Operator

'----------------------------------------------------------------
Operator <> (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 Operator

'----------------------------------------------------------------
Operator < (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 -1 Else Return 0
End Operator

'----------------------------------------------------------------
Operator > (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 Operator

'----------------------------------------------------------------
Operator <= (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 Operator

'----------------------------------------------------------------
Operator >= (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 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

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

'===================================================================
' Conversion Cast to Integer, LongInteger and Double
'===================================================================
Function BigInt_2_Integer(Byref a As BigInt) As Integer
Dim As String s = a.s
If Len(s) <> 4 Then
Print " Overflow in BigInt to Integer conversion. "
Print a
Sleep : Stop
End If
Return *Cptr(Integer Ptr, Strptr(s))
End Function

Operator bigint.Cast() As Integer   ' CInt
Return BigInt_2_Integer(this)
End Operator

'----------------------------------------------------------------
Function BigInt_2_LongInt(Byref a As BigInt) As Longint
Dim As String s = a.s
If Len(s) > 8 Then
Print " Overflow in BigInt to LongInteger conversion. "
Print a
Sleep : Stop
End If
If Len(s) = 4 Then  ' sign extend 4 byte integer to 8 byte LongInt
If Bit(s[3], 7) Then
s = s + String(4, Chr(255))
Else
s = s + String(4, Chr(0))
End If
End If
Return *Cptr(Longint Ptr, Strptr(s))
End Function

Operator bigint.Cast() As LongInt   ' CLongInt
Return BigInt_2_LongInt(this)
End Operator

'----------------------------------------------------------------
Function BigInt_2_Double ( Byref a As BigInt ) As Double
Dim As Bigint b = a
Dim As Ulongint uli, Sign_Bit
If Bit( b.s[Len(b.s) - 1], 7) Then  ' extract the sign bit
Sign_Bit = Culngint(1) Shl 63   ' remember the sign
b = -b  ' rectify
End If  ' b is now a positive BigInt
' overflows double? if mag > 1.797693134862310e308 = signed infinity
If Len(b.s) > 128 Then ' 128 bytes * 8 bits per byte = 1024 bits
If Len(b.s) = 132 Then  ' special case of sign bit = entire block
If ( b.s[131] Or b.s[130] Or b.s[129] Or b.s[128] ) <> 0 Then
uli = Sign_Bit Or &h7FF0000000000000 ' all ones exponent
Return *Cast(Double Ptr, @uli)  ' bit pattern is a double
End If
End If
End If
If Len(b.s) = 4 Then    ' test for simple zero
If ( b.s[3] Or b.s[2] Or b.s[1] Or b.s[0] ) = 0 Then Return 0
End If
Dim As Longint expo = 8 * Len(b.s) + 1022 ' = bits + expo_bias - 1
' if needed for the conversion, extend tail with two LS blocks of zero
If Len(b.s) < 12 Then b.s = Chr(0,0,0,0, 0,0,0,0) + b.s
' the ms block still contains the data, so no change to expo
Dim As Ubyte Ptr ubp = Strptr(b.s) + Len(b.s) - 1 ' point to the MSbyte
Dim As Integer i
For i = 0 To 4  ' find the leading non-zero byte, MS block may be zero
If *ubp > 0 Then Exit For
ubp = ubp - 1
expo = expo - 8 ' expo reduction of 8 bits per zero byte skipped
Next i  ' ubp now points to the MS non-zero byte
uli = *Cast(Ulongint Ptr, ubp - 7)  ' normalize bytes, left justify
For i = 63 To 56 Step -1    ' find the MS set bit
If Bit(uli, i) Then Exit For
expo = expo - 1
Next i  ' i now points to MSbit
uli = uli Shr (i - 52)  ' shift right to put MSbit in bit 52
uli = Bitreset(uli, 52) ' kill only the implicit bit now in bit 52
uli = Sign_Bit Or (expo Shl 52) Or uli  ' build the double
Return *Cast(Double Ptr, @uli)  ' return the bit pattern as a double
End Function

Operator bigint.Cast() As Double    ' Cdbl
Return BigInt_2_Double(this)
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

'===============================================================
'       E n d    o f    B i g    I n t e g e r    C o d e
'===============================================================
counting_pine
Posts: 6221
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs
Wow, you've put a lot of work into this, maybe it should be in Projects now.
stephanbrunker
Posts: 62
Joined: Nov 02, 2013 14:57

### Re: Arbitrarily Big Integer Routines

I found one bug in the file, but otherwise it works great - and the're a few additions:

- shift_right was only commented

----------

EDIT:
I needed to understand the differences between two's complement and unsigned first and how this program works.

So my first try was slow and buggy.

For my purpose (a Carter-Wegman MAC - requires 128-bit operation), I need the conversion functions:

- bigint to unsigned binary string and vice versa
- unsigned hex to bigint for easier input

The binary string is little endian, the hex is big endian.
Last edited by stephanbrunker on Nov 03, 2013 22:01, edited 1 time in total.
stephanbrunker
Posts: 62
Joined: Nov 02, 2013 14:57

### Re: Arbitrarily Big Integer Routines

I found a major bug in addition:

Try this:

Code: Select all

#include "Big_Integer_2.bas"

dim as bigint a,b,c
dim as uinteger bi

a.s=chr(0,0,0,0,183,38,83,152)
bi=3036098375
b=bi

print "     a:", show_hex(a)
print "   + b:", "        ";show_hex(b)
print "==============================="

c=a+b

print "   = c:", show_hex(c)

sleep

I found the bug in the operator bigint.let for uinteger:

Code: Select all

If (128 And rhs) Then this.s += Chr(0,0,0,0) ' make it positive

but 128 is not the sign bit for uinteger ...

I added my conversion functions, fixed the bug at uinteger, added the constructor for ulong.
The new conversion functions are (signed and unsigned):
hex_bigint
uhex_bigint
bin_bigint
ubin_bigint
bigint_bin
bigint_ubin

So, here the code:

Code: Select all

'version 1.3 4 November 2013
'=============================
'bugfix constructor and operator uinteger, added ulong
'conversion functions hex_bigint, uhex_bigint,
'bin_bigint, bigint_bin;  ubin_bigint, bigint_ubin

' version 1.2   3 July 2011
'================================================================
' Big_Integer_2.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
' 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 Integer)
Declare Constructor (Byref As Longint)
Declare Constructor (Byref As String)
Declare Constructor (Byref As Ulong)
' 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 Integer)
Declare Operator Let(Byval rhs As Uinteger)
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)
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 Operator < (Byref As bigint, Byref As bigint) As Integer
Declare Function is_LT(Byref As bigint, Byref As bigint) As Integer

'================================================================
Constructor bigint  ' default constructor
this.s = Chr(0,0,0,0)   ' zero
End Constructor

Constructor bigint (Byref rhs As Integer)
this = rhs
End Constructor

Constructor bigint (Byref rhs As ulong)
this = rhs
End Constructor

Constructor bigint (Byref rhs As Longint)
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 Integer)
this.s = Chr(0,0,0,0)
Dim As Integer Ptr bip = Cptr(Integer Ptr, Strptr(this.s))
Dim As Integer Ptr  ip = Cptr(Integer Ptr, @rhs)
*bip = *ip
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Uinteger)
this.s = Chr(0,0,0,0)
Dim As Uinteger Ptr bip = Cptr(Uinteger Ptr, Strptr(this.s))
Dim As Uinteger Ptr uip = Cptr(Uinteger 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 Ulong)
this.s = Chr(0,0,0,0)
Dim As Ulong Ptr bip = Cptr(Ulong Ptr, Strptr(this.s))
Dim As Ulong Ptr uip = Cptr(Ulong 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

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

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

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

'================================================================
' 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 Ptr uipd, uipc, uips
uipd = Cast( Uinteger Ptr, Strptr(s.s))' the Uinteger data in bigint
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.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

'================================================================
' pack ascii numeric string to a straight binary number in a bigint
' the output bigint will have a length that is a multiple of 4 bytes
Function pack( Byref ascii As String) As bigint
Dim As String a  ' a is the ascii input
Dim As bigint b  ' 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.s = String(blocks * 4, Chr(0) ) ' binary destination string
'------------------------------------------------------------
Dim As Uinteger Ptr bp, bpz, bpcarry, bpdata
bpz = Cast(Uinteger Ptr, Strptr(b.s)) ' binary output string[0]
Dim As Ulongint product, carry, multiplier = 1e9
bpdata = Cast(Uinteger Ptr, @product) ' bottom half of product
bpcarry = bpdata + 1                ' 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 = Valulng(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 ' an exact count of the blocks used
End If
Next i
b.s = Left(b.s, blocks * 4) ' keep only used blocks
'-------------------------------------------------------------
If b.s[Len(b.s)-1] And 128 Then b.s = b.s + bigint_0.s ' 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 bb As bigint) As String
Dim As bigint b
Dim As String d
b = bb
d = Chr(0)   ' initial decimal output string
Dim As Integer i, j, product, carry, sign
' if negative then negate, append the sign later
If b.s[Len(b.s)-1] And 128 Then   ' negative
sign = -1
b = negate(b)
End If
' change from base 2 to base 10
For j = Len(b.s)-1 To 0 Step -1 ' each byte in string is base 256
carry = b.s[j]   ' the next byte to add after multiply
For i = Len(d)-1 To 0 Step -1
product = 256 * d[i] + carry
d[i] = product Mod 10
carry = product \ 10
Next i
Do While carry > 0  ' output string overflows
d = Chr(carry Mod 10) + d   ' extend output string
carry = carry \ 10            '  as needed
Loop
Next j
' change from Ubyte to ASCII
For i = 0 To Len(d) - 1
d[i] = d[i] + Asc("0")
Next i
If sign Then d = "-" + d Else d = "+" + d
Return d
End Function

'===============================================================
Function addition(Byref aa As bigint, Byref bb As bigint) As bigint
Dim As bigint a, b
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 Ptr uia, uib, low, carry
uia = Cast(Uinteger Ptr, Strptr(a.s) )
uib = Cast(Uinteger Ptr, Strptr(b.s) )
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 bigint, Byref bb As bigint) As bigint
Dim As bigint cc = Negate(bb)
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 Ptr ia, ib, ic, iaz, icz
iaz = Cast(Uinteger Ptr, Strptr(a.s) )  ' the base addresses of bigints
icz = Cast(Uinteger Ptr, Strptr(c.s) )
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.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 ' number of blocks in a
bsize = bsize Shr 2 ' Shr 2 is faster than \4
'------------------------------------------------------------
' pointers into all the bigints
Dim As Uinteger Ptr ia, ib, ic, iaz, ibz, icz
iaz = Cast(Uinteger Ptr, Strptr(a.s) )
ibz = Cast(Uinteger Ptr, Strptr(b.s) )
icz = Cast(Uinteger Ptr, Strptr(c.s) )
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 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

'=======================================================================
Function Factorial Overload(Byref a As bigint) As bigint
Dim As bigint f, n
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

'=======================================================================
' bit functions
'=======================================================================
' 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

'================================================================
' convert a bigint to binary
Function show_bin(Byref s As bigint) As String
Dim As Integer i
Dim As String h     ' lsb is string[0] = little endian
For i = Len(s.s)-1 To 0 Step -1
h = h + Bin(s.s[i], 8)
Next i
Return h
End Function

'================================================================
' convert a bigint to hexadecimal
Function show_hex(Byref s As bigint) As String
Dim As Integer i
Dim As String h     ' lsb is string[0] = little endian
For i = Len(s.s)-1 To 0 Step -1
h = h + Hex(s.s[i], 2)
Next i
Return h
End Function

'================================================================
' convert a bigint to octal
Function show_oct(Byref a As bigint) As String
Dim As String b, c
Dim As bigint s
s = a
If 128 And a.s[Len(a.s)-1] Then ' extend the sign
s.s = s.s + bigint_m1.s
Else
s.s = s.s + bigint_0.s
End If
Dim As Integer i
Dim As Ulongint u
For i = 1 To Len(a.s) Step 3
b = Mid(s.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 bigint) As Integer
If Len(Trim(a.s, Chr(0))) Then Return 0 Else Return -1
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
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 = "0" Then Return 0 ' is zero
If is_zero(a) Then Return 0 ' is zero
Return 1 ' is positive
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 = "0" then return 0    ' 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 bigint, Byref bb As bigint,_
Byref intdiv As bigint, Byref intmod As bigint)
Dim As bigint a = aa, b = bb, c
If b = "0" Then 'If is_zero(b) Then
Print " Division by zero. "
Stop
End If
Dim As Integer lena = Len(a.s), lenb = Len(b.s)
If (lena <= 8) And (lenb <= 8) Then ' arguments are one or two blocks ++++++++++++++++++++++
intdiv.s = Chr(0,0,0,0,0,0,0,0)   ' use 64 bit long integers      signed or unsigned?
intmod.s = Chr(0,0,0,0,0,0,0,0)   ' allocate space for results    is this buggy because of Longint ?
If lena = 4 Then a.s = a.s + String(4, Bit(a.s[lena - 1], 7))
If lenb = 4 Then b.s = b.s + String(4, Bit(b.s[lenb - 1], 7))
Dim As Longint Ptr pa, pb, pc, pd
pa = Cast(Longint Ptr, Strptr(a.s))
pb = Cast(Longint Ptr, Strptr(b.s))
pc = Cast(Longint Ptr, Strptr(intdiv.s))
pd = Cast(Longint Ptr, Strptr(intmod.s))
*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.s[Len(a.s)-1]
sb = 128 And b.s[Len(b.s)-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 = bigint_0
intmod = aa
Else    ' proceed with the long slow division
b = shift_Left(b, bitshifts)    ' align the first bits
lena = Len(a.s)
lenb = Len(b.s)
If lena <> lenb Then
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(b.s[lenb - 1], 7) )
lena = Len(a.s)   ' update length of lena only
End If
End If
'------------------------------------------------------------
Dim As bigint qu
qu.s = String(lena, Chr(0))  ' quotient
Dim As Integer i, j, qbit, blocks = lena Shr 2
c.s = 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.s))
pb = Cast(Uinteger Ptr, Strptr(b.s))  ' strings can move
pc = Cast(Uinteger Ptr, Strptr(c.s))
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.s))
pc = Cast(Uinteger Ptr, Strptr(c.s))  ' 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.s) - 1
sum16 = qu.s[j] + qu.s[j] + qbit
qbit = Hibyte(sum16)
qu.s[j] = Lobyte(sum16)
Next j
'------------------------------------
' shift b one bit right
For j = 0 To Len(a.s)-2   ' all except the top byte
b.s[j] = ( b.s[j] Shr 1 ) + (128 * (b.s[j+1] And 1))
Next j
b.s[lena - 1] = b.s[lena - 1] Shr 1 ' top byte
'------------------------------------
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 bigint, Byref b As bigint) As bigint
Dim As bigint 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 bigint, Byref b As bigint) As bigint
Dim As bigint c, d
divmod(a, b, c, d)
Return d
End Function

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

Operator + (Byref x As bigint, Byref y As bigint) As bigint
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
Return divide(x, y)
End Operator

Operator / (Byref x As bigint, Byref y As bigint) As bigint
Return divide(x, y) ' should this WARN or be prohibited ?
End Operator

'----------------------------------------------------------------
' operate and assign
Operator bigint.+= (Byref rhs As bigint)
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)
this = divide(this, rhs)
End Operator

'=================================================================
' relational operators, return a boolean as an integer, 0 is false.
'=================================================================
Operator = (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 Operator

'----------------------------------------------------------------
Operator <> (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 Operator

'----------------------------------------------------------------
Operator < (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 -1 Else Return 0
End Operator

'----------------------------------------------------------------
Operator > (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 Operator

'----------------------------------------------------------------
Operator <= (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 Operator

'----------------------------------------------------------------
Operator >= (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 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

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

'===================================================================
' Conversion Cast to Integer, LongInteger and Double
'===================================================================
Function BigInt_2_Integer(Byref a As BigInt) As Integer
Dim As String s = a.s
If Len(s) <> 4 Then
Print " Overflow in BigInt to Integer conversion. "
Print a
Sleep : Stop
End If
Return *Cptr(Integer Ptr, Strptr(s))
End Function

Operator bigint.Cast() As Integer   ' CInt
Return BigInt_2_Integer(this)
End Operator

'----------------------------------------------------------------
Function BigInt_2_LongInt(Byref a As BigInt) As Longint
Dim As String s = a.s
If Len(s) > 8 Then
Print " Overflow in BigInt to LongInteger conversion. "
Print a
Sleep : Stop
End If
If Len(s) = 4 Then  ' sign extend 4 byte integer to 8 byte LongInt
If Bit(s[3], 7) Then
s = s + String(4, Chr(255))
Else
s = s + String(4, Chr(0))
End If
End If
Return *Cptr(Longint Ptr, Strptr(s))
End Function

Operator bigint.Cast() As LongInt   ' CLongInt
Return BigInt_2_LongInt(this)
End Operator

'----------------------------------------------------------------
Function BigInt_2_Double ( Byref a As BigInt ) As Double
Dim As Bigint b = a
Dim As Ulongint uli, Sign_Bit
If Bit( b.s[Len(b.s) - 1], 7) Then  ' extract the sign bit
Sign_Bit = Culngint(1) Shl 63   ' remember the sign
b = -b  ' rectify
End If  ' b is now a positive BigInt
' overflows double? if mag > 1.797693134862310e308 = signed infinity
If Len(b.s) > 128 Then ' 128 bytes * 8 bits per byte = 1024 bits
If Len(b.s) = 132 Then  ' special case of sign bit = entire block
If ( b.s[131] Or b.s[130] Or b.s[129] Or b.s[128] ) <> 0 Then
uli = Sign_Bit Or &h7FF0000000000000 ' all ones exponent
Return *Cast(Double Ptr, @uli)  ' bit pattern is a double
End If
End If
End If
If Len(b.s) = 4 Then    ' test for simple zero
If ( b.s[3] Or b.s[2] Or b.s[1] Or b.s[0] ) = 0 Then Return 0
End If
Dim As Longint expo = 8 * Len(b.s) + 1022 ' = bits + expo_bias - 1
' if needed for the conversion, extend tail with two LS blocks of zero
If Len(b.s) < 12 Then b.s = Chr(0,0,0,0, 0,0,0,0) + b.s
' the ms block still contains the data, so no change to expo
Dim As Ubyte Ptr ubp = Strptr(b.s) + Len(b.s) - 1 ' point to the MSbyte
Dim As Integer i
For i = 0 To 4  ' find the leading non-zero byte, MS block may be zero
If *ubp > 0 Then Exit For
ubp = ubp - 1
expo = expo - 8 ' expo reduction of 8 bits per zero byte skipped
Next i  ' ubp now points to the MS non-zero byte
uli = *Cast(Ulongint Ptr, ubp - 7)  ' normalize bytes, left justify
For i = 63 To 56 Step -1    ' find the MS set bit
If Bit(uli, i) Then Exit For
expo = expo - 1
Next i  ' i now points to MSbit
uli = uli Shr (i - 52)  ' shift right to put MSbit in bit 52
uli = Bitreset(uli, 52) ' kill only the implicit bit now in bit 52
uli = Sign_Bit Or (expo Shl 52) Or uli  ' build the double
Return *Cast(Double Ptr, @uli)  ' return the bit pattern as a double
End Function

Operator bigint.Cast() As Double    ' Cdbl
Return BigInt_2_Double(this)
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

'----------------------------------------------------------------
' convert a bigint to unsigned hexadecimal
Function show_uhex(Byref s As bigint) As String
Dim As Integer i
dim as bigint a = s
If 128 And a.s[Len(a.s)-1] Then     'bigint is negative
print "cannot convert negative to uniform"
Sleep : Stop
End If
Dim As String h     ' hex is big-endian
For i = Len(s.s)-1 To 0 Step -1
h = h + Hex(s.s[i], 2)
Next i
if len(h) <> 8 then h=ltrim(h,"00000000")
Return h
End Function

'two's compliment hex to bigint
Function hex_bigint(byval a as string) as bigint
dim c as bigint
dim b as integer
dim d as ulong
dim i as ushort
if (len(a) mod 8) <> 0 then
for i=1 to (8-(len(a) mod 8))
a="0" & a
next i
end if
b=valint("&h" & mid(a,1,8))     'can be positive or negative
c=b
if b >= 0 then
for i=9 to len(a) step 8
c=shift_left(c,32)
d=valuint("&h" & mid(a,i,8))
c+=d
next i
else
for  i=9 to len(a) step 8
c=shift_left(c,32)
d=valuint("&h" & mid(a,i,8))
c-=d
next i
end if
return c
End Function

'unsigned hex to bigint
Function uhex_bigint(byval a as string) as bigint
dim c as bigint
dim b as ulong
dim i as ulong
if (len(a) mod 8) <> 0 then
for i=1 to (8-(len(a) mod 8))
a="0" & a
next i
end if
for i=1 to len(a) step 8
c=shift_left(c,32)
b=valuint("&h" & mid(a,i,8))
c+=b
next i
return c
End Function

'bigint to twos compliment binary
Function bigint_bin(byref a as bigint) as string
Dim As String s = a.s
return s
End Function

'twos compliment binary to bigint
Function bin_bigint(byref a as string) as bigint
dim as bigint b
b.s = a
return b
End Function

'bigint to unsigned binary
Function bigint_ubin(byref a as bigint) as string
Dim As String s
If 128 And a.s[Len(a.s)-1] Then     'bigint is negative
print "cannot convert negative to unsigned"
Sleep : Stop
Else
s = a.s
End If
s=rtrim(s,chr(0,0,0,0))
return s
End Function

'unsigned binary to bigint
Function ubin_bigint(byref a as string) as bigint
Dim as bigint b
b.s=a
If (128 And b.s[Len(b.s)-1]) Then b.s += Chr(0,0,0,0) ' make it positive
return b
End Function

'===============================================================
'       E n d    o f    B i g    I n t e g e r    C o d e
'===============================================================