Arbitrarily Big Integer Routines

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

Arbitrarily Big Integer Routines

Postby Richard » Mar 17, 2010 2:36

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)
    cc = addition(aa, cc)
    Return cc
End Function

'================================================================
' multiply
Function multiply(Byref aa As String, Byref bb As String) As String
    '------------------------------------------------------------
    ' sort out the signs and rectify the inputs
    Dim As String a = aa, b = bb
    Dim As Integer sign_a, sign_b, sign_c
    sign_a = a[Len(a)-1] And 128
    sign_b = b[Len(b)-1] And 128
    sign_c = sign_a Xor sign_b
    If sign_a Then a = negate(a)
    If sign_b Then b = negate(b)
    '------------------------------------------------------------
    ' find the dimensions of the problem
    Dim As Integer i, j, asize, bsize
    asize = Len(a) ' number of bytes in a
    bsize = Len(b) ' number of bytes in b
    Dim As String c = String(asize + bsize, Chr(0)) ' initialise accumulator
    asize = asize Shr 2 ' number of blocks in a
    bsize = bsize Shr 2 ' Shr 2 is faster than \4
    '------------------------------------------------------------
    ' pointers into all the strings
    Dim As Uinteger Ptr ia, ib, ic, iaz, ibz, icz
    iaz = Cast(Uinteger Ptr, Strptr(a) )
    ibz = Cast(Uinteger Ptr, Strptr(b) )
    icz = Cast(Uinteger Ptr, Strptr(c) )
    Dim As Ulongint product
    Dim As Uinteger Ptr uipcarry, uipdata = Cast(Uinteger Ptr, @product)
    uipcarry = uipdata + 1  ' top half of product is carry
    '------------------------------------------------------------
    ia = iaz
    For i = 0 To asize - 1
        ib = ibz
        product = 0 ' clear carry
        For j = 0 To bsize - 1
            ic = icz + i + j
            product = Culngint(*ia) * Culngint(*ib) + Culngint(*ic) + Culngint(*uipcarry)
            *ic = *uipdata
            ib += 1
        Next j
        *(ic+1) = *uipcarry
        ia += 1
    Next i
    '------------------------------------------------------------
    If sign_c = 128 Then c = negate(c)
    prune(c)
    Return c
End Function   

'   digits    time
'      10    3.9 usec
'     100   11.5 usec
'      1k   640. usec
'     10k    64. msec
'    100k   6.35 sec

'=======================================================================
' square a number, approaches twice the speed of multiply for big numbers
Function square(Byref aa As String) As String
    '------------------------------------------------------------
    Dim As String a = aa
    If (128 And a[Len(a)-1]) Then a = negate(a)
    '------------------------------------------------------------
    ' find the dimension of the problem
    Dim As Integer i, j, asize = Len(a) ' number of bytes in a
    Dim As String c = String(2 * asize, Chr(0) ) ' initialise accumulator
    asize = (asize Shr 2) - 1   ' highest block number in a
    '------------------------------------------------------------
    ' pointers into all the strings
    Dim As Uinteger Ptr ia, ib, ic, iaz, icz
    iaz = Cast(Uinteger Ptr, Strptr(a) )  ' the base addresses of strings
    icz = Cast(Uinteger Ptr, Strptr(c) )
    Dim As Ulongint product ' bottom half is data, top half will be carry
    Dim As Uinteger Ptr uipcarry, uipdata = Cast(Uinteger Ptr, @product)
    uipcarry = uipdata + 1  ' top half of product is carry
    '------------------------------------------------------------
    ' multiply one triangular corner only
    ia = iaz
    For i = 1 To asize
        ia += 1 ' ia starts at 1 not zero because 0,0 is on the diagonal
        ib = iaz    ' ib is second pointer into a, ib starts at zero
        ic = icz + i        ' ic is the accumulator ic = icz + i + j
        product = 0     ' clear carry
        For j = 0 To i - 1
            product = Culngint(*ia) * Culngint(*ib) + Culngint(*ic) + Culngint(*uipcarry)
            *ic = *uipdata
            ib += 1
            ic += 1
        Next j
        *ic = *uipcarry     ' line of blocks gets one longer each loop
    Next i
    '------------------------------------------------------------
    ' double the triangle, cannot do it at same time as squaring diagonal
    ic = icz        ' because it can cause the product to overflow
    product = 0 ' clear carry
    For i = 0 To (2 * asize) + 1
        product = Culngint(*ic) + Culngint(*ic) + Culngint(*uipcarry)
        *ic = *uipdata
        ic += 1
    Next i
    '------------------------------------------------------------
    ' square and accumulate the diagonal elements
    ia = iaz
    ic = icz
    product = 0     ' clear carry
    For i = 0 To asize
        ' square the diagonal entry, while propagating carry
        product = Culngint(*ia) * Culngint(*ia) + Culngint(*ic) + Culngint(*uipcarry)
        *ic = *uipdata
        ic += 1
        ' propagate carry through the following odd block of C
        product = Culngint(*ic) + Culngint(*uipcarry)
        *ic = *uipdata
        ic += 1
        ia += 1
    Next i
    '------------------------------------------------------------
    If 128 And c[Len(c)-1] Then c = c + Chr(0,0,0,0)   ' sign is positive
    prune(c)
    Return c
