]> Cypherpunks repositories - gostls13.git/commitdiff
- factored out 128-bit muladd and div into arith.go
authorRobert Griesemer <gri@golang.org>
Tue, 11 Aug 2009 00:29:55 +0000 (17:29 -0700)
committerRobert Griesemer <gri@golang.org>
Tue, 11 Aug 2009 00:29:55 +0000 (17:29 -0700)
- wrote corresponding fast versions in fast.arith.s
- implemented in-place operations for some routines
- updated existing code to be compatible with in-place
  routines

These changes allow the pidigits benchmark to run
approx. 30% faster. Enabling the assembly routines
in fast.arith.s will give another approx. 3%.

R=r
DELTA=486  (252 added, 68 deleted, 166 changed)
OCL=32980
CL=33003

src/pkg/bignum/Makefile
src/pkg/bignum/arith.go [new file with mode: 0644]
src/pkg/bignum/bignum.go
src/pkg/bignum/fast.arith.s [new file with mode: 0644]
src/pkg/bignum/integer.go
src/pkg/bignum/rational.go

index 12d1f3193d52a4bc1f92d8709b337bd14c9233c8..20f75bc279ed92a3fe864eeb09c167ed0a896fd1 100644 (file)
@@ -4,7 +4,7 @@
 
 
 # DO NOT EDIT.  Automatically generated by gobuild.
-# gobuild -m bignum.go integer.go rational.go >Makefile
+# gobuild -m arith.go bignum.go integer.go rational.go >Makefile
 
 D=
 
@@ -33,30 +33,37 @@ coverage: packages
        $(AS) $*.s
 
 O1=\
-       bignum.$O\
+       arith.$O\
 
 O2=\
-       integer.$O\
+       bignum.$O\
 
 O3=\
+       integer.$O\
+
+O4=\
        rational.$O\
 
 
-phases: a1 a2 a3
+phases: a1 a2 a3 a4
 _obj$D/bignum.a: phases
 
 a1: $(O1)
-       $(AR) grc _obj$D/bignum.a bignum.$O
+       $(AR) grc _obj$D/bignum.a arith.$O
        rm -f $(O1)
 
 a2: $(O2)
-       $(AR) grc _obj$D/bignum.a integer.$O
+       $(AR) grc _obj$D/bignum.a bignum.$O
        rm -f $(O2)
 
 a3: $(O3)
-       $(AR) grc _obj$D/bignum.a rational.$O
+       $(AR) grc _obj$D/bignum.a integer.$O
        rm -f $(O3)
 
+a4: $(O4)
+       $(AR) grc _obj$D/bignum.a rational.$O
+       rm -f $(O4)
+
 
 newpkg: clean
        mkdir -p _obj$D
@@ -66,6 +73,7 @@ $(O1): newpkg
 $(O2): a1
 $(O3): a2
 $(O4): a3
+$(O5): a4
 
 nuke: clean
        rm -f $(GOROOT)/pkg/$(GOOS)_$(GOARCH)$D/bignum.a
