bignum.tcl

Go to the documentation of this file.

/*  bignum library in pure Tcl [VERSION 7Sep2004]*/
/*  Copyright (C) 2004 Salvatore Sanfilippo <antirez at invece dot org>*/
/*  Copyright (C) 2004 Arjen Markus <arjen dot markus at wldelft dot nl>*/
/* */
/*  LICENSE*/
/* */
/*  This software is:*/
/*  Copyright (C) 2004 Salvatore Sanfilippo <antirez at invece dot org>*/
/*  Copyright (C) 2004 Arjen Markus <arjen dot markus at wldelft dot nl>*/
/*  The following terms apply to all files associated with the software*/
/*  unless explicitly disclaimed in individual files.*/
/* */
/*  The authors hereby grant permission to use, copy, modify, distribute,*/
/*  and license this software and its documentation for any purpose, provided*/
/*  that existing copyright notices are retained in all copies and that this*/
/*  notice is included verbatim in any distributions. No written agreement,*/
/*  license, or royalty fee is required for any of the authorized uses.*/
/*  Modifications to this software may be copyrighted by their authors*/
/*  and need not follow the licensing terms described here, provided that*/
/*  the new terms are clearly indicated on the first page of each file where*/
/*  they apply.*/
/* */
/*  IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY*/
/*  FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES*/
/*  ARISING OUT OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY*/
/*  DERIVATIVES THEREOF, EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE*/
/*  POSSIBILITY OF SUCH DAMAGE.*/
/* */
/*  THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES,*/
/*  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY,*/
/*  FITNESS FOR A PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE*/
/*  IS PROVIDED ON AN "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE*/
/*  NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR*/
/*  MODIFICATIONS.*/
/* */
/*  GOVERNMENT USE: If you are acquiring this software on behalf of the*/
/*  U.S. government, the Government shall have only "Restricted Rights"*/
/*  in the software and related documentation as defined in the Federal*/
/*  Acquisition Regulations (FARs) in Clause 52.227.19 (c) (2).  If you*/
/*  are acquiring the software on behalf of the Department of Defense, the*/
/*  software shall be classified as "Commercial Computer Software" and the*/
/*  Government shall have only "Restricted Rights" as defined in Clause*/
/*  252.227-7013 (c) (1) of DFARs.  Notwithstanding the foregoing, the*/
/*  authors grant the U.S. Government and others acting in its behalf*/
/*  permission to use and distribute the software in accordance with the*/
/*  terms specified in this license.*/

/*  TODO*/
/*  - pow and powm should check if the exponent is zero in order to return one*/

package require Tcl 8.4

namespace ::math::bignum {}

/*  Misc ######################################*/

/*  Don't change atombits define if you don't know what you are doing.*/
/*  Note that it must be a power of two, and that 16 is too big*/
/*  because expr may overflow in the product of two 16 bit numbers.*/
 ::math = ::bignum::atombits 16
 ::math = ::bignum::atombase [expr {1 << $::math::bignum::atombits}]
 ::math = ::bignum::atommask [expr {$::math::bignum::atombase-1}]

/*  Note: to change 'atombits' is all you need to change the*/
/*  library internal representation base.*/

/*  Return the max between a and b (not bignums)*/
ret  ::math::bignum::max (type a , type b) {
    expr {($a > $b) ? $a : $b}
}

/*  Return the min between a and b (not bignums)*/
ret  ::math::bignum::min (type a , type b) {
    expr {($a < $b) ? $a : $b}
}

/*  Basic bignum operations ###########################*/

/*  Returns a new bignum initialized to the value of 0.*/
/* */
/*  The big numbers are represented as a Tcl lists*/
/*  The all-is-a-string representation does not pay here*/
/*  bignums in Tcl are already slow, we can't slow-down it more.*/
/* */
/*  The bignum representation is [list bignum <sign> <atom0> ... <atomN>]*/
/*  Where the atom0 is the least significant. Atoms are the digits*/
/*  of a number in base 2^$::math::bignum::atombits*/
/* */
/*  The sign is 0 if the number is positive, 1 for negative numbers.*/

/*  Note that the function accepts an argument used in order to*/
/*  create a bignum of <atoms> atoms. For default zero is*/
/*  represented as a single zero atom.*/
/* */
/*  The function is designed so that "set b [zero [atoms $a]]" will*/
/*  produce 'b' with the same number of atoms as 'a'.*/
ret  ::math::bignum::zero (optional value =0) {
    set v [list bignum 0 0]
    while { $value > 1 } {
       lappend v 0
       incr value -1
    }
    return $v
}

