// Copyright 2010 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 multi-precision rational numbers.

package big

import (
	
	
)

// A Rat represents a quotient a/b of arbitrary precision.
// The zero value for a Rat represents the value 0.
//
// Operations always take pointer arguments (*Rat) rather
// than Rat values, and each unique Rat value requires
// its own unique *Rat pointer. To "copy" a Rat value,
// an existing (or newly allocated) Rat must be set to
// a new value using the Rat.Set method; shallow copies
// of Rats are not supported and may lead to errors.
type Rat struct {
	// To make zero values for Rat work w/o initialization,
	// a zero value of b (len(b) == 0) acts like b == 1. At
	// the earliest opportunity (when an assignment to the Rat
	// is made), such uninitialized denominators are set to 1.
	// a.neg determines the sign of the Rat, b.neg is ignored.
	a, b Int
}

// NewRat creates a new Rat with numerator a and denominator b.
func (,  int64) *Rat {
	return new(Rat).SetFrac64(, )
}

// SetFloat64 sets z to exactly f and returns z.
// If f is not finite, SetFloat returns nil.
func ( *Rat) ( float64) *Rat {
	const  = 1<<11 - 1
	 := math.Float64bits()
	 :=  & (1<<52 - 1)
	 := int(( >> 52) & )
	switch  {
	case : // non-finite
		return nil
	case 0: // denormal
		 -= 1022
	default: // normal
		 |= 1 << 52
		 -= 1023
	}

	 := 52 - 

	// Optimization (?): partially pre-normalise.
	for &1 == 0 &&  > 0 {
		 >>= 1
		--
	}

	.a.SetUint64()
	.a.neg =  < 0
	.b.Set(intOne)
	if  > 0 {
		.b.Lsh(&.b, uint())
	} else {
		.a.Lsh(&.a, uint(-))
	}
	return .norm()
}

// quotToFloat32 returns the non-negative float32 value
// nearest to the quotient a/b, using round-to-even in
// halfway cases. It does not mutate its arguments.
// Preconditions: b is non-zero; a and b have no common factors.
func quotToFloat32(,  nat) ( float32,  bool) {
	const (
		// float size in bits
		 = 32

		// mantissa
		  = 23
		 =  + 1 // incl. implicit 1
		 =  + 1

		// exponent
		 =  - 
		 = 1<<(-1) - 1
		  = 1 - 
		  = 
	)

	// TODO(adonovan): specialize common degenerate cases: 1.0, integers.
	 := .bitLen()
	if  == 0 {
		return 0, true
	}
	 := .bitLen()
	if  == 0 {
		panic("division by zero")
	}

	// 1. Left-shift A or B such that quotient A/B is in [1<<Msize1, 1<<(Msize2+1)
	// (Msize2 bits if A < B when they are left-aligned, Msize2+1 bits if A >= B).
	// This is 2 or 3 more than the float32 mantissa field width of Msize:
	// - the optional extra bit is shifted away in step 3 below.
	// - the high-order 1 is omitted in "normal" representation;
	// - the low-order 1 will be used during rounding then discarded.
	 :=  - 
	var ,  nat
	 = .set()
	 = .set()
	if  :=  - ;  > 0 {
		 = .shl(, uint())
	} else if  < 0 {
		 = .shl(, uint(-))
	}

	// 2. Compute quotient and remainder (q, r).  NB: due to the
	// extra shift, the low-order bit of q is logically the
	// high-order bit of r.
	var  nat
	,  := .div(, , ) // (recycle a2)
	 := low32()
	 := len() > 0 // mantissa&1 && !haveRem => remainder is exactly half

	// 3. If quotient didn't fit in Msize2 bits, redo division by b2<<1
	// (in effect---we accomplish this incrementally).
	if >> == 1 {
		if &1 == 1 {
			 = true
		}
		 >>= 1
		++
	}
	if >> != 1 {
		panic(fmt.Sprintf("expected exactly %d bits of result", ))
	}

	// 4. Rounding.
	if - <=  &&  <=  {
		// Denormal case; lose 'shift' bits of precision.
		 := uint( - ( - 1)) // [1..Esize1)
		 :=  & (1<< - 1)
		 =  ||  != 0
		 >>= 
		 = 2 -  // == exp + shift
	}
	// Round q using round-half-to-even.
	 = !
	if &1 != 0 {
		 = false
		if  || &2 != 0 {
			if ++;  >= 1<< {
				// Complete rollover 11...1 => 100...0, so shift is safe
				 >>= 1
				++
			}
		}
	}
	 >>= 1 // discard rounding bit.  Mantissa now scaled by 1<<Msize1.

	 = float32(math.Ldexp(float64(), -))
	if math.IsInf(float64(), 0) {
		 = false
	}
	return
}