End Function
'   digits   multiply(x,x)  square(x)
'      10      3.9 usec     3.9 usec
'     100     11.5 usec     9.0 usec
'      1k      640 usec     340 usec
'     10k       64 msec      33 msec
'    100k     6.35 sec     3.21 sec

'=======================================================================
' shift up one bit, low towards high
Function mul2(Byref a As String) As String
    Dim As Integer i, sign, sum, carry = 0
    Dim As String b = a
    sign = b[Len(b) - 1] And 128    ' snag the msb of highest byte
    For i = 0 To Len(b) - 1
        sum = b[i] + b[i] + carry
        carry = Hibyte(sum)
        b[i] = Lobyte(sum)
    Next i
    If ( b[Len(b) - 1] And 128 ) <> sign Then
        carry = carry * 255
        b = b + Chr(carry,carry,carry,carry)    ' sign extend four bytes
    End If
    Return b
End Function

'=======================================================================
' shift down one bit, high towards low
Function div2(Byref a As String) As String
    Dim As Integer i, carry = 0
    Dim As String b = a
    For i = 0 To Len(a)-2   ' all except the top byte of four
        b[i] = ( b[i] Shr 1 ) + (128 * (b[i+1] And 1))
    Next i
    i = Len(b) - 1
    b[i] = b[i] Shr 1
    b[i] = b[i] Or (2 * ( b[i] And 64)) ' sign extend the msb
    prune(b)
    Return b
End Function

'=======================================================================
' raise a number to a positive power
Function power(Byref x As String, Byval n As Ulongint) As String
    Dim As Integer i = 64
    Do  ' find first set bit
        i = i - 1
    Loop Until Bit(n, i) Or (i = 0)
    i = i + 1
    Dim As String pwr = Chr(1,0,0,0)    ' one
    Do
        i = i - 1
        pwr = square(pwr)   ' this was a multiply but square is faster
        If Bit(n, i) Then pwr = multiply(pwr, x)
    Loop Until i = 0
    Return pwr  ' pwr was pruned by square and by multiply
End Function
' times with     multiply    square     prune&square
'  2 ^ 3          25 usec   13.5 usec   27.4 usec
'  2 ^ 35        300 usec    147 usec   46.4 usec
'123456789^127  1.43 msec    780 usec    260 usec
'  10 ^ 1000      82 msec     41 msec    186 usec
'  10 ^ 10000     14 sec       7 sec    10.9 msec
'  10 ^ 100000                          1.07 sec

'=======================================================================
Declare Function is_LT(Byref As String, Byref As String) As Integer

'=======================================================================
Function Factorial(Byref a As String) As String
    Dim As String f = Chr(1,0,0,0), n = a
    Do Until is_LT(n, Chr(2,0,0,0))
        f = multiply(f, n)
        n = subtract(n, Chr(1,0,0,0))
    Loop
    Return f
End Function

'=======================================================================
' bit functions
'=======================================================================
' NOT. Invert all the bits in a string
Function complement(Byref aa As String) As String
    Dim As String a = aa
    For i As Integer = 0 To Len(a)-1
        a[i] = 255 - a[i]
    Next i
    Return a
End Function