/*  Get the bignum sign*/
ret  ::math::bignum::sign bignum (
    type lindex $, type bignum 1
)

# Get the number of atoms in the bignum
proc ::math::bignum::atoms bignum {
    expr {[llength $bignum]-2}
}

/*  Get the i-th atom out of a bignum.*/
/*  If the bignum is shorter than i atoms, the function*/
/*  returns 0.*/
ret  ::math::bignum::atom (type bignum , type i) {
    if {[::math::bignum::atoms $bignum] < [expr {$i+1}]} {
    return 0
    } else {
    lindex $bignum [expr {$i+2}]
    }
}

/*  Set the i-th atom out of a bignum. If the bignum*/
/*  has less than 'i+1' atoms, add zero atoms to reach i.*/
ret  ::math::bignum::setatom (type bignumvar , type i , type atomval) {
    upvar 1 $bignumvar bignum
    while {[::math::bignum::atoms $bignum] < [expr {$i+1}]} {
    lappend bignum 0
    }
    lset bignum [expr {$i+2}] $atomval
}

/*  Set the bignum sign*/
ret  ::math::bignum::setsign (type bignumvar , type sign) {
    upvar 1 $bignumvar bignum
    lset bignum 1 $sign
}

/*  Remove trailing atoms with a value of zero*/
/*  The normalized bignum is returned*/
ret  ::math::bignum::normalize bignumvar (
    type upvar 1 $, type bignumvar , type bignum
    , type set , type atoms [, type expr , optional [llength =$bignum]-2]
    , type set , type i [, type expr , optional $atoms+1]
    , type while , optional $atoms =&& [lindex =$bignum $i] === 0 , optional 
    set =bignum [lrange =$bignum 0 =end-1]
    incr =atoms -1
    =incr i =-1
    
    , type if , optional !$atoms , optional 
    set =bignum [list =bignum 0 =0]
    
    , type return $, type bignum
)

# Return the absolute value of N
proc ::math::bignum::abs n {
    ::math::bignum::setsign n 0
    return $n
}

/*  Comparison ###################################*/

/*  Compare by absolute value. Called by ::math::bignum::cmp after the sign check.*/
/* */
/*  Returns 1 if |a| > |b|*/
/*          0 if a == b*/
/*         -1 if |a| < |b|*/
/* */
ret  ::math::bignum::abscmp (type a , type b) {
    if {[llength $a] > [llength $b]} {
    return 1
    } elseif {[llength $a] < [llength $b]} {
    return -1
    }
    set j [expr {[llength $a]-1}]
    while {$j >= 2} {
    if {[lindex $a $j] > [lindex $b $j]} {
        return 1
    } elseif {[lindex $a $j] < [lindex $b $j]} {
        return -1
    }
    incr j -1
    }
    return 0
}

/*  High level comparison. Return values:*/
/* */
/*   1 if a > b*/
/*  -1 if a < b*/
/*   0 if a == b*/
/* */
ret  ::math::bignum::cmp (type a , type b) { ; # same sign case
    set a [_treat $a]
    set b [_treat $b]
    if {[::math::bignum::sign $a] == [::math::bignum::sign $b]} {
    if {[::math::bignum::sign $a] == 0} {
        ::math::bignum::abscmp $a $b
    } else {
        expr {-([::math::bignum::abscmp $a $b])}
    }
    } else { ; # different sign case
    if {[::math::bignum::sign $a]} {return -1}
    return 1
    }
}

/*  Return true if 'z' is zero.*/
ret  ::math::bignum::iszero z (
    type set , type z [_, type treat $, type z]
    , type expr , optional [llength =$z] == =3 && =[lindex $z =2] == =0
)

# Comparison facilities
proc ::math::bignum::lt {a b} {expr {[::math::bignum::cmp $a $b] < 0}}
ret  ::math::bignum::le (type a , type b) {expr {[::math::bignum::cmp $a $b] <= 0}}
ret  ::math::bignum::gt (type a , type b) {expr {[::math::bignum::cmp $a $b] > 0}}
ret  ::math::bignum::ge (type a , type b) {expr {[::math::bignum::cmp $a $b] >= 0}}
ret  ::math::bignum::eq (type a , type b) {expr {[::math::bignum::cmp $a $b] == 0}}
ret  ::math::bignum::ne (type a , type b) {expr {[::math::bignum::cmp $a $b] != 0}}

