Updated 2013-01-20 08:49:27 by pooryorick

I speed up more extensively adda and multa in doing computation in base 10000 instead of 10. Here is code
```global base
set base 10000

global base
set rs {}
set carry 0
if {[llength \$xs]>[llength \$ys]} {
set len [llength \$ys]
set flag 0
} else {
set len [llength \$xs]
set flag 1
}
set pos 0
foreach x \$xs y \$ys {
if {\$pos>=\$len} {
if {\$flag} {set x 0} else {set y 0}
}
set n [expr {\$x+\$y+\$carry}]
set carry [expr {\$n / \$base}]
lappend rs [expr {\$n % \$base}]
incr pos ;# who stole this line? restored 11/00 MM
}
if {\$carry} {lappend rs \$carry}
return \$rs
}

proc multi {xs a} {
global base
set rs {}
#puts "multi : xs=\$xs"
if { \$a == \$base } {
set rs \$xs
set rs [ linsert \$rs 0 0 ]
#puts "multi: rs=\$rs"
return \$rs
}
set carry 0
foreach x \$xs {
set n [expr {\$x*\$a+\$carry}]
set carry [expr {\$n / \$base}]
lappend rs [expr {\$n % \$base}]
}
if {\$carry} {lappend rs \$carry}
return \$rs
}

proc multa {xs ys} {
global base
set rs {}
set d \$ys
foreach x \$xs {
set n [multi \$d \$x]
set rs [adda \$rs \$n ]
set d [multi \$d \$base ]
}
return \$rs
}

# Convert an integer to a (list-formatted) number.
proc tolist {n} {
global base
set rs {}
if {[catch {incr n 0}]} {
set k 0
set nn {}
set rs {}
set l [string length \$n ]
incr l -1
for { set i \$l } { \$i >= 0 } { incr i -1 } {
set digit [ string index \$n \$i ]
incr k
set nn "\$digit\$nn"
if { \$k == 4} {
set k 0
lappend rs  \$nn
#puts -nonewline " \$nn"
set nn {}
}
}
if { \$nn != {} } {
lappend rs  \$nn
#puts -nonewline " \$nn"
}
#puts ""
} else {
while {\$n > 0} {
lappend rs [expr {\$n%\$base}]
set n [expr {\$n/\$base}]
}
}
#puts "tolist rs= \$rs"
return \$rs
}

# Convert a (list-formatted) number to a human-readable form.
# Will be an integer (modulo type-shimmering) if representable.
proc fromlist {xs} {
set r {}
foreach x \$xs {
if { \$x == 0 } {
set r "0000\$r"
} else {
set r \$x\$r
}
}
if {![string length \$r]} {set r 0}
return \$r
}

# DEMO CODE: Factorial.  Try calculating [fact 40] for a 48-digit number
proc fact {n} {
set this 1
for {set i 1} {\$i <= \$n} {incr i} {
set this [multa \$this [tolist \$i]]
}
return [fromlist \$this]
}

# DEMO CODE: Fibonacci Sequence.  Doesn't grow anything like as fast as factorial (in fact, two-term recurrences such as Fibonacci grow at quadratic rates),
# but still grows beyond 32-bit integers quite easily.
proc fib {n} {
set this 1
set last 0
for {set i 1} {\$i<\$n} {incr i} {
set last \$this
set this \$new
}
return [fromlist \$this]
}```

JJD

Here are a few procedures I put together to perform arbitrary precision arithmetic. If you are using these a lot, perhaps you should consider the mpexpr extension [1] instead.

The basic representation I use is a list of decimal digits (easy to convert to and from human-readable form) stored least-significant digit first (so the code handles the magnitude of the numbers implicitly.)

This all comes from a discussion I was having in email relating to some messages online about calculating factorials. I've decided to put the procs here so I can find them again. Feel free to optimise them, add in relevant links, more examples, other representations, etc.

DKF

WHOA! Donal, am I missing something here?
``` % tolist 75676576576576576567575
can't use floating-point value as operand of "%"
%
-PSE```

That's not an integer, is it? It's just a string of digits... Oh well, I've fixed the code below now, so your example should now work. - DKF

HOW is that not an integer?
```inýý·teýý·ger (nt-jr)
n. Mathematics

1.A member of the set of positive whole numbers (1, 2, 3, . . .),
negative whole numbers (-1, -2, -3, . . .),
and zero (0).
2.A complete unit or entity.```

BTW, I fixed your missing bracket problem down there... things are quite spiffy now. How would we ever get along without our own knight-of-the-lambda-calculus? Thanks Donal! -PSE

Hah! Dictionary definitions are for wimps! It cannot be parsed by Tcl_GetInt(), so obviously it cannot be an integer. This is confirmed by the use of the [string is integer 75676576576576576567575] test, which naturally fails.