'=======================================================================
' get the value of a specified bit in a big integer
Function Bit_Value(Byref v As String, Byval b As Ulongint) As Integer
    Dim As Integer bitval, by = b \ 8
    If by < Len(v) Then
        bitval = Bit(v[by], b Mod 8)
    Else
        If v[Len(v)-1] And 128 Then bitval = -1 ' the extended sign bit
    End If
    Return bitval
End Function

'================================================================
' set a specified bit in a big integer
Function Bit_Set(Byref vv As String, Byval b As Ulongint) As String
    Dim As String v = vv
    Dim As Integer by, bi, delta, sign
    by = b \ 8      ' byte number
    bi = b Mod 8    ' bit number
    delta = by - Len(v) + 1
    If bi = 7 Then delta = delta + 1    ' protect the sign bit
    If delta > 0 Then    ' lengthen the number
        delta = ((delta + 3)\ 4) * 4
        If v[Len(v)-1] And 128 Then sign = -1 ' the extended sign bit
        v =  v + String(delta, sign)
    End If
    v[by] = Bitset(v[by], bi)
    Return v
End Function

'================================================================
' clear a specified bit in a big integer
Function Bit_Reset(Byref vv As String, Byval b As Ulongint) As String
    Dim As String v = vv
    Dim As Integer by, bi, delta, sign
    by = b \ 8      ' byte number
    bi = b Mod 8    ' bit number
    delta = by - Len(v) + 1
    If bi = 7 Then delta = delta + 1    ' protect the sign bit
    If delta > 0 Then    ' lengthen the number
        delta = ((delta + 3)\ 4) * 4
        If v[Len(v)-1] And 128 Then sign = -1 ' the extended sign bit
        v =  v + String(delta, sign)
    End If
    v[by] = Bitreset(v[by], bi)
    Return v
End Function

'================================================================
' shift string n bits left
Function shift_left(Byref a As String, Byval n As Uinteger) As String
    If n = 0 Then Return a
    Dim As Integer nblocks = n \ 32, nbits = n Mod 32
    Dim As String s = String(nblocks * 4, Chr(0)) + a ' put zeros on the rhs
    Dim As String m = bit_set(Chr(0,0,0,0), nbits)
    s = multiply(m, s)
    Return s
End Function

'================================================================
' shift string n bits right, by shifting left nbits and right nblocks
Function shift_right(Byref a As String, Byval n As Uinteger) As String
    If n = 0 Then Return a
    If n > (8*Len(a)) Then Return Chr(0,0,0,0)
    Dim As Integer nblocks = n \ 32, nbits = n Mod 32
    Dim As String s = a
    Dim As String m = bit_set(Chr(0,0,0,0), 32 - nbits )
    s = multiply(m, s)  ' move bits left
    s = Right(s, Len(s) - (nblocks+1)*4 )
    Return s
End Function

'================================================================
' bitwise AND
Function Bit_And(Byref aa As String, Byref bb As String) As String
    Dim As String a = aa, b = bb
    Dim As Integer lena, lenb, i
    lena = Len(a)
    lenb = Len(b)
    If lena > lenb Then
        b = b + String(lena - lenb, Bit(b[lenb - 1], 7))
    Else
        a = a + String(lenb - lena, Bit(a[lena - 1], 7))
    End If
    Dim As String c = b
    For i = 0 To Len(c) - 1
        c[i] = c[i] And a[i]
    Next i
    Return c
End Function

'================================================================
' bitwise Or
Function Bit_Or(Byref aa As String, Byref bb As String) As String
    Dim As String a = aa, b = bb
    Dim As Integer lena, lenb, i
    lena = Len(a)
    lenb = Len(b)
    If lena > lenb Then
        b = b + String(lena - lenb, Bit(b[lenb - 1], 7))
    Else
        a = a + String(lenb - lena, Bit(a[lena - 1], 7))
    End If
    Dim As String c = b
    For i = 0 To Len(c) - 1
        c[i] = c[i] Or a[i]
    Next i
    Return c
End Function

'================================================================
' bitwise Xor
Function Bit_Xor(Byref aa As String, Byref bb As String) As String
    Dim As String a = aa, b = bb
    Dim As Integer lena, lenb, i
    lena = Len(a)
    lenb = Len(b)
    If lena > lenb Then
        b = b + String(lena - lenb, Bit(b[lenb - 1], 7))
    Else
        a = a + String(lenb - lena, Bit(a[lena - 1], 7))
    End If
    Dim As String c = b
    For i = 0 To Len(c) - 1
        c[i] = c[i] Xor a[i]
    Next i
    Return c
