// x = x[n-1]*B^(n-1) + x[n-2]*B^(n-2) + ... + x[1]*B + x[0]
//
// with 0 <= x[i] < B and 0 <= i < n is stored in an array of length n,
-// with the digits x[i] as the array elements. 0 is represented as an
-// empty array (length == 0).
+// with the digits x[i] as the array elements.
//
// A natural number is normalized if the array contains no leading 0 digits.
// During arithmetic operations, denormalized values may occur which are
-// always normalized before returning the final result.
+// always normalized before returning the final result. The normalized
+// representation of 0 is the empty array (length = 0).
//
// The base B is chosen as large as possible on a given platform but there
-// are a few constraints besides the largest unsigned integer type available.
+// are a few constraints besides the size of the largest unsigned integer
+// type available.
// TODO describe the constraints.
-type Word uint64;
const LogW = 64;
-
const LogH = 4; // bits for a hex digit (= "small" number)
-const H = 1 << LogH;
-
const LogB = LogW - LogH;
-const L = LogB;
-const B = 1 << LogB;
-const M = B - 1;
-// For division
-
const (
- L3 = L / 3;
+ L3 = LogB / 3;
B3 = 1 << L3;
M3 = B3 - 1;
+
+ L2 = L3 * 2;
+ B2 = 1 << L2;
+ M2 = B2 - 1;
+
+ L = L3 * 3;
+ B = 1 << L;
+ M = B - 1;
)
type (
- Word3 uint32;
- Natural3 [] Word3;
+ Digit3 uint32;
+ Digit uint64;
)
}
-func IsSmall(x Word) bool {
- return x < H;
+func IsSmall(x Digit) bool {
+ return x < 1<<LogH;
}
-func Split(x Word) (Word, Word) {
+func Split(x Digit) (Digit, Digit) {
return x>>L, x&M;
}
-export func Dump(x *[]Word) {
+export func Dump(x *[]Digit) {
+ print("[", len(x), "]");
+ for i := len(x) - 1; i >= 0; i-- {
+ print(" ", x[i]);
+ }
+ println();
+}
+
+
+export func Dump3(x *[]Digit3) {
print("[", len(x), "]");
for i := len(x) - 1; i >= 0; i-- {
print(" ", x[i]);
// ----------------------------------------------------------------------------
// Natural numbers
-export type Natural []Word;
+export type Natural []Digit;
export var NatZero *Natural = new(Natural, 0);
-export func NewNat(x Word) *Natural {
+export func Nat(x Digit) *Natural {
var z *Natural;
switch {
case x == 0:
}
-func Normalize3(x *Natural3) *Natural3 {
+func Normalize3(x *[]Digit3) *[]Digit3 {
n := len(x);
for n > 0 && x[n - 1] == 0 { n-- }
if n < len(x) {
assert(n >= m);
z := new(Natural, n + 1);
- c := Word(0);
+ c := Digit(0);
for i := 0; i < m; i++ { c, z[i] = Split(x[i] + y[i] + c); }
for i := m; i < n; i++ { c, z[i] = Split(x[i] + c); }
z[n] = c;
assert(n >= m);
z := new(Natural, n);
- c := Word(0);
- for i := 0; i < m; i++ { c, z[i] = Split(x[i] - y[i] + c); }
+ c := Digit(0);
+ for i := 0; i < m; i++ { c, z[i] = Split(x[i] - y[i] + c); } // TODO verify asr!!!
for i := m; i < n; i++ { c, z[i] = Split(x[i] + c); }
assert(c == 0); // x.Sub(y) must be called with x >= y
// Computes x = x*a + c (in place) for "small" a's.
-func (x* Natural) MulAdd1(a, c Word) *Natural {
+func (x* Natural) MulAdd1(a, c Digit) *Natural {
assert(IsSmall(a-1) && IsSmall(c));
n := len(x);
z := new(Natural, n + 1);
// Returns c = x*y div B, z = x*y mod B.
-func Mul1(x, y Word) (Word, Word) {
+func Mul1(x, y Digit) (Digit, Digit) {
// Split x and y into 2 sub-digits each (in base sqrt(B)),
// multiply the digits separately while avoiding overflow,
// and return the product as two separate digits.
for j := 0; j < m; j++ {
d := y[j];
if d != 0 {
- c := Word(0);
+ c := Digit(0);
for i := 0; i < n; i++ {
// z[i+j] += x[i]*d + c;
z1, z0 := Mul1(x[i], d);
}
-func Shl1(x, c Word, s uint) (Word, Word) {
+func Shl1(x, c Digit, s uint) (Digit, Digit) {
assert(s <= LogB);
return x >> (LogB - s), x << s & M | c
}
-func Shr1(x, c Word, s uint) (Word, Word) {
+func Shr1(x, c Digit, s uint) (Digit, Digit) {
assert(s <= LogB);
return x << (LogB - s) & M, x >> s | c
}
s = s % LogB;
z := new(Natural, n + si + 1);
- c := Word(0);
+ c := Digit(0);
for i := 0; i < n; i++ { c, z[i+si] = Shl1(x[i], c, s); }
z[n+si] = c;
assert(si <= n);
z := new(Natural, n - si);
- c := Word(0);
+ c := Digit(0);
for i := n - 1; i >= si; i-- { c, z[i-si] = Shr1(x[i], c, s); }
return Normalize(z);
}
-func SplitBase(x *Natural) *Natural3 {
- xl := len(x);
- z := new(Natural3, xl * 3);
- for i, j := 0, 0; i < xl; i, j = i + 1, j + 3 {
+// DivMod needs multi-precision division which is not available if Digit
+// is already using the largest uint size. Split base before division,
+// and merge again after. Each Digit is split into 3 Digit3's.
+
+func SplitBase(x *Natural) *[]Digit3 {
+ // TODO Use Log() for better result - don't need Normalize3 at the end!
+ n := len(x);
+ z := new([]Digit3, n*3 + 1); // add space for extra digit (used by DivMod)
+ for i, j := 0, 0; i < n; i, j = i+1, j+3 {
t := x[i];
- z[j] = Word3(t & M3); t >>= L3; j++;
- z[j] = Word3(t & M3); t >>= L3; j++;
- z[j] = Word3(t & M3); t >>= L3; j++;
+ z[j+0] = Digit3(t >> (L3*0) & M3);
+ z[j+1] = Digit3(t >> (L3*1) & M3);
+ z[j+2] = Digit3(t >> (L3*2) & M3);
}
return Normalize3(z);
}
-func Scale(x *Natural, f Word) *Natural3 {
- return nil;
+func MergeBase(x *[]Digit3) *Natural {
+ i := len(x);
+ j := (i+2)/3;
+ z := new(Natural, j);
+
+ switch i%3 {
+ case 1: z[j-1] = Digit(x[i-1]); i--; j--;
+ case 2: z[j-1] = Digit(x[i-1])<<L3 | Digit(x[i-2]); i -= 2; j--;
+ case 0:
+ }
+
+ for i >= 3 {
+ z[j-1] = ((Digit(x[i-1])<<L3) | Digit(x[i-2]))<<L3 | Digit(x[i-3]);
+ i -= 3;
+ j--;
+ }
+ assert(j == 0);
+
+ return Normalize(z);
+}
+
+
+func Split3(x Digit) (Digit, Digit3) {
+ return uint64(int64(x)>>L3), Digit3(x&M3)
+}
+
+
+func Product(x *[]Digit3, y Digit) {
+ n := len(x);
+ c := Digit(0);
+ for i := 0; i < n; i++ { c, x[i] = Split3(Digit(x[i])*y + c) }
+ assert(c == 0);
}
-func TrialDigit(r, d *Natural3, k, m int) Word {
- km := k + m;
- assert(2 <= m && m <= km);
- r3 := (Word(r[km]) << L3 + Word(r[km - 1])) << L3 + Word(r[km - 2]);
- d2 := Word(d[m - 1]) << L3 + Word(d[m - 2]);
- qt := r3 / d2;
- if qt >= B {
- qt = B - 1;
+func Quotient(x *[]Digit3, y Digit) {
+ n := len(x);
+ c := Digit(0);
+ for i := n-1; i >= 0; i-- {
+ t := c*B3 + Digit(x[i]);
+ c, x[i] = t%y, Digit3(t/y);
}
- return qt;
+ assert(c == 0);
}
-func DivMod(x, y *Natural) {
- xl := len(x);
- yl := len(y);
- assert(2 <= yl && yl <= xl); // use special-case algorithm otherwise
+// Division and modulo computation - destroys x and y. Based on the
+// algorithms described in:
+//
+// 1) D. Knuth, "The Art of Computer Programming. Volume 2. Seminumerical
+// Algorithms." Addison-Wesley, Reading, 1969.
+//
+// 2) P. Brinch Hansen, Multiple-length division revisited: A tour of the
+// minefield. "Software - Practice and Experience 24", (June 1994),
+// 579-601. John Wiley & Sons, Ltd.
+//
+// Specifically, the inplace computation of quotient and remainder
+// is described in 1), while 2) provides the background for a more
+// accurate initial guess of the trial digit.
+
+func DivMod(x, y *[]Digit3) (*[]Digit3, *[]Digit3) {
+ const b = B3;
- f := B / (y[yl - 1] + 1);
- r := Scale(x, f);
- d := Scale(y, f);
- n := len(r);
- m := len(d);
+ n := len(x);
+ m := len(y);
+ assert(m > 0); // division by zero
+ assert(n+1 <= cap(x)); // space for one extra digit (should it be == ?)
+ x = x[0 : n + 1];
- for k := n - m; k >= 0; k-- {
- qt := TrialDigit(r, d, k, m);
+ if m == 1 {
+ // division by single digit
+ d := Digit(y[0]);
+ c := Digit(0);
+ for i := n; i > 0; i-- {
+ t := c*b + Digit(x[i-1]);
+ c, x[i] = t%d, Digit3(t/d);
+ }
+ x[0] = Digit3(c);
+
+ } else if m > n {
+ // quotient = 0, remainder = x
+ // TODO in this case we shouldn't even split base - FIX THIS
+ m = n;
+
+ } else {
+ // general case
+ assert(2 <= m && m <= n);
+ assert(x[n] == 0);
+
+ // normalize x and y
+ f := b/(Digit(y[m-1]) + 1);
+ Product(x, f);
+ Product(y, f);
+ assert(b/2 <= y[m-1] && y[m-1] < b); // incorrect scaling
+
+ d2 := Digit(y[m-1])*b + Digit(y[m-2]);
+ for i := n-m; i >= 0; i-- {
+ k := i+m;
+
+ // compute trial digit
+ r3 := (Digit(x[k])*b + Digit(x[k-1]))*b + Digit(x[k-2]);
+ q := r3/d2;
+ if q >= b { q = b-1 }
+
+ // subtract y*q
+ c := Digit(0);
+ for j := 0; j < m; j++ {
+ c, x[i+j] = Split3(c + Digit(x[i+j]) - Digit(y[j])*q);
+ }
+
+ // correct if trial digit was too large
+ if c + Digit(x[k]) != 0 {
+ // add y
+ c := Digit(0);
+ for j := 0; j < m; j++ {
+ c, x[i+j] = Split3(c + Digit(x[i+j]) + Digit(y[j]));
+ }
+ // correct trial digit
+ q--;
+ }
+
+ x[k] = Digit3(q);
+ }
+ // undo normalization for remainder
+ Quotient(x[0 : m], f);
}
+
+ return x[m : n+1], x[0 : m];
}
func (x *Natural) Div(y *Natural) *Natural {
- panic("UNIMPLEMENTED");
- return nil;
+ q, r := DivMod(SplitBase(x), SplitBase(y));
+ return MergeBase(q);
}
func (x *Natural) Mod(y *Natural) *Natural {
- panic("UNIMPLEMENTED");
- return nil;
+ q, r := DivMod(SplitBase(x), SplitBase(y));
+ return MergeBase(r);
}
func (x *Natural) Cmp(y *Natural) int {
- xl := len(x);
- yl := len(y);
+ n := len(x);
+ m := len(y);
- if xl != yl || xl == 0 {
- return xl - yl;
+ if n != m || n == 0 {
+ return n - m;
}
- i := xl - 1;
+ i := n - 1;
for i > 0 && x[i] == y[i] { i--; }
d := 0;
}
-func Log1(x Word) int {
+func Log1(x Digit) int {
n := -1;
for x != 0 { x >>= 1; n++; }
return n;
// Computes x = x div d (in place - the recv maybe modified) for "small" d's.
// Returns updated x and x mod d.
-func (x *Natural) DivMod1(d Word) (*Natural, Word) {
+func (x *Natural) DivMod1(d Digit) (*Natural, Digit) {
assert(0 < d && IsSmall(d - 1));
- c := Word(0);
+ c := Digit(0);
for i := len(x) - 1; i >= 0; i-- {
c = c<<L + x[i];
x[i] = c/d;
}
-func (x *Natural) String(base Word) string {
+func (x *Natural) String(base Digit) string {
if x.IsZero() {
return "0";
}
x = Copy(x); // don't destroy recv
for !x.IsZero() {
i--;
- var d Word;
+ var d Digit;
x, d = x.DivMod1(base);
s[i] = hex[d];
};
}
-func MulRange(a, b Word) *Natural {
+export func MulRange(a, b Digit) *Natural {
switch {
- case a > b: return NewNat(1);
- case a == b: return NewNat(a);
- case a + 1 == b: return NewNat(a).Mul(NewNat(b));
+ case a > b: return Nat(1);
+ case a == b: return Nat(a);
+ case a + 1 == b: return Nat(a).Mul(Nat(b));
}
m := (a + b)>>1;
assert(a <= m && m < b);
}
-export func Fact(n Word) *Natural {
+export func Fact(n Digit) *Natural {
// Using MulRange() instead of the basic for-loop
// lead to faster factorial computation.
return MulRange(2, n);
}
-func HexValue(ch byte) Word {
- d := Word(H);
+func HexValue(ch byte) Digit {
+ d := Digit(1 << LogH);
switch {
- case '0' <= ch && ch <= '9': d = Word(ch - '0');
- case 'a' <= ch && ch <= 'f': d = Word(ch - 'a') + 10;
- case 'A' <= ch && ch <= 'F': d = Word(ch - 'A') + 10;
+ case '0' <= ch && ch <= '9': d = Digit(ch - '0');
+ case 'a' <= ch && ch <= 'f': d = Digit(ch - 'a') + 10;
+ case 'A' <= ch && ch <= 'F': d = Digit(ch - 'A') + 10;
}
return d;
}
// TODO auto-detect base if base argument is 0
-export func NatFromString(s string, base Word) *Natural {
+export func NatFromString(s string, base Digit) *Natural {
x := NatZero;
for i := 0; i < len(s); i++ {
d := HexValue(s[i]);
}
+export func Int(x int64) *Integer {
+ return nil;
+}
+
+
func (x *Integer) Add(y *Integer) *Integer {
var z *Integer;
if x.sign == y.sign {
}
-func (x *Integer) String(base Word) string {
+func (x *Integer) String(base Digit) string {
if x.mant.IsZero() {
return "0";
}
}
-export func IntFromString(s string, base Word) *Integer {
+export func IntFromString(s string, base Digit) *Integer {
// get sign, if any
sign := false;
if len(s) > 0 && (s[0] == '-' || s[0] == '+') {
package main
-import Bignum "bignum"
+import Big "bignum"
const (
sa = "991";
var (
- a = Bignum.NatFromString(sa, 10);
- b = Bignum.NatFromString(sb, 10);
- c = Bignum.NatFromString(sc, 10);
+ a = Big.NatFromString(sa, 10);
+ b = Big.NatFromString(sb, 10);
+ c = Big.NatFromString(sc, 10);
)
var test_msg string;
-func TEST(n int, b bool) {
+func TEST(n uint, b bool) {
if !b {
panic("TEST failed: ", test_msg, "(", n, ")\n");
}
}
+func TEST_EQ(n uint, x, y *Big.Natural) {
+ if x.Cmp(y) != 0 {
+ println("TEST failed: ", test_msg, "(", n, ")\n");
+ println("x = ", x.String(10));
+ println("y = ", y.String(10));
+ panic();
+ }
+}
+
+
func TestConv() {
test_msg = "TestConv";
- TEST(0, a.Cmp(Bignum.NewNat(991)) == 0);
- TEST(1, b.Cmp(Bignum.Fact(20)) == 0);
- TEST(2, c.Cmp(Bignum.Fact(100)) == 0);
+ TEST(0, a.Cmp(Big.Nat(991)) == 0);
+ TEST(1, b.Cmp(Big.Fact(20)) == 0);
+ TEST(2, c.Cmp(Big.Fact(100)) == 0);
TEST(3, a.String(10) == sa);
TEST(4, b.String(10) == sb);
TEST(5, c.String(10) == sc);
TEST(1, c.Shr(1).Cmp(c) < 0);
test_msg = "TestShift2";
- for i := 0; i < 100; i++ {
- TEST(i, c.Shl(uint(i)).Shr(uint(i)).Cmp(c) == 0);
+ for i := uint(0); i < 100; i++ {
+ TEST(i, c.Shl(i).Shr(i).Cmp(c) == 0);
}
test_msg = "TestShift3L";
{ const m = 3;
p := b;
- f := Bignum.NewNat(1<<m);
- for i := 0; i < 100; i++ {
- TEST(i, b.Shl(uint(i*m)).Cmp(p) == 0);
+ f := Big.Nat(1<<m);
+ for i := uint(0); i < 100; i++ {
+ TEST_EQ(i, b.Shl(i*m), p);
p = p.Mul(f);
}
}
test_msg = "TestShift3R";
{ p := c;
- for i := 0; c.Cmp(Bignum.NatZero) == 0; i++ {
- TEST(i, c.Shr(uint(i)).Cmp(p) == 0);
+ for i := uint(0); c.Cmp(Big.NatZero) == 0; i++ {
+ TEST_EQ(i, c.Shr(i), p);
p = p.Shr(1);
}
}
}
+func TestMul() {
+ test_msg = "TestMulA";
+ TEST_EQ(0, b.Mul(Big.MulRange(0, 100)), Big.Nat(0));
+ TEST_EQ(0, b.Mul(Big.MulRange(21, 100)), c);
+
+ test_msg = "TestMulB";
+ const n = 100;
+ p := b.Mul(c).Shl(n);
+ for i := uint(0); i < n; i++ {
+ TEST_EQ(i, b.Shl(i).Mul(c.Shl(n-i)), p);
+ }
+}
+
+
+func TestDiv() {
+ test_msg = "TestDivA";
+ TEST_EQ(0, c.Div(Big.Nat(1)), c);
+ TEST_EQ(1, c.Div(Big.Nat(100)), Big.Fact(99));
+ TEST_EQ(2, b.Div(c), Big.Nat(0));
+ TEST_EQ(4, Big.Nat(1).Shl(100).Div(Big.Nat(1).Shl(90)), Big.Nat(1).Shl(10));
+ TEST_EQ(5, c.Div(b), Big.MulRange(21, 100));
+
+ test_msg = "TestDivB";
+ const n = 100;
+ p := Big.Fact(n);
+ for i := uint(0); i < n; i++ {
+ TEST_EQ(i, p.Div(Big.MulRange(1, uint64(i))), Big.MulRange(uint64(i+1), n));
+ }
+}
+
+
+func TestMod() {
+ test_msg = "TestModA";
+ for i := uint(0); ; i++ {
+ d := Big.Nat(1).Shl(i);
+ if d.Cmp(c) < 0 {
+ TEST_EQ(i, c.Add(d).Mod(c), d);
+ } else {
+ TEST_EQ(i, c.Add(d).Div(c), Big.Nat(2));
+ //TEST_EQ(i, c.Add(d).Mod(c), d.Sub(c));
+ break;
+ }
+ }
+}
+
+
func main() {
TestConv();
TestShift();
+ TestMul();
+ TestDiv();
+ TestMod();
print("PASSED\n");
}