Commit 496a9353 authored by Robert Griesemer's avatar Robert Griesemer

bignum: delete package - functionality subsumed by package big

R=rsc
CC=golang-dev
https://golang.org/cl/1742045
parent b2a919fc
......@@ -63,7 +63,6 @@ DIRS=\
encoding/hex\
encoding/pem\
exec\
exp/bignum\
exp/datafmt\
exp/draw\
exp/draw/x11\
......
# 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.
include ../../../Make.$(GOARCH)
TARG=exp/bignum
GOFILES=\
arith.go\
bignum.go\
integer.go\
rational.go\
include ../../../Make.pkg
// 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")
}
// 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 ·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 ·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 ·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
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
// 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 implements Newton-Raphson division and uses
// it as an additional test case for bignum.
//
// Division of x/y is achieved by computing r = 1/y to
// obtain the quotient q = x*r = x*(1/y) = x/y. The
// reciprocal r is the solution for f(x) = 1/x - y and
// the solution is approximated through iteration. The
// iteration does not require division.
package bignum
import "testing"
// An fpNat is a Natural scaled by a power of two
// (an unsigned floating point representation). The
// value of an fpNat x is x.m * 2^x.e .
//
type fpNat struct {
m Natural
e int
}
// sub computes x - y.
func (x fpNat) sub(y fpNat) fpNat {
switch d := x.e - y.e; {
case d < 0:
return fpNat{x.m.Sub(y.m.Shl(uint(-d))), x.e}
case d > 0:
return fpNat{x.m.Shl(uint(d)).Sub(y.m), y.e}
}
return fpNat{x.m.Sub(y.m), x.e}
}
// mul2 computes x*2.
func (x fpNat) mul2() fpNat { return fpNat{x.m, x.e + 1} }
// mul computes x*y.
func (x fpNat) mul(y fpNat) fpNat { return fpNat{x.m.Mul(y.m), x.e + y.e} }
// mant computes the (possibly truncated) Natural representation
// of an fpNat x.
//
func (x fpNat) mant() Natural {
switch {
case x.e > 0:
return x.m.Shl(uint(x.e))
case x.e < 0:
return x.m.Shr(uint(-x.e))
}
return x.m
}
// nrDivEst computes an estimate of the quotient q = x0/y0 and returns q.
// q may be too small (usually by 1).
//
func nrDivEst(x0, y0 Natural) Natural {
if y0.IsZero() {
panic("division by zero")
return nil
}
// y0 > 0
if y0.Cmp(Nat(1)) == 0 {
return x0
}
// y0 > 1
switch d := x0.Cmp(y0); {
case d < 0:
return Nat(0)
case d == 0:
return Nat(1)
}
// x0 > y0 > 1
// Determine maximum result length.
maxLen := int(x0.Log2() - y0.Log2() + 1)
// In the following, each number x is represented
// as a mantissa x.m and an exponent x.e such that
// x = xm * 2^x.e.
x := fpNat{x0, 0}
y := fpNat{y0, 0}
// Determine a scale factor f = 2^e such that
// 0.5 <= y/f == y*(2^-e) < 1.0
// and scale y accordingly.
e := int(y.m.Log2()) + 1
y.e -= e
// t1
var c = 2.9142
const n = 14
t1 := fpNat{Nat(uint64(c * (1 << n))), -n}
// Compute initial value r0 for the reciprocal of y/f.
// r0 = t1 - 2*y
r := t1.sub(y.mul2())
two := fpNat{Nat(2), 0}
// Newton-Raphson iteration
p := Nat(0)
for i := 0; ; i++ {
// check if we are done
// TODO: Need to come up with a better test here
// as it will reduce computation time significantly.
// q = x*r/f
q := x.mul(r)
q.e -= e
res := q.mant()
if res.Cmp(p) == 0 {
return res
}
p = res
// r' = r*(2 - y*r)
r = r.mul(two.sub(y.mul(r)))
// reduce mantissa size
// TODO: Find smaller bound as it will reduce
// computation time massively.
d := int(r.m.Log2()+1) - maxLen
if d > 0 {
r = fpNat{r.m.Shr(uint(d)), r.e + d}
}
}
panic("unreachable")
return nil
}
func nrdiv(x, y Natural) (q, r Natural) {
q = nrDivEst(x, y)
r = x.Sub(y.Mul(q))
// if r is too large, correct q and r
// (usually one iteration)
for r.Cmp(y) >= 0 {
q = q.Add(Nat(1))
r = r.Sub(y)
}
return
}
func div(t *testing.T, x, y Natural) {
q, r := nrdiv(x, y)
qx, rx := x.DivMod(y)
if q.Cmp(qx) != 0 {
t.Errorf("x = %s, y = %s, got q = %s, want q = %s", x, y, q, qx)
}
if r.Cmp(rx) != 0 {
t.Errorf("x = %s, y = %s, got r = %s, want r = %s", x, y, r, rx)
}
}
func idiv(t *testing.T, x0, y0 uint64) { div(t, Nat(x0), Nat(y0)) }
func TestNRDiv(t *testing.T) {
idiv(t, 17, 18)
idiv(t, 17, 17)
idiv(t, 17, 1)
idiv(t, 17, 16)
idiv(t, 17, 10)
idiv(t, 17, 9)
idiv(t, 17, 8)
idiv(t, 17, 5)
idiv(t, 17, 3)
idiv(t, 1025, 512)
idiv(t, 7489595, 2)
idiv(t, 5404679459, 78495)
idiv(t, 7484890589595, 7484890589594)
div(t, Fact(100), Fact(91))
div(t, Fact(1000), Fact(991))
//div(t, Fact(10000), Fact(9991)); // takes too long - disabled for now
}
// 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.
// Rational numbers
package bignum
import "fmt"
// Rational represents a quotient a/b of arbitrary precision.
//
type Rational struct {
a *Integer // numerator
b Natural // denominator
}
// MakeRat makes a rational number given a numerator a and a denominator b.
//
func MakeRat(a *Integer, b Natural) *Rational {
f := a.mant.Gcd(b) // f > 0
if f.Cmp(Nat(1)) != 0 {
a = MakeInt(a.sign, a.mant.Div(f))
b = b.Div(f)
}
return &Rational{a, b}
}
// Rat creates a small rational number with value a0/b0.
//
func Rat(a0 int64, b0 int64) *Rational {
a, b := Int(a0), Int(b0)
if b.sign {
a = a.Neg()
}
return MakeRat(a, b.mant)
}
// Value returns the numerator and denominator of x.
//
func (x *Rational) Value() (numerator *Integer, denominator Natural) {
return x.a, x.b
}
// Predicates
// IsZero returns true iff x == 0.
//
func (x *Rational) IsZero() bool { return x.a.IsZero() }
// IsNeg returns true iff x < 0.
//
func (x *Rational) IsNeg() bool { return x.a.IsNeg() }
// IsPos returns true iff x > 0.
//
func (x *Rational) IsPos() bool { return x.a.IsPos() }
// IsInt returns true iff x can be written with a denominator 1
// 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 }
// Operations
// Neg returns the negated value of x.
//
func (x *Rational) Neg() *Rational { return MakeRat(x.a.Neg(), x.b) }
// Add returns the sum x + y.
//
func (x *Rational) Add(y *Rational) *Rational {
return MakeRat((x.a.MulNat(y.b)).Add(y.a.MulNat(x.b)), x.b.Mul(y.b))
}
// Sub returns the difference x - y.
//
func (x *Rational) Sub(y *Rational) *Rational {
return MakeRat((x.a.MulNat(y.b)).Sub(y.a.MulNat(x.b)), x.b.Mul(y.b))
}
// Mul returns the product x * y.
//
func (x *Rational) Mul(y *Rational) *Rational { return MakeRat(x.a.Mul(y.a), x.b.Mul(y.b)) }
// Quo returns the quotient x / y for y != 0.
// If y == 0, a division-by-zero run-time error occurs.
//
func (x *Rational) Quo(y *Rational) *Rational {
a := x.a.MulNat(y.b)
b := y.a.MulNat(x.b)
if b.IsNeg() {
a = a.Neg()
}
return MakeRat(a, b.mant)
}
// Cmp compares x and y. The result is an int value
//
// < 0 if x < y
// == 0 if x == y
// > 0 if x > y
//
func (x *Rational) Cmp(y *Rational) int { return (x.a.MulNat(y.b)).Cmp(y.a.MulNat(x.b)) }
// ToString converts x to a string for a given base, with 2 <= base <= 16.
// The string representation is of the form "n" if x is an integer; otherwise
// it is of form "n/d".
//
func (x *Rational) ToString(base uint) string {
s := x.a.ToString(base)
if !x.IsInt() {
s += "/" + x.b.ToString(base)
}
return s
}
// String converts x to its decimal string representation.
// x.String() is the same as x.ToString(10).
//
func (x *Rational) String() string { return x.ToString(10) }
// Format is a support routine for fmt.Formatter. It accepts
// the formats 'b' (binary), 'o' (octal), and 'x' (hexadecimal).
//
func (x *Rational) Format(h fmt.State, c int) { fmt.Fprintf(h, "%s", x.ToString(fmtbase(c))) }
// RatFromString returns the rational number corresponding to the
// longest possible prefix of s representing a rational number in a
// given conversion base, the actual conversion base used, and the
// prefix length. The syntax of a rational number is:
//
// rational = mantissa [exponent] .
// mantissa = integer ('/' natural | '.' natural) .
// exponent = ('e'|'E') integer .
//
// If the base argument is 0, the string prefix determines the actual
// conversion base for the mantissa. A prefix of ``0x'' or ``0X'' selects
// base 16; the ``0'' prefix selects base 8. Otherwise the selected base is 10.
// If the mantissa is represented via a division, both the numerator and
// denominator may have different base prefixes; in that case the base of
// of the numerator is returned. If the mantissa contains a decimal point,
// the base for the fractional part is the same as for the part before the
// decimal point and the fractional part does not accept a base prefix.
// The base for the exponent is always 10.
//
func RatFromString(s string, base uint) (*Rational, uint, int) {
// read numerator
a, abase, alen := IntFromString(s, base)
b := Nat(1)
// read denominator or fraction, if any
var blen int
if alen < len(s) {
ch := s[alen]
if ch == '/' {
alen++
b, base, blen = NatFromString(s[alen:], base)
} else if ch == '.' {
alen++
b, base, blen = NatFromString(s[alen:], abase)
assert(base == abase)
f := Nat(uint64(base)).Pow(uint(blen))
a = MakeInt(a.sign, a.mant.Mul(f).Add(b))
b = f
}
}
// read exponent, if any
rlen := alen + blen
if rlen < len(s) {
ch := s[rlen]
if ch == 'e' || ch == 'E' {
rlen++
e, _, elen := IntFromString(s[rlen:], 10)
rlen += elen
m := Nat(10).Pow(uint(e.mant.Value()))
if e.sign {
b = b.Mul(m)
} else {
a = a.MulNat(m)
}
}
}
return MakeRat(a, b), base, rlen
}
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment