Arbitrarily Big Integer Routines

User projects written in or related to FreeBASIC.
stephanbrunker
Posts: 62
Joined: Nov 02, 2013 14:57

Re: Arbitrarily Big Integer Routines

Post by stephanbrunker »

I liked the code because it's usage is very easy. Unfortunately, it was much, much slower than other bigint implementations. So I've rewritten the division and comparation, optimized some other functions, and the speed is the same now than for example the largeint.bas. I couldn't fix the constructors for integer completely, but maybe someone else can fix that. In the meantime, just stick on the other datatypes and avoid "integer" and "uinteger".

Code: Select all

' version 2.0 25 November 2013
'=============================
' by Stephan Brunker (stephanbrunker at web punkt de)
'=============================
' The easy way to use - analog to the "hello world" the 1+1:
' #include "big_integer.bas"
' Dim as bigint a,b,c
' a = "1"
' b = "1"
' c = a + b
' print c
' a and b can be so long as you want - up to 1.7E10 bits
' all the operators: + - * \ mod ^ += -= *= \= mod= are overloaded
' also all the comparators: = < > <> >= <= 
' for additional functions look in the code, it should be
' well commented to see for what the functions are for
' you can also input and output in various formats, see section
' "Conversion Functions"
'===============================================================
'known bug: Please dont input "integer" or "uinteger". 
'Integer<32> is internally handled as long, and Integer<64> as longint. 
'But simply Integer can be both, and I couldn't get the constructors 
'to work properly - you can try the commented sections below.
'same with undefined assignments: 
'(bigint) = (not bigint) {operator} (not bigint) doesn't  work too.
'please use a cuint<32> in the syntax, then the conversion is clear.
'maybe someone else can fix that ..!
'===============================================================
' Changelog
'-----------
'-removed comments @shift_right
'-bugfix constructor and operator uinteger, added ulong
'-conversion functions hex_bigint, uhex_bigint,
'    bin_bigint, bigint_bin;  ubin_bigint, bigint_ubin
'-added mod_power, removed Divide and Remainder
'-integer and uinteger defined for 32-bit/64-bit 
'-bug exploitation in negate: 
'    culngint(not(uinteger)) = cuint(not(uinteger) fixed
'-rewritten divmod(30x faster now)
'-rewritten comparation - 20x faster
'-general speed optimisation: Variables are faster than pointers, 
'    pointers faster than string operations: 
'    Square is 10% faster now, Multiply 5% faster
'-using the internal functions is slightly faster than the overloaded
'    operators because of copying the input twice
'
' version 1.2   3 July 2011
'================================================================
' Big_Integer_2.bas  Arbitrary Precision, Two's Complement Integer.
' by Richard @ freebasic.net/forum
'================================================================
' supported functions are;
' (short constants -1, 0, 1, 2 ) Bigint_m1, Bigint_0, Bigint_1, Bigint_2
' Msbit, Prune, Neg, Absolute, pack, unpack, Addition, Subtract
' Multiply, Square, Mul2, Div2, Power, Factorial
' Shift_Left, Shift_Right, Complement, Bit_value, Bit_set, Bit_reset
' Bit_And, Bit_Or, Bit_Xor, Bit_Eqv, Bit_Imp, Show_Bin, Show_Hex, Show_Oct
' is_ZERO, is_NZ, is_POS, is_NEG, is_EQ, is_NEQ
' is_LT, is_GT, is_LT_EQ, is_GT_EQ, BigInt_Sgn, DivMod, Divide, Remainder
' CInt, CLongInt, CDbl
'----------------------------------------------------------------
' This package only handles integers encoded in a two's complement format.
' The first bit in the two's complement format is always the sign which
' takes the value of 1 = negative, 0 = positive or zero. Byte examples are;
' +127 = 0111 1111  most positive
'   +8 = 0000 1000
'   +7 = 0000 0111
'   +4 = 0000 0100
'   +2 = 0000 0010
'   +1 = 0000 0001
'    0 = 0000 0000  zero
'   -1 = 1111 1111
'   -2 = 1111 1110
'   -4 = 1111 1100
'   -7 = 1111 1001
'   -8 = 1111 1000
' -127 = 1000 0001
' -128 = 1000 0000  most negative
'----------------------------------------------------------------
' Each successive 32 bits are packed into a 4 byte block.
' Each block is stored in 4 successive bytes of a string.
'----------------------------------------------------------------
' Strings in FreeBASIC are indexed and printed from left to right. Big
' integers appear to be stored in strings backwards which can be confusing.
' The Left side of the string is the Right side of the number and vice versa.
' Beware: Shift_Left moves a string to the right, Shift_Right moves it left.
'----------------------------------------------------------------
' The first byte in the string is the least significant byte of the number.
' The last block in the string is the most significant block of the number.
' String s indexing has s[0] as the LS byte and s[len(s)-1] as the MS byte.
' These big integer strings are always multiples of 4 bytes long.
' The msb of a number is sign extended to the MS bit of the MS block.
'----------------------------------------------------------------
' A number is always stored in a way that correctly represents that number.
' Where an operation would overflow into the MSB, a sign block is
' appended to the number so as to prevent overflow or sign change.
' Unnecessary leading zero or all ones blocks are removed by prune.
'----------------------------------------------------------------
' String pointers may change if a string is created or any length is changed.

'================================================================
Type bigint     ' a little endian, two's complement binary number
    s As String ' packed into a string, in blocks four bytes long
    ' constructor
    Declare Constructor ' default constructor 
    Declare Constructor (Byref As String)
    Declare Constructor (Byref As Long)
    Declare Constructor (Byref As Ulong)
'    Declare Constructor (Byref As Integer)
'    Declare Constructor (Byref As Uinteger)
    Declare Constructor (Byref As Longint)    
    Declare Constructor (Byref As Ulongint)
    ' Let
    Declare Operator Let(Byref rhs As String)
    Declare Operator Let(Byval rhs As Byte)
    Declare Operator Let(Byval rhs As Ubyte)
    Declare Operator Let(Byval rhs As Short)
    Declare Operator Let(Byval rhs As Ushort)
    Declare Operator Let(Byval rhs As Long)
    Declare Operator Let(Byval rhs As Ulong)
    Declare Operator Let(Byval rhs As Longint)
    Declare Operator Let(Byval rhs As Ulongint)
    Declare Operator Let(Byval rhs As Double)
    ' cast
    Declare Operator Cast() As String
    Declare Operator Cast() As Integer ' CInt
    Declare Operator Cast() As LongInt ' CLongInt
    Declare Operator Cast() As Double  ' Cdbl
    ' implicit step versions
    Declare Operator For ()
    Declare Operator Step()
    Declare Operator Next(Byref end_cond As bigint) As Integer
    ' explicit step versions
    Declare Operator For (Byref step_var As bigint)
    Declare Operator Step(Byref step_var As bigint)
    Declare Operator Next(Byref end_cond As bigint, Byref step_var As bigint ) As Integer
    ' operate and assign
    Declare Operator += (Byref rhs As bigint)
    Declare Operator -= (Byref rhs As bigint)
    Declare Operator *= (Byref rhs As bigint)
    Declare Operator \= (Byref rhs As bigint)
    Declare Operator mod= (Byref rhs As bigint)