diff --git a/src/pkg/bignum/arith.go b/src/pkg/bignum/arith.go
new file mode 100644 (file)
index 0000000..4ae7e3e
--- /dev/null
@@ -0,0 +1,121 @@
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Fast versions of the routines in this file are in fast.arith.s.
+// Simply replace this file with arith.s (renamed from fast.arith.s)
+// and the bignum package will build and run on a platform that
+// supports the assembly routines.
+
+package bignum
+
+import "unsafe"
+
+// z1<<64 + z0 = x*y
+func Mul128(x, y uint64) (z1, z0 uint64) {
+       // Split x and y into 2 halfwords each, multiply
+       // the halfwords separately while avoiding overflow,
+       // and return the product as 2 words.
+
+       const (
+               W = uint(unsafe.Sizeof(x))*8;
+               W2 = W/2;
+               B2 = 1<<W2;
+               M2 = B2-1;
+       )
+
+       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)
+               // z = z[0] = x*y
+               z0 = x*y;
+               return;
+       }
+
+       if y < B2 {
+               // sub-digits of x and y are (x1, x0) and (0, y)
+               // x = (x1*B2 + x0)
+               // y = (y1*B2 + y0)
+               x1, x0 := x>>W2, x&M2;
+       
+               // x*y = t2*B2*B2 + t1*B2 + t0
+               t0 := x0*y;
+               t1 := x1*y;
+
+               // compute result digits but avoid overflow
+               // z = z[1]*B + z[0] = x*y
+               z0 = t1<<W2 + t0;
+               z1 = (t1 + t0>>W2) >> 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;
+       y1, y0 := y>>W2, y&M2;
+
+       // x*y = t2*B2*B2 + t1*B2 + t0
+       t0 := x0*y0;
+       t1 := x1*y0 + x0*y1;
+       t2 := x1*y1;
+
+       // compute result digits but avoid overflow
+       // z = z[1]*B + z[0] = x*y
+       z0 = t1<<W2 + t0;
+       z1 = t2 + (t1 + t0>>W2) >> W2;
+       return;
+}
+
+
+// z1<<64 + z0 = x*y + c
+func MulAdd128(x, y, c uint64) (z1, z0 uint64) {
+       // Split x and y into 2 halfwords each, multiply
+       // the halfwords separately while avoiding overflow,
+       // and return the product as 2 words.
+
+       const (
+               W = uint(unsafe.Sizeof(x))*8;
+               W2 = W/2;
+               B2 = 1<<W2;
+               M2 = B2-1;
+       )
+
+       // TODO(gri) Should implement special cases for faster execution.
+
+       // general case
+       // sub-digits of x, y, and c are (x1, x0), (y1, y0), (c1, c0)
+       // x = (x1*B2 + x0)
+       // y = (y1*B2 + y0)
+       x1, x0 := x>>W2, x&M2;
+       y1, y0 := y>>W2, y&M2;
+       c1, c0 := c>>W2, c&M2;
+
+       // x*y + c = t2*B2*B2 + t1*B2 + t0
+       t0 := x0*y0 + c0;
+       t1 := x1*y0 + x0*y1 + c1;
+       t2 := x1*y1;
+
+       // compute result digits but avoid overflow
+       // z = z[1]*B + z[0] = x*y
+       z0 = t1<<W2 + t0;
+       z1 = t2 + (t1 + t0>>W2) >> W2;
+       return;
+}
+
+
+// q = (x1<<64 + x0)/y + r
+func Div128(x1, x0, y uint64) (q, r uint64) {
+       if x1 == 0 {
+               q, r = x0/y, x0%y;
+               return;
+       }
+
+       // TODO(gri) Implement general case.
+       panic("Div128 not implemented for x > 1<<64-1");
+}
index 63f385fca4b6a72a56adc4129fc764a833c07813..2b36fbd75b2d7d4b0e5e07a2a3db93225c412595 100755 (executable)
 //
 package bignum
 
-import "fmt"
+import (
+       "bignum";
+       "fmt";
+)
 
+// TODO(gri) Complete the set of in-place operations.
 
 // ----------------------------------------------------------------------------
 // Internal representation