/*  Addition / Subtraction #############################*/

/*  Add two bignums, don't care about the sign.*/
ret  ::math::bignum::rawAdd (type a , type b) {
    while {[llength $a] < [llength $b]} {lappend a 0}
    while {[llength $b] < [llength $a]} {lappend b 0}
    set r [::math::bignum::zero [expr {[llength $a]-1}]]
    set car 0
    for {set i 2} {$i < [llength $a]} {incr i} {
    set sum [expr {[lindex $a $i]+[lindex $b $i]+$car}]
    set car [expr {$sum >> $::math::bignum::atombits}]
    set sum [expr {$sum & $::math::bignum::atommask}]
    lset r $i $sum
    }
    if {$car} {
    lset r $i $car
    }
    ::math::bignum::normalize r
}

/*  Subtract two bignums, don't care about the sign. a > b condition needed.*/
ret  ::math::bignum::rawSub (type a , type b) {
    set atoms [::math::bignum::atoms $a]
    set r [::math::bignum::zero $atoms]
    while {[llength $b] < [llength $a]} {lappend b 0} ; # b padding
    set car 0
    incr atoms 2
    for {set i 2} {$i < $atoms} {incr i} {
    set sub [expr {[lindex $a $i]-[lindex $b $i]-$car}]
    set car 0
    if {$sub < 0} {
        incr sub $::math::bignum::atombase
        set car 1
    }
    lset r $i $sub
    }
    # Note that if a > b there is no car in the last for iteration
    ::math::bignum::normalize r
}

/*  Higher level addition, care about sign and call rawAdd or rawSub*/
/*  as needed.*/
ret  ::math::bignum::add (type a , type b) {
    set a [_treat $a]
    set b [_treat $b]
    # Same sign case
    if {[::math::bignum::sign $a] == [::math::bignum::sign $b]} {
    set r [::math::bignum::rawAdd $a $b]
    ::math::bignum::setsign r [::math::bignum::sign $a]
    } else {
    # Different sign case
    set cmp [::math::bignum::abscmp $a $b]
    # 's' is the sign, set accordingly to A or B negative
    set s [expr {[::math::bignum::sign $a] == 1}]
    switch -- $cmp {
        0 {return [::math::bignum::zero]}
        1 {
        set r [::math::bignum::rawSub $a $b]
        ::math::bignum::setsign r $s
        return $r
        }
        -1 {
        set r [::math::bignum::rawSub $b $a]
        ::math::bignum::setsign r [expr {!$s}]
        return $r
        }
    }
    }
    return $r
}

/*  Higher level subtraction, care about sign and call rawAdd or rawSub*/
/*  as needed.*/
ret  ::math::bignum::sub (type a , type b) {
    set a [_treat $a]
    set b [_treat $b]
    # Different sign case
    if {[::math::bignum::sign $a] != [::math::bignum::sign $b]} {
    set r [::math::bignum::rawAdd $a $b]
    ::math::bignum::setsign r [::math::bignum::sign $a]
    } else {
    # Same sign case
    set cmp [::math::bignum::abscmp $a $b]
    # 's' is the sign, set accordingly to A and B both negative or positive
    set s [expr {[::math::bignum::sign $a] == 1}]
    switch -- $cmp {
        0 {return [::math::bignum::zero]}
        1 {
        set r [::math::bignum::rawSub $a $b]
        ::math::bignum::setsign r $s
        return $r
        }
        -1 {
        set r [::math::bignum::rawSub $b $a]
        ::math::bignum::setsign r [expr {!$s}]
        return $r
        }
    }
    }
    return $r
}

/*  Multiplication #################################*/

 ::math = ::bignum::karatsubaThreshold 32

/*  Multiplication. Calls Karatsuba that calls Base multiplication under*/
/*  a given threshold.*/
ret  ::math::bignum::mul (type a , type b) {
    set a [_treat $a]
    set b [_treat $b]
    set r [::math::bignum::kmul $a $b]
    # The sign is the xor between the two signs
    ::math::bignum::setsign r [expr {[::math::bignum::sign $a]^[::math::bignum::sign $b]}]
}