End Type

'================================================================
' only the necessary declarations are here
Declare Function pack(Byref As String) As bigint
Declare Function unpack(Byref As bigint) As String
Declare Function negate(Byref As bigint) As bigint
Declare Function bigint_compare(byref a as bigint, byref b as bigint) as integer
Declare Operator < (Byref As bigint, Byref As bigint) As Integer
Declare Operator > (Byref As bigint, Byref As bigint) As Integer
Declare Operator = (Byref As bigint, Byref As bigint) As Integer
Declare Operator <> (Byref a As bigint, Byref b As bigint) As Integer
Declare Operator <= (Byref a As bigint, Byref b As bigint) As Integer
Declare Operator >= (Byref a As bigint, Byref b As bigint) As Integer
Declare Operator + (Byref x As bigint) As bigint
Declare Operator - (Byref x As bigint) As bigint
Declare Operator + (Byref x As bigint, Byref y As bigint) As bigint
Declare Operator - (Byref x As bigint, Byref y As bigint) As bigint
Declare Operator * (Byref x As bigint, Byref y As bigint) As bigint
Declare Operator \ (Byref x As bigint, Byref y As bigint) As bigint
Declare Operator mod (Byref x As bigint, Byref y As bigint) as bigint
Declare Operator ^ (Byref x As bigint, Byref y As Longint) As bigint
Declare Function is_LT(Byref As bigint, Byref As bigint) As Integer
Declare Function show_hex(Byref s As bigint) As String
Declare Function msbit(Byref a As bigint) As Integer
Declare Function BigInt_2_Double ( Byref a As BigInt ) As Double
Declare Function Bit_Value(Byref v As bigint, Byval b As Ulongint) As Integer

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

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

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

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

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

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

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

'================================================================
' assignment operator

Operator bigint.let (Byref rhs As String)
this = pack(rhs)
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Byte)
If (128 And rhs) Then
    this.s = Chr(rhs,-1,-1,-1) 
Else 
    this.s = Chr(rhs,0,0,0)
End If
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Ubyte)
this.s = Chr(rhs,0,0,0)
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Short)
If (32768 And rhs) Then
    this.s = Chr(Lobyte(rhs), Hibyte(rhs), -1, -1 )
Else
    this.s = Chr(Lobyte(rhs), Hibyte(rhs), 0, 0 )
End If
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Ushort)
this.s = Chr(Lobyte(rhs), Hibyte(rhs), 0, 0 )
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Long)
this.s = Chr(0,0,0,0)
Dim As Integer<32> Ptr bip = Cptr(Integer<32> Ptr, Strptr(this.s))
Dim As Integer<32> Ptr  ip = Cptr(Integer<32> Ptr, @rhs)
*bip = *ip
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval rhs As ULong)
this.s = Chr(0,0,0,0)
Dim As Uinteger<32> Ptr bip = Cptr(Uinteger<32> Ptr, Strptr(this.s))
Dim As Uinteger<32> Ptr uip = Cptr(Uinteger<32> Ptr, @rhs)
*bip = *uip
If (128 And this.s[3]) Then this.s &= Chr(0,0,0,0) ' make it positive
End Operator
''----------------------------------------------------------------
'Operator bigint.let (Byval rhs As Integer)
'this.s = Chr(0,0,0,0)
'if rhs > 2147483467 then print "64-bit Integers not yet supported" : sleep : end
'Dim As Integer<32> Ptr bip = Cptr(Integer<32> Ptr, Strptr(this.s))
'Dim As Integer<32> Ptr  ip = Cptr(Integer<32> Ptr, @rhs)
'*bip = *ip
'End Operator
'
''----------------------------------------------------------------
'Operator bigint.let (Byval rhs As UInteger)
'this.s = Chr(0,0,0,0)
'if rhs > 4294967295 then print "64-bit Integers not yet supported" : sleep : end
'Dim As Uinteger<32> Ptr bip = Cptr(Uinteger<32> Ptr, Strptr(this.s))
'Dim As Uinteger<32> Ptr uip = Cptr(Uinteger<32> Ptr, @rhs)
'*bip = *uip
'If (128 And this.s[3]) Then this.s &= Chr(0,0,0,0) ' make it positive
'End Operator
'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Longint)
this.s = String(8, Chr(0) )
Dim As Longint Ptr bip = Cptr(Longint Ptr, Strptr(this.s))
Dim As Longint Ptr lip = Cptr(Longint Ptr, @rhs)
*bip = *lip
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval rhs As Ulongint)
this.s = String(8, Chr(0) )
Dim As Ulongint Ptr  bip = Cptr(Ulongint Ptr, Strptr(this.s))
Dim As Ulongint Ptr ulip = Cptr(Ulongint Ptr, @rhs)
*bip = *ulip
If (128 And this.s[7]) Then this.s &= Chr(0,0,0,0) ' make it positive
End Operator

'----------------------------------------------------------------
Operator bigint.let (Byval d As Double)
    Const As Ulongint implicit_bit = 2^52          ' 4503599627370496
    Const As Ulongint mant_mask = implicit_bit - 1 ' 4503599627370495
    Dim As Ulongint u, mant
    Dim As Integer negative, sign, expo
    Dim As bigint a
    '----------------------------------------------------
    If d < 0 Then negative = -1 ' remember sign
    d = Int(Abs(d) + 0.5) ' rectify and round to closest integer
    '----------------------------------------------------
    u = *(Cptr(Ulongint Ptr, @d))   ' copy Double into a Ulongint frame
    'sign = (u Shr 63)   ' extract the sign bit
    expo = (u Shr 52) And 2047  ' the 11 bit exponent
    mant = (u And mant_mask )   ' 52 mantissa bits
    '----------------------------------------------------
    If expo = 0 Then    ' the double has zero value or is de-normalized
        this.s = Chr(0,0,0,0)   ' de-normalized is very very small
    Else
        mant = mant + implicit_bit ' insert the missing implicit bit
        expo = expo - 1023  ' remove the excess exponent 
        If expo < 53 Then
            mant = mant Shr (52 - expo) ' reduce it to integer
            a.s = String(8, Chr(0))
            *(Cptr(Ulongint Ptr, Strptr(a.s))) = mant   ' make bigint
            If negative Then a = negate(a)
        Else
            Print "WARNING: Double argument was unreliable."
            Stop
        End If
        this = a
    End If
End Operator

'================================================================
' create a string for printing a bigint
Operator bigint.cast() As String
Dim As String t
t = unpack(this)
Return t
End Operator

Operator bigint.Cast() As Double    ' Cdbl
    Return BigInt_2_Double(this)
