## Arbitrarily Big Integer Routines

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

### Arbitrarily Big Integer Routines

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

Save this file as “Big_Integer_0.bas”

Code: Select all

`' version 1.0   17 March 2010'================================================================' Big_Integer_0.bas  Arbitrary Precision, Two's Complement Integer.'================================================================' supported functions are;' (short constants -1, 0, 1, 2 ) Bigint_m1, Bigint_0, Bigint_1, Bigint_2' Msbit, Prune, Neg, Absolute, pack, unpack, Addition, Subtract' Multiply, Square, Mul2, Div2, Power, Factorial' Shift_Left, Shift_Right, Complement, Bit_value, Bit_set, Bit_reset' Bit_And, Bit_Or, Bit_Xor, Bit_Eqv, Bit_Imp, Show_Bin, Show_Hex, Show_Oct' is_ZERO, is_NZ, is_POS, is_NEG' is_EQ, is_NEQ, is_LT, is_GT, is_LT_EQ, is_GT_EQ' BigInt_Sgn, DivMod, Divide, Remainder'----------------------------------------------------------------' This package only handles integers encoded in a two's complement format.' The first bit in the two's complement format is always the sign which' takes the value of 1 = negative, 0 = positive or zero. Byte examples are;' +127 = 0111 1111  most positive'   +8 = 0000 1000'   +7 = 0000 0111'   +4 = 0000 0100'   +2 = 0000 0010'   +1 = 0000 0001'    0 = 0000 0000  zero'   -1 = 1111 1111'   -2 = 1111 1110'   -4 = 1111 1100'   -7 = 1111 1001'   -8 = 1111 1000' -127 = 1000 0001' -128 = 1000 0000  most negative'----------------------------------------------------------------' Each successive 32 bits are packed into a 4 byte block.' Each block is stored in 4 successive bytes of a string.'----------------------------------------------------------------' Strings in FreeBASIC are indexed and printed from left to right. Big' integers appear to be stored in strings backwards which can be confusing.' The Left side of the string is the right side of the number and vice versa.' Beware: Shift_Left moves a string to the right, Shift_Right moves it left.'----------------------------------------------------------------' The first byte in the string is the least significant byte of the number.' The last block in the string is the most significant block of the number.' String s indexing has s[0] as the LS byte and s[len(s)-1] as the MS byte.' These big integer strings are always multiples of 4 bytes long.' The msb of a number is sign extended to the MS bit of the MS block.'----------------------------------------------------------------' A number is always stored in a way that correctly represents that number.' Where an operation would overflow into the MSB, a sign block is' appended to the number so as to prevent overflow or sign change.' Unnecessary leading zero or all ones blocks are removed by prune.'----------------------------------------------------------------' String pointers can change if a string is created or any length is changed.'================================================================'Type bigint     ' this will be used to overload FreeBASIC operators later'    s As String ' a little endian, two's complement number binary'End Type        ' number packed into a string in blocks of four bytes'================================================================' common small constantsDim Shared As String bigint_m1, bigint_0, bigint_1, bigint_2bigint_m1 = Chr(255,255,255,255)' minus 1bigint_0 = Chr(0,0,0,0) ' zerobigint_1 = Chr(1,0,0,0) ' onebigint_2 = Chr(2,0,0,0) ' two'=======================================================================' find the bit position of the first bit that differs from the sign bitFunction 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 kEnd Function'-----------------------------------------------------------------------' remove unnecessary leading blocks, prune time easily earns it's keepSub 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 stringFunction 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 sEnd 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 positiveFunction absolute(Byref a As String) As String    Dim As String b = a    If 128 And b[Len(b)-1] Then b = negate(b)    Return bEnd 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 bytesFunction 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 bEnd 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 stringFunction 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 dsEnd 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 aEnd Function' 100 digits in 4.2 usec'===============================================================' subtractFunction subtract(Byref aa As String, Byref bb As String) As String    Dim As String cc = Negate(bb)    cc = addition(aa, cc)    Return ccEnd Function'================================================================' multiplyFunction 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 cEnd 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 numbersFunction 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 cEnd 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 highFunction 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 bEnd Function'=======================================================================' shift down one bit, high towards lowFunction 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 bEnd Function'=======================================================================' raise a number to a positive powerFunction 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 multiplyEnd 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 fEnd Function'=======================================================================' bit functions'=======================================================================' NOT. Invert all the bits in a stringFunction 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 aEnd Function'=======================================================================' get the value of a specified bit in a big integerFunction 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 bitvalEnd Function'================================================================' set a specified bit in a big integerFunction 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 vEnd Function'================================================================' clear a specified bit in a big integerFunction 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 vEnd Function'================================================================' shift string n bits leftFunction 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 sEnd Function'================================================================' shift string n bits right, by shifting left nbits and right nblocksFunction 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 sEnd Function'================================================================' bitwise ANDFunction 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 cEnd Function'================================================================' bitwise OrFunction 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 cEnd Function'================================================================' bitwise XorFunction 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 cEnd Function'================================================================' bitwise Eqv,     equivalence is the complement of XorFunction 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 cEnd Function'================================================================' bitwise Imp,     implicationFunction 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 cEnd Function'================================================================' convert a string to binaryFunction 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 hEnd Function'================================================================' convert a string to hexadecimalFunction 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 hEnd Function'================================================================' convert a string to octalFunction 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 cEnd Function'================================================================' Relational functions'----------------------------------------------------------------Function is_ZERO(Byref a As String) As Integer    If Len(Trim(a, Chr(0))) Then Return 0 Else Return -1End Function'----------------------------------------------------------------Function is_NZ(Byref a As String) As Integer    If Len(Trim(a, Chr(0))) Then Return -1 Else Return 0End 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 -1End Function'----------------------------------------------------------------Function is_NEG(Byref a As String) As Integer    If (128 And a[Len(a)-1]) Then Return -1 Else Return 0End 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 -1End 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 0End 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 0End 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 0End 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 -1End 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 -1End 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 positiveEnd Function'=======================================================================' integer divide, a \ b, return div and modSub 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 IfEnd Sub' 100 digits divided by 1 takes 775. usec'==============================================================='' divide a by b, return the integer resultFunction divide(Byref a As String, Byref b As String) As String    Dim As String c, d    divmod(a, b, c, d)    Return cEnd Function'===============================================================' divide a by b and return only the integer remainderFunction remainder(Byref a As String, Byref b As String) As String    Dim As String c, d    divmod(a, b, c, d)    Return dEnd 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 filedim as string a, b, c, da = pack("123456789012345678901234567890")print "  input ="; unpack(a)b = square(a)print " square ="; unpack(b)c = divide(b, a)print " divide ="; unpack(c)d = subtract(c, a)print "subtract="; unpack(d)sleep`

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