// quotToFloat64 returns the non-negative float64 value
// nearest to the quotient a/b, using round-to-even in
// halfway cases. It does not mutate its arguments.
// Preconditions: b is non-zero; a and b have no common factors.
func quotToFloat64(,  nat) ( float64,  bool) {
	const (
		// float size in bits
		 = 64

		// mantissa
		  = 52
		 =  + 1 // incl. implicit 1
		 =  + 1

		// exponent
		 =  - 
		 = 1<<(-1) - 1
		  = 1 - 
		  = 
	)

	// TODO(adonovan): specialize common degenerate cases: 1.0, integers.
	 := .bitLen()
	if  == 0 {
		return 0, true
	}
	 := .bitLen()
	if  == 0 {
		panic("division by zero")
	}

	// 1. Left-shift A or B such that quotient A/B is in [1<<Msize1, 1<<(Msize2+1)
	// (Msize2 bits if A < B when they are left-aligned, Msize2+1 bits if A >= B).
	// This is 2 or 3 more than the float64 mantissa field width of Msize:
	// - the optional extra bit is shifted away in step 3 below.
	// - the high-order 1 is omitted in "normal" representation;
	// - the low-order 1 will be used during rounding then discarded.
	 :=  - 
	var ,  nat
	 = .set()
	 = .set()
	if  :=  - ;  > 0 {
		 = .shl(, uint())
	} else if  < 0 {
		 = .shl(, uint(-))
	}

	// 2. Compute quotient and remainder (q, r).  NB: due to the
	// extra shift, the low-order bit of q is logically the
	// high-order bit of r.
	var  nat
	,  := .div(, , ) // (recycle a2)
	 := low64()
	 := len() > 0 // mantissa&1 && !haveRem => remainder is exactly half

	// 3. If quotient didn't fit in Msize2 bits, redo division by b2<<1
	// (in effect---we accomplish this incrementally).
	if >> == 1 {
		if &1 == 1 {
			 = true
		}
		 >>= 1
		++
	}
	if >> != 1 {
		panic(fmt.Sprintf("expected exactly %d bits of result", ))
	}

	// 4. Rounding.
	if - <=  &&  <=  {
		// Denormal case; lose 'shift' bits of precision.
		 := uint( - ( - 1)) // [1..Esize1)
		 :=  & (1<< - 1)
		 =  ||  != 0
		 >>= 
		 = 2 -  // == exp + shift
	}
	// Round q using round-half-to-even.
	 = !
	if &1 != 0 {
		 = false
		if  || &2 != 0 {
			if ++;  >= 1<< {
				// Complete rollover 11...1 => 100...0, so shift is safe
				 >>= 1
				++
			}
		}
	}
	 >>= 1 // discard rounding bit.  Mantissa now scaled by 1<<Msize1.

	 = math.Ldexp(float64(), -)
	if math.IsInf(, 0) {
		 = false
	}
	return
}

// Float32 returns the nearest float32 value for x and a bool indicating
// whether f represents x exactly. If the magnitude of x is too large to
// be represented by a float32, f is an infinity and exact is false.
// The sign of f always matches the sign of x, even if f == 0.
func ( *Rat) () ( float32,  bool) {
	 := .b.abs
	if len() == 0 {
		 = natOne
	}
	,  = quotToFloat32(.a.abs, )
	if .a.neg {
		 = -
	}
	return
}