End Operator

'================================================================
' common small constants
'Dim Shared As bigint bigint_m1, bigint_0, bigint_1, bigint_2
'bigint_m1.s = Chr(255,255,255,255)' minus 1
'bigint_0.s = Chr(0,0,0,0) ' zero
'bigint_1.s = Chr(1,0,0,0) ' one
'bigint_2.s = Chr(2,0,0,0) ' two

Dim Shared As bigint bigint_m1=-1
Dim Shared As bigint bigint_0=0
Dim Shared As bigint bigint_1=1
Dim Shared As bigint bigint_2=2

'-----------------------------------------------------------------------
' remove unnecessary leading blocks, prune time easily earns it's keep
Sub prune(Byref a As bigint)
    a.s = Left(a.s, (((msbit(a) + 1) Shr 5) + 1 ) Shl 2)
End Sub
' Times and space are improved through the sensible use of prune.
' Addition of opposite signs or subtraction can cancel many blocks. 
' A redundant block can appear during multiplication. Square or div2. 
' Mul2, Complement, Negate and Absolute do not generate unnecessary blocks.
' Power is pruned internally by the prune in multiply and square.

'================================================================
'                  ARITHMETRIC OPERATIONS
'================================================================
' Negate the twos complement binary number in a bigint
Function negate(Byref a As bigint) As bigint
    Dim As bigint s = a
    Dim As Integer blocks = Len(s.s) Shr 2
    Dim As Ulongint sum
    Dim As Uinteger<32> carry
    Dim As Uinteger<32> Ptr ps
    ps = Cast( Uinteger<32> Ptr, Strptr(s.s))' the Uinteger data in bigint
    carry = 1       ' set carry
    Do  ' slow ahead until clear of the carry
        sum = cuint(Not *ps) + carry
        *ps = sum
        carry = sum shr 32
        ps += 1
        blocks -= 1
    Loop Until (carry = 0) Or (blocks = 0)
    If blocks > 0 Then
        Do  ' no carry, so full speed ahead
            *ps = Not *ps
            ps +=1
            blocks -= 1
        Loop Until blocks = 0
    End If
    ' Negating the most negative integer is a problem because carry propagation
    ' flips the sign which should have become positive. But negation of zero
    ' does not change the sign either, so we need to differentiate between zero
    ' and one by quickly examining the final carry bit from the two's complement.
    If carry = 0 Then ' carry was not generated by the most negative number
        If (128 And a.s[Len(a.s)-1]) = (128 And s.s[Len(s.s)-1]) Then s.s = s.s + bigint_0.s
    End If  ' this prevents a negated zero being extended by an extra zero block
    Return s
End Function

'   digits    time on 2 GHz Pentium
'     10    2.2 usec
'    100    2.5 usec
'     1k    4.2 usec
'    10k    16. usec
'   100k   200. usec

'================================================================
' absolute value is positive
Function absolute(Byref a As bigint) As bigint
    Dim As bigint b = a
    If 128 And b.s[Len(b.s)-1] Then b = negate(b)
    Return b
End Function
'================================================================
' addition
Function addition(Byref aa As bigint, Byref bb As bigint) As bigint
    Dim As bigint a = aa, b = bb
    Dim As Integer blocks, i, j, lena, lenb, sa, sb, delta
    '------------------------------------------------------------
    ' test to see if the two most significant digits differ which 
    lena = Len(a.s) - 1   ' might change the sign without carry
    If a.s[lena] And 128 Then sa = 255  ' sign as a byte
    i = a.s[lena] And 192 ' if MSBs differ then extend the sign
    If (i = 64) Or (i = 128) Then a.s = a.s + String(4, Chr(sa) )
    '------------------------------------------------------------
    lenb = Len(b.s) - 1
    If b.s[lenb] And 128 Then sb = 255
    i = b.s[lenb] And 192
    If (i = 64) Or (i = 128) Then b.s = b.s + String(4, Chr(sb) )
    '------------------------------------------------------------
    ' align the two bigints to be added
    delta = Len(a.s) - Len(b.s) 'new values
    If delta <> 0 Then  ' sign extend the shorter
        If delta > 0 Then
            ' a = a
            If b.s[Len(b.s)-1] And 128 Then i = 255 Else i = 0
            b.s = b.s + String(delta, Chr(i) )  ' extend b
        Else
            If aa.s[Len(aa.s)-1] And 128 Then i = 255 Else i = 0
            a.s = a.s + String(-delta, Chr(i) )  ' extend a
            ' b = b
        End If
    End If  ' a and b are now the same length
    '------------------------------------------------------------
    ' accumulate b into a
    blocks = Len(a.s) Shr 2
    Dim As Ulongint sum = 0 ' clear carry
    Dim As Uinteger<32> carry
    Dim As Uinteger<32> Ptr pa, pb
    pa = Cast(Uinteger<32> Ptr, Strptr(a.s) )
    pb = Cast(Uinteger<32> Ptr, Strptr(b.s) )
    For i = 0 To blocks-1
        sum = Clngint(*(pa+i)) + Clngint(*(pb+i)) + carry
        *(pa+i) = sum
        carry = sum shr 32
    Next i
    prune(a)
    Return a
End Function
' 100 digits in 4.2 usec

'===============================================================
' subtract
Function subtract(Byref aa As bigint, Byref bb As bigint) As bigint
    Dim As bigint cc = Negate(bb)
    cc = addition(aa, cc)
    Return cc
End Function

'=======================================================================
' square a number, approaches twice the speed of multiply for big numbers
Function square(Byref aa As bigint) As bigint
    '------------------------------------------------------------
    Dim As bigint a = aa, c
    If (128 And a.s[Len(a.s)-1]) Then a = negate(a)
    '------------------------------------------------------------
    ' find the dimension of the problem
    Dim As Integer i, j, asize = Len(a.s) ' number of bytes in a
    c.s = String(2 * asize, Chr(0) ) ' initialise accumulator
    asize = (asize Shr 2) - 1   ' highest block number in a
    '------------------------------------------------------------
    ' pointers into all the bigints
    Dim As Uinteger<32> Ptr pa, pc
    pa = Cast(Uinteger<32> Ptr, Strptr(a.s) )  ' the base addresses of bigints
    pc = Cast(Uinteger<32> Ptr, Strptr(c.s) )
    Dim As Ulongint product ' bottom half is data, top half will be carry
    Dim As Uinteger<32> carry, sum
    '------------------------------------------------------------
    ' multiply one triangular corner only
    For i = 1 To asize
            ' pa starts at 1 not zero because 0,0 is on the diagonal
            ' the second element in a starts at zero
            ' pc is the accumulator ic = icz + i + j
        carry = 0     ' clear carry
        For j = 0 To i - 1
            product = Culngint(*(pa+i)) * Culngint(*(pa+j)) + *(pc+j+i) + carry 
            *(pc+j+i) = product
            carry = product shr 32
        Next j
        *(pc+j+i) = carry     ' line of blocks gets one longer each loop
    Next i
    '------------------------------------------------------------
    ' double the triangle, cannot do it at same time as squaring diagonal
                    ' because it can cause the product to overflow
    carry = 0 ' clear carry
    For i = 0 To (2 * asize) + 1
        sum = *(pc+i)
        product = Culngint(sum) + sum + carry
        *(pc+i) = product
        carry = product shr 32
    Next i
    '------------------------------------------------------------
    ' square and accumulate the diagonal elements
    carry = 0     ' clear carry
    For i = 0 To asize
        ' square the diagonal entry, while propagating carry
        sum = *(pa+i)
        product = Culngint(sum) * Culngint(sum) + *(pc+i+i) + carry
        *(pc+i+i) = product
        carry = product shr 32
        ' propagate carry through the following odd block of C
        product = Culngint(*(pc+i+i+1)) + carry
        *(pc+i+i+1) = product
        carry = product shr 32
    Next i
    '------------------------------------------------------------
    If 128 And c.s[Len(c.s)-1] Then c.s = c.s + bigint_0.s   ' sign is positive
    prune(c)
    Return c