### This is freaking sweet

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

### Cannot grasp what is wrong inside the function

I am trying to use a power-mod function:

Code: Select all

`#include "big_integer_0.bas"'----------------------------------------------------------------------' raise a number to a positive power; take the number mod m'   result = b^e mod m         [base^exponent : modulus]'Function powermod(Byval Number As String, Byval Exponent As string, byval Modulus as string) As String        Dim As String pwr, t, b, e, m        b = Number    e = Exponent    m = Modulus           t = multiply(m,b)            print "inside function powermod()"    print "    b= "; unpack(b)    print "    e= "; unpack(e)    print "    m= "; unpack(m)    print "t=b*m= "; unpack(t)        t = multiply(m,e)    print " b=m*e ";unpack(t)    sleep        pwr = pack("1")        Do while is_NZ( e )        if Bit_Value(e,0) then             pwr = multiply(pwr,b)            pwr = remainder(pwr,m)        endif        e = div2(e)        b = square(b)        b = remainder(b,m)    Loop    Return pwr    End Function'=========================================='sample trial programdim 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= +2047s=b*m= +108491-------------------------inside function powermod()    b= +53    e= +1023    m= +2047t=b*m= +0    t = multiply(m,e)b=m*e  +0`

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

Code: Select all

`' version 1.1   6 June 2011'================================================================' Big_Integer_1.bas  Arbitrary Precision, Two's Complement Integer.'================================================================' supported functions are;' (short constants -1, 0, 1, 2 ) Bigint_m1, Bigint_0, Bigint_1, Bigint_2' Msbit, Prune, Neg, Absolute, pack, unpack, Addition, Subtract' Multiply, Square, Mul2, Div2, Power, Factorial' Shift_Left, Shift_Right, Complement, Bit_value, Bit_set, Bit_reset' Bit_And, Bit_Or, Bit_Xor, Bit_Eqv, Bit_Imp, Show_Bin, Show_Hex, Show_Oct' is_ZERO, is_NZ, is_POS, is_NEG, is_EQ, is_NEQ' is_LT, is_GT, is_LT_EQ, is_GT_EQ, BigInt_Sgn, DivMod, Divide, Remainder'----------------------------------------------------------------' This package only handles integers encoded in a two's complement format.' The first bit in the two's complement format is always the sign which' takes the value of 1 = negative, 0 = positive or zero. Byte examples are;' +127 = 0111 1111  most positive'   +8 = 0000 1000'   +7 = 0000 0111'   +4 = 0000 0100'   +2 = 0000 0010'   +1 = 0000 0001'    0 = 0000 0000  zero'   -1 = 1111 1111'   -2 = 1111 1110'   -4 = 1111 1100'   -7 = 1111 1001'   -8 = 1111 1000' -127 = 1000 0001' -128 = 1000 0000  most negative'----------------------------------------------------------------' Each successive 32 bits are packed into a 4 byte block.' Each block is stored in 4 successive bytes of a string.'----------------------------------------------------------------' Strings in FreeBASIC are indexed and printed from left to right. Big' integers appear to be stored in strings backwards which can be confusing.' The Left side of the string is the Right side of the number and vice versa.' Beware: Shift_Left moves a string to the right, Shift_Right moves it left.'----------------------------------------------------------------' The first byte in the string is the least significant byte of the number.' The last block in the string is the most significant block of the number.' String s indexing has s[0] as the LS byte and s[len(s)-1] as the MS byte.' These big integer strings are always multiples of 4 bytes long.' The msb of a number is sign extended to the MS bit of the MS block.'----------------------------------------------------------------' A number is always stored in a way that correctly represents that number.' Where an operation would overflow into the MSB, a sign block is' appended to the number so as to prevent overflow or sign change.' Unnecessary leading zero or all ones blocks are removed by prune.'----------------------------------------------------------------' String pointers may change if a string is created or any length is changed.'================================================================Type bigint     ' a little endian, two's complement binary number    s As String ' packed into a string, in blocks four bytes long    ' constructor    Declare Constructor ' default constructor    Declare Constructor (Byref As Integer)    Declare Constructor (Byref As Longint)    Declare Constructor (Byref As String)    ' Let    Declare Operator Let(Byref rhs As String)    Declare Operator Let(Byval rhs As Byte)    Declare Operator Let(Byval rhs As Ubyte)    Declare Operator Let(Byval rhs As Short)    Declare Operator Let(Byval rhs As Ushort)    Declare Operator Let(Byval rhs As Integer)    Declare Operator Let(Byval rhs As Uinteger)    Declare Operator Let(Byval rhs As Longint)    Declare Operator Let(Byval rhs As Ulongint)    Declare Operator Let(Byval rhs As Double)    ' cast    Declare Operator Cast() As String    ' implicit step versions    Declare Operator For ()    Declare Operator Step()    Declare Operator Next(Byref end_cond As bigint) As Integer    ' explicit step versions    Declare Operator For (Byref step_var As bigint)    Declare Operator Step(Byref step_var As bigint)    Declare Operator Next(Byref end_cond As bigint, Byref step_var As bigint ) As Integer    ' operate and assign    Declare Operator += (Byref rhs As bigint)    Declare Operator -= (Byref rhs As bigint)    Declare Operator *= (Byref rhs As bigint)    Declare Operator \= (Byref rhs As bigint)End Type'================================================================' only the necessary declarations are hereDeclare Function pack(Byref As String) As bigintDeclare Function unpack(Byref As bigint) As StringDeclare Function negate(Byref As bigint) As bigintDeclare Operator < (Byref As bigint, Byref As bigint) As IntegerDeclare Function is_LT(Byref As bigint, Byref As bigint) As Integer'================================================================Constructor bigint  ' default constructorthis.s = Chr(0,0,0,0)   ' zeroEnd ConstructorConstructor bigint (Byref rhs As Integer)this = rhsEnd ConstructorConstructor bigint (Byref rhs As Longint)this = rhsEnd ConstructorConstructor bigint (Byref rhs As String)this = rhsEnd Constructor'================================================================' assignment operatorOperator 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 IfEnd 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 IfEnd 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 = *ipEnd 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 = *uipIf (128 And rhs) Then this.s += Chr(0,0,0,0) ' make it positiveEnd 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 = *lipEnd 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 = *ulipIf (128 And this.s[7]) Then this.s += Chr(0,0,0,0) ' make it positiveEnd Operator'----------------------------------------------------------------Operator bigint.let (Byval d As Double)Const As Ulongint implicit_bit = 2^52          ' 4503599627370496Const As Ulongint mant_mask = implicit_bit - 1 ' 4503599627370495Dim As Ulongint u, mantDim As Integer negative, sign, expoDim As bigint a'----------------------------------------------------If d < 0 Then negative = -1 ' remember signd = 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 bitexpo = (u Shr 52) And 2047  ' the 11 bit exponentmant = (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 smallElse    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 = aEnd IfEnd Operator'================================================================' create a string for printing a bigintOperator bigint.cast() As StringDim As String tt = Unpack(this)Return tEnd Operator'================================================================' common small constantsDim Shared As bigint bigint_m1, bigint_0, bigint_1, bigint_2bigint_m1.s = Chr(255,255,255,255)' minus 1bigint_0.s = Chr(0,0,0,0) ' zerobigint_1.s = Chr(1,0,0,0) ' onebigint_2.s = Chr(2,0,0,0) ' two'================================================================' find the bit position of the first bit that differs from the sign bitFunction 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 kEnd Function'-----------------------------------------------------------------------' remove unnecessary leading blocks, prune time easily earns it's keepSub 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 bigintFunction 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 sEnd 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 positiveFunction 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 bEnd 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 bytesFunction 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 bEnd 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 stringFunction 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 dEnd 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 aEnd Function' 100 digits in 4.2 usec'===============================================================' subtractFunction subtract(Byref aa As bigint, Byref bb As bigint) As bigint    Dim As bigint cc = Negate(bb)    cc = addition(aa, cc)    Return ccEnd Function'=======================================================================' square a number, approaches twice the speed of multiply for big numbersFunction 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 cEnd 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'================================================================' multiplyFunction 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 IfEnd 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 highFunction 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 bEnd Function'=======================================================================' shift down one bit, high towards lowFunction 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 bEnd Function'=======================================================================' raise a number to a positive powerFunction 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 multiplyEnd 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 fEnd Function'=======================================================================' bit functions'=======================================================================' NOT. Invert all the bits in a bigintFunction 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 aEnd Function'=======================================================================' get the value of a specified bit in a big integerFunction 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 bitvalEnd Function'================================================================' set a specified bit in a big integerFunction 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 vEnd Function'================================================================' clear a specified bit in a big integerFunction 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 vEnd Function'================================================================' shift bigint n bits leftFunction 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 sEnd 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 ANDFunction 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 cEnd Function'================================================================' bitwise OrFunction 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 cEnd Function'================================================================' bitwise XorFunction 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 cEnd Function'================================================================' bitwise Eqv,     equivalence is the complement of XorFunction 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 cEnd Function'================================================================' bitwise Imp,     implicationFunction 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 cEnd Function'================================================================' convert a bigint to binaryFunction 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 hEnd Function'================================================================' convert a bigint to hexadecimalFunction 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 hEnd Function'================================================================' convert a bigint to octalFunction 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 cEnd Function'================================================================' Relational functions'----------------------------------------------------------------Function is_ZERO(Byref a As bigint) As Integer    If Len(Trim(a.s, Chr(0))) Then Return 0 Else Return -1End 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 0End 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 positiveEnd 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 modSub 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 IfEnd Sub' 100 digits divided by 1 takes 775. usec'===============================================================' divide a by b, return the integer resultFunction divide(Byref a As bigint, Byref b As bigint) As bigint    Dim As bigint c, d    divmod(a, b, c, d)    Return cEnd Function'----------------------------------------------------------------' divide a by b and return only the integer remainderFunction remainder(Byref a As bigint, Byref b As bigint) As bigint    Dim As bigint c, d    divmod(a, b, c, d)    Return dEnd Function'=====================================================================' Overload arithmetic operators'=====================================================================' unary plus +Operator + (Byref x As bigint) As bigintReturn xEnd Operator' unary minus -Operator - (Byref x As bigint) As bigintReturn negate(x)End Operator' additionOperator + (Byref x As bigint, Byref y As bigint) As bigintReturn addition(x, y)End Operator' subtractionOperator - (Byref x As bigint, Byref y As bigint) As bigintReturn subtract(x, y)End OperatorOperator * (Byref x As bigint, Byref y As bigint) As bigintReturn multiply(x, y)End OperatorOperator \ (Byref x As bigint, Byref y As bigint) As bigintReturn divide(x, y)End OperatorOperator / (Byref x As bigint, Byref y As bigint) As bigintReturn divide(x, y) ' should this be prohibited ?End Operator'----------------------------------------------------------------' operate and assignOperator bigint.+= (Byref rhs As bigint)this = addition(this, rhs)End OperatorOperator bigint.-= (Byref rhs As bigint)this = subtract(this, rhs)End OperatorOperator bigint.*= (Byref rhs As bigint)this = multiply(this, rhs)End OperatorOperator 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 IntegerDim As bigint cc = subtract(a, b)If Len(Trim(c.s, Chr(0))) Then Return 0 Else Return -1End Operator'----------------------------------------------------------------Operator <> (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint cc = subtract(a, b)If Len(Trim(c.s, Chr(0))) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------Operator < (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint cc = subtract(a, b)If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------Operator > (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint c = subtract(b, a)If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------Operator <= (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint cc = subtract(b, a)If (128 And c.s[Len(c.s)-1]) Then Return 0 Else Return -1End Operator'----------------------------------------------------------------Operator >= (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint c = subtract(b, a)If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------' raise a bigint to the power of a long integer, result = x^y Operator ^ (Byref x As bigint, Byref y As Longint) As bigintReturn power(x, y)End Operator'=================================================================' more intrinsic numerical functions'=================================================================Function Abs_(Byref x As Double) As Double    Return Abs(x)End Function#undef AbsFunction Abs Overload (Byref x As Double) As Double    Return abs_(x)End FunctionFunction 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 SgnFunction Sgn Overload (Byref x As Double) As Double    Return Sgn_(x)End FunctionFunction 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 1Operator bigint.for( )End OperatorOperator bigint.step( )this += bigint_1End OperatorOperator bigint.next( Byref end_cond As bigint ) As IntegerReturn this <= end_condEnd Operator'----------------------------------------------------------------' explicit step versionsOperator bigint.for( Byref step_var As bigint )End OperatorOperator bigint.step( Byref step_var As bigint )this += step_var   End OperatorOperator bigint.next( Byref end_cond As bigint, Byref step_var As bigint ) As IntegerIf step_var < bigint_0 Then    Return this >= end_condElse    Return this <= end_condEnd IfEnd 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, stoptimeDim As bigint a, b = 2, cDim As String s1, s2Printdim as double bits, blocks, bytes, dig = 10000 ' decimal digitsprint dig; " decimal digits"bits = dig * 3.3219281blocks = int(bits / 32)bytes = blocks * 4a.s = string(bytes, chr(0)) + chr(1,0,0,0) print "length "; len(a.s); " bytes"startime = timerb = a * astoptime = timerprint " squareing ", 1e3*(stoptime - startime); " millisec"b = astartime = timerb = b * astoptime = timerprint "multiplying", 1e3*(stoptime - startime); " millisec"print "length "; len(b.s); " bytes"print "done"Sleep`
Richard
Posts: 2984
Joined: Jan 15, 2007 20:44
Location: Australia
I think due to Byval using Zstrings and finding a zero byte ?
Sorry about the delay in posting the latest Bigint, I am still not satisfied.
Bigint Overload seems to handle the Byval mode without the Zstring problem.
@ Integer. Here is your new code, does it now work as you would expect ?

