# bignum library in pure Tcl [VERSION 7Sep2004] # Copyright (C) 2004 Salvatore Sanfilippo # Copyright (C) 2004 Arjen Markus # # LICENSE # # This software is: # Copyright (C) 2004 Salvatore Sanfilippo # Copyright (C) 2004 Arjen Markus # 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 ... ] # 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. 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