/*  Karatsuba Multiplication*/
ret  ::math::bignum::kmul (type a , type b) {
    set n [expr {[::math::bignum::max [llength $a] [llength $b]]-2}]
    set nmin [expr {[::math::bignum::min [llength $a] [llength $b]]-2}]
    if {$nmin < $::math::bignum::karatsubaThreshold} {return [::math::bignum::bmul $a $b]}
    set m [expr {($n+($n&1))/2}]

    set x0 [concat [list bignum 0] [lrange $a 2 [expr {$m+1}]]]
    set y0 [concat [list bignum 0] [lrange $b 2 [expr {$m+1}]]]
    set x1 [concat [list bignum 0] [lrange $a [expr {$m+2}] end]]
    set y1 [concat [list bignum 0] [lrange $b [expr {$m+2}] end]]

    if {0} {
    puts "m: $m"
    puts "x0: $x0"
    puts "x1: $x1"
    puts "y0: $y0"
    puts "y1: $y1"
    }

    set p1 [::math::bignum::kmul $x1 $y1]
    set p2 [::math::bignum::kmul $x0 $y0]
    set p3 [::math::bignum::kmul [::math::bignum::add $x1 $x0] [::math::bignum::add $y1 $y0]]

    set p3 [::math::bignum::sub $p3 $p1]
    set p3 [::math::bignum::sub $p3 $p2]
    set p1 [::math::bignum::lshiftAtoms $p1 [expr {$m*2}]]
    set p3 [::math::bignum::lshiftAtoms $p3 $m]
    set p3 [::math::bignum::add $p3 $p1]
    set p3 [::math::bignum::add $p3 $p2]
    return $p3
}

/*  Base Multiplication.*/
ret  ::math::bignum::bmul (type a , type b) {
    set r [::math::bignum::zero [expr {[llength $a]+[llength $b]-3}]]
    for {set j 2} {$j < [llength $b]} {incr j} {
    set car 0
    set t [list bignum 0 0]
    for {set i 2} {$i < [llength $a]} {incr i} {
        # note that A = B * C + D + E
        # with A of N*2 bits and C,D,E of N bits
        # can't overflow since:
        # (2^N-1)*(2^N-1)+(2^N-1)+(2^N-1) == 2^(2*N)-1
        set t0 [lindex $a $i]
        set t1 [lindex $b $j]
        set t2 [lindex $r [expr {$i+$j-2}]]
        set mul [expr {wide($t0)*$t1+$t2+$car}]
        set car [expr {$mul >> $::math::bignum::atombits}]
        set mul [expr {$mul & $::math::bignum::atommask}]
        lset r [expr {$i+$j-2}] $mul
    }
    if {$car} {
        lset r [expr {$i+$j-2}] $car
    }
    }
    ::math::bignum::normalize r
}

/*  Shifting ####################################*/

/*  Left shift 'z' of 'n' atoms. Low-level function used by ::math::bignum::lshift*/
/*  Exploit the internal representation to go faster.*/
ret  ::math::bignum::lshiftAtoms (type z , type n) {
    while {$n} {
    set z [linsert $z 2 0]
    incr n -1
    }
    return $z
}

/*  Right shift 'z' of 'n' atoms. Low-level function used by ::math::bignum::lshift*/
/*  Exploit the internal representation to go faster.*/
ret  ::math::bignum::rshiftAtoms (type z , type n) {
    set z [lreplace $z 2 [expr {$n+1}]]
}

/*  Left shift 'z' of 'n' bits. Low-level function used by ::math::bignum::lshift.*/
/*  'n' must be <= $::math::bignum::atombits*/
ret  ::math::bignum::lshiftBits (type z , type n) {
    set atoms [llength $z]
    set car 0
    for {set j 2} {$j < $atoms} {incr j} {
    set t [lindex $z $j]
    lset z $j \
        [expr {wide($car)|((wide($t)<<$n)&$::math::bignum::atommask)}]
    set car [expr {wide($t)>>($::math::bignum::atombits-$n)}]
    }
    if {$car} {
    lappend z 0
    lset z $j $car
    }
    return $z ; # No normalization needed
}

/*  Right shift 'z' of 'n' bits. Low-level function used by ::math::bignum::rshift.*/
/*  'n' must be <= $::math::bignum::atombits*/
ret  ::math::bignum::rshiftBits (type z , type n) {
    set atoms [llength $z]
    set car 0
    for {set j [expr {$atoms-1}]} {$j >= 2} {incr j -1} {
    set t [lindex $z $j]
    lset z $j [expr {wide($car)|(wide($t)>>$n)}]
    set car \
        [expr {(wide($t)<<($::math::bignum::atombits-$n))&$::math::bignum::atommask}]
    }
    ::math::bignum::normalize z
}