Code: Select all

`#include "big_integer_1.bas"'----------------------------------------------------------------------' raise a number to a positive power; take the number mod m'   result = b^e mod m         [base^exponent : modulus]'----------------------------------------------------------------------' convert Byval to Byreffunction 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 programDim As Bigint b,e,m,r,sb = "53"e = "1023"m = "2047"s = b * mPrint "-------------------------"Print "before the function call"Print "    b= "; bPrint "    e= "; ePrint "    m= "; mPrint "s=b*m= "; sPrint "-------------------------"r = powermod( b, e, m )Print "b^e mod m= "; rSleep`
integer
Posts: 386
Joined: Feb 01, 2007 16:54
Location: usa
That is SUPER!
It works as expected.

Thank you, Richard.

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

P.S.

WOW!

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

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

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

This is my routine for that:

Code: Select all

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

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

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

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

Code: Select all

`#include "Big_Integer_1.bas"Function BigInt_2_LongInt ( Byref a As BigInt ) As Longint    Dim As String s = a.s    If Len(s) > 8 Then        Print " Overflow in BigInt to LongInteger conversion. "        Print a        Sleep : Stop    End If    If Len(s) = 4 Then  ' sign extend 4 byte integer to 8 byte LongInt        If Bit(s[3], 7) Then            s = s + String(4, Chr(255))        Else            s = s + String(4, Chr(0))        End If    End If    Return *Cptr(Longint Ptr, Strptr(s))End FunctionFunction 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 FunctionDim As BigInt xx = 3For i As Integer = 1 To 20    Print BigInt_2_Integer( x ),    Print x    x = x * 10Next iSleep`
Adulation is always welcome from an Integer.
It still surprises me that I wrote BigInt and that it actually works.
Richard
Posts: 2984
Joined: Jan 15, 2007 20:44
Location: Australia
Here is a fast BigInt to Double converter with some test code.

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

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