End Function

'   digits   multiply(x,x)  square(x)
'      10      3.9 usec     3.9 usec
'     100     11.5 usec     9.0 usec
'      1k      640 usec     340 usec
'     10k       64 msec      33 msec
'    100k     6.35 sec     3.21 sec
'      1M                   320 sec

'================================================================
' multiply
Function multiply(Byref aa As bigint, Byref bb As bigint) As bigint
'    If @aa = @bb Then   ' detect if square, half time if very big
'        Return square(aa)
'    Else
        '------------------------------------------------------------
        ' sort out the signs and rectify the inputs
        Dim As bigint a = aa, b = bb, c
        Dim As Integer sign_a, sign_b, sign_c
        sign_a = a.s[Len(a.s)-1] And 128
        sign_b = b.s[Len(b.s)-1] And 128
        sign_c = sign_a Xor sign_b 
        If sign_a Then a = negate(a)
        If sign_b Then b = negate(b)
        '------------------------------------------------------------
        ' find the dimensions of the problem
        Dim As Integer i, j, asize, bsize
        asize = Len(a.s) ' number of bytes in a
        bsize = Len(b.s) ' number of bytes in b
        c.s = String(asize + bsize, Chr(0)) ' initialise accumulator
        asize = asize Shr 2 - 1 ' number of blocks in a 
        bsize = bsize Shr 2 - 1 ' Shr 2 is faster than \4
        '------------------------------------------------------------
        ' pointers into all the bigints
        Dim As Uinteger<32> Ptr ia, ib, ic
        ia = Cast(Uinteger<32> Ptr, Strptr(a.s) )
        ib = Cast(Uinteger<32> Ptr, Strptr(b.s) )
        ic = Cast(Uinteger<32> Ptr, Strptr(c.s) )
        Dim As Ulongint product
        Dim As Uinteger<32> carry
        '------------------------------------------------------------
        For i = 0 To asize
            carry = 0 ' clear carry
            For j = 0 To bsize
                product = Culngint(*(ia+i)) * Culngint(*(ib+j)) + *(ic+i+j) + carry
                *(ic+i+j) = product
                carry = product shr 32
            Next j
            *(ic+i+j) = carry
        Next i
        '------------------------------------------------------------
        If sign_c = 128 Then c = negate(c)
        prune(c)
        Return c
'    End If
End Function    

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

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

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

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

'=======================================================================
' integer divide, a \ b, return div and mod
Sub divmod(_
    Byref aa As bigint, Byref bb As bigint,_
    Byref q As bigint, Byref r As bigint)
    Dim As Integer lenaa = Len(aa.s), lenbb = Len(bb.s)
    '------------------------------------------------------------       
    'Test if Division by zero
    Dim as uinteger<32> ptr pbb
    pbb = Cast(Uinteger<32> Ptr, Strptr(bb.s))
    If *pbb = 0 and lenbb = 4 then
        Print " Division by zero. "
        sleep: End
    End If 
    '------------------------------------------------------------       
    'Test if Longint Division possible
    If (lenaa <= 8) And (lenbb <= 8) Then ' arguments are one or two blocks
        q.s = Chr(0,0,0,0,0,0,0,0)   ' use 64 bit long integers
        r.s = Chr(0,0,0,0,0,0,0,0)   ' allocate space for results
        If lenaa = 4 Then aa.s = aa.s + String(4, Bit(aa.s[lenaa - 1], 7))
        If lenbb = 4 Then bb.s = bb.s + String(4, Bit(bb.s[lenbb - 1], 7))
        Dim As Longint Ptr pa1, pb1, pq1, pr1
        pa1 = Cast(Longint Ptr, Strptr(aa.s))    '2 block bigint -> Longint
        pb1 = Cast(Longint Ptr, Strptr(bb.s))
        pq1 = Cast(Longint Ptr, Strptr(q.s))
        pr1 = Cast(Longint Ptr, Strptr(r.s))
        If *pb1 = 0 then print " Divison by zero. " : sleep : end
        *pq1 = *pa1 \ *pb1
        *pr1 = *pa1 Mod *pb1
        Exit Sub
    End If
    '------------------------------------------------------------       
    'Test if divisor is bigger than dividend
    If bb > aa then
        q = bigint_0
        r = aa
        exit sub
    End if
    '------------------------------------------------------------       
    'Read Signs
    Dim As Integer sa, sb, sq
    sa = 128 And aa.s[lenaa-1]
    sb = 128 And bb.s[lenbb-1]
    sq = sa Xor sb  ' sign of the result
    '---------------------------------------------------------------------
    ' Setup variables and pointers
    ' r=dividend and remainder
    ' b=divisor
    ' q=quotient
    ' sum=interim result for subtraction
    Dim As bigint b = bb
    If sb Then b = negate(b)    
    r = aa
    If sa Then r = negate(r)  
    Dim As Ushort Ptr pr, pb, pq
    Dim As Uinteger<32> sum, high48b, offset=&hFFFF0000
    Dim as Uinteger rblocks,bblocks,blockshift,rounds,blocks,i,j
    Dim as Uinteger lenr=len(r.s), lenb=len(b.s)
    Dim as Ushort qi,carry
    Dim as Integer substract
    Dim As Ulongint high48r=0
    pr = Cast(Ushort Ptr, Strptr(r.s))
    pb = Cast(Ushort Ptr, Strptr(b.s)) 
    rblocks=(lenr \ 2)-1
    bblocks=(lenb \ 2)-1
    '------------------------------------------------------------       
    'convert to 16bit blocklength for effective testdivison
    'append exactly two zero blocks and ignore the second leading zeroblock
    if *(pr+rblocks) <> 0 then
        r.s = r.s & chr(0,0,0,0)
        rblocks += 1
    elseif *(pr+rblocks)=0 and *(pr+rblocks-1)=0 then
        rblocks -= 1
    end if
    if *(pb+bblocks) <> 0 then
        b.s = b.s & chr(0,0,0,0)
        bblocks += 1
    elseif *(pb+bblocks)=0 and *(pb+bblocks-1)=0 then
        bblocks -= 1
    end if      
    'new pointer after changing string
    pr = Cast(Ushort Ptr, Strptr(r.s))
    pb = Cast(Ushort Ptr, Strptr(b.s)) 
    blockshift=rblocks-bblocks
    rounds=blockshift + 1
    '------------------------------------------------------------       
    ' setup quotient
    if bit(rounds,0)=0 then
        q.s = string(rounds*2,0) 'quotient 
         pq = Cast(Ushort Ptr, Strptr(q.s))
         pq += (len(q.s) \ 2) - 1 'start at msb
    else
        q.s = string((rounds+1)*2,0)
        pq = Cast(Ushort Ptr, Strptr(q.s)) 
        pq += (len(q.s) \ 2) - 2 'start at msb 
    end if
    '------------------------------------------------------------ 
    ' start calculation
    'most significant 48bits of divisor are constant
    if bblocks = 1 then
        high48b = Cuint(*(pb+bblocks-1) shl 16) 'second block is empty
    else
        high48b = Cuint(*(pb+bblocks-1) shl 16) + *(pb+bblocks-2) + 1 
    end if
    for blocks = 1 to rounds
        'testdivision: calculate highest three blocks and divide
        high48r = Culngint(*(pr+rblocks)) shl 32 + _
            Cuint(*(pr+rblocks-1) shl 16) + *(pr+rblocks-2)
        qi=(high48r \ high48b )            'estimate quotient, max error 1

        if qi > 0 then
            'r= r - q*b for every block, begin at lsblock, with carry
            carry=0
            for i=blockshift to rblocks      
                if carry=0 then     'don't apply the offset!
                    sum = *(pr+i) - (*(pb+i-blockshift) * qi) + carry
                else          'offset because the overflow 16->32 bit
                    sum = *(pr+i) - (*(pb+i-blockshift) * qi) + carry + offset
                end if
                *(pr+i)=sum 
                carry=sum shr 16
            next i
        end if
        '------------------------------------------------------------ 
        'if c>b, additional substraction needed
        substract=0
        for i=rblocks to blockshift step -1
            substract= *(pr+i) - *(pb+i-blockshift)
            if substract <> 0 then exit for
        next i
        if substract > 0 then
            carry=1
            qi += 1     
            for i=blockshift to rblocks
                'substract from remainder from first step
                sum = *(pr+i) + Cushort(Not(*(pb+i-blockshift))) + carry
                *(pr+i)=sum
                carry=sum shr 16
            next i              
        End if 
        '------------------------------------------------------------ 
        rblocks -= 1
        blockshift -= 1
        *pq=qi      'write quotient with pointer      
        pq -= 1     'next block of quotient
    next blocks
    '------------------------------------------------------------ 
    'finalisation
    prune(r)    ' trim result
    prune(q)
    If sq Then q = negate(q)  ' Xor of input signs
    If sa Then r = negate(r)    ' sign of original input A
    '------------------------------------------------------------