// Float64 returns the nearest float64 value for x and a bool indicating
// whether f represents x exactly. If the magnitude of x is too large to
// be represented by a float64, f is an infinity and exact is false.
// The sign of f always matches the sign of x, even if f == 0.
func ( *Rat) () ( float64,  bool) {
	 := .b.abs
	if len() == 0 {
		 = natOne
	}
	,  = quotToFloat64(.a.abs, )
	if .a.neg {
		 = -
	}
	return
}

// SetFrac sets z to a/b and returns z.
// If b == 0, SetFrac panics.
func ( *Rat) (,  *Int) *Rat {
	.a.neg = .neg != .neg
	 := .abs
	if len() == 0 {
		panic("division by zero")
	}
	if &.a ==  || alias(.a.abs, ) {
		 = nat(nil).set() // make a copy
	}
	.a.abs = .a.abs.set(.abs)
	.b.abs = .b.abs.set()
	return .norm()
}

// SetFrac64 sets z to a/b and returns z.
// If b == 0, SetFrac64 panics.
func ( *Rat) (,  int64) *Rat {
	if  == 0 {
		panic("division by zero")
	}
	.a.SetInt64()
	if  < 0 {
		 = -
		.a.neg = !.a.neg
	}
	.b.abs = .b.abs.setUint64(uint64())
	return .norm()
}

// SetInt sets z to x (by making a copy of x) and returns z.
func ( *Rat) ( *Int) *Rat {
	.a.Set()
	.b.abs = .b.abs.setWord(1)
	return 
}

// SetInt64 sets z to x and returns z.
func ( *Rat) ( int64) *Rat {
	.a.SetInt64()
	.b.abs = .b.abs.setWord(1)
	return 
}

// SetUint64 sets z to x and returns z.
func ( *Rat) ( uint64) *Rat {
	.a.SetUint64()
	.b.abs = .b.abs.setWord(1)
	return 
}

// Set sets z to x (by making a copy of x) and returns z.
func ( *Rat) ( *Rat) *Rat {
	if  !=  {
		.a.Set(&.a)
		.b.Set(&.b)
	}
	if len(.b.abs) == 0 {
		.b.abs = .b.abs.setWord(1)
	}
	return 
}

// Abs sets z to |x| (the absolute value of x) and returns z.
func ( *Rat) ( *Rat) *Rat {
	.Set()
	.a.neg = false
	return 
}

// Neg sets z to -x and returns z.
func ( *Rat) ( *Rat) *Rat {
	.Set()
	.a.neg = len(.a.abs) > 0 && !.a.neg // 0 has no sign
	return 
}

// Inv sets z to 1/x and returns z.
// If x == 0, Inv panics.
func ( *Rat) ( *Rat) *Rat {
	if len(.a.abs) == 0 {
		panic("division by zero")
	}
	.Set()
	.a.abs, .b.abs = .b.abs, .a.abs
	return 
}

// Sign returns:
//
//	-1 if x <  0
//	 0 if x == 0
//	+1 if x >  0
//
func ( *Rat) () int {
	return .a.Sign()
}

// IsInt reports whether the denominator of x is 1.
func ( *Rat) () bool {
	return len(.b.abs) == 0 || .b.abs.cmp(natOne) == 0
}

// Num returns the numerator of x; it may be <= 0.
// The result is a reference to x's numerator; it
// may change if a new value is assigned to x, and vice versa.
// The sign of the numerator corresponds to the sign of x.
func ( *Rat) () *Int {
	return &.a
}