Code: Select all

`'===================================================================' Convert a BigInt to a double with 1.#INF overflow if too big'===================================================================#include "Big_Integer_1.bas"Function BigInt_2_Double ( Byref a As BigInt ) As Double    Dim As Bigint b = a    Dim As Ulongint uli, Sign_Bit    Dim As Integer i    If Bit( b.s[Len(b.s) - 1], 7) Then  ' extract the sign bit        Sign_Bit = Culngint(1) Shl 63   ' remember the sign         b = -b  ' rectify    End If  ' b is now a positive BigInt    ' overflows double? if mag > 1.797693134862310e308 = signed infinity    If Len(b.s) > 128 Then ' 128 bytes * 8 bits per byte = 1024 bits        If Len(b.s) = 132 Then  ' special case of sign bit = entire block            If ( b.s[131] Or b.s[130] Or b.s[129] Or b.s[128] ) <> 0 Then                uli = Sign_Bit Or &h7FF0000000000000 ' all ones exponent                Return *Cast(Double Ptr, @uli)  ' return the bit pattern            End If        End If    End If    If Len(b.s) = 4 Then    ' test for simple zero        If ( b.s[3] Or b.s[2] Or b.s[1] Or b.s[0] ) = 0 Then Return 0    End If    Dim As Longint expo = 8 * Len(b.s) + 1022 ' = bits - 1 + bias of expo    ' if needed for the conversion, extend tail with two zero LS blocks    If Len(b.s) < 12 Then b.s = Chr(0,0,0,0, 0,0,0,0) + b.s    ' the ms block still contains the data, so no change to expo    Dim As Ubyte Ptr ubp = Strptr(b.s) + Len(b.s) - 1 ' point to the MSbyte    For i = 0 To 4  ' find the leading non-zero byte, MS block may be zero        If *ubp > 0 Then Exit For        ubp = ubp - 1        expo = expo - 8 ' expo reduction of 8 bits per zero byte skipped    Next i  ' ubp now points to the MS non-zero byte    uli = *Cast(Ulongint Ptr, ubp - 7)  ' normalize bytes, left justify    For i = 63 To 56 Step -1    ' find the MS set bit        If Bit(uli, i) Then Exit For        expo = expo - 1    Next i  ' i now points to MSbit    uli = uli Shr (i - 52)  ' move MSbit into bit 52    uli = Bitreset(uli, 52) ' kill the MSbit implied bit in 52    uli = Sign_Bit Or (expo Shl 52) Or uli  ' build the double    Return *Cast(Double Ptr, @uli)End Function'===================================================================' test the BigInt_2_Double functionDim As BigInt xx = -8Dim As Double dDim As Integer iFor i = 0 To 17    Print x; "  ";     Print Tab(8)    Print BigInt_2_Double( x )    x = x +1Next iPrintFor 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 iPrintFor 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 iPrintx.s = String(128, Chr(255)) + Chr(0,0,0,0) ' most positive doublePrint xPrint Len(x.s),Print BigInt_2_Double(x)Printx = x + 1Print xPrint Len(x.s),Print BigInt_2_Double(x)Printx.s = Chr(1) + String(127, Chr(0)) + Chr(255,255,255,255) ' most negativePrint xPrint Len(x.s),Print BigInt_2_Double(x)Printx = x - 1Print xPrint 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: 2984
Joined: Jan 15, 2007 20:44
Location: Australia
CInt(), CLongInt() and CDbl() are now supported in this latest version called Big_Integer_2.bas

Code: Select all

`' version 1.2   3 July 2011'================================================================' Big_Integer_2.bas  Arbitrary Precision, Two's Complement Integer.'================================================================' supported functions are;' (short constants -1, 0, 1, 2 ) Bigint_m1, Bigint_0, Bigint_1, Bigint_2' Msbit, Prune, Neg, Absolute, pack, unpack, Addition, Subtract' Multiply, Square, Mul2, Div2, Power, Factorial' Shift_Left, Shift_Right, Complement, Bit_value, Bit_set, Bit_reset' Bit_And, Bit_Or, Bit_Xor, Bit_Eqv, Bit_Imp, Show_Bin, Show_Hex, Show_Oct' is_ZERO, is_NZ, is_POS, is_NEG, is_EQ, is_NEQ' is_LT, is_GT, is_LT_EQ, is_GT_EQ, BigInt_Sgn, DivMod, Divide, Remainder' CInt, CLongInt, CDbl'----------------------------------------------------------------' This package only handles integers encoded in a two's complement format.' The first bit in the two's complement format is always the sign which' takes the value of 1 = negative, 0 = positive or zero. Byte examples are;' +127 = 0111 1111  most positive'   +8 = 0000 1000'   +7 = 0000 0111'   +4 = 0000 0100'   +2 = 0000 0010'   +1 = 0000 0001'    0 = 0000 0000  zero'   -1 = 1111 1111'   -2 = 1111 1110'   -4 = 1111 1100'   -7 = 1111 1001'   -8 = 1111 1000' -127 = 1000 0001' -128 = 1000 0000  most negative'----------------------------------------------------------------' Each successive 32 bits are packed into a 4 byte block.' Each block is stored in 4 successive bytes of a string.'----------------------------------------------------------------' Strings in FreeBASIC are indexed and printed from left to right. Big' integers appear to be stored in strings backwards which can be confusing.' The Left side of the string is the Right side of the number and vice versa.' Beware: Shift_Left moves a string to the right, Shift_Right moves it left.'----------------------------------------------------------------' The first byte in the string is the least significant byte of the number.' The last block in the string is the most significant block of the number.' String s indexing has s[0] as the LS byte and s[len(s)-1] as the MS byte.' These big integer strings are always multiples of 4 bytes long.' The msb of a number is sign extended to the MS bit of the MS block.'----------------------------------------------------------------' A number is always stored in a way that correctly represents that number.' Where an operation would overflow into the MSB, a sign block is' appended to the number so as to prevent overflow or sign change.' Unnecessary leading zero or all ones blocks are removed by prune.'----------------------------------------------------------------' String pointers may change if a string is created or any length is changed.'================================================================Type bigint     ' a little endian, two's complement binary number    s As String ' packed into a string, in blocks four bytes long    ' constructor    Declare Constructor ' default constructor    Declare Constructor (Byref As Integer)    Declare Constructor (Byref As Longint)    Declare Constructor (Byref As String)    ' Let    Declare Operator Let(Byref rhs As String)    Declare Operator Let(Byval rhs As Byte)    Declare Operator Let(Byval rhs As Ubyte)    Declare Operator Let(Byval rhs As Short)    Declare Operator Let(Byval rhs As Ushort)    Declare Operator Let(Byval rhs As Integer)    Declare Operator Let(Byval rhs As Uinteger)    Declare Operator Let(Byval rhs As Longint)    Declare Operator Let(Byval rhs As Ulongint)    Declare Operator Let(Byval rhs As Double)    ' cast    Declare Operator Cast() As String    Declare Operator Cast() As Integer ' CInt    Declare Operator Cast() As LongInt ' CLongInt    Declare Operator Cast() As Double  ' Cdbl    ' implicit step versions    Declare Operator For ()    Declare Operator Step()    Declare Operator Next(Byref end_cond As bigint) As Integer    ' explicit step versions    Declare Operator For (Byref step_var As bigint)    Declare Operator Step(Byref step_var As bigint)    Declare Operator Next(Byref end_cond As bigint, Byref step_var As bigint ) As Integer    ' operate and assign    Declare Operator += (Byref rhs As bigint)    Declare Operator -= (Byref rhs As bigint)    Declare Operator *= (Byref rhs As bigint)    Declare Operator \= (Byref rhs As bigint)End Type'================================================================' only the necessary declarations are hereDeclare Function pack(Byref As String) As bigintDeclare Function unpack(Byref As bigint) As StringDeclare Function negate(Byref As bigint) As bigintDeclare Operator < (Byref As bigint, Byref As bigint) As IntegerDeclare Function is_LT(Byref As bigint, Byref As bigint) As Integer'================================================================Constructor bigint  ' default constructorthis.s = Chr(0,0,0,0)   ' zeroEnd ConstructorConstructor bigint (Byref rhs As Integer)this = rhsEnd ConstructorConstructor bigint (Byref rhs As Longint)this = rhsEnd ConstructorConstructor bigint (Byref rhs As String)this = rhsEnd Constructor'================================================================' assignment operatorOperator 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 IfEnd 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 IfEnd 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 = *ipEnd 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 = *uipIf (128 And rhs) Then this.s += Chr(0,0,0,0) ' make it positiveEnd 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 = *lipEnd 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 = *ulipIf (128 And this.s[7]) Then this.s += Chr(0,0,0,0) ' make it positiveEnd Operator'----------------------------------------------------------------Operator bigint.let (Byval d As Double)Const As Ulongint implicit_bit = 2^52          ' 4503599627370496Const As Ulongint mant_mask = implicit_bit - 1 ' 4503599627370495Dim As Ulongint u, mantDim As Integer negative, sign, expoDim As bigint a'----------------------------------------------------If d < 0 Then negative = -1 ' remember signd = 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 bitexpo = (u Shr 52) And 2047  ' the 11 bit exponentmant = (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 smallElse    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 = aEnd IfEnd Operator'================================================================' create a string for printing a bigintOperator bigint.cast() As StringDim As String tt = Unpack(this)Return tEnd Operator'================================================================' common small constantsDim Shared As bigint bigint_m1, bigint_0, bigint_1, bigint_2bigint_m1.s = Chr(255,255,255,255)' minus 1bigint_0.s = Chr(0,0,0,0) ' zerobigint_1.s = Chr(1,0,0,0) ' onebigint_2.s = Chr(2,0,0,0) ' two'================================================================' find the bit position of the first bit that differs from the sign bitFunction 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 kEnd Function'-----------------------------------------------------------------------' remove unnecessary leading blocks, prune time easily earns it's keepSub 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 bigintFunction 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 sEnd 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 positiveFunction 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 bEnd 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 bytesFunction 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 bEnd 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 stringFunction 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 dEnd 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 aEnd Function' 100 digits in 4.2 usec'===============================================================' subtractFunction subtract(Byref aa As bigint, Byref bb As bigint) As bigint    Dim As bigint cc = Negate(bb)    cc = addition(aa, cc)    Return ccEnd Function'=======================================================================' square a number, approaches twice the speed of multiply for big numbersFunction 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 cEnd 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'================================================================' multiplyFunction 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 IfEnd 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 highFunction 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 bEnd Function'=======================================================================' shift down one bit, high towards lowFunction 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 bEnd Function'=======================================================================' raise a number to a positive powerFunction 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 multiplyEnd 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 fEnd Function'=======================================================================' bit functions'=======================================================================' NOT. Invert all the bits in a bigintFunction 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 aEnd Function'=======================================================================' get the value of a specified bit in a big integerFunction 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 bitvalEnd Function'================================================================' set a specified bit in a big integerFunction 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 vEnd Function'================================================================' clear a specified bit in a big integerFunction 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 vEnd Function'================================================================' shift bigint n bits leftFunction 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 sEnd 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 ANDFunction 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 cEnd Function'================================================================' bitwise OrFunction 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 cEnd Function'================================================================' bitwise XorFunction 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 cEnd Function'================================================================' bitwise Eqv,     equivalence is the complement of XorFunction 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 cEnd Function'================================================================' bitwise Imp,     implicationFunction 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 cEnd Function'================================================================' convert a bigint to binaryFunction 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 hEnd Function'================================================================' convert a bigint to hexadecimalFunction 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 hEnd Function'================================================================' convert a bigint to octalFunction 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 cEnd Function'================================================================' Relational functions'----------------------------------------------------------------Function is_ZERO(Byref a As bigint) As Integer    If Len(Trim(a.s, Chr(0))) Then Return 0 Else Return -1End 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 0End 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 positiveEnd 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 modSub 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 IfEnd Sub' 100 digits divided by 1 takes 775. usec'===============================================================' divide a by b, return the integer resultFunction divide(Byref a As bigint, Byref b As bigint) As bigint    Dim As bigint c, d    divmod(a, b, c, d)    Return cEnd Function'----------------------------------------------------------------' divide a by b and return only the integer remainderFunction remainder(Byref a As bigint, Byref b As bigint) As bigint    Dim As bigint c, d    divmod(a, b, c, d)    Return dEnd Function'=====================================================================' Overload arithmetic operators'=====================================================================' unary plus +Operator + (Byref x As bigint) As bigintReturn xEnd Operator' unary minus -Operator - (Byref x As bigint) As bigintReturn negate(x)End Operator' additionOperator + (Byref x As bigint, Byref y As bigint) As bigintReturn addition(x, y)End Operator' subtractionOperator - (Byref x As bigint, Byref y As bigint) As bigintReturn subtract(x, y)End OperatorOperator * (Byref x As bigint, Byref y As bigint) As bigintReturn multiply(x, y)End OperatorOperator \ (Byref x As bigint, Byref y As bigint) As bigintReturn divide(x, y)End OperatorOperator / (Byref x As bigint, Byref y As bigint) As bigintReturn divide(x, y) ' should this WARN or be prohibited ?End Operator'----------------------------------------------------------------' operate and assignOperator bigint.+= (Byref rhs As bigint)this = addition(this, rhs)End OperatorOperator bigint.-= (Byref rhs As bigint)this = subtract(this, rhs)End OperatorOperator bigint.*= (Byref rhs As bigint)this = multiply(this, rhs)End OperatorOperator 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 IntegerDim As bigint cc = subtract(a, b)If Len(Trim(c.s, Chr(0))) Then Return 0 Else Return -1End Operator'----------------------------------------------------------------Operator <> (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint cc = subtract(a, b)If Len(Trim(c.s, Chr(0))) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------Operator < (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint cc = subtract(a, b)If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------Operator > (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint c = subtract(b, a)If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------Operator <= (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint cc = subtract(b, a)If (128 And c.s[Len(c.s)-1]) Then Return 0 Else Return -1End Operator'----------------------------------------------------------------Operator >= (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint c = subtract(b, a)If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------' raise a bigint to the power of a long integer, result = x^y Operator ^ (Byref x As bigint, Byref y As Longint) As bigintReturn power(x, y)End Operator'=================================================================' more intrinsic numerical functions'=================================================================Function Abs_(Byref x As Double) As Double    Return Abs(x)End Function#undef AbsFunction Abs Overload (Byref x As Double) As Double    Return abs_(x)End FunctionFunction 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 SgnFunction Sgn Overload (Byref x As Double) As Double    Return Sgn_(x)End FunctionFunction 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 FunctionOperator 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 FunctionOperator 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 doubleEnd FunctionOperator 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 1Operator bigint.for( )End OperatorOperator bigint.step( )this += bigint_1End OperatorOperator bigint.next( Byref end_cond As bigint ) As IntegerReturn this <= end_condEnd Operator'----------------------------------------------------------------' explicit step versionsOperator bigint.for( Byref step_var As bigint )End OperatorOperator bigint.step( Byref step_var As bigint )this += step_var   End OperatorOperator bigint.next( Byref end_cond As bigint, Byref step_var As bigint ) As IntegerIf step_var < bigint_0 Then    Return this >= end_condElse    Return this <= end_condEnd IfEnd Operator'==============================================================='       E n d    o f    B i g    I n t e g e r    C o d e'===============================================================`
counting_pine
Posts: 6178
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs
Wow, you've put a lot of work into this, maybe it should be in Projects now.
stephanbrunker
Posts: 62
Joined: Nov 02, 2013 14:57

### Re: Arbitrarily Big Integer Routines

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

- shift_right was only commented

----------

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

So my first try was slow and buggy.

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

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

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

### Re: Arbitrarily Big Integer Routines

I found a major bug in addition:

Try this:

Code: Select all

`#include "Big_Integer_2.bas"dim as bigint a,b,cdim as uinteger bia.s=chr(0,0,0,0,183,38,83,152)bi=3036098375b=biprint "     a:", show_hex(a)print "   + b:", "        ";show_hex(b)print "==============================="c=a+bprint "   = 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 hereDeclare Function pack(Byref As String) As bigintDeclare Function unpack(Byref As bigint) As StringDeclare Function negate(Byref As bigint) As bigintDeclare Operator < (Byref As bigint, Byref As bigint) As IntegerDeclare Function is_LT(Byref As bigint, Byref As bigint) As Integer'================================================================Constructor bigint  ' default constructorthis.s = Chr(0,0,0,0)   ' zeroEnd ConstructorConstructor bigint (Byref rhs As Integer)this = rhsEnd ConstructorConstructor bigint (Byref rhs As ulong)this = rhsEnd ConstructorConstructor bigint (Byref rhs As Longint)this = rhsEnd ConstructorConstructor bigint (Byref rhs As String)this = rhsEnd Constructor'================================================================' assignment operatorOperator 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 IfEnd 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 IfEnd 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 = *ipEnd 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 = *uipIf (128 And this.s[3]) Then this.s += Chr(0,0,0,0) ' make it positiveEnd 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 = *uipIf (128 And this.s[3]) Then this.s += Chr(0,0,0,0) ' make it positiveEnd 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 = *lipEnd 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 = *ulipIf (128 And this.s[7]) Then this.s += Chr(0,0,0,0) ' make it positiveEnd Operator'----------------------------------------------------------------Operator bigint.let (Byval d As Double)Const As Ulongint implicit_bit = 2^52          ' 4503599627370496Const As Ulongint mant_mask = implicit_bit - 1 ' 4503599627370495Dim As Ulongint u, mantDim As Integer negative, sign, expoDim As bigint a'----------------------------------------------------If d < 0 Then negative = -1 ' remember signd = 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 bitexpo = (u Shr 52) And 2047  ' the 11 bit exponentmant = (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 smallElse    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 = aEnd IfEnd Operator'================================================================' create a string for printing a bigintOperator bigint.cast() As StringDim As String tt = unpack(this)Return tEnd Operator'================================================================' common small constantsDim Shared As bigint bigint_m1, bigint_0, bigint_1, bigint_2bigint_m1.s = Chr(255,255,255,255)' minus 1bigint_0.s = Chr(0,0,0,0) ' zerobigint_1.s = Chr(1,0,0,0) ' onebigint_2.s = Chr(2,0,0,0) ' two'================================================================' find the bit position of the first bit that differs from the sign bitFunction 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 kEnd Function'-----------------------------------------------------------------------' remove unnecessary leading blocks, prune time easily earns it's keepSub 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 bigintFunction 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 sEnd 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 positiveFunction 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 bEnd 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 bytesFunction 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 bEnd 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 stringFunction 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 dEnd 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 aEnd Function' 100 digits in 4.2 usec'===============================================================' subtractFunction subtract(Byref aa As bigint, Byref bb As bigint) As bigint    Dim As bigint cc = Negate(bb)    cc = addition(aa, cc)    Return ccEnd Function'=======================================================================' square a number, approaches twice the speed of multiply for big numbersFunction 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 cEnd 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'================================================================' multiplyFunction 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 IfEnd 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 highFunction 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 bEnd Function'=======================================================================' shift down one bit, high towards lowFunction 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 bEnd Function'=======================================================================' raise a number to a positive powerFunction 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 multiplyEnd 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 fEnd Function'=======================================================================' bit functions'=======================================================================' NOT. Invert all the bits in a bigintFunction 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 aEnd Function'=======================================================================' get the value of a specified bit in a big integerFunction 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 bitvalEnd Function'================================================================' set a specified bit in a big integerFunction 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 vEnd Function'================================================================' clear a specified bit in a big integerFunction 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 vEnd Function'================================================================' shift bigint n bits leftFunction 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 sEnd Function'================================================================' shift bigint n bits right, by shifting left nbits and right nblocksFunction 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 sEnd Function'================================================================' bitwise ANDFunction 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 cEnd Function'================================================================' bitwise OrFunction 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 cEnd Function'================================================================' bitwise XorFunction 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 cEnd Function'================================================================' bitwise Eqv,     equivalence is the complement of XorFunction 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 cEnd Function'================================================================' bitwise Imp,     implicationFunction 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 cEnd Function'================================================================' convert a bigint to binaryFunction 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 hEnd Function'================================================================' convert a bigint to hexadecimalFunction 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 hEnd Function'================================================================' convert a bigint to octalFunction 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 cEnd Function'================================================================' Relational functions'----------------------------------------------------------------Function is_ZERO(Byref a As bigint) As Integer    If Len(Trim(a.s, Chr(0))) Then Return 0 Else Return -1End 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 0End 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 positiveEnd 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 modSub 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 IfEnd Sub' 100 digits divided by 1 takes 775. usec'===============================================================' divide a by b, return the integer resultFunction divide(Byref a As bigint, Byref b As bigint) As bigint    Dim As bigint c, d    divmod(a, b, c, d)    Return cEnd Function'----------------------------------------------------------------' divide a by b and return only the integer remainderFunction remainder(Byref a As bigint, Byref b As bigint) As bigint    Dim As bigint c, d    divmod(a, b, c, d)    Return dEnd Function'=====================================================================' Overload arithmetic operators'=====================================================================' unary plus +Operator + (Byref x As bigint) As bigintReturn xEnd Operator' unary minus -Operator - (Byref x As bigint) As bigintReturn negate(x)End Operator' additionOperator + (Byref x As bigint, Byref y As bigint) As bigintReturn addition(x, y)End Operator' subtractionOperator - (Byref x As bigint, Byref y As bigint) As bigintReturn subtract(x, y)End OperatorOperator * (Byref x As bigint, Byref y As bigint) As bigintReturn multiply(x, y)End OperatorOperator \ (Byref x As bigint, Byref y As bigint) As bigintReturn divide(x, y)End OperatorOperator / (Byref x As bigint, Byref y As bigint) As bigintReturn divide(x, y) ' should this WARN or be prohibited ?End Operator'----------------------------------------------------------------' operate and assignOperator bigint.+= (Byref rhs As bigint)this = addition(this, rhs)End OperatorOperator bigint.-= (Byref rhs As bigint)this = subtract(this, rhs)End OperatorOperator bigint.*= (Byref rhs As bigint)this = multiply(this, rhs)End OperatorOperator 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 IntegerDim As bigint cc = subtract(a, b)If Len(Trim(c.s, Chr(0))) Then Return 0 Else Return -1End Operator'----------------------------------------------------------------Operator <> (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint cc = subtract(a, b)If Len(Trim(c.s, Chr(0))) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------Operator < (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint cc = subtract(a, b)If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------Operator > (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint c = subtract(b, a)If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------Operator <= (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint cc = subtract(b, a)If (128 And c.s[Len(c.s)-1]) Then Return 0 Else Return -1End Operator'----------------------------------------------------------------Operator >= (Byref a As bigint, Byref b As bigint) As IntegerDim As bigint c = subtract(b, a)If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0End Operator'----------------------------------------------------------------' raise a bigint to the power of a long integer, result = x^y Operator ^ (Byref x As bigint, Byref y As Longint) As bigintReturn power(x, y)End Operator'=================================================================' more intrinsic numerical functions'=================================================================Function Abs_(Byref x As Double) As Double    Return Abs(x)End Function#undef AbsFunction Abs Overload (Byref x As Double) As Double    Return abs_(x)End FunctionFunction 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 SgnFunction Sgn Overload (Byref x As Double) As Double    Return Sgn_(x)End FunctionFunction 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 FunctionOperator 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 FunctionOperator 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 doubleEnd FunctionOperator 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 1Operator bigint.for( )End OperatorOperator bigint.step( )this += bigint_1End OperatorOperator bigint.next( Byref end_cond As bigint ) As IntegerReturn this <= end_condEnd Operator'----------------------------------------------------------------' explicit step versionsOperator bigint.for( Byref step_var As bigint )End OperatorOperator bigint.step( Byref step_var As bigint )this += step_var   End OperatorOperator bigint.next( Byref end_cond As bigint, Byref step_var As bigint ) As IntegerIf step_var < bigint_0 Then    Return this >= end_condElse    Return this <= end_condEnd IfEnd Operator'----------------------------------------------------------------' convert a bigint to unsigned hexadecimalFunction 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 hEnd Function'two's compliment hex to bigintFunction 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 cEnd Function'unsigned hex to bigintFunction 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 cEnd Function'bigint to twos compliment binaryFunction bigint_bin(byref a as bigint) as string    Dim As String s = a.s    return sEnd Function'twos compliment binary to bigintFunction bin_bigint(byref a as string) as bigint    dim as bigint b    b.s = a    return bEnd Function'bigint to unsigned binaryFunction 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 sEnd Function'unsigned binary to bigintFunction 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 bEnd Function'==============================================================='       E n d    o f    B i g    I n t e g e r    C o d e'===============================================================`