## 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

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 hereDeclare Function pack(Byref As String) As bigintDeclare Function unpack(Byref As bigint) As StringDeclare Function negate(Byref As bigint) As bigintDeclare Function bigint_compare(byref a as bigint, byref b as bigint) as integerDeclare Operator < (Byref As bigint, Byref As bigint) As IntegerDeclare Operator > (Byref As bigint, Byref As bigint) As IntegerDeclare Operator = (Byref As bigint, Byref As bigint) As IntegerDeclare Operator <> (Byref a As bigint, Byref b As bigint) As IntegerDeclare Operator <= (Byref a As bigint, Byref b As bigint) As IntegerDeclare Operator >= (Byref a As bigint, Byref b As bigint) As IntegerDeclare Operator + (Byref x As bigint) As bigintDeclare Operator - (Byref x As bigint) As bigintDeclare Operator + (Byref x As bigint, Byref y As bigint) As bigintDeclare Operator - (Byref x As bigint, Byref y As bigint) As bigintDeclare Operator * (Byref x As bigint, Byref y As bigint) As bigintDeclare Operator \ (Byref x As bigint, Byref y As bigint) As bigintDeclare Operator mod (Byref x As bigint, Byref y As bigint) as bigintDeclare Operator ^ (Byref x As bigint, Byref y As Longint) As bigintDeclare Function is_LT(Byref As bigint, Byref As bigint) As IntegerDeclare Function show_hex(Byref s As bigint) As StringDeclare Function msbit(Byref a As bigint) As IntegerDeclare Function BigInt_2_Double ( Byref a As BigInt ) As DoubleDeclare Function Bit_Value(Byref v As bigint, Byval b As Ulongint) As Integer'================================================================Constructor bigint  ' default constructorthis.s = chr(0,0,0,0)   ' zeroEnd ConstructorConstructor bigint (Byref rhs As Long)this = rhsEnd ConstructorConstructor bigint (Byref rhs As ULong)this = rhsEnd Constructor'Constructor bigint (Byref rhs As Integer)'this = rhs'End Constructor''Constructor bigint (Byref rhs As Uinteger)'this = rhs'End ConstructorConstructor bigint (Byref rhs As Longint)this = rhsEnd ConstructorConstructor bigint (Byref rhs As ULongint)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 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 = *ipEnd 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 = *uipIf (128 And this.s[3]) Then this.s &= Chr(0,0,0,0) ' make it positiveEnd 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 = *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          ' 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 IfEnd Operator'================================================================' create a string for printing a bigintOperator bigint.cast() As StringDim As String tt = unpack(this)Return tEnd OperatorOperator 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) ' twoDim Shared As bigint bigint_m1=-1Dim Shared As bigint bigint_0=0Dim Shared As bigint bigint_1=1Dim Shared As bigint bigint_2=2'-----------------------------------------------------------------------' 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.'================================================================'                  ARITHMETRIC OPERATIONS'================================================================' 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<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 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'================================================================' additionFunction 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 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<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 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 - 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 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'=======================================================================' integer divide, a \ b, return div and modSub 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 fEnd Function'======================================================================='exponentiation modulusFunction 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 rEnd Function'======================================================================='               BIT FUNCTIONS'=======================================================================' 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'================================================================' 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'================================================================' 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 cyclesEnd 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 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=bigint_0 Then Return 0 ' is zero    Return 1 ' is positiveEnd 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 bigintdim as bigint a,bdivmod(x,y,a,b)Return aEnd OperatorOperator / (Byref x As bigint, Byref y As bigint) As bigintprint "no floating point division in Big_Integer"sleep: endreturn 0End OperatorOperator mod (Byref x As bigint, Byref y As bigint) as bigintdim as bigint c,ddivmod(x,y,c,d)return dEnd 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)dim as bigint ddivmod(this,rhs,this,d)End OperatorOperator bigint.mod= (Byref rhs as bigint)dim as bigint cdivmod(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 0end function '----------------------------------------------------------------------Operator = (Byref a As bigint, Byref b As bigint) As IntegerDim As Integer c = bigint_compare(a, b)if c=0 then return -1 else return 0End Operator'----------------------------------------------------------------Operator <> (Byref a As bigint, Byref b As bigint) As IntegerDim As Integer c = bigint_compare(a, b)if c=0 then return 0 else return -1End Operator'----------------------------------------------------------------Operator < (Byref a As bigint, Byref b As bigint) As IntegerDim As Integer c = bigint_compare(a,b)if c = 1 then return -1 else return 0End Operator'----------------------------------------------------------------Operator > (Byref a As bigint, Byref b As bigint) As IntegerDim As Integer c = bigint_compare(a,b)if c = -1 then return -1 else return 0End Operator'----------------------------------------------------------------Operator <= (Byref a As bigint, Byref b As bigint) As IntegerDim As Integer c = bigint_compare(a,b)if c = 1 or c = 0 then return -1 else return 0End Operator'----------------------------------------------------------------Operator >= (Byref a As bigint, Byref b As bigint) As IntegerDim As Integer c = bigint_compare(a, b)if c = -1 or c = 0 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'================================================================' 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'=================================================================' 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`
stephanbrunker
Posts: 62
Joined: Nov 02, 2013 14:57

### Re: Arbitrarily Big Integer Routines

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 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<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 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'===============================================================' 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 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 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 hEnd Function'================================================================' convert a bigint to octalFunction 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 cEnd Function'===================================================================' convert a bigint to integerFunction 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'===================================================================' convert a bigint to LongintFunction 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'===================================================================' convert a bigint to doubleFunction 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 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    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 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    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 bEnd 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

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-beispiele/mathematik/einfache-schnelle-bigint-erweiterung-284.html
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

### Re: Arbitrarily Big Integer Routines

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[i] + carry FreeBasic has only very limited knowledge of optimization, the constant part needs to be after the * to get optimized.product = d[i] * 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[i] 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 file=replacement1.bas]
'----------------------------------------------------------------
' 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[/code]
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: 184
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

### Re: Arbitrarily Big Integer Routines

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=0Dim 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 IfNextSleep`
srvaldez
Posts: 2482
Joined: Sep 25, 2005 21:54

### Re: Arbitrarily Big Integer Routines

hello Makoto WATANABE
the version you mention does not work for me either, use the version posted about 3 posts above viewtopic.php?p=192857#p192857
that version works for me.
Makoto WATANABE
Posts: 184
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

### Re: Arbitrarily Big Integer Routines

Dear srvaldez

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

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: 184
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

### Re: Arbitrarily Big Integer Routines

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/viewtopic.php?f=3&t=23225&p=226650#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

@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: 48
Joined: Mar 02, 2016 12:44

### Re: Arbitrarily Big Integer Routines

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

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: 184
Joined: Apr 10, 2010 11:41
Location: Japan
Contact:

### Re: Arbitrarily Big Integer Routines

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

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 * bprint cend`

then you can call your executable:
bigintcmd.exe "5555555555555555" "88888888888888888888888888"

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

### Re: Arbitrarily Big Integer Routines

Dear stephanbrunker

I'd like to enter numbers using "Input" command instead of the command line.
Please teach me a sample program.
ur_naz
Posts: 48
Joined: Mar 02, 2016 12:44

### Re: Arbitrarily Big Integer Routines

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