// Denom returns the denominator of x; it is always > 0.
// The result is a reference to x's denominator, unless
// x is an uninitialized (zero value) Rat, in which case
// the result is a new Int of value 1. (To initialize x,
// any operation that sets x will do, including x.Set(x).)
// If the result is a reference to x's denominator it
// may change if a new value is assigned to x, and vice versa.
func ( *Rat) () *Int {
	.b.neg = false // the result is always >= 0
	if len(.b.abs) == 0 {
		// Note: If this proves problematic, we could
		//       panic instead and require the Rat to
		//       be explicitly initialized.
		return &Int{abs: nat{1}}
	}
	return &.b
}

func ( *Rat) () *Rat {
	switch {
	case len(.a.abs) == 0:
		// z == 0; normalize sign and denominator
		.a.neg = false
		fallthrough
	case len(.b.abs) == 0:
		// z is integer; normalize denominator
		.b.abs = .b.abs.setWord(1)
	default:
		// z is fraction; normalize numerator and denominator
		 := .a.neg
		.a.neg = false
		.b.neg = false
		if  := NewInt(0).lehmerGCD(nil, nil, &.a, &.b); .Cmp(intOne) != 0 {
			.a.abs, _ = .a.abs.div(nil, .a.abs, .abs)
			.b.abs, _ = .b.abs.div(nil, .b.abs, .abs)
		}
		.a.neg = 
	}
	return 
}

// mulDenom sets z to the denominator product x*y (by taking into
// account that 0 values for x or y must be interpreted as 1) and
// returns z.
func mulDenom(, ,  nat) nat {
	switch {
	case len() == 0 && len() == 0:
		return .setWord(1)
	case len() == 0:
		return .set()
	case len() == 0:
		return .set()
	}
	return .mul(, )
}

// scaleDenom sets z to the product x*f.
// If f == 0 (zero value of denominator), z is set to (a copy of) x.
func ( *Int) ( *Int,  nat) {
	if len() == 0 {
		.Set()
		return
	}
	.abs = .abs.mul(.abs, )
	.neg = .neg
}

// Cmp compares x and y and returns:
//
//   -1 if x <  y
//    0 if x == y
//   +1 if x >  y
//
func ( *Rat) ( *Rat) int {
	var ,  Int
	.scaleDenom(&.a, .b.abs)
	.scaleDenom(&.a, .b.abs)
	return .Cmp(&)
}

// Add sets z to the sum x+y and returns z.
func ( *Rat) (,  *Rat) *Rat {
	var ,  Int
	.scaleDenom(&.a, .b.abs)
	.scaleDenom(&.a, .b.abs)
	.a.Add(&, &)
	.b.abs = mulDenom(.b.abs, .b.abs, .b.abs)
	return .norm()
}

// Sub sets z to the difference x-y and returns z.
func ( *Rat) (,  *Rat) *Rat {
	var ,  Int
	.scaleDenom(&.a, .b.abs)
	.scaleDenom(&.a, .b.abs)
	.a.Sub(&, &)
	.b.abs = mulDenom(.b.abs, .b.abs, .b.abs)
	return .norm()
}

// Mul sets z to the product x*y and returns z.
func ( *Rat) (,  *Rat) *Rat {
	if  ==  {
		// a squared Rat is positive and can't be reduced (no need to call norm())
		.a.neg = false
		.a.abs = .a.abs.sqr(.a.abs)
		if len(.b.abs) == 0 {
			.b.abs = .b.abs.setWord(1)
		} else {
			.b.abs = .b.abs.sqr(.b.abs)
		}
		return 
	}
	.a.Mul(&.a, &.a)
	.b.abs = mulDenom(.b.abs, .b.abs, .b.abs)
	return .norm()
}

// Quo sets z to the quotient x/y and returns z.
// If y == 0, Quo panics.
func ( *Rat) (,  *Rat) *Rat {
	if len(.a.abs) == 0 {
		panic("division by zero")
	}
	var ,  Int
	.scaleDenom(&.a, .b.abs)
	.scaleDenom(&.a, .b.abs)
	.a.abs = .abs
	.b.abs = .abs
	.a.neg = .neg != .neg
	return .norm()
}