End Sub

'=======================================================================
Function Factorial Overload(Byref a As bigint) As bigint
    Dim As bigint f=bigint_1, n=a
    Do Until is_LT(n, bigint_2)
        f = multiply(f, n)
        n = subtract(n, bigint_1)
    Loop 
    Return f
End Function

'=======================================================================
'exponentiation modulus
Function mod_power(byref bb as bigint, byref ee as bigint, byref m as bigint) as bigint
    if (ee.s[Len(ee.s)-1] And 128) then print "Cannot raise a bigint to a negative power"
    if ee = bigint_0 then return 1
    dim as bigint r= bigint_1, b=bb,x 'x is dump for the quotient
    dim as integer bitlen,i,z
    Dim as Uinteger<32> ptr pee
    pee=cast(Uinteger<32> ptr, strptr(ee.s))
    Dim as Uinteger<32> spee
    spee = *pee         'load first block of exponent in variable
    
    bitlen=msbit(ee)    'the highest set bit
    divmod(b,m,x,b)     'initial reduction

    'i counts from the lsb to msb, z counts the 32 bits in a block
    for i=0 to bitlen-1
        if Bit(spee,z) then     'if bit is set then multiply
            r *= b
            divmod(r,m,x,r)
        end if
        b=square(b)
        divmod(b,m,x,b)
        if z=31 then            'reset z, next block
            z=0 
            pee += 1
            spee = *pee
        else 
            z += 1
        end if
    next i
    r *= b  'bitvalue for highest bit=1
    divmod(r,m,x,r)
    return r
End Function

'=======================================================================
'               BIT FUNCTIONS
'=======================================================================
' find the bit position of the first bit that differs from the sign bit
Function msbit(Byref a As bigint) As Integer
    Dim As Integer i, j, k = 0
    i = Len(a.s)
    If 128 And a.s[Len(a.s)-1] Then ' negative
        Do  ' find the highest non-255 byte in the string
            i = i - 1
            j = a.s[i]
        Loop Until (j < 255) Or (i = 0)
        j = 255 - j
    Else                ' positive
        Do  ' find the highest non-zero byte in the string
            i = i - 1
            j = a.s[i]
        Loop Until (j > 0) Or (i = 0)
    End If
    ' find the highest non-sign bit in the byte
    If j And   1 Then k = 1 ' straight code is faster than a loop
    If j And   2 Then k = 2
    If j And   4 Then k = 3
    If j And   8 Then k = 4
    If j And  16 Then k = 5
    If j And  32 Then k = 6
    If j And  64 Then k = 7
    If j And 128 Then k = 8
    k = k + (i * 8) - 1 ' if no bits differ (-1 or 0) then return -1
    Return k
End Function

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

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

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

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

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

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

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

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

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

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

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

'================================================================
' Relational functions - Don't know why these are outcommented
'----------------------------------------------------------------
Function is_ZERO(Byref a As bigint) As Integer
    If Len(Trim(a.s, Chr(0))) Then Return 0 Else Return -1
    'if a=bigint_0 would be faster 942:1764 cycles
End Function

