# 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 eval ::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.
set ::math::bignum::atombits 16
set ::math::bignum::atombase [expr {1 << $::math::bignum::atombits}]
set ::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)
proc ::math::bignum::max {a b} {
expr {($a > $b) ? $a : $b}
}
# Return the min between a and b (not bignums)
proc ::math::bignum::min {a 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'.
proc ::math::bignum::zero {{value 0}} {
set v [list bignum 0 0]
while { $value > 1 } {
lappend v 0
incr value -1
}
return $v
}
# Get the bignum sign
proc ::math::bignum::sign bignum {
lindex $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.
proc ::math::bignum::atom {bignum 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.
proc ::math::bignum::setatom {bignumvar i 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
proc ::math::bignum::setsign {bignumvar sign} {
upvar 1 $bignumvar bignum
lset bignum 1 $sign
}
# Remove trailing atoms with a value of zero
# The normalized bignum is returned
proc ::math::bignum::normalize bignumvar {
upvar 1 $bignumvar bignum
set atoms [expr {[llength $bignum]-2}]
set i [expr {$atoms+1}]
while {$atoms && [lindex $bignum $i] == 0} {
set bignum [lrange $bignum 0 end-1]
incr atoms -1
incr i -1
}
if {!$atoms} {
set bignum [list bignum 0 0]
}
return $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|
#
proc ::math::bignum::abscmp {a 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
#
proc ::math::bignum::cmp {a 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.
proc ::math::bignum::iszero z {
set z [_treat $z]
expr {[llength $z] == 3 && [lindex $z 2] == 0}
}
# Comparison facilities
proc ::math::bignum::lt {a b} {expr {[::math::bignum::cmp $a $b] < 0}}
proc ::math::bignum::le {a b} {expr {[::math::bignum::cmp $a $b] <= 0}}
proc ::math::bignum::gt {a b} {expr {[::math::bignum::cmp $a $b] > 0}}
proc ::math::bignum::ge {a b} {expr {[::math::bignum::cmp $a $b] >= 0}}
proc ::math::bignum::eq {a b} {expr {[::math::bignum::cmp $a $b] == 0}}
proc ::math::bignum::ne {a b} {expr {[::math::bignum::cmp $a $b] != 0}}
########################### Addition / Subtraction #############################
# Add two bignums, don't care about the sign.
proc ::math::bignum::rawAdd {a 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.
proc ::math::bignum::rawSub {a 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.
proc ::math::bignum::add {a 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.
proc ::math::bignum::sub {a 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 #################################
set ::math::bignum::karatsubaThreshold 32
# Multiplication. Calls Karatsuba that calls Base multiplication under
# a given threshold.
proc ::math::bignum::mul {a 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
proc ::math::bignum::kmul {a 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.
proc ::math::bignum::bmul {a 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.
proc ::math::bignum::lshiftAtoms {z 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.
proc ::math::bignum::rshiftAtoms {z 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
proc ::math::bignum::lshiftBits {z 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
proc ::math::bignum::rshiftBits {z 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.
proc ::math::bignum::lshift {z 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.
proc ::math::bignum::rshift {z 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'
proc ::math::bignum::setbit {bignumvar 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'
proc ::math::bignum::clearbit {bignumvar 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.
proc ::math::bignum::testbit {z 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
proc ::math::bignum::bitand {a 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
proc ::math::bignum::bitxor {a 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
proc ::math::bignum::bitor {a 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'.
proc ::math::bignum::bits z {
set atoms [::math::bignum::atoms $z]
set bits [expr {($atoms-1)*$::math::bignum::atombits}]
set atom [lindex $z [expr {$atoms+1}]]
while {$atom} {
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}]
set r [list bignum 0 0]
set q [::math::bignum::zero [expr {[llength $n]-2}]]
while {$bit >= 0} {
set b_atom [expr {($bit / $::math::bignum::atombits) + 2}]
set b_bit [expr {1 << ($bit & ($::math::bignum::atombits-1))}]
set r [::math::bignum::lshiftBits $r 1]
if {[lindex $n $b_atom]&$b_bit} {
lset r 2 [expr {[lindex $r 2] | 1}]
}
if {[::math::bignum::abscmp $r $d] >= 0} {
set r [::math::bignum::rawSub $r $d]
lset 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).
proc ::math::bignum::rawDivByAtom {n 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.
proc ::math::bignum::divqr {n 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.
proc ::math::bignum::div {n d} {
lindex [::math::bignum::divqr $n $d] 0
}
# Like divqr, but only the remainder is returned.
proc ::math::bignum::rem {n d} {
lindex [::math::bignum::divqr $n $d] 1
}
# Modular reduction. Returns N modulo M
proc ::math::bignum::mod {n 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
proc ::math::bignum::isodd n {
expr {[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
proc ::math::bignum::pow {b 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
proc ::math::bignum::powm {b e 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.
proc ::math::bignum::sqrt n {
if {[lindex $n 1]} {
error "Square root of a negative number"
}
set i [expr {(([::math::bignum::bits $n]-1)/2)+1}]
set b [expr {$i*2}] ; # Bit to set to get 2^i*2^i
set r [::math::bignum::zero] ; # guess
set x [::math::bignum::zero] ; # guess^2
set s [::math::bignum::zero] ; # guess^2 backup
set t [::math::bignum::zero] ; # intermediate result
for {} {$i >= 0} {incr i -1; incr b -2} {
::math::bignum::setbit t $b
set x [::math::bignum::rawAdd $s $t]
::math::bignum::clearbit t $b
if {[::math::bignum::abscmp $x $n] <= 0} {
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
}
set r [::math::bignum::rshiftBits $r $shift]
return $r
}
############################ Convertion to/from string #########################
# The string representation charset. Max base is 36
set ::math::bignum::cset "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].
proc ::math::bignum::tostr {z {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'.
proc ::math::bignum::fromstr {str {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
#
proc ::math::bignum::_treat {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 eval ::math::bignum {
namespace export *
}
# Announce the package
package provide math::bignum 3.1.1