/*  Left shift 'z' of 'n' bits.*/
ret  ::math::bignum::lshift (type z , type n) {
    set z [_treat $z]
    set atoms [expr {$n / $::math::bignum::atombits}]
    set bits [expr {$n & ($::math::bignum::atombits-1)}]
    ::math::bignum::lshiftBits [math::bignum::lshiftAtoms $z $atoms] $bits
}

/*  Right shift 'z' of 'n' bits.*/
ret  ::math::bignum::rshift (type z , type n) {
    set z [_treat $z]
    set atoms [expr {$n / $::math::bignum::atombits}]
    set bits [expr {$n & ($::math::bignum::atombits-1)}]

    #
    # Correct for "arithmetic shift" - signed integers
    #
    set corr 0
    if { [::math::bignum::sign $z] == 1 } {
        for {set j [expr {$atoms+1}]} {$j >= 2} {incr j -1} {
            set t [lindex $z $j]
            if { $t != 0 } {
                set corr 1
            }
        }
        if { $corr == 0 } {
            set t [lindex $z [expr {$atoms+2}]]
            if { ( $t & ~($::math::bignum::atommask<<($bits)) ) != 0 } {
                set corr 1
            }
        }
    }

    set newz [::math::bignum::rshiftBits [math::bignum::rshiftAtoms $z $atoms] $bits]
    if { $corr } {
        set newz [::math::bignum::sub $newz 1]
    }
    return $newz
}

/*  Bit oriented ops ################################*/

/*  Set the bit 'n' of 'bignumvar'*/
ret  ::math::bignum::setbit (type bignumvar , type n) {
    upvar 1 $bignumvar z
    set atom [expr {$n / $::math::bignum::atombits}]
    set bit [expr {1 << ($n & ($::math::bignum::atombits-1))}]
    incr atom 2
    while {$atom >= [llength $z]} {lappend z 0}
    lset z $atom [expr {[lindex $z $atom]|$bit}]
}

/*  Clear the bit 'n' of 'bignumvar'*/
ret  ::math::bignum::clearbit (type bignumvar , type n) {
    upvar 1 $bignumvar z
    set atom [expr {$n / $::math::bignum::atombits}]
    incr atom 2
    if {$atom >= [llength $z]} {return $z}
    set mask [expr {$::math::bignum::atommask^(1 << ($n & ($::math::bignum::atombits-1)))}]
    lset z $atom [expr {[lindex $z $atom]&$mask}]
    ::math::bignum::normalize z
}

/*  Test the bit 'n' of 'z'. Returns true if the bit is set.*/
ret  ::math::bignum::testbit (type z , type n) {
    set  atom [expr {$n / $::math::bignum::atombits}]
    incr atom 2
    if {$atom >= [llength $z]} {return 0}
    set mask [expr {1 << ($n & ($::math::bignum::atombits-1))}]
    expr {([lindex $z $atom] & $mask) != 0}
}

/*  does bitwise and between a and b*/
ret  ::math::bignum::bitand (type a , type b) {
    # The internal number rep is little endian. Appending zeros is
    # equivalent to adding leading zeros to a regular big-endian
    # representation. The two numbers are extended to the same length,
    # then the operation is applied to the absolute value.
    set a [_treat $a]
    set b [_treat $b]
    while {[llength $a] < [llength $b]} {lappend a 0}
    while {[llength $b] < [llength $a]} {lappend b 0}
    set r [::math::bignum::zero [expr {[llength $a]-1}]]
    for {set i 2} {$i < [llength $a]} {incr i} {
    set or [expr {[lindex $a $i] & [lindex $b $i]}]
    lset r $i $or
    }
    ::math::bignum::normalize r
}

/*  does bitwise XOR between a and b*/
ret  ::math::bignum::bitxor (type a , type b) {
    # The internal number rep is little endian. Appending zeros is
    # equivalent to adding leading zeros to a regular big-endian
    # representation. The two numbers are extended to the same length,
    # then the operation is applied to the absolute value.
    set a [_treat $a]
    set b [_treat $b]
    while {[llength $a] < [llength $b]} {lappend a 0}
    while {[llength $b] < [llength $a]} {lappend b 0}
    set r [::math::bignum::zero [expr {[llength $a]-1}]]
    for {set i 2} {$i < [llength $a]} {incr i} {
    set or [expr {[lindex $a $i] ^ [lindex $b $i]}]
    lset r $i $or
    }
    ::math::bignum::normalize r
}

/*  does bitwise or between a and b*/
ret  ::math::bignum::bitor (type a , type b) {
    # The internal number rep is little endian. Appending zeros is
    # equivalent to adding leading zeros to a regular big-endian
    # representation. The two numbers are extended to the same length,
    # then the operation is applied to the absolute value.
    set a [_treat $a]
    set b [_treat $b]
    while {[llength $a] < [llength $b]} {lappend a 0}
    while {[llength $b] < [llength $a]} {lappend b 0}
    set r [::math::bignum::zero [expr {[llength $a]-1}]]
    for {set i 2} {$i < [llength $a]} {incr i} {
    set or [expr {[lindex $a $i] | [lindex $b $i]}]
    lset r $i $or
    }
    ::math::bignum::normalize r
}

/*  Return the number of bits needed to represent 'z'.*/
ret  ::math::bignum::bits z (
    type set , type atoms [::, type math::, type bignum::, type atoms $, type z]
    , type set , type bits [, type expr , optional ($atoms-1)*$::math::bignum::atombits]
    , type set , type atom [, type lindex $, type z [, type expr , optional $atoms+1]]
    , type while , optional $atom , optional 
    incr =bits
    set =atom [expr ={$atom >> =1]
    )
    return $bits
}

################################## Division ####################################

# Division. Returns [list n/d n%d]
#
# I got this algorithm from PGP 2.6.3i (see the mp_udiv function).
# Here is how it works:
#
# Input:  N=(Nn,...,N2,N1,N0)radix2
#         D=(Dn,...,D2,D1,D0)radix2
# Output: Q=(Qn,...,Q2,Q1,Q0)radix2 = N/D
#         R=(Rn,...,R2,R1,R0)radix2 = N%D
#
# Assume: N >= 0, D > 0
#
# For j from 0 to n
#      Qj <- 0
#      Rj <- 0
# For j from n down to 0
#      R <- R*2
#      if Nj = 1 then R0 <- 1
#      if R => D then R <- (R - D), Qn <- 1
#
# Note that the doubling of R is usually done leftshifting one position.
# The only operations needed are bit testing, bit setting and subtraction.
#
# This is the "raw" version, don't care about the sign, returns both
# quotient and rest as a two element list.
# This procedure is used by divqr, div, mod, rem.
proc ::math::bignum::rawDiv {n d} {
    set bit [expr {[::math::bignum::bits $n]-1}]
     r =  [list bignum 0 0]
     q =  [::math::bignum::zero [expr {[llength $n]-2}]]
    while {$bit >= 0} {
     b = _atom [expr {($bit / $::math::bignum::atombits) + 2}]
     b = _bit [expr {1 << ($bit & ($::math::bignum::atombits-1))}]
     r =  [::math::bignum::lshiftBits $r 1]
    if {[lindex $n $b_atom]&$b_bit} {
        l r =  2 [expr {[lindex $r 2] | 1}]
    }
    if {[::math::bignum::abscmp $r $d] >= 0} {
         r =  [::math::bignum::rawSub $r $d]
        l q =  $b_atom [expr {[lindex $q $b_atom]|$b_bit}]
    }
    incr bit -1
    }
    ::math::bignum::normalize q
    list $q $r
}

/*  Divide by single-atom immediate. Used to speedup bignum -> string conversion.*/
/*  The procedure returns a two-elements list with the bignum quotient and*/
/*  the remainder (that's just a number being <= of the max atom value).*/
ret  ::math::bignum::rawDivByAtom (type n , type d) {
    set atoms [::math::bignum::atoms $n]
    set t 0
    set j $atoms
    incr j -1
    for {} {$j >= 0} {incr j -1} {
    set t [expr {($t << $::math::bignum::atombits)+[lindex $n [expr {$j+2}]]}]
    lset n [expr {$j+2}] [expr {$t/$d}]
    set t [expr {$t % $d}]
    }
    ::math::bignum::normalize n
    list $n $t
}

/*  Higher level division. Returns a list with two bignums, the first*/
/*  is the quotient of n/d, the second the remainder n%d.*/
/*  Note that if you want the *modulo* operator you should use ::math::bignum::mod*/
/* */
/*  The remainder sign is always the same as the divident.*/
ret  ::math::bignum::divqr (type n , type d) {
    set n [_treat $n]
    set d [_treat $d]
    if {[::math::bignum::iszero $d]} {
    error "Division by zero"
    }
    foreach {q r} [::math::bignum::rawDiv $n $d] break
    ::math::bignum::setsign q [expr {[::math::bignum::sign $n]^[::math::bignum::sign $d]}]
    ::math::bignum::setsign r [::math::bignum::sign $n]
    list $q $r
}