'----------------------------------------------------------------
'Function is_NZ(Byref a As bigint) As Integer
'    If Len(Trim(a.s, Chr(0))) Then Return -1 Else Return 0
'End Function
'
''----------------------------------------------------------------
'Function is_POS(Byref a As bigint) As Integer
'    ' note that zero is deemed to be positive
'    If (128 And a.s[Len(a.s)-1]) Then Return 0 Else Return -1
'End Function
'
''----------------------------------------------------------------
'Function is_NEG(Byref a As bigint) As Integer
'    If (128 And a.s[Len(a.s)-1]) Then Return -1 Else Return 0
'End Function
'
''----------------------------------------------------------------
'Function is_EQ(Byref a As bigint, Byref b As bigint) As Integer
'    Dim As bigint c
'    c = subtract(a, b)
'    If Len(Trim(c.s, Chr(0))) Then Return 0 Else Return -1
'End Function
'
''----------------------------------------------------------------
'Function is_NEQ(Byref a As bigint, Byref b As bigint) As Integer
'    Dim As bigint c
'    c = subtract(a, b)
'    If Len(Trim(c.s, Chr(0))) Then Return -1 Else Return 0
'End Function
'
'----------------------------------------------------------------
Function is_LT(Byref a As bigint, Byref b As bigint) As Integer
    Dim As bigint c = subtract(a,b)
    If (128 And c.s[Len(c.s)-1]) Then Return -1 Else Return 0
End Function

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

''----------------------------------------------------------------
'Function is_GT_EQ(Byref a As bigint, Byref b As bigint) As Integer
'    Dim As bigint c
'    c = subtract(a, b)
'    If (128 And c.s[Len(c.s)-1]) Then Return 0 Else Return -1
'End Function
'
'----------------------------------------------------------------
Function BigInt_Sgn(Byref a As bigint) As Integer
    Dim As Integer i = 128 And a.s[Len(a.s)-1]
    If i Then Return -1 ' is negative
    If a=bigint_0 Then Return 0 ' is zero
    Return 1 ' is positive
End Function

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

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

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

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

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

Operator \ (Byref x As bigint, Byref y As bigint) As bigint
dim as bigint a,b
divmod(x,y,a,b)
Return a
End Operator

Operator / (Byref x As bigint, Byref y As bigint) As bigint
print "no floating point division in Big_Integer"
sleep: end
return 0
End Operator

Operator mod (Byref x As bigint, Byref y As bigint) as bigint
dim as bigint c,d
divmod(x,y,c,d)
return d
End Operator

'----------------------------------------------------------------
' operate and assign
Operator bigint.+= (Byref rhs As bigint)
this = addition(this, rhs)
End Operator

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

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

Operator bigint.\= (Byref rhs As bigint)
dim as bigint d
divmod(this,rhs,this,d)
End Operator

Operator bigint.mod= (Byref rhs as bigint)
dim as bigint c
divmod(this,rhs,c,this)
End Operator

'===================================================================
' relational operators, return a boolean as an integer, 0 is false. 
'===================================================================
'Main Comparation function - faster than substract ...
Function bigint_compare(byref a as bigint, byref b as bigint) as integer
    'return -1 for a>b; 1 for a<b and 0 for equal
    dim as integer i,z
    dim as ubyte vala, valb
    dim as byte signa, signb
    signa=128 And a.s[Len(a.s)-1]   '-1 (true) negative, 0 (false) positive
    signb=128 and b.s[len(b.s)-1]
    '-------------------------------------------
    'sign is different - easy:
    If signa = 0 and signb = -128 Then        
            return -1
            exit function
    elseif signa = -128 and signb = 0 then
            return 1
            exit function
    end if
    '-------------------------------------------
    'len is different - easy:
    if len(a.s) > len(b.s) then
        if signa = 0 then 
            return -1
        else
            return 1
        end if
        exit function
    elseif len(a.s) < len(b.s) then
        if signa = 0 then 
            return 1
        else
            return -1
        end if
        exit function
    end if
    '-------------------------------------------
    'compare block for block:
    for i=len(a.s)-1 to 0 step -1
        z=a.s[i] - b.s[i]
        if z > 0 then
            return -1
            exit function
        elseif z < 0 then
            return 1
            exit function
        end if
    next i
    return 0
end function 

'----------------------------------------------------------------------

Operator = (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a, b)
if c=0 then return -1 else return 0
End Operator

'----------------------------------------------------------------
Operator <> (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a, b)
if c=0 then return 0 else return -1
End Operator

'----------------------------------------------------------------
Operator < (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a,b)
if c = 1 then return -1 else return 0
End Operator

'----------------------------------------------------------------
Operator > (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a,b)
if c = -1 then return -1 else return 0
End Operator

'----------------------------------------------------------------
Operator <= (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a,b)
if c = 1 or c = 0 then return -1 else return 0
End Operator

'----------------------------------------------------------------
Operator >= (Byref a As bigint, Byref b As bigint) As Integer
Dim As Integer c = bigint_compare(a, b)
if c = -1 or c = 0 then return -1 else return 0
End Operator

'----------------------------------------------------------------
' raise a bigint to the power of a long integer, result = x^y 
Operator ^ (Byref x As bigint, Byref y As Longint) As bigint
Return power(x, y)
End Operator

'================================================================
' FOR bigint, NEXT bigint and STEP bigint 
'=====================================================================
' implicit step versions. Implicit step is 1
Operator bigint.for( )
End Operator

Operator bigint.step( )
this += bigint_1
End Operator

Operator bigint.next( Byref end_cond As bigint ) As Integer
Return this <= end_cond
End Operator

'----------------------------------------------------------------
' explicit step versions
Operator bigint.for( Byref step_var As bigint )
End Operator

Operator bigint.step( Byref step_var As bigint )
this += step_var   
End Operator

Operator bigint.next( Byref end_cond As bigint, Byref step_var As bigint ) As Integer
If step_var < bigint_0 Then
    Return this >= end_cond
Else
    Return this <= end_cond
End If
End Operator

'=================================================================
' more intrinsic numerical functions
'=================================================================
Function Abs_(Byref x As Double) As Double
    Return Abs(x)
End Function
#undef Abs
Function Abs Overload (Byref x As Double) As Double
    Return abs_(x)
End Function

Function Abs(Byref x As bigint) As bigint
    Return Absolute(x)
End Function

'----------------------------------------------------------------
Function Sgn_(Byref x As Double) As Double
    Return Sgn(x)
End Function
#undef Sgn
Function Sgn Overload (Byref x As Double) As Double
    Return Sgn_(x)
End Function

Function Sgn(Byref x As bigint) As Integer
    Return BigInt_Sgn(x)
End Function

stephanbrunker
Posts: 62
Joined: Nov 02, 2013 14:57

Re: Arbitrarily Big Integer Routines

Post by stephanbrunker »

Wow, it's too big now ... so here the last part:

Code: Select all


