From 081bc69cde2d1919e1844f46998c6d545ca9921c Mon Sep 17 00:00:00 2001 From: Robert Griesemer <gri@golang.org> Date: Thu, 6 Aug 2009 18:16:51 -0700 Subject: [PATCH] - initial version of pidigits.go benchmark - extra bignum.go functionality for pidigits.go - tuned bignum multiplication R=r DELTA=193 (186 added, 0 deleted, 7 changed) OCL=32852 CL=32856 --- src/pkg/bignum/bignum.go | 84 +++++++++++++++++++++++-- src/pkg/bignum/integer.go | 15 +++++ test/bench/pidigits.go | 128 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 221 insertions(+), 6 deletions(-) create mode 100644 test/bench/pidigits.go diff --git a/src/pkg/bignum/bignum.go b/src/pkg/bignum/bignum.go index dd2200b6aa..63f385fca4 100755 --- a/src/pkg/bignum/bignum.go +++ b/src/pkg/bignum/bignum.go @@ -283,7 +283,7 @@ func (x Natural) Sub(y Natural) Natural { // Returns c = x*y div B, z = x*y mod B. // -func mul11(x, y digit) (digit, digit) { +func mul11(x, y digit) (z1, z0 digit) { // Split x and y into 2 sub-digits each, // multiply the digits separately while avoiding overflow, // and return the product as two separate digits. @@ -296,7 +296,44 @@ func mul11(x, y digit) (digit, digit) { const B2 = 1<<W2; const M2 = _B2 - 1; - // split x and y into sub-digits + if x < y { + x, y = y, x; + } + + if x < _B2 { + // y < _B2 because y <= x + // sub-digits of x and y are (0, x) and (0, y) + // x = x + // y = y + t0 := x*y; + + // compute result digits but avoid overflow + // z = z1*B + z0 = x*y + z0 = t0 & _M; + z1 = (t0>>W2) >> (_W-W2); + return; + } + + if y < _B2 { + // split x and y into sub-digits + // sub-digits of y are (x1, x0) and (0, y) + // x = (x1*B2 + x0) + // y = y + x1, x0 := x>>W2, x&M2; + + // x*y = t1*B2 + t0 + t0 := x0*y; + t1 := x1*y; + + // compute result digits but avoid overflow + // z = z1*B + z0 = x*y + z0 = (t1<<W2 + t0)&_M; + z1 = (t1 + t0>>W2) >> (_W-W2); + return; + } + + // general case + // sub-digits of x and y are (x1, x0) and (y1, y0) // x = (x1*B2 + x0) // y = (y1*B2 + y0) x1, x0 := x>>W2, x&M2; @@ -307,12 +344,40 @@ func mul11(x, y digit) (digit, digit) { t1 := x1*y0 + x0*y1; t2 := x1*y1; - // compute the result digits but avoid overflow + // compute result digits but avoid overflow // z = z1*B + z0 = x*y - z0 := (t1<<W2 + t0)&_M; - z1 := t2<<DW + (t1 + t0>>W2)>>(_W-W2); + z0 = (t1<<W2 + t0)&_M; + z1 = t2<<DW + (t1 + t0>>W2) >> (_W-W2); + return; +} + + +func (x Natural) Mul(y Natural) Natural + +// Mul1 returns the product x * d. +// +func (x Natural) Mul1(d uint64) Natural { + switch { + case d == 0: return nat[0]; + case d == 1: return x; + case d >= _B: return x.Mul(Nat(d)); + } - return z1, z0; + n := len(x); + z := make(Natural, n + 1); + if d != 0 { + c := digit(0); + for i := 0; i < n; i++ { + // z[i] += c + x[i]*d; + z1, z0 := mul11(x[i], digit(d)); + t := c + z[i] + z0; + c, z[i] = t>>_W, t&_M; + c += z1; + } + z[n] = c; + } + + return normalize(z); } @@ -321,6 +386,13 @@ func mul11(x, y digit) (digit, digit) { func (x Natural) Mul(y Natural) Natural { n := len(x); m := len(y); + if n < m { + return y.Mul(x); + } + + if m == 1 && y[0] < _B { + return x.Mul1(uint64(y[0])); + } z := make(Natural, n + m); for j := 0; j < m; j++ { diff --git a/src/pkg/bignum/integer.go b/src/pkg/bignum/integer.go index f3c111d0b4..bb1429aee1 100644 --- a/src/pkg/bignum/integer.go +++ b/src/pkg/bignum/integer.go @@ -161,6 +161,21 @@ func (x *Integer) Sub(y *Integer) *Integer { } +// Mul1 returns the product x * d. +// +func (x *Integer) Mul1(d int64) *Integer { + // x * y == x * y + // x * (-y) == -(x * y) + // (-x) * y == -(x * y) + // (-x) * (-y) == x * y + f := uint64(d); + if d < 0 { + f = uint64(-d); + } + return MakeInt(x.sign != (d < 0), x.mant.Mul1(f)); +} + + // Mul returns the product x * y. // func (x *Integer) Mul(y *Integer) *Integer { diff --git a/test/bench/pidigits.go b/test/bench/pidigits.go new file mode 100644 index 0000000000..6e1e7e0534 --- /dev/null +++ b/test/bench/pidigits.go @@ -0,0 +1,128 @@ +/* +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + * Neither the name of "The Computer Language Benchmarks Game" nor the + name of "The Computer Language Shootout Benchmarks" nor the names of + its contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. +*/ + +/* The Computer Language Benchmarks Game + * http://shootout.alioth.debian.org/ + * + * contributed by The Go Authors. + * based on pidigits.c (by Paolo Bonzini & Sean Bartlett, + * modified by Michael Mellor) + */ + +package main + +import ( + "bignum"; + "flag"; + "fmt"; +) + +var n = flag.Int("n", 27, "number of digits"); + +// TODO for easier profiling, remove eventually +var silent = flag.Bool("s", false, "don't print result"); + +var ( + tmp1 *bignum.Integer; + tmp2 *bignum.Integer; + numer = bignum.Int(1); + accum = bignum.Int(0); + denom = bignum.Int(1); +) + +func extract_digit() int64 { + if numer.Cmp(accum) > 0 { + return -1; + } + + /* Compute (numer * 3 + accum) / denom */ + tmp1, tmp2 = numer.MulNat(bignum.Nat(3)).Add(accum).QuoRem(denom); + + /* Now, if (numer * 4 + accum) % denom... */ + tmp2 = tmp2.Add(numer); + + /* ... is normalized, then the two divisions have the same result. */ + if tmp2.Cmp(denom) >= 0 { + return -1; + } + + return tmp1.Value(); +} + +func next_term(k int64) { + y2 := k*2 + 1; + + tmp1 = numer.Shl(1); + accum = accum.Add(tmp1); + accum = accum.Mul1(y2); + numer = numer.Mul1(k); + denom = denom.Mul1(y2); +} + +func eliminate_digit(d int64) { + accum = accum.Sub(denom.Mul1(d)); + accum = accum.Mul1(10); + numer = numer.Mul1(10); +} + +func printf(s string, arg ...) { + if !*silent { + fmt.Printf(s, arg); + } +} + +func main() { + flag.Parse(); + + var m int; // 0 <= m < 10 + for i, k := 0, int64(0); ; { + d := int64(-1); + for d < 0 { + k++; + next_term(k); + d = extract_digit(); + } + + printf("%c", d + '0'); + + i++; + m = i%10; + if m == 0 { + printf("\t:%d\n", i); + } + if i >= *n { + break; + } + eliminate_digit(d); + } + + if m > 0 { + printf("%s\t:%d\n", " "[m : 10], *n); + } +} -- 2.30.9