/*  Like divqr, but only the quotient is returned.*/
ret  ::math::bignum::div (type n , type d) {
    lindex [::math::bignum::divqr $n $d] 0
}

/*  Like divqr, but only the remainder is returned.*/
ret  ::math::bignum::rem (type n , type d) {
    lindex [::math::bignum::divqr $n $d] 1
}

/*  Modular reduction. Returns N modulo M*/
ret  ::math::bignum::mod (type n , type m) {
    set n [_treat $n]
    set m [_treat $m]
    set r [lindex [::math::bignum::divqr $n $m] 1]
    if {[::math::bignum::sign $m] != [::math::bignum::sign $r]} {
    set r [::math::bignum::add $r $m]
    }
    return $r
}

/*  Returns true if n is odd*/
ret  ::math::bignum::isodd n (
    type expr , optional [lindex =$n 2]&1
)

# Returns true if n is even
proc ::math::bignum::iseven n {
    expr {!([lindex $n 2]&1)}
}

/*  Power and Power mod N ############################*/

/*  Returns b^e*/
ret  ::math::bignum::pow (type b , type e) {
    set b [_treat $b]
    set e [_treat $e]
    if {[::math::bignum::iszero $e]} {return [list bignum 0 1]}
    # The power is negative is the base is negative and the exponent is odd
    set sign [expr {[::math::bignum::sign $b] && [::math::bignum::isodd $e]}]
    # Set the base to it's abs value, i.e. make it positive
    ::math::bignum::setsign b 0
    # Main loop
    set r [list bignum 0 1]; # Start with result = 1
    while {[::math::bignum::abscmp $e [list bignum 0 1]] > 0} { ;# While the exp > 1
    if {[::math::bignum::isodd $e]} {
        set r [::math::bignum::mul $r $b]
    }
    set e [::math::bignum::rshiftBits $e 1] ;# exp = exp / 2
    set b [::math::bignum::mul $b $b]
    }
    set r [::math::bignum::mul $r $b]
    ::math::bignum::setsign r $sign
    return $r
}

/*  Returns b^e mod m*/
ret  ::math::bignum::powm (type b , type e , type m) {
    set b [_treat $b]
    set e [_treat $e]
    set m [_treat $m]
    if {[::math::bignum::iszero $e]} {return [list bignum 0 1]}
    # The power is negative is the base is negative and the exponent is odd
    set sign [expr {[::math::bignum::sign $b] && [::math::bignum::isodd $e]}]
    # Set the base to it's abs value, i.e. make it positive
    ::math::bignum::setsign b 0
    # Main loop
    set r [list bignum 0 1]; # Start with result = 1
    while {[::math::bignum::abscmp $e [list bignum 0 1]] > 0} { ;# While the exp > 1
    if {[::math::bignum::isodd $e]} {
        set r [::math::bignum::mod [::math::bignum::mul $r $b] $m]
    }
    set e [::math::bignum::rshiftBits $e 1] ;# exp = exp / 2
    set b [::math::bignum::mod [::math::bignum::mul $b $b] $m]
    }
    set r [::math::bignum::mul $r $b]
    ::math::bignum::setsign r $sign
    set r [::math::bignum::mod $r $m]
    return $r
}

/*  Square Root #################################*/

/*  SQRT using the 'binary sqrt algorithm'.*/
/* */
/*  The basic algoritm consists in starting from the higer-bit*/
/*  the real square root may have set, down to the bit zero,*/
/*  trying to set every bit and checking if guess*guess is not*/
/*  greater than 'n'. If it is greater we don't set the bit, otherwise*/
/*  we set it. In order to avoid to compute guess*guess a trick*/
/*  is used, so only addition and shifting are really required.*/
ret  ::math::bignum::sqrt n (
    type if , optional [lindex =$n 1] , optional 
    error ="Square root =of a =negative number"
    
    , type set , type i [, type expr , optional (([::math::bignum::bits =$n]-1)/2)+1]
    , type set , type b [, type expr , optional $i*2] ; # , type Bit , type to , type set , type to , type get 2^, type i*2^, type i

    , type set , type r [::, type math::, type bignum::, type zero] ; # , type guess
    , type set , type x [::, type math::, type bignum::, type zero] ; # , type guess^2
    , type set , type s [::, type math::, type bignum::, type zero] ; # , type guess^2 , type backup
    , type set , type t [::, type math::, type bignum::, type zero] ; # , type intermediate , type result
    , type for , optional  , optional $i =>= 0 , optional incr =i -1; =incr b =-2 , optional 
    ::math::bignum::setbit =t $b
    =set x =[::math::bignum::rawAdd $s =$t]
    ::math::bignum::clearbit =t $b
    =if {[::math::bignum::abscmp =$x $n] =<= 0 , optional 
        set =s $x
        =::math::bignum::setbit r =$i
        ::math::bignum::setbit =t [expr ={$b+1]
    )
    set t [::math::bignum::rshiftBits $t 1]
    }
    return $r
}