'================================================================
'                  CONVERSION FUNCTIONS
'================================================================
' pack ascii numeric string to a straight binary number in a bigint
' the output bigint will have a length that is a multiple of 4 bytes
Function pack( Byref ascii As String) As bigint
    Dim As String a  ' a is the ascii input
    Dim As bigint b  ' b is the binary output
    Dim As Integer p, i, j, ch, blocks, sign
    a = ascii
    If Instr(a, "-") <> 0 Then sign = -1
    '------------------------------------------------------------
    ' extract numerics from input string
    j = 0   '  in-place write pointer
    For i = 0 To Len(a) - 1
        ch = a[i]
        If (ch >= Asc("0")) And (ch <= Asc("9")) Then
            a[j] = ch
            j += 1
        End If
    Next i
    a = Left(a, j)  ' j is one ahead from string index = length of a
    '------------------------------------------------------------
    ' extend to next multiple of 9 digits
    i = Len(a)
    blocks = Int(0.99 + i / 9)  ' number of 9 digit blocks needed
    p = 9 * blocks
    a = String(p - i, "0") + a  ' pad to next multiple of 9 digits
    '------------------------------------------------------------
    ' decimal to binary conversion
    i = ( 8 + Len(a) * 3.32192809488) \ 8   ' bytes needed for binary
    blocks = 1 + (i \ 4)                    ' adjust to multiple of 4
    b.s = String(blocks * 4, Chr(0) ) ' binary destination string
    '------------------------------------------------------------
    Dim As Uinteger<32> Ptr bp, bpz, bpcarry, bpdata
    bpz = Cast(Uinteger<32> Ptr, Strptr(b.s)) ' binary output string[0]
    Dim As Ulongint product, carry, multiplier = 1e9
    bpdata = Cast(Uinteger<32> 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<32>(*bpdata)
            carry = Culngint(*bpcarry)
            bp += 1
        Next j
        ' advancing blocks only as needed doubles the speed of conversion
        If Carry Then
            *bp = carry
            blocks += 1 ' an exact count of the blocks used
        End If
    Next i
    b.s = Left(b.s, blocks * 4) ' keep only used blocks
    '-------------------------------------------------------------
    If b.s[Len(b.s)-1] And 128 Then b.s = b.s + bigint_0.s ' MSB must be 0
    If sign Then b = negate(b)
    '-------------------------------------------------------------
    Return b
End Function

    ' on 2 GHz Pentium      original time
    ' digits    new time    without blocks 
    '    100     45. usec      46 us
    '   1000    520. usec     800 us
    '    10k     29. msec      58 ms
    '   100k     2.7 sec      5.7 seconds
    '     1M     4.5 min     9.42 minutes

'===============================================================
' unpack a straight binary string to a decimal ascii string
Function unpack(Byref bb As bigint) As String
    Dim As bigint b
    Dim As String d
    b = bb
    d = Chr(0)   ' initial decimal output string
    Dim As Integer i, j, product, carry, sign
    ' if negative then negate, append the sign later
    If b.s[Len(b.s)-1] And 128 Then   ' negative
        sign = -1
        b = negate(b)
    End If
    ' change from base 2 to base 10
    For j = Len(b.s)-1 To 0 Step -1 ' each byte in string is base 256
        carry = b.s[j]   ' the next byte to add after multiply
        For i = Len(d)-1 To 0 Step -1
            product = 256 * d[i] + carry
            d[i] = product Mod 10
            carry = product \ 10
        Next i
        Do While carry > 0  ' output string overflows
            d = Chr(carry Mod 10) + d   ' extend output string
            carry = carry \ 10            '  as needed
        Loop
    Next j
    ' change from Ubyte to ASCII
    For i = 0 To Len(d) - 1
        d[i] = d[i] + Asc("0")
    Next i
    If sign Then d = "-" + d Else d = "+" + d
    Return d
End Function

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

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

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

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

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

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

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

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

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

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

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

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

'===================================================================
'twos compliment binary to bigint
Function bin_bigint(byref a as string) as bigint
    dim as bigint b
    if (128 And b.s[Len(b.s)-1]) then           'negative, pad to blocklen with FF
        b.s = a & string((4-len(a) mod 4),chr(255))
    else
        b.s = a & string((4-len(a) mod 4),chr(0))     'positive, pad to blocklen with 00
    end if
    return b
End Function

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

'===================================================================
'unsigned binary to bigint
Function ubin_bigint(byref a as string) as bigint
    Dim as bigint b
    Dim as integer pad
    pad = 4-(len(a) mod 4)
    if pad = 4 then pad = 0
    b.s=a & string(pad,0)    'Pad to blocklen
    If (128 And b.s[Len(b.s)-1]) Then b.s &= Chr(0,0,0,0) ' make it positive
    return b
End Function

'===============================================================
'       E n d    o f    B i g    I n t e g e r    C o d e
'===============================================================
stephanbrunker
Posts: 62
Joined: Nov 02, 2013 14:57

Re: Arbitrarily Big Integer Routines

Post by stephanbrunker »

For better Integration, I renamed all the functions, fixed the bug in the constructors and so on. The new Version is displayed here: http://www.freebasic-portal.de/code-bei ... g-284.html
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

Re: Arbitrarily Big Integer Routines

Post by frisian »

stephanbrunker

First some rather small things not very important.
Dim As Integer i, sign, sum, carry = 0, dim sets a variable automatic to 0 so carry = 0 is not necessary, you can drop the = 0 part.

Constructor bigint (ByRef a As String)
Dim As Integer p, i, j, ch, blocks, sign
If InStr(a, "-") <> 0 Then sign = -1
'------------------------------------------------------------
' extract numerics from input string
j = 0 ' in-place write pointer

When the program reaches the statement j = 0 “j” hasn't been used so it has the value of “0”, you can drop this line completely.

In Operator bigint.cast() As String you have this statement:
product = 256 * d + carry FreeBasic has only very limited knowledge of optimization, the constant part needs to be after the * to get optimized.product = d * 256 + carry get optimized.

Possible error.
'twos compliment binary to bigint
Function bin_bigint(ByRef a As String) As bigint
Dim As bigint b
If (128 And b.s[Len(b.s)-1]) Then 'negative, pad to blocklen with FF

I think b.s in the previous line should be a.s since b is just created and it can not be negative, a on the other hand can be negative.

I took a look at big_integer.bas just out of curiosity, normally I use GMP for big numbers. I converted a small program I made to use big_integer.bas. It did take a reasonable time for raising a number to a power and convert a big string into a bigint but it took a very long time to convert a bigint to a string. Ok it was a big number (500000 digits) and I had compiled the program using -gen GAS.

Looking at Operator bigint.cast() As String I realized that it had potential for some speed up's. I replaced the variable product by carry that made some optimization possible. I also made carry a Uinteger. (That makes “mod” and “\” use DIV instead of IDIV).I also replaced d by a pointer k.(Did only give a very small gain).
Tested it with a big number (string(300000.”2”)). It went from 523 seconds to 420 seconds (FB). For GCC it went from 57.3 second to 49 seconds. I could not get it faster. Ok replacing the “mod” and “\” with assembler code would make it about 50 % faster for FB, but assembler code can not be used with GCC. I went on doing some other things. I then realized that a byte can hold 2 digits, that would mean that the string d would be half as long, meaning that the total number of divisions needed would be halved. I then needed a extra string to hold the answer, splitting the byte into 2 digits and adjust.

Replacement code for Operator bigint.cast() As String

Code: Select all

'----------------------------------------------------------------
' unpack a straight binary string to a decimal ascii string
Operator bigint.cast() As String
  Dim As bigint b
  Dim As String d
  b = This
  d = Chr(0)                                     ' initial decimal output string
  Dim As Integer i, j, sign
  dim as uinteger carry
  ' if negative then negate, append the sign later
  If b.s[Len(b.s) - 1] And 128 Then              ' negative
    sign = - 1
    b = bigint_negate(b)
  End If
  ' change from base 2 to base 100
  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
    dim as ubyte ptr k = strptr(d) + len(d) - 1  ' added pointer
    For i = Len(d) - 1 To 0 Step - 1
      ' product = 256 * d[i] + carry
      carry += *k * 256
      *k = carry Mod 100
      carry \= 100
      k -= 1
    Next i
    Do While carry > 0                           ' output string overflows
      d = Chr(carry Mod 100) + d                 ' extend output string
      carry = carry \ 100                        '  as needed
    Loop
  Next j
  ' need to split the value of the byte into 2 digits and to convert into ASCII
  dim as string d2 = chr(d[0] \ 10 + asc( "0")) + chr(d[0] mod 10 + asc( "0"))
  ' remove a possible 0 in front of the number
  d2 = ltrim(d2, "0")
  For i = 1 To Len(d) - 1
    d2 += chr(d[i] \ 10 + asc( "0")) + chr(d[i] mod 10 + asc( "0"))
  Next i
  If sign Then d2 = "-" + d2 Else d2 = "+" + d2
  Return d2
End Operator
The timing is 230 seconds for FB and for GCC 30 seconds.
FB is now 2 times faster and GCC nearly 2 times.

Regards.

PS
GCC means compiled with -gen gcc -O max
Makoto WATANABE
Posts: 231
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

Re: Arbitrarily Big Integer Routines

Post by Makoto WATANABE »

Mysterious result of Mod operator
I tried the following program with "version 2.4 17 October 2015".
It is seemed the same value , but one is always bigger when compared.
Please teach me how to fix my program.

Code: Select all

#Include "big_integer.bas"

Dim As Bigint a=12,b=0,c=0,d=0
Dim As LongInt i,e,f

      Print "i", "c", "d", "comparison", "e", "f"
      
For i=1 To 10
   b=i
   c=a Mod b

   If c = d Then
      Print "Equal"
   ElseIf c>d Then
      e=c
      f=d
      Print i, c, d, "c greater", e, f
   ElseIf c<d Then
      Print "d greater"
   End If
Next
Sleep
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: Arbitrarily Big Integer Routines

Post by srvaldez »

hello Makoto WATANABE
the version you mention does not work for me either, use the version posted about 3 posts above http://www.freebasic.net/forum/viewtopi ... 57#p192857
that version works for me.
Makoto WATANABE
Posts: 231
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

Re: Arbitrarily Big Integer Routines

Post by Makoto WATANABE »

Dear srvaldez

Thanks for your quick reply.
I was able to correctly compare the result of the remainder, using "version 1.3 4 November 2013".
Thank you for teaching me how to fix my program.
******************

Dear Richard
Dear stephanbrunker

I would like to translate your "Big_Integer.bas" into Japanese and introduce it to Japanese people on my website.
Please consent to this.
*****************

P.S.
Using GNU with FreeBASIC is difficult for beginners.
http://www.freebasic.net/wiki/wikka.php?wakka=ExtLibgmp

I would appreciate it if you could increase available functions such as square root with "Big Integer Routines".
Makoto WATANABE
Posts: 231
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

Re: Arbitrarily Big Integer Routines

Post by Makoto WATANABE »

Dear Richard
Dear stephanbrunker

I found out that I could easily treat Big Integer by user definition functions of GMP.
http://www.freebasic.net/forum/viewtopi ... 50#p226650

I will not treat this Big Integer contents in my website.
Sorry, I beg you to understand my situation.
stephanbrunker
Posts: 62
Joined: Nov 02, 2013 14:57

Re: Arbitrarily Big Integer Routines

Post by stephanbrunker »

@frisian:
Unfortunately, I didn't got an notification four your post. In the meantime, Hartmut Ruske fixed some bugs and I am going to apply your recommendations in the next commit at www.freebasic-portal.de

So, that library is not dead, even if it is kind of exotic … but still I think there are advantages because of the overloaded operators which makes it simple to use and because it is not an compiled library, you have the source at hand and it is only 66kB of additional code.
ur_naz
Posts: 49
Joined: Mar 02, 2016 12:44

Re: Arbitrarily Big Integer Routines

Post by ur_naz »

Why don't you use private fields and private methods? Why don't you use public methods instead standalone functions?
stephanbrunker
Posts: 62
Joined: Nov 02, 2013 14:57

Re: Arbitrarily Big Integer Routines

Post by stephanbrunker »

Simply said - and its also written as comment in the code: Because these operations tend to be extensive computationally, you have the possibility to access them directly. The "private" functions are compare, square, mul2, div2, div and prune and if I understand my comments correctly, they also cannot be made private because the operators access them. Also if you want the modulus and the result of the division at the same time, you call div directly.

The only field is the data string, and this too you want to access directly in some cases.

Be careful to look at the actual version: https://www.freebasic-portal.de/code-be ... g-284.html

There all the functions are (static) member functions, not standalone, means bigint.compare, bigint.div and so on. This extension has evolved a lot from the first version at the beginning of this thread … I think I am moving / copying the code to Github, it is better to manage.
Makoto WATANABE
Posts: 231
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

Re: Arbitrarily Big Integer Routines

Post by Makoto WATANABE »

Dear stephanbrunker

Thank you for publishing Big_Integer version 2.8.
I confirmed that my usage applications running in version 1.3 works with 2.8.

I have a question.
I want to use programs by entering numerical values of Bigint from a console screen.
I think that it may be input by strings and cast, please teach me a concrete program example.
stephanbrunker
Posts: 62
Joined: Nov 02, 2013 14:57

Re: Arbitrarily Big Integer Routines

Post by stephanbrunker »

That is very, very easy, given FB's command():

bigintcmd.exe

Code: Select all

#include "big_integer.bi"

dim a as bigint = Command(1)
dim b as bigint = Command(2)

dim c as bigint = a * b
print c
end
then you can call your executable:
bigintcmd.exe "5555555555555555" "88888888888888888888888888"

Prints:
+493827160493827111111111106172839506172840
Makoto WATANABE
Posts: 231
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

Re: Arbitrarily Big Integer Routines

Post by Makoto WATANABE »

Dear stephanbrunker

Thanks for your quick reply.
Please excuse my inadequate expressions.
I'd like to enter numbers using "Input" command instead of the command line.
Please teach me a sample program.
ur_naz
Posts: 49
Joined: Mar 02, 2016 12:44

Re: Arbitrarily Big Integer Routines

Post by ur_naz »

Not bad but very slow against SBCL or even python 3.7 built-in bigint types.
tested on factorial(50 000)
Post Reply