End Function

'================================================================
' bitwise Eqv,     equivalence is the complement of Xor
Function Bit_Eqv(Byref aa As String, Byref bb As String) As String
    Dim As String a = aa, b = bb
    Dim As Integer lena, lenb, i
    lena = Len(a)
    lenb = Len(b)
    If lena > lenb Then
        b = b + String(lena - lenb, Bit(b[lenb - 1], 7))
    Else
        a = a + String(lenb - lena, Bit(a[lena - 1], 7))
    End If
    Dim As String c = b
    For i = 0 To Len(c) - 1
        c[i] = (c[i] Eqv a[i])
    Next i
    Return c
End Function

'================================================================
' bitwise Imp,     implication
Function Bit_Imp(Byref aa As String, Byref bb As String) As String
    Dim As String a = aa, b = bb
    Dim As Integer lena, lenb, i
    lena = Len(a)
    lenb = Len(b)
    If lena > lenb Then
        b = b + String(lena - lenb, Bit(b[lenb - 1], 7))
    Else
        a = a + String(lenb - lena, Bit(a[lena - 1], 7))
    End If
    Dim As String c = b
    For i = 0 To Len(c) - 1
        c[i] = (c[i] Imp a[i])
    Next i
    Return c
End Function

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

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

'================================================================
' convert a string to octal
Function show_oct(Byref a As String) As String
    Dim As String b, c, s = a
    If 128 And a[Len(a)-1] Then ' extend the sign
        s = s + Chr(255,255,255,255)
    Else
        s = s + Chr(0,0,0,0)
    End If
    Dim As Integer i
    Dim As Ulongint u
    For i = 1 To Len(a) Step 3
        b = Mid(s, i, 3)    ' take three bytes = 24 bits
        u = b[0] + 256 * (b[1] + 256 * b[2])
        c = Oct(u, 8) + c ' eight symbols per three bytes
    Next i
    Return c
End Function

'================================================================
' Relational functions
'----------------------------------------------------------------
Function is_ZERO(Byref a As String) As Integer
    If Len(Trim(a, Chr(0))) Then Return 0 Else Return -1
End Function

'----------------------------------------------------------------
Function is_NZ(Byref a As String) As Integer
    If Len(Trim(a, Chr(0))) Then Return -1 Else Return 0
End Function

'----------------------------------------------------------------
Function is_POS(Byref a As String) As Integer
    ' note that zero is deemed to be positive
    If (128 And a[Len(a)-1]) Then Return 0 Else Return -1
End Function

'----------------------------------------------------------------
Function is_NEG(Byref a As String) As Integer
    If (128 And a[Len(a)-1]) Then Return -1 Else Return 0
End Function

'----------------------------------------------------------------
Function is_EQ(Byref a As String, Byref b As String) As Integer
    If Len(Trim(subtract(a, b), Chr(0))) Then Return 0 Else Return -1
End Function

'----------------------------------------------------------------
Function is_NEQ(Byref a As String, Byref b As String) As Integer
    If Len(Trim(subtract(a, b), Chr(0))) Then Return -1 Else Return 0
End Function

'----------------------------------------------------------------
Function is_LT(Byref a As String, Byref b As String) As Integer
    Dim As String c = subtract(a, b)
    If (128 And c[Len(c)-1]) Then Return -1 Else Return 0
End Function

'----------------------------------------------------------------
Function is_GT(Byref a As String, Byref b As String) As Integer
    Dim As String c = subtract(b, a)
    If (128 And c[Len(c)-1]) Then Return -1 Else Return 0
End Function

'----------------------------------------------------------------
Function is_LT_EQ(Byref a As String, Byref b As String) As Integer
    Dim As String c = subtract(b, a)
    If (128 And c[Len(c)-1]) Then Return 0 Else Return -1
End Function

'----------------------------------------------------------------
Function is_GT_EQ(Byref a As String, Byref b As String) As Integer
    Dim As String c = subtract(a, b)
    If (128 And c[Len(c)-1]) Then Return 0 Else Return -1
End Function

'=======================================================================
Function BigInt_Sgn(Byref a As String) As Integer
    Dim As Integer i = 128 And a[Len(a)-1]
    If i Then Return -1 ' is negative
    If is_zero(a) Then Return 0 ' is zero
    Return 1 ' is positive