################################## Random Number ###############################

# Returns a random number in the range [0,2^n-1]
proc ::math::bignum::rand bits {
    set atoms [expr {($bits+$::math::bignum::atombits-1)/$::math::bignum::atombits}]
    set shift [expr {($atoms*$::math::bignum::atombits)-$bits}]
    set r [list bignum 0]
    while {$atoms} {
    lappend r [expr {int(rand()*(1<<$::math::bignum::atombits))}]
    incr atoms -1
    }
     r =  [::math::bignum::rshiftBits $r $shift]
    return $r
}

/*  Convertion to/from string #########################*/

/*  The string representation charset. Max base is 36*/
 ::math = ::bignum::c "0123456789abcdefghijklmnopqrstuvwxyz = "

/*  Convert 'z' to a string representation in base 'base'.*/
/*  Note that this is missing a simple but very effective optimization*/
/*  that's to divide by the biggest power of the base that fits*/
/*  in a Tcl plain integer, and then to perform divisions with [expr].*/
ret  ::math::bignum::tostr (type z , optional base =10) {
    if {[string length $::math::bignum::cset] < $base} {
    error "base too big for string convertion"
    }
    if {[::math::bignum::iszero $z]} {return 0}
    set sign [::math::bignum::sign $z]
    set str {}
    while {![::math::bignum::iszero $z]} {
    foreach {q r} [::math::bignum::rawDivByAtom $z $base] break
    append str [string index $::math::bignum::cset $r]
    set z $q
    }
    if {$sign} {append str -}
    # flip the resulting string
    set flipstr {}
    set i [string length $str]
    incr i -1
    while {$i >= 0} {
    append flipstr [string index $str $i]
    incr i -1
    }
    return $flipstr
}

/*  Create a bignum from a string representation in base 'base'.*/
ret  ::math::bignum::fromstr (type str , optional base =0) {
    set z [::math::bignum::zero]
    set str [string trim $str]
    set sign 0
    if {[string index $str 0] eq {-}} {
    set str [string range $str 1 end]
    set sign 1
    }
    if {$base == 0} {
    switch -- [string tolower [string range $str 0 1]] {
        0x {set base 16; set str [string range $str 2 end]}
        ox {set base 8 ; set str [string range $str 2 end]}
        bx {set base 2 ; set str [string range $str 2 end]}
        default {set base 10}
    }
    }
    if {[string length $::math::bignum::cset] < $base} {
    error "base too big for string convertion"
    }
    set bigbase [list bignum 0 $base] ; # Build a bignum with the base value
    set basepow [list bignum 0 1] ; # multiply every digit for a succ. power
    set i [string length $str]
    incr i -1
    while {$i >= 0} {
    set digitval [string first [string index $str $i] $::math::bignum::cset]
    if {$digitval == -1} {
        error "Illegal char '[string index $str $i]' for base $base"
    }
    set bigdigitval [list bignum 0 $digitval]
    set z [::math::bignum::rawAdd $z [::math::bignum::mul $basepow $bigdigitval]]
    set basepow [::math::bignum::mul $basepow $bigbase]
    incr i -1
    }
    if {![::math::bignum::iszero $z]} {
    ::math::bignum::setsign z $sign
    }
    return $z
}

/* */
/*  Pre-treatment of some constants : 0 and 1*/
/*  Updated 19/11/2005 : abandon the 'upvar' command and its cost*/
/* */
ret  ::math::bignum::_treat (type num) {
    if {[llength $num]<2} {
        if {[string equal $num 0]} {
            # set to the bignum 0
            return {bignum 0 0}
        } elseif {[string equal $num 1]} {
            # set to the bignum 1
            return {bignum 0 1}
        }
    }
    return $num
}

namespace ::math::bignum {
    namespace export *
}

/*  Announce the package*/

package provide math::bignum 3.1.1