@@ -59,7 +63,7 @@ type (
 
 
 const (
-       logW = 64;
+       logW = 64;  // word width
        logH = 4;  // bits for a hex digit (= small number)
        logB = logW - logH;  // largest bit-width available
 
@@ -92,7 +96,7 @@ func isSmall(x digit) bool {
 
 // For debugging. Keep around.
 /*
-func dump(x []digit) {
+func dump(x Natural) {
        print("[", len(x), "]");
        for i := len(x) - 1; i >= 0; i-- {
                print(" ", x[i]);
@@ -110,26 +114,16 @@ func dump(x []digit) {
 type Natural []digit;
 
 
-// Common small values - allocate once.
-var nat [16]Natural;
-
-func init() {
-       nat[0] = Natural{};  // zero has no digits
-       for i := 1; i < len(nat); i++ {
-               nat[i] = Natural{digit(i)};
-       }
-}
-
-
 // Nat creates a small natural number with value x.
 //
 func Nat(x uint64) Natural {
-       // avoid allocation for common small values
-       if x < uint64(len(nat)) {
-               return nat[x];
+       if x == 0 {
+               return nil;  // len == 0
        }
 
        // single-digit values
+       // (note: cannot re-use preallocated values because
+       //        the in-place operations may overwrite them)
        if x < _B {
                return Natural{digit(x)};
        }
@@ -211,25 +205,40 @@ func (x Natural) IsZero() bool {
 
 func normalize(x Natural) Natural {
        n := len(x);
-       for n > 0 && x[n - 1] == 0 { n-- }
-       if n < len(x) {
-               x = x[0 : n];  // trim leading 0's
+       for n > 0 && x[n-1] == 0 { n-- }
+       return x[0 : n];
+}
+
+
+// nalloc returns a Natural of n digits. If z is large
+// enough, z is resized and returned. Otherwise, a new
+// Natural is allocated.
+//
+func nalloc(z Natural, n int) Natural {
+       size := n;
+       if size <= 0 {
+               size = 4;
+       }
+       if size <= cap(z) {
+               return z[0 : n];
        }
-       return x;
+       return make(Natural, n, size);
 }
 
 
-// Add returns the sum x + y.
+// Nadd sets *zp to the sum x + y.
+// *zp may be x or y.
 //
-func (x Natural) Add(y Natural) Natural {
+func Nadd(zp *Natural, x, y Natural) {
        n := len(x);
        m := len(y);
        if n < m {
-               return y.Add(x);
+               Nadd(zp, y, x);
+               return;
        }
 
+       z := nalloc(*zp, n+1);
        c := digit(0);
-       z := make(Natural, n + 1);
        i := 0;
        for i < m {
                t := c + x[i] + y[i];
@@ -245,23 +254,32 @@ func (x Natural) Add(y Natural) Natural {
                z[i] = c;
                i++;
        }
+       *zp = z[0 : i]
+}
+
 
-       return z[0 : i];
+// Add returns the sum z = x + y.
+//
+func (x Natural) Add(y Natural) Natural {
+       var z Natural;
+       Nadd(&z, x, y);
+       return z;
 }
 
 
-// Sub returns the difference x - y for x >= y.
+// Nsub sets *zp to the difference x - y for x >= y.
 // If x < y, an underflow run-time error occurs (use Cmp to test if x >= y).
+// *zp may be x or y.
 //
-func (x Natural) Sub(y Natural) Natural {
+func Nsub(zp *Natural, x, y Natural) {
        n := len(x);
        m := len(y);
        if n < m {
                panic("underflow")
        }
 
+       z := nalloc(*zp, n);
        c := digit(0);
-       z := make(Natural, n);
        i := 0;
        for i < m {
                t := c + x[i] - y[i];
@@ -273,110 +291,104 @@ func (x Natural) Sub(y Natural) Natural {
                c, z[i] = digit(int64(t)>>_W), t&_M;  // requires arithmetic shift!
                i++;
        }
-       for i > 0 && z[i - 1] == 0 {  // normalize
-               i--;
+       if int64(c) < 0 {
+               panic("underflow");
        }
+       *zp = normalize(z);
+}
+
 
-       return z[0 : i];
+// Sub returns the difference x - y for x >= y.
+// If x < y, an underflow run-time error occurs (use Cmp to test if x >= y).
+//
+func (x Natural) Sub(y Natural) Natural {
+       var z Natural;
+       Nsub(&z, x, y);
+       return z;
 }
 
 
-// Returns c = x*y div B, z = x*y mod B.
+// MulAdd128 is defined in arith.go and arith.s .
+func MulAdd128(x, y, c uint64) (z1, z0 uint64)
+
+// Returns z1 = (x*y + c) div B, z0 = (x*y + c) mod B.
 //
-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.
+func muladd11(x, y, c digit) (digit, digit) {
+       z1, z0 := MulAdd128(uint64(x), uint64(y), uint64(c));
+       return digit(z1<<(64 - logB) | z0>>logB), digit(z0&_M);
+}
 
-       // This code also works for non-even bit widths W
-       // which is why there are separate constants below
-       // for half-digits.
-       const W2 = (_W + 1)/2;
-       const DW = W2*2 - _W;  // 0 or 1
-       const B2  = 1<<W2;
-       const M2  = _B2 - 1;
 
-       if x < y {
-               x, y = y, x;
+func mul1(z, x Natural, y digit) (c digit) {
+       for i := 0; i < len(x); i++ {
+               c, z[i] = muladd11(x[i], y, c);
        }
+       return;
+}
 
-       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);
+// Nscale sets *z to the scaled value (*z) * d.
+//
+func Nscale(z *Natural, d uint64) {
+       switch {
+       case d == 0:
+               *z = Nat(0);
+               return
+       case d == 1:
+               return;
+       case d >= _B:
+               *z = z.Mul1(d);
                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;
+       c := mul1(*z, *z, digit(d));
 
-               // compute result digits but avoid overflow
-               // z = z1*B + z0 = x*y
-               z0 = (t1<<W2 + t0)&_M;
-               z1 = (t1 + t0>>W2) >> (_W-W2);
-               return;
+       if c != 0 {
+               n := len(*z);
+               if n >= cap(*z) {
+                       zz := make(Natural, n+1);
+                       for i, d := range *z {
+                               zz[i] = d;
+                       }
+                       *z = zz
+               } else {
+                       *z = (*z)[0 : n+1];
+               }
+               (*z)[n] = c;
        }
+}
 
-       // 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;
-       y1, y0 := y>>W2, y&M2;
 
-       // x*y = t2*B2^2 + t1*B2 + t0
-       t0 := x0*y0;
-       t1 := x1*y0 + x0*y1;
-       t2 := x1*y1;
+// Computes x = x*d + c for small d's.
+//
+func muladd1(x Natural, d, c digit) Natural {
+       assert(isSmall(d-1) && isSmall(c));
+       n := len(x);
+       z := make(Natural, n + 1);
 
-       // 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);
-       return;
-}
+       for i := 0; i < n; i++ {
+               t := c + x[i]*d;
+               c, z[i] = t>>_W, t&_M;
+       }
+       z[n] = c;
 
+       return normalize(z);
+}
 
-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 == 0: return Nat(0);
        case d == 1: return x;
+       case isSmall(digit(d)): muladd1(x, digit(d), 0);
        case d >= _B: return x.Mul(Nat(d));
        }
 
-       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;
-       }
-
+       z := make(Natural, len(x) + 1);
+       c := mul1(z, x, digit(d));
+       z[len(x)] = c;
        return normalize(z);
 }
 
@@ -390,6 +402,10 @@ func (x Natural) Mul(y Natural) Natural {
                return y.Mul(x);
        }
 
+       if m == 0 {
+               return Nat(0);
+       }
+
        if m == 1 && y[0] < _B {
                return x.Mul1(uint64(y[0]));
        }
@@ -400,11 +416,7 @@ func (x Natural) Mul(y Natural) Natural {
                if d != 0 {
                        c := digit(0);
                        for i := 0; i < n; i++ {
-                               // z[i+j] += c + x[i]*d;
-                               z1, z0 := mul11(x[i], d);
-                               t := c + z[i+j] + z0;
-                               c, z[i+j] = t>>_W, t&_M;
-                               c += z1;
+                               c, z[i+j] = muladd11(x[i], d, z[i+j] + c);
                        }
                        z[n+j] = c;
                }
@@ -450,11 +462,10 @@ func pack(x []digit2) Natural {
 }
 
 
-func mul1(z, x []digit2, y digit2) digit2 {
-       n := len(x);
+func mul21(z, x []digit2, y digit2) digit2 {
        c := digit(0);
        f := digit(y);
-       for i := 0; i < n; i++ {
+       for i := 0; i < len(x); i++ {
                t := c + digit(x[i])*f;
                c, z[i] = t>>_W2, digit2(t&_M2);
        }
@@ -462,12 +473,11 @@ func mul1(z, x []digit2, y digit2) digit2 {
 }
 
 
-func div1(z, x []digit2, y digit2) digit2 {
-       n := len(x);
+func div21(z, x []digit2, y digit2) digit2 {
        c := digit(0);
        d := digit(y);
-       for i := n-1; i >= 0; i-- {
-               t := c*_B2 + digit(x[i]);
+       for i := len(x)-1; i >= 0; i-- {
+               t := c<<_W2 + digit(x[i]);
                c, z[i] = t%d, digit2(t/d);
        }
        return digit2(c);
@@ -501,13 +511,13 @@ func divmod(x, y []digit2) ([]digit2, []digit2) {
                panic("division by zero");
        }
        assert(n+1 <= cap(x));  // space for one extra digit
-       x = x[0 : n + 1];
+       x = x[0 : n+1];
        assert(x[n] == 0);
 
        if m == 1 {
                // division by single digit
                // result is shifted left by 1 in place!
-               x[0] = div1(x[1 : n+1], x[0 : n], y[0]);
+               x[0] = div21(x[1 : n+1], x[0 : n], y[0]);
 
        } else if m > n {
                // y > x => quotient = 0, remainder = x
@@ -524,13 +534,12 @@ func divmod(x, y []digit2) ([]digit2, []digit2) {
                //      satisfied (as done in Hacker's Delight).
                f := _B2 / (digit(y[m-1]) + 1);
                if f != 1 {
-                       mul1(x, x, digit2(f));
-                       mul1(y, y, digit2(f));
+                       mul21(x, x, digit2(f));
+                       mul21(y, y, digit2(f));
                }
                assert(_B2/2 <= y[m-1] && y[m-1] < _B2);  // incorrect scaling
 
                y1, y2 := digit(y[m-1]), digit(y[m-2]);
-               d2 := digit(y1)<<_W2 + digit(y2);
                for i := n-m; i >= 0; i-- {
                        k := i+m;
 
@@ -540,7 +549,7 @@ func divmod(x, y []digit2) ([]digit2, []digit2) {
                                if x0 != y1 {
                                        q = (x0<<_W2 + x1)/y1;
                                } else {
-                                       q = _B2 - 1;
+                                       q = _B2-1;
                                }
                                for y2*q > (x0<<_W2 + x1 - y1*q)<<_W2 + x2 {
                                        q--
@@ -572,7 +581,7 @@ func divmod(x, y []digit2) ([]digit2, []digit2) {
 
                // undo normalization for remainder
                if f != 1 {
-                       c := div1(x[0 : m], x[0 : m], digit2(f));
+                       c := div21(x[0 : m], x[0 : m], digit2(f));
                        assert(c == 0);
                }
        }
@@ -610,7 +619,7 @@ func (x Natural) DivMod(y Natural) (Natural, Natural) {
 }
 
 
-func shl(z, x []digit, s uint) digit {
+func shl(z, x Natural, s uint) digit {
        assert(s <= _W);
        n := len(x);
        c := digit(0);
@@ -634,7 +643,7 @@ func (x Natural) Shl(s uint) Natural {
 }
 
 
-func shr(z, x []digit, s uint) digit {
+func shr(z, x Natural, s uint) digit {
        assert(s <= _W);
        n := len(x);
        c := digit(0);
@@ -680,7 +689,7 @@ func (x Natural) And(y Natural) Natural {
 }
 
 
-func copy(z, x []digit) {
+func copy(z, x Natural) {
        for i, e := range x {
                z[i] = e
        }
@@ -881,23 +890,6 @@ func hexvalue(ch byte) uint {
 }
 
 
-// Computes x = x*d + c for small d's.
-//
-func muladd1(x Natural, d, c digit) Natural {
-       assert(isSmall(d-1) && isSmall(c));
-       n := len(x);
-       z := make(Natural, n + 1);
-
-       for i := 0; i < n; i++ {
-               t := c + x[i]*d;
-               c, z[i] = t>>_W, t&_M;
-       }
-       z[n] = c;
-
-       return normalize(z);
-}
-
-
 // NatFromString returns the natural number corresponding to the
 // longest possible prefix of s representing a natural number in a
 // given conversion base, the actual conversion base used, and the
@@ -924,7 +916,7 @@ func NatFromString(s string, base uint) (Natural, uint, int) {
 
        // convert string
        assert(2 <= base && base <= 16);
-       x := nat[0];
+       x := Nat(0);
        for ; i < n; i++ {
                d := hexvalue(s[i]);
                if d < base {
@@ -964,7 +956,7 @@ func (x Natural) Pop() uint {
 // Pow computes x to the power of n.
 //
 func (xp Natural) Pow(n uint) Natural {
-       z := nat[1];
+       z := Nat(1);
        x := xp;
        for n > 0 {
                // z * x^n == x^n0
@@ -982,7 +974,7 @@ func (xp Natural) Pow(n uint) Natural {
 //
 func MulRange(a, b uint) Natural {
        switch {
-       case a > b: return nat[1];
+       case a > b: return Nat(1);
        case a == b: return Nat(uint64(a));
        case a + 1 == b: return Nat(uint64(a)).Mul(Nat(uint64(b)));
        }
diff --git a/src/pkg/bignum/fast.arith.s b/src/pkg/bignum/fast.arith.s
new file mode 100644 (file)
index 0000000..4474e38
--- /dev/null
@@ -0,0 +1,41 @@
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// This file provides fast assembly versions
+// of the routines in arith.go.
+
+// func Mul128(x, y uint64) (z1, z0 uint64)
+// z1<<64 + z0 = x*y
+//
+TEXT bignum·Mul128(SB),7,$0
+       MOVQ a+0(FP), AX
+       MULQ a+8(FP)
+       MOVQ DX, a+16(FP)
+       MOVQ AX, a+24(FP)
+       RET
+
+
+// func MulAdd128(x, y, c uint64) (z1, z0 uint64)
+// z1<<64 + z0 = x*y + c
+//
+TEXT bignum·MulAdd128(SB),7,$0
+       MOVQ a+0(FP), AX
+       MULQ a+8(FP)
+       ADDQ a+16(FP), AX
+       ADCQ $0, DX
+       MOVQ DX, a+24(FP)
+       MOVQ AX, a+32(FP)
+       RET
+
+
+// func Div128(x1, x0, y uint64) (q, r uint64)
+// q = (x1<<64 + x0)/y + r
+//
+TEXT bignum·Div128(SB),7,$0
+       MOVQ a+0(FP), DX
+       MOVQ a+8(FP), AX
+       DIVQ a+16(FP)
+       MOVQ AX, a+24(FP)
+       MOVQ DX, a+32(FP)
+       RET
index bb1429aee13de2fa571b8405dfb4c4f2ea8d7932..30e13092fac6e7993e3d03e9f3d7ddf47723f834 100644 (file)
@@ -9,9 +9,12 @@
 
 package bignum
 
-import "bignum"
-import "fmt"
+import (
+       "bignum";
+       "fmt";
+)
 
+// TODO(gri) Complete the set of in-place operations.
 
 // Integer represents a signed integer value of arbitrary precision.
 //
@@ -113,61 +116,82 @@ func (x *Integer) Neg() *Integer {
 }
 
 
-// Add returns the sum x + y.
+// Iadd sets z to the sum x + y.
+// z must exist and may be x or y.
 //
-func (x *Integer) Add(y *Integer) *Integer {
-       var z *Integer;
+func Iadd(z, x, y *Integer) {
        if x.sign == y.sign {
                // x + y == x + y
                // (-x) + (-y) == -(x + y)
-               z = MakeInt(x.sign, x.mant.Add(y.mant));
+               z.sign = x.sign;
+               Nadd(&z.mant, x.mant, y.mant);
        } else {
                // x + (-y) == x - y == -(y - x)
                // (-x) + y == y - x == -(x - y)
                if x.mant.Cmp(y.mant) >= 0 {
-                       z = MakeInt(false, x.mant.Sub(y.mant));
+                       z.sign = x.sign;
+                       Nsub(&z.mant, x.mant, y.mant);
                } else {
-                       z = MakeInt(true, y.mant.Sub(x.mant));
+                       z.sign = !x.sign;
+                       Nsub(&z.mant, y.mant, x.mant);
                }
        }
-       if x.sign {
-               z.sign = !z.sign;
-       }
-       return z;
 }
 
 
-// Sub returns the difference x - y.
+// Add returns the sum x + y.
 //
-func (x *Integer) Sub(y *Integer) *Integer {
-       var z *Integer;
+func (x *Integer) Add(y *Integer) *Integer {
+       var z Integer;
+       Iadd(&z, x, y);
+       return &z;
+}
+
+
+func Isub(z, x, y *Integer) {
        if x.sign != y.sign {
                // x - (-y) == x + y
                // (-x) - y == -(x + y)
-               z = MakeInt(false, x.mant.Add(y.mant));
+               z.sign = x.sign;
+               Nadd(&z.mant, x.mant, y.mant);
        } else {
                // x - y == x - y == -(y - x)
                // (-x) - (-y) == y - x == -(x - y)
                if x.mant.Cmp(y.mant) >= 0 {
-                       z = MakeInt(false, x.mant.Sub(y.mant));
+                       z.sign = x.sign;
+                       Nsub(&z.mant, x.mant, y.mant);
                } else {
-                       z = MakeInt(true, y.mant.Sub(x.mant));
+                       z.sign = !x.sign;
+                       Nsub(&z.mant, y.mant, x.mant);
                }
        }
-       if x.sign {
-               z.sign = !z.sign;
+}
+
+
+// Sub returns the difference x - y.
+//
+func (x *Integer) Sub(y *Integer) *Integer {
+       var z Integer;
+       Isub(&z, x, y);
+       return &z;
+}
+
+
+// Nscale sets *z to the scaled value (*z) * d.
+//
+func Iscale(z *Integer, d int64) {
+       f := uint64(d);
+       if d < 0 {
+               f = uint64(-d);
        }
-       return z;
+       z.sign = z.sign != (d < 0);
+       Nscale(&z.mant, f);
 }
 
 
 // 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);
@@ -326,7 +350,7 @@ func (x *Integer) Shl(s uint) *Integer {
 func (x *Integer) Shr(s uint) *Integer {
        if x.sign {
                // (-x) >> s == ^(x-1) >> s == ^((x-1) >> s) == -(((x-1) >> s) + 1)
-               return MakeInt(true, x.mant.Sub(nat[1]).Shr(s).Add(nat[1]));
+               return MakeInt(true, x.mant.Sub(Nat(1)).Shr(s).Add(Nat(1)));
        }
        
        return MakeInt(false, x.mant.Shr(s));
@@ -337,11 +361,11 @@ func (x *Integer) Shr(s uint) *Integer {
 func (x *Integer) Not() *Integer {
        if x.sign {
                // ^(-x) == ^(^(x-1)) == x-1
-               return MakeInt(false, x.mant.Sub(nat[1]));
+               return MakeInt(false, x.mant.Sub(Nat(1)));
        }
 
        // ^x == -x-1 == -(x+1)
-       return MakeInt(true, x.mant.Add(nat[1]));
+       return MakeInt(true, x.mant.Add(Nat(1)));
 }
 
 
@@ -351,7 +375,7 @@ func (x *Integer) And(y *Integer) *Integer {
        if x.sign == y.sign {
                if x.sign {
                        // (-x) & (-y) == ^(x-1) & ^(y-1) == ^((x-1) | (y-1)) == -(((x-1) | (y-1)) + 1)
-                       return MakeInt(true, x.mant.Sub(nat[1]).Or(y.mant.Sub(nat[1])).Add(nat[1]));
+                       return MakeInt(true, x.mant.Sub(Nat(1)).Or(y.mant.Sub(Nat(1))).Add(Nat(1)));
                }
 
                // x & y == x & y
@@ -364,7 +388,7 @@ func (x *Integer) And(y *Integer) *Integer {
        }
 
        // x & (-y) == x & ^(y-1) == x &^ (y-1)
-       return MakeInt(false, x.mant.AndNot(y.mant.Sub(nat[1])));
+       return MakeInt(false, x.mant.AndNot(y.mant.Sub(Nat(1))));
 }
 
 
@@ -374,7 +398,7 @@ func (x *Integer) AndNot(y *Integer) *Integer {
        if x.sign == y.sign {
                if x.sign {
                        // (-x) &^ (-y) == ^(x-1) &^ ^(y-1) == ^(x-1) & (y-1) == (y-1) &^ (x-1)
-                       return MakeInt(false, y.mant.Sub(nat[1]).AndNot(x.mant.Sub(nat[1])));
+                       return MakeInt(false, y.mant.Sub(Nat(1)).AndNot(x.mant.Sub(Nat(1))));
                }
 
                // x &^ y == x &^ y
@@ -383,11 +407,11 @@ func (x *Integer) AndNot(y *Integer) *Integer {
 
        if x.sign {
                // (-x) &^ y == ^(x-1) &^ y == ^(x-1) & ^y == ^((x-1) | y) == -(((x-1) | y) + 1)
-               return MakeInt(true, x.mant.Sub(nat[1]).Or(y.mant).Add(nat[1]));
+               return MakeInt(true, x.mant.Sub(Nat(1)).Or(y.mant).Add(Nat(1)));
        }
 
        // x &^ (-y) == x &^ ^(y-1) == x & (y-1)
-       return MakeInt(false, x.mant.And(y.mant.Sub(nat[1])));
+       return MakeInt(false, x.mant.And(y.mant.Sub(Nat(1))));
 }
 
 
@@ -397,7 +421,7 @@ func (x *Integer) Or(y *Integer) *Integer {
        if x.sign == y.sign {
                if x.sign {
                        // (-x) | (-y) == ^(x-1) | ^(y-1) == ^((x-1) & (y-1)) == -(((x-1) & (y-1)) + 1)
-                       return MakeInt(true, x.mant.Sub(nat[1]).And(y.mant.Sub(nat[1])).Add(nat[1]));
+                       return MakeInt(true, x.mant.Sub(Nat(1)).And(y.mant.Sub(Nat(1))).Add(Nat(1)));
                }
 
                // x | y == x | y
@@ -410,7 +434,7 @@ func (x *Integer) Or(y *Integer) *Integer {
        }
 
        // x | (-y) == x | ^(y-1) == ^((y-1) &^ x) == -(^((y-1) &^ x) + 1)
-       return MakeInt(true, y.mant.Sub(nat[1]).AndNot(x.mant).Add(nat[1]));
+       return MakeInt(true, y.mant.Sub(Nat(1)).AndNot(x.mant).Add(Nat(1)));
 }
 
 
@@ -420,7 +444,7 @@ func (x *Integer) Xor(y *Integer) *Integer {
        if x.sign == y.sign {
                if x.sign {
                        // (-x) ^ (-y) == ^(x-1) ^ ^(y-1) == (x-1) ^ (y-1)
-                       return MakeInt(false, x.mant.Sub(nat[1]).Xor(y.mant.Sub(nat[1])));
+                       return MakeInt(false, x.mant.Sub(Nat(1)).Xor(y.mant.Sub(Nat(1))));
                }
 
                // x ^ y == x ^ y
@@ -433,7 +457,7 @@ func (x *Integer) Xor(y *Integer) *Integer {
        }
 
        // x ^ (-y) == x ^ ^(y-1) == ^(x ^ (y-1)) == -((x ^ (y-1)) + 1)
-       return MakeInt(true, x.mant.Xor(y.mant.Sub(nat[1])).Add(nat[1]));
+       return MakeInt(true, x.mant.Xor(y.mant.Sub(Nat(1))).Add(Nat(1)));
 }
 
 
index c00022b296b8094785f8130a3f43dc96b3abb7c7..92b5d8883bfdc426acf16f132017baffa7e4e044 100644 (file)
@@ -22,7 +22,7 @@ type Rational struct {
 //
 func MakeRat(a *Integer, b Natural) *Rational {
        f := a.mant.Gcd(b);  // f > 0
-       if f.Cmp(nat[1]) != 0 {
+       if f.Cmp(Nat(1)) != 0 {
                a = MakeInt(a.sign, a.mant.Div(f));
                b = b.Div(f);
        }
@@ -75,7 +75,7 @@ func (x *Rational) IsPos() bool {
 // in the form x == x'/1; i.e., if x is an integer value.
 //
 func (x *Rational) IsInt() bool {
-       return x.b.Cmp(nat[1]) == 0;
+       return x.b.Cmp(Nat(1)) == 0;
 }
 
 
@@ -184,7 +184,7 @@ func (x *Rational) Format(h fmt.State, c int) {
 func RatFromString(s string, base uint) (*Rational, uint, int) {
        // read numerator
        a, abase, alen := IntFromString(s, base);
-       b := nat[1];
+       b := Nat(1);
 
        // read denominator or fraction, if any
        var blen int;
@@ -211,7 +211,7 @@ func RatFromString(s string, base uint) (*Rational, uint, int) {
                        rlen++;
                        e, _, elen := IntFromString(s[rlen : len(s)], 10);
                        rlen += elen;
-                       m := nat[10].Pow(uint(e.mant.Value()));
+                       m := Nat(10).Pow(uint(e.mant.Value()));
                        if e.sign {
                                b = b.Mul(m);
                        } else {