End Function

'=======================================================================
' integer divide, a \ b, return div and mod
Sub divmod(_
    Byref aa As String, Byref bb As String,_
    Byref intdiv As String, Byref intmod As String)
    Dim As String a = aa, b = bb, c
    If is_zero(b) Then
        Print " Division by zero. "
        Stop
    End If
    Dim As Integer lena = Len(a), lenb = Len(b)
    If (lena <= 8) And (lenb <= 8) Then ' arguments are one or two blocks
        intdiv = Chr(0,0,0,0,0,0,0,0)   ' use 64 bit long integers
        intmod = Chr(0,0,0,0,0,0,0,0)   ' allocate space for results
        If lena = 4 Then a = a + String(4, Bit(a[lena - 1], 7))
        If lenb = 4 Then b = b + String(4, Bit(b[lenb - 1], 7))
        Dim As Longint Ptr pa, pb, pc, pd
        pa = Cast(Longint Ptr, Strptr(a))
        pb = Cast(Longint Ptr, Strptr(b))
        pc = Cast(Longint Ptr, Strptr(intdiv))
        pd = Cast(Longint Ptr, Strptr(intmod))
        *pc = *pa \ *pb
        *pd = *pa Mod *pb
        '------------------------------------------------------------
    Else    ' more than two blocks so do a long division
        '------------------------------------------------------------
        ' sign of inputs and output
        Dim As Integer sa, sb, sq
        sa = 128 And a[Len(a)-1]
        sb = 128 And b[Len(b)-1]
        sq = sa Xor sb  ' sign of the result
        ' rectify the inputs
        If sa Then a = negate(a)
        If sb Then b = negate(b)
        '------------------------------------------------------------
        ' prepare to normalise the divisor
        Dim As Integer ja = msbit(a), jb = msbit(b), bitshifts =  ja - jb
        If ja < jb Then ' abandon the division as result is zero
            intdiv = Chr(0,0,0,0)
            intmod = aa
        Else    ' proceed with the long slow division
            b = shift_Left(b, bitshifts)    ' align the first bits
            lena = Len(a)
            lenb = Len(b)
            If lena <> lenb Then
                If lena > lenb Then
                    b = b + String(lena - lenb, Bit(b[lenb - 1], 7) )
                Else
                    a = a + String(lenb - lena, Bit(b[lenb - 1], 7) )
                    lena = Len(a)   ' update length of lena only
                End If
            End If
            '------------------------------------------------------------
            Dim As String qu = String(lena, Chr(0))  ' quotient
            Dim As Integer i, j, qbit, blocks = lena Shr 2
            c = String(lena, Chr(0) )   ' prepare to receive the result
            '------------------------------------------------------------
            ' setup all pointers for 32 bit block arithmetic
            Dim As Uinteger Ptr pa, pb, pc, pcarry, pdata
            Dim As Ulongint sum
            pdata = Cast(Uinteger Ptr, @sum)  ' bottom half of sum is data
            pcarry = pdata + 1   ' point to carry in top half of sum
            For i = 0 To bitshifts
                '------------------------------------
                ' setup and perform subtraction
                pa = Cast(Uinteger Ptr, Strptr(a))
                pb = Cast(Uinteger Ptr, Strptr(b))  ' strings can move
                pc = Cast(Uinteger Ptr, Strptr(c))
                j = blocks
                *pcarry = 1
                Do  ' inline subtract, c = a - b
                    sum = Culngint(*pa) + Culngint(Not(*pb)) + Culngint(*pcarry)
                    *pc = *pdata ' save this result
                    pa += 1 ' advance all three pointers
                    pb += 1
                    pc += 1
                    j = j - 1  ' count the 32 bit blocks
                Loop Until j = 0
                '------------------------------------
                ' test sign of result
                If *pcarry Then
                    ' result of subtraction was positive
                    pa = Cast(Uinteger Ptr, Strptr(a))
                    pc = Cast(Uinteger Ptr, Strptr(c))  ' strings can move
                    For j = 1 To blocks ' copy c to a
                        *pa = *pc
                        pa += 1
                        pc += 1
                    Next j
                    qbit = 1
                Else ' result of subtraction was negative
                    qbit = 0
                End If
                '------------------------------------
                ' shift quotient one bit left, insert qbit on right
                Dim As Uinteger sum16
                For j = 0 To Len(b) - 1
                    sum16 = qu[j] + qu[j] + qbit
                    qbit = Hibyte(sum16)
                    qu[j] = Lobyte(sum16)
                Next j
                '------------------------------------
                ' shift b one bit right
                For j = 0 To Len(a)-2   ' all except the top byte
                    b[j] = ( b[j] Shr 1 ) + (128 * (b[j+1] And 1))
                Next j
                b[lena - 1] = b[lena - 1] Shr 1
                '------------------------------------
            Next i  ' do next bit
            prune(a)    ' finished division
            prune(qu)
            If sq Then qu = negate(qu)  ' Xor of input signs
            If sa Then a = negate(a)    ' sign of original input A
            intdiv = qu ' the quotient
            intmod = a  ' the remainder
            '------------------------------------------------------------
        End If
    End If