(And cheers for fixing any missing bracket problems.)

DKF

All of the words in the dictionary cannot express my gratitude that I do not have to sit across the table from you at the Monday morning status meeting... -PSE

My status meetings are usually on Friday afternoons, and I have the advantage of working with people with heavy admin and/or teaching loads, so any work I choose to do is definitely appreciated. And I drink coffee that is substantially stronger than Starbucks makes theirs...

Now, I've got a compiler spec. that I really ought to get written for as soon as possible. So I do hope you'll forgive me if I leave this pleasant little discussion for now and get on with a little work. -DKF

So should I, but I like this place ;-) how 'bout
`proc string_is_dictionary_integer {s} {regexp {^-?[0-9]+\$} \$s} ;#RS`

```proc K {x y} {set x};# For a speedup

set rs {}
set carry 0
if {[llength \$xs]>[llength \$ys]} {
set len [llength \$ys]
set flag 0
} else {
set len [llength \$xs]
set flag 1
}
set pos 0
foreach x \$xs y \$ys {
if {\$pos>=\$len} {
if {\$flag} {set x 0} else {set y 0}
}
set n [expr {\$x+\$y+\$carry}]
set carry [expr {\$n / 10}]
lappend rs [expr {\$n % 10}]
incr pos ;# who stole this line? restored 11/00 MM
}
if {\$carry} {lappend rs \$carry}
return \$rs
}

# Multiply two numbers.
proc multa {xs ys} {
set rs 0
foreach x \$xs {
for {} {\$x>0} {incr x -1} {
set rs [adda [K \$rs [set rs {}]] \$ys]
}
set ys [linsert [K \$ys [set ys {}]] 0 0]
}
return \$rs
}

# Convert an integer to a (list-formatted) number.
proc tolist {n} {
set rs {}
if {[catch {incr n 0}]} {
foreach digit [split \$n {}] {
set rs [linsert [K \$rs [set rs {}]] 0 \$digit]
}
} else {
while {\$n > 0} {
lappend rs [expr {\$n%10}]
set n [expr {\$n/10}]
}
}
return \$rs
}

# Convert a (list-formatted) number to a human-readable form.
# Will be an integer (modulo type-shimmering) if representable.
proc fromlist {xs} {
set r {}
foreach x \$xs {
set r \$x\$r
}
if {![string length \$r]} {set r 0}
return \$r
}

# DEMO CODE: Factorial.  Try calculating [fact 40] for a 48-digit number
proc fact {n} {
set this 1
for {set i 1} {\$i <= \$n} {incr i} {
set this [multa [K \$this [set this {}]] [tolist \$i]]
}
return [fromlist \$this]
}

# DEMO CODE: Fibonacci Sequence.  Doesn't grow anything like as fast as factorial,
# but still grows beyond 32-bit integers quite easily.
proc fib {n} {
set this 1
set last 0
for {set i 1} {\$i<\$n} {incr i} {
set last \$this
set this \$new
}
return [fromlist \$this]
}```

See math::bignum from tcllib. It is pure-Tcl. And in Tcl 8.5, default integers are now Bignums as described in TIP #237.

See mkTulip at http://mkextensions.sourceforge.net/ .

Also perhaps checkout Decimal Arithmetic Package for tcl 8.5

Check out [Bignums Fast and Dirty ] and arbint and mpexpr.

[Maybe write here a bit about arithmetic based on Schönhage and Strassen, and also Priest's '91 paper on "Algorithms for Arbitrary Precision Floating Point Arithmetic". Is this the right place?]

• Also Cody's Software for the Elementary Function
• R. P. Brent. Fast multiple-precision evaluation of elementary functions. Journal of the Association for Computing Machinery, 23(2):242-251, April 1976.
• R.P. Brent. A FORTRAN multiple-precision arithmetic package ACM Trans. on Math. Software, 4, no 1, 57-70, 1978.
• Douglas M. Priest, "Algorithms for arbitrary precision floating point arithmetic," 10th IEEE Symposium on Computer Arithmetic, 1991, pp. 132-143
• D. M. Smith, "A FORTRAN package for floating-point multiple-precision arithmetic," ACM transactions on mathematical software 17 (2) (Jun. 1991) 273-283.
• David H. Bailey, "Multiprecision Translation and Execution of Fortran Programs", ACM Transactions on Mathematical Software, vol. 19, no. 3 (Sept. 1993), pg. 288--319
• Smith, D.M. 1998. Multiple Precision Complex Arithmetic and Functions. ACM Trans. Math. Softw. 24, 4 (December), 359--367
• http://www.netlib.org