PageRenderTime 51ms CodeModel.GetById 23ms app.highlight 24ms RepoModel.GetById 2ms app.codeStats 0ms

/big/rat.go

https://code.google.com/p/algogo/
Go | 371 lines | 267 code | 44 blank | 60 comment | 48 complexity | 7bce3e25a89f6bfdff156dcd9f748b05 MD5 | raw file
  1// Copyright 2010 The Go Authors. All rights reserved.
  2// Use of this source code is governed by a BSD-style
  3// license that can be found in the LICENSE file.
  4
  5// This file implements multi-precision rational numbers.
  6
  7package big
  8
  9import (
 10	"encoding/binary"
 11	"fmt"
 12	"os"
 13	"strings"
 14)
 15
 16// A Rat represents a quotient a/b of arbitrary precision. The zero value for
 17// a Rat, 0/0, is not a legal Rat.
 18type Rat struct {
 19	a Int
 20	b nat
 21}
 22
 23// NewRat creates a new Rat with numerator a and denominator b.
 24func NewRat(a, b int64) *Rat {
 25	return new(Rat).SetFrac64(a, b)
 26}
 27
 28// SetFrac sets z to a/b and returns z.
 29func (z *Rat) SetFrac(a, b *Int) *Rat {
 30	z.a.Set(a)
 31	z.a.neg = a.neg != b.neg
 32	z.b = z.b.set(b.abs)
 33	return z.norm()
 34}
 35
 36// SetFrac64 sets z to a/b and returns z.
 37func (z *Rat) SetFrac64(a, b int64) *Rat {
 38	z.a.SetInt64(a)
 39	if b < 0 {
 40		b = -b
 41		z.a.neg = !z.a.neg
 42	}
 43	z.b = z.b.setUint64(uint64(b))
 44	return z.norm()
 45}
 46
 47// SetInt sets z to x (by making a copy of x) and returns z.
 48func (z *Rat) SetInt(x *Int) *Rat {
 49	z.a.Set(x)
 50	z.b = z.b.setWord(1)
 51	return z
 52}
 53
 54// SetInt64 sets z to x and returns z.
 55func (z *Rat) SetInt64(x int64) *Rat {
 56	z.a.SetInt64(x)
 57	z.b = z.b.setWord(1)
 58	return z
 59}
 60
 61// Sign returns:
 62//
 63//	-1 if x <  0
 64//	 0 if x == 0
 65//	+1 if x >  0
 66//
 67func (x *Rat) Sign() int {
 68	return x.a.Sign()
 69}
 70
 71// IsInt returns true if the denominator of x is 1.
 72func (x *Rat) IsInt() bool {
 73	return len(x.b) == 1 && x.b[0] == 1
 74}
 75
 76// Num returns the numerator of z; it may be <= 0.
 77// The result is a reference to z's numerator; it
 78// may change if a new value is assigned to z.
 79func (z *Rat) Num() *Int {
 80	return &z.a
 81}
 82
 83// Denom returns the denominator of z; it is always > 0.
 84// The result is a reference to z's denominator; it
 85// may change if a new value is assigned to z.
 86func (z *Rat) Denom() *Int {
 87	return &Int{false, z.b}
 88}
 89
 90func gcd(x, y nat) nat {
 91	// Euclidean algorithm.
 92	var a, b nat
 93	a = a.set(x)
 94	b = b.set(y)
 95	for len(b) != 0 {
 96		var q, r nat
 97		_, r = q.div(r, a, b)
 98		a = b
 99		b = r
100	}
101	return a
102}
103
104func (z *Rat) norm() *Rat {
105	f := gcd(z.a.abs, z.b)
106	if len(z.a.abs) == 0 {
107		// z == 0
108		z.a.neg = false // normalize sign
109		z.b = z.b.setWord(1)
110		return z
111	}
112	if f.cmp(natOne) != 0 {
113		z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f)
114		z.b, _ = z.b.div(nil, z.b, f)
115	}
116	return z
117}
118
119func mulNat(x *Int, y nat) *Int {
120	var z Int
121	z.abs = z.abs.mul(x.abs, y)
122	z.neg = len(z.abs) > 0 && x.neg
123	return &z
124}
125
126// Cmp compares x and y and returns:
127//
128//   -1 if x <  y
129//    0 if x == y
130//   +1 if x >  y
131//
132func (x *Rat) Cmp(y *Rat) (r int) {
133	return mulNat(&x.a, y.b).Cmp(mulNat(&y.a, x.b))
134}
135
136// Abs sets z to |x| (the absolute value of x) and returns z.
137func (z *Rat) Abs(x *Rat) *Rat {
138	z.a.Abs(&x.a)
139	z.b = z.b.set(x.b)
140	return z
141}
142
143// Add sets z to the sum x+y and returns z.
144func (z *Rat) Add(x, y *Rat) *Rat {
145	a1 := mulNat(&x.a, y.b)
146	a2 := mulNat(&y.a, x.b)
147	z.a.Add(a1, a2)
148	z.b = z.b.mul(x.b, y.b)
149	return z.norm()
150}
151
152// Sub sets z to the difference x-y and returns z.
153func (z *Rat) Sub(x, y *Rat) *Rat {
154	a1 := mulNat(&x.a, y.b)
155	a2 := mulNat(&y.a, x.b)
156	z.a.Sub(a1, a2)
157	z.b = z.b.mul(x.b, y.b)
158	return z.norm()
159}
160
161// Mul sets z to the product x*y and returns z.
162func (z *Rat) Mul(x, y *Rat) *Rat {
163	z.a.Mul(&x.a, &y.a)
164	z.b = z.b.mul(x.b, y.b)
165	return z.norm()
166}
167
168// Quo sets z to the quotient x/y and returns z.
169// If y == 0, a division-by-zero run-time panic occurs.
170func (z *Rat) Quo(x, y *Rat) *Rat {
171	if len(y.a.abs) == 0 {
172		panic("division by zero")
173	}
174	a := mulNat(&x.a, y.b)
175	b := mulNat(&y.a, x.b)
176	z.a.abs = a.abs
177	z.b = b.abs
178	z.a.neg = a.neg != b.neg
179	return z.norm()
180}
181
182// Neg sets z to -x (by making a copy of x if necessary) and returns z.
183func (z *Rat) Neg(x *Rat) *Rat {
184	z.a.Neg(&x.a)
185	z.b = z.b.set(x.b)
186	return z
187}
188
189// Set sets z to x (by making a copy of x if necessary) and returns z.
190func (z *Rat) Set(x *Rat) *Rat {
191	z.a.Set(&x.a)
192	z.b = z.b.set(x.b)
193	return z
194}
195
196func ratTok(ch int) bool {
197	return strings.IndexRune("+-/0123456789.eE", ch) >= 0
198}
199
200// Scan is a support routine for fmt.Scanner. It accepts the formats
201// 'e', 'E', 'f', 'F', 'g', 'G', and 'v'. All formats are equivalent.
202func (z *Rat) Scan(s fmt.ScanState, ch int) os.Error {
203	tok, err := s.Token(true, ratTok)
204	if err != nil {
205		return err
206	}
207	if strings.IndexRune("efgEFGv", ch) < 0 {
208		return os.NewError("Rat.Scan: invalid verb")
209	}
210	if _, ok := z.SetString(string(tok)); !ok {
211		return os.NewError("Rat.Scan: invalid syntax")
212	}
213	return nil
214}
215
216// SetString sets z to the value of s and returns z and a boolean indicating
217// success. s can be given as a fraction "a/b" or as a floating-point number
218// optionally followed by an exponent. If the operation failed, the value of z
219// is undefined.
220func (z *Rat) SetString(s string) (*Rat, bool) {
221	if len(s) == 0 {
222		return z, false
223	}
224
225	// check for a quotient
226	sep := strings.Index(s, "/")
227	if sep >= 0 {
228		if _, ok := z.a.SetString(s[0:sep], 10); !ok {
229			return z, false
230		}
231		s = s[sep+1:]
232		var err os.Error
233		if z.b, _, err = z.b.scan(strings.NewReader(s), 10); err != nil {
234			return z, false
235		}
236		return z.norm(), true
237	}
238
239	// check for a decimal point
240	sep = strings.Index(s, ".")
241	// check for an exponent
242	e := strings.IndexAny(s, "eE")
243	var exp Int
244	if e >= 0 {
245		if e < sep {
246			// The E must come after the decimal point.
247			return z, false
248		}
249		if _, ok := exp.SetString(s[e+1:], 10); !ok {
250			return z, false
251		}
252		s = s[0:e]
253	}
254	if sep >= 0 {
255		s = s[0:sep] + s[sep+1:]
256		exp.Sub(&exp, NewInt(int64(len(s)-sep)))
257	}
258
259	if _, ok := z.a.SetString(s, 10); !ok {
260		return z, false
261	}
262	powTen := nat{}.expNN(natTen, exp.abs, nil)
263	if exp.neg {
264		z.b = powTen
265		z.norm()
266	} else {
267		z.a.abs = z.a.abs.mul(z.a.abs, powTen)
268		z.b = z.b.setWord(1)
269	}
270
271	return z, true
272}
273
274// String returns a string representation of z in the form "a/b" (even if b == 1).
275func (z *Rat) String() string {
276	return z.a.String() + "/" + z.b.decimalString()
277}
278
279// RatString returns a string representation of z in the form "a/b" if b != 1,
280// and in the form "a" if b == 1.
281func (z *Rat) RatString() string {
282	if z.IsInt() {
283		return z.a.String()
284	}
285	return z.String()
286}
287
288// FloatString returns a string representation of z in decimal form with prec
289// digits of precision after the decimal point and the last digit rounded.
290func (z *Rat) FloatString(prec int) string {
291	if z.IsInt() {
292		s := z.a.String()
293		if prec > 0 {
294			s += "." + strings.Repeat("0", prec)
295		}
296		return s
297	}
298
299	q, r := nat{}.div(nat{}, z.a.abs, z.b)
300
301	p := natOne
302	if prec > 0 {
303		p = nat{}.expNN(natTen, nat{}.setUint64(uint64(prec)), nil)
304	}
305
306	r = r.mul(r, p)
307	r, r2 := r.div(nat{}, r, z.b)
308
309	// see if we need to round up
310	r2 = r2.add(r2, r2)
311	if z.b.cmp(r2) <= 0 {
312		r = r.add(r, natOne)
313		if r.cmp(p) >= 0 {
314			q = nat{}.add(q, natOne)
315			r = nat{}.sub(r, p)
316		}
317	}
318
319	s := q.decimalString()
320	if z.a.neg {
321		s = "-" + s
322	}
323
324	if prec > 0 {
325		rs := r.decimalString()
326		leadingZeros := prec - len(rs)
327		s += "." + strings.Repeat("0", leadingZeros) + rs
328	}
329
330	return s
331}
332
333// Gob codec version. Permits backward-compatible changes to the encoding.
334const ratGobVersion byte = 1
335
336// GobEncode implements the gob.GobEncoder interface.
337func (z *Rat) GobEncode() ([]byte, os.Error) {
338	buf := make([]byte, 1+4+(len(z.a.abs)+len(z.b))*_S) // extra bytes for version and sign bit (1), and numerator length (4)
339	i := z.b.bytes(buf)
340	j := z.a.abs.bytes(buf[0:i])
341	n := i - j
342	if int(uint32(n)) != n {
343		// this should never happen
344		return nil, os.NewError("Rat.GobEncode: numerator too large")
345	}
346	binary.BigEndian.PutUint32(buf[j-4:j], uint32(n))
347	j -= 1 + 4
348	b := ratGobVersion << 1 // make space for sign bit
349	if z.a.neg {
350		b |= 1
351	}
352	buf[j] = b
353	return buf[j:], nil
354}
355
356// GobDecode implements the gob.GobDecoder interface.
357func (z *Rat) GobDecode(buf []byte) os.Error {
358	if len(buf) == 0 {
359		return os.NewError("Rat.GobDecode: no data")
360	}
361	b := buf[0]
362	if b>>1 != ratGobVersion {
363		return os.NewError(fmt.Sprintf("Rat.GobDecode: encoding version %d not supported", b>>1))
364	}
365	const j = 1 + 4
366	i := j + binary.BigEndian.Uint32(buf[j-4:j])
367	z.a.neg = b&1 != 0
368	z.a.abs = z.a.abs.setBytes(buf[j:i])
369	z.b = z.b.setBytes(buf[i:])
370	return nil
371}