End Sub
' 100 digits divided by 1 takes 775. usec

'===============================================================
'' divide a by b, return the integer result
Function divide(Byref a As String, Byref b As String) As String
    Dim As String c, d
    divmod(a, b, c, d)
    Return c
End Function

'===============================================================
' divide a by b and return only the integer remainder
Function remainder(Byref a As String, Byref b As String) As String
    Dim As String c, d
    divmod(a, b, c, d)
    Return d
End Function

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

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

Postby Ophelius » Mar 17, 2010 3:09

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

Postby PyRoe » Mar 17, 2010 9:01

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:

Postby duke4e » Mar 17, 2010 11:18

Great! Can't wait for overloaded version
integer
Posts: 376
Joined: Feb 01, 2007 16:54
Location: usa

Cannot grasp what is wrong inside the function

Postby integer » Jun 05, 2011 23:05

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: 2907
Joined: Jan 15, 2007 20:44
Location: Australia

Postby Richard » Jun 06, 2011 1:14

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)
    cc = addition(aa, cc)
    Return cc
End Function

'=======================================================================
' square a number, approaches twice the speed of multiply for big numbers
Function square(Byref aa As bigint) As bigint
    '------------------------------------------------------------
    Dim As bigint a = aa, c
    If (128 And a.s[Len(a.s)-1]) Then a = negate(a)
    '------------------------------------------------------------
    ' find the dimension of the problem
    Dim As Integer i, j, asize = Len(a.s) ' number of bytes in a
    c.s = String(2 * asize, Chr(0) ) ' initialise accumulator
    asize = (asize Shr 2) - 1   ' highest block number in a
    '------------------------------------------------------------
    ' pointers into all the bigints
    Dim As Uinteger 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

'=====================================================================
' Overload arithmetic operators
'=====================================================================
' unary plus +
Operator + (Byref x As bigint) As bigint
Return x
End Operator

' unary minus -
Operator - (Byref x As bigint) As bigint
Return negate(x)
End Operator

' addition
Operator + (Byref x As bigint, Byref y As bigint) As bigint
Return addition(x, y)
End Operator

' subtraction
Operator - (Byref x As bigint, Byref y As bigint) As bigint
Return subtract(x, y)
End Operator

Operator * (Byref x As bigint, Byref y As bigint) As bigint
Return multiply(x, y)
End Operator

Operator \ (Byref x As bigint, Byref y As bigint) As bigint
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)
this = addition(this, rhs)
End Operator

Operator bigint.-= (Byref rhs As bigint)
this = subtract(this, rhs)
End Operator

Operator bigint.*= (Byref rhs As bigint)
this = multiply(this, rhs)
End Operator

Operator bigint.\= (Byref rhs As bigint)
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: 2907
Joined: Jan 15, 2007 20:44
Location: Australia

Postby Richard » Jun 06, 2011 1:58

Replacing your Byval with Byref fixes the problem with your code;
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: 376
Joined: Feb 01, 2007 16:54
Location: usa

Postby integer » Jun 06, 2011 4:54

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!
I like the overloaded version.

Thanks.
integer
Posts: 376
Joined: Feb 01, 2007 16:54
Location: usa

Postby integer » Jun 25, 2011 3:17

@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: 2907
Joined: Jan 15, 2007 20:44
Location: Australia

Postby Richard » Jun 25, 2011 6:32

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: 2907
Joined: Jan 15, 2007 20:44
Location: Australia

Postby Richard » Jul 01, 2011 3:49

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: 2907
Joined: Jan 15, 2007 20:44
Location: Australia

Postby Richard » Jul 03, 2011 3:45

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)
    cc = addition(aa, cc)
    Return cc
End Function

'=======================================================================
' square a number, approaches twice the speed of multiply for big numbers
Function square(Byref aa As bigint) As bigint
    '------------------------------------------------------------
    Dim As bigint a = aa, c
    If (128 And a.s[Len(a.s)-1]) Then a = negate(a)
    '------------------------------------------------------------
    ' find the dimension of the problem
    Dim As Integer i, j, asize = Len(a.s) ' number of bytes in a
    c.s = String(2 * asize, Chr(0) ) ' initialise accumulator
    asize = (asize Shr 2) - 1   ' highest block number in a
    '------------------------------------------------------------
    ' pointers into all the bigints
    Dim As Uinteger 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

'=====================================================================
' Overload arithmetic operators
'=====================================================================
' unary plus +
Operator + (Byref x As bigint) As bigint
Return x
End Operator

' unary minus -
Operator - (Byref x As bigint) As bigint
Return negate(x)
End Operator

' addition
Operator + (Byref x As bigint, Byref y As bigint) As bigint
Return addition(x, y)
End Operator

' subtraction
Operator - (Byref x As bigint, Byref y As bigint) As bigint
Return subtract(x, y)
End Operator

Operator * (Byref x As bigint, Byref y As bigint) As bigint
Return multiply(x, y)
End Operator

Operator \ (Byref x As bigint, Byref y As bigint) As bigint
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)
this = addition(this, rhs)
End Operator

Operator bigint.-= (Byref rhs As bigint)
this = subtract(this, rhs)
End Operator

Operator bigint.*= (Byref rhs As bigint)
this = multiply(this, rhs)
End Operator

Operator bigint.\= (Byref rhs As bigint)
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
Site Admin
Posts: 6166
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs

Postby counting_pine » Jul 05, 2011 16:25

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

Postby stephanbrunker » Nov 02, 2013 15:26

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

Postby stephanbrunker » Nov 03, 2013 13:07

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
'=============================
'removed comments @shift_right
'bugfix constructor and operator uinteger, added ulong
'conversion functions hex_bigint, uhex_bigint,
'bin_bigint, bigint_bin;  ubin_bigint, bigint_ubin

' 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)
    cc = addition(aa, cc)
    Return cc
End Function

'=======================================================================
' square a number, approaches twice the speed of multiply for big numbers
Function square(Byref aa As bigint) As bigint
    '------------------------------------------------------------
    Dim As bigint a = aa, c
    If (128 And a.s[Len(a.s)-1]) Then a = negate(a)
    '------------------------------------------------------------
    ' find the dimension of the problem
    Dim As Integer i, j, asize = Len(a.s) ' number of bytes in a
    c.s = String(2 * asize, Chr(0) ) ' initialise accumulator
    asize = (asize Shr 2) - 1   ' highest block number in a
    '------------------------------------------------------------
    ' pointers into all the bigints
    Dim As Uinteger 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

'=====================================================================
' Overload arithmetic operators
'=====================================================================
' unary plus +
Operator + (Byref x As bigint) As bigint
Return x
End Operator

' unary minus -
Operator - (Byref x As bigint) As bigint
Return negate(x)
End Operator

' addition
Operator + (Byref x As bigint, Byref y As bigint) As bigint
Return addition(x, y)
End Operator

' subtraction
Operator - (Byref x As bigint, Byref y As bigint) As bigint
Return subtract(x, y)
End Operator

Operator * (Byref x As bigint, Byref y As bigint) As bigint
Return multiply(x, y)
End Operator

Operator \ (Byref x As bigint, Byref y As bigint) As bigint
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)
this = addition(this, rhs)
End Operator

Operator bigint.-= (Byref rhs As bigint)
this = subtract(this, rhs)
End Operator

Operator bigint.*= (Byref rhs As bigint)
this = multiply(this, rhs)
End Operator

Operator bigint.\= (Byref rhs As bigint)
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
    'pad leading zeroes for positive hex
    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
    'pad leading zeroes
    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
'===============================================================

Return to “Projects”

Who is online

Users browsing this forum: No registered users and 5 guests