// 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 unsigned multi-precision integers (natural
// numbers). They are the building blocks for the implementation
// of signed integers, rationals, and floating-point numbers.
//
// Caution: This implementation relies on the function "alias"
//          which assumes that (nat) slice capacities are never
//          changed (no 3-operand slice expressions). If that
//          changes, alias needs to be updated for correctness.

package big

import (
	
	
	
	
	
)

// An unsigned integer x of the form
//
//	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 a slice of length n,
// with the digits x[i] as the slice elements.
//
// A number is normalized if the slice contains no leading 0 digits.
// During arithmetic operations, denormalized values may occur but are
// always normalized before returning the final result. The normalized
// representation of 0 is the empty or nil slice (length = 0).
type nat []Word

var (
	natOne  = nat{1}
	natTwo  = nat{2}
	natFive = nat{5}
	natTen  = nat{10}
)

func ( nat) () string {
	return "0x" + string(.itoa(false, 16))
}

func ( nat) () nat {
	 := len()
	for  > 0 && [-1] == 0 {
		--
	}
	return [0:]
}

func ( nat) ( int) nat {
	if  <= cap() {
		return [:] // reuse z
	}
	if  == 1 {
		// Most nats start small and stay that way; don't over-allocate.
		return make(nat, 1)
	}
	// Choosing a good value for e has significant performance impact
	// because it increases the chance that a value can be reused.
	const  = 4 // extra capacity
	return make(nat, , +)
}

func ( nat) ( Word) nat {
	if  == 0 {
		return [:0]
	}
	 = .make(1)
	[0] = 
	return 
}

func ( nat) ( uint64) nat {
	// single-word value
	if  := Word(); uint64() ==  {
		return .setWord()
	}
	// 2-word value
	 = .make(2)
	[1] = Word( >> 32)
	[0] = Word()
	return 
}

func ( nat) ( nat) nat {
	 = .make(len())
	copy(, )
	return 
}

func ( nat) (,  nat) nat {
	 := len()
	 := len()

	switch {
	case  < :
		return .(, )
	case  == 0:
		// n == 0 because m >= n; result is 0
		return [:0]
	case  == 0:
		// result is x
		return .set()
	}
	// m > 0

	 = .make( + 1)
	 := addVV([:], [:], [:])
	if  >  {
		 = addVW([:], [:], )
	}
	[] = 

	return .norm()
}

func ( nat) (,  nat) nat {
	 := len()
	 := len()

	switch {
	case  < :
		panic("underflow")
	case  == 0:
		// n == 0 because m >= n; result is 0
		return [:0]
	case  == 0:
		// result is x
		return .set()
	}
	// m > 0

	 = .make()
	 := subVV([:], [:], [:])
	if  >  {
		 = subVW([:], [:], )
	}
	if  != 0 {
		panic("underflow")
	}

	return .norm()
}

func ( nat) ( nat) ( int) {
	 := len()
	 := len()
	if  !=  ||  == 0 {
		switch {
		case  < :
			 = -1
		case  > :
			 = 1
		}
		return
	}

	 :=  - 1
	for  > 0 && [] == [] {
		--
	}

	switch {
	case [] < []:
		 = -1
	case [] > []:
		 = 1
	}
	return
}

// montgomery computes z mod m = x*y*2**(-n*_W) mod m,
// assuming k = -1/m mod 2**_W.
// z is used for storing the result which is returned;
// z must not alias x, y or m.
// See Gueron, "Efficient Software Implementations of Modular Exponentiation".
// https://eprint.iacr.org/2011/239.pdf
// In the terminology of that paper, this is an "Almost Montgomery Multiplication":
// x and y are required to satisfy 0 <= z < 2**(n*_W) and then the result
// z is guaranteed to satisfy 0 <= z < 2**(n*_W), but it may not be < m.
func ( nat) (, ,  nat,  Word,  int) nat {
	// This code assumes x, y, m are all the same length, n.
	// (required by addMulVVW and the for loop).
	// It also assumes that x, y are already reduced mod m,
	// or else the result will not be properly reduced.
	if len() !=  || len() !=  || len() !=  {
		panic("math/big: mismatched montgomery number lengths")
	}
	 = .make( * 2)
	clear()
	var  Word
	for  := 0;  < ; ++ {
		 := []
		 := addMulVVWW([:+], [:+], , , 0)
		 := [] * 
		 := addMulVVWW([:+], [:+], , , 0)
		 :=  + 
		 :=  + 
		[+] = 
		if  <  ||  <  {
			 = 1
		} else {
			 = 0
		}
	}
	if  != 0 {
		subVV([:], [:], )
	} else {
		copy([:], [:])
	}
	return [:]
}

// alias reports whether x and y share the same base array.
//
// Note: alias assumes that the capacity of underlying arrays
// is never changed for nat values; i.e. that there are
// no 3-operand slice expressions in this code (or worse,
// reflect-based operations to the same effect).
func alias(,  nat) bool {
	return cap() > 0 && cap() > 0 && &[0:cap()][cap()-1] == &[0:cap()][cap()-1]
}

// addTo implements z += x; z must be long enough.
// (we don't use nat.add because we need z to stay the same
// slice, and we don't need to normalize z after each addition)
func addTo(,  nat) {
	if  := len();  > 0 {
		if  := addVV([:], [:], [:]);  != 0 {
			if  < len() {
				addVW([:], [:], )
			}
		}
	}
}

// mulRange computes the product of all the unsigned integers in the
// range [a, b] inclusively. If a > b (empty range), the result is 1.
// The caller may pass stk == nil to request that mulRange obtain and release one itself.
func ( nat) ( *stack, ,  uint64) nat {
	switch {
	case  == 0:
		// cut long ranges short (optimization)
		return .setUint64(0)
	case  > :
		return .setUint64(1)
	case  == :
		return .setUint64()
	case +1 == :
		return .mul(, nat(nil).setUint64(), nat(nil).setUint64())
	}

	if  == nil {
		 = getStack()
		defer .free()
	}

	 :=  + (-)/2 // avoid overflow
	return .mul(, nat(nil).(, , ), nat(nil).(, +1, ))
}

// A stack provides temporary storage for complex calculations
// such as multiplication and division.
// The stack is a simple slice of words, extended as needed
// to hold all the temporary storage for a calculation.
// In general, if a function takes a *stack, it expects a non-nil *stack.
// However, certain functions may allow passing a nil *stack instead,
// so that they can handle trivial stack-free cases without forcing the
// caller to obtain and free a stack that will be unused. These functions
// document that they accept a nil *stack in their doc comments.
type stack struct {
	w []Word
}

var stackPool sync.Pool

// getStack returns a temporary stack.
// The caller must call [stack.free] to give up use of the stack when finished.
func getStack() *stack {
	,  := stackPool.Get().(*stack)
	if  == nil {
		 = new(stack)
	}
	return 
}

// free returns the stack for use by another calculation.
func ( *stack) () {
	.w = .w[:0]
	stackPool.Put()
}

// save returns the current stack pointer.
// A future call to restore with the same value
// frees any temporaries allocated on the stack after the call to save.
func ( *stack) () int {
	return len(.w)
}

// restore restores the stack pointer to n.
// It is almost always invoked as
//
//	defer stk.restore(stk.save())
//
// which makes sure to pop any temporaries allocated in the current function
// from the stack before returning.
func ( *stack) ( int) {
	.w = .w[:]
}

// nat returns a nat of n words, allocated on the stack.
func ( *stack) ( int) nat {
	 := ( + 3) &^ 3 // round up to multiple of 4
	 := len(.w)
	.w = slices.Grow(.w, )
	.w = .w[:+]
	 := .w[ : + : +]
	if  > 0 {
		[0] = 0xfedcb
	}
	return 
}

// bitLen returns the length of x in bits.
// Unlike most methods, it works even if x is not normalized.
func ( nat) () int {
	// This function is used in cryptographic operations. It must not leak
	// anything but the Int's sign and bit size through side-channels. Any
	// changes must be reviewed by a security expert.
	if  := len() - 1;  >= 0 {
		// bits.Len uses a lookup table for the low-order bits on some
		// architectures. Neutralize any input-dependent behavior by setting all
		// bits after the first one bit.
		 := uint([])
		 |=  >> 1
		 |=  >> 2
		 |=  >> 4
		 |=  >> 8
		 |=  >> 16
		 |=  >> 16 >> 16 // ">> 32" doesn't compile on 32-bit architectures
		return *_W + bits.Len()
	}
	return 0
}

// trailingZeroBits returns the number of consecutive least significant zero
// bits of x.
func ( nat) () uint {
	if len() == 0 {
		return 0
	}
	var  uint
	for [] == 0 {
		++
	}
	// x[i] != 0
	return *_W + uint(bits.TrailingZeros(uint([])))
}

// isPow2 returns i, true when x == 2**i and 0, false otherwise.
func ( nat) () (uint, bool) {
	var  uint
	for [] == 0 {
		++
	}
	if  == uint(len())-1 && []&([]-1) == 0 {
		return *_W + uint(bits.TrailingZeros(uint([]))), true
	}
	return 0, false
}

func same(,  nat) bool {
	return len() == len() && len() > 0 && &[0] == &[0]
}

// z = x << s
func ( nat) ( nat,  uint) nat {
	if  == 0 {
		if same(, ) {
			return 
		}
		if !alias(, ) {
			return .set()
		}
	}

	 := len()
	if  == 0 {
		return [:0]
	}
	// m > 0

	 :=  + int(/_W)
	 = .make( + 1)
	if  %= _W;  == 0 {
		copy([-:], )
		[] = 0
	} else {
		[] = lshVU([-:], , )
	}
	clear([0 : -])

	return .norm()
}

// z = x >> s
func ( nat) ( nat,  uint) nat {
	if  == 0 {
		if same(, ) {
			return 
		}
		if !alias(, ) {
			return .set()
		}
	}

	 := len()
	 :=  - int(/_W)
	if  <= 0 {
		return [:0]
	}
	// n > 0

	 = .make()
	if  %= _W;  == 0 {
		copy(, [-:])
	} else {
		rshVU(, [-:], )
	}

	return .norm()
}

func ( nat) ( nat,  uint,  uint) nat {
	 := int( / _W)
	 := Word(1) << ( % _W)
	 := len()
	switch  {
	case 0:
		 = .make()
		copy(, )
		if  >=  {
			// no need to grow
			return 
		}
		[] &^= 
		return .norm()
	case 1:
		if  >=  {
			 = .make( + 1)
			clear([:])
		} else {
			 = .make()
		}
		copy(, )
		[] |= 
		// no need to normalize
		return 
	}
	panic("set bit is not 0 or 1")
}

// bit returns the value of the i'th bit, with lsb == bit 0.
func ( nat) ( uint) uint {
	 :=  / _W
	if  >= uint(len()) {
		return 0
	}
	// 0 <= j < len(x)
	return uint([] >> ( % _W) & 1)
}

// sticky returns 1 if there's a 1 bit within the
// i least significant bits, otherwise it returns 0.
func ( nat) ( uint) uint {
	 :=  / _W
	if  >= uint(len()) {
		if len() == 0 {
			return 0
		}
		return 1
	}
	// 0 <= j < len(x)
	for ,  := range [:] {
		if  != 0 {
			return 1
		}
	}
	if []<<(_W-%_W) != 0 {
		return 1
	}
	return 0
}

func ( nat) (,  nat) nat {
	 := len()
	 := len()
	if  >  {
		 = 
	}
	// m <= n

	 = .make()
	for  := 0;  < ; ++ {
		[] = [] & []
	}

	return .norm()
}

// trunc returns z = x mod 2ⁿ.
func ( nat) ( nat,  uint) nat {
	 := ( + _W - 1) / _W
	if uint(len()) <  {
		return .set()
	}
	 = .make(int())
	copy(, )
	if %_W != 0 {
		[len()-1] &= 1<<(%_W) - 1
	}
	return .norm()
}

func ( nat) (,  nat) nat {
	 := len()
	 := len()
	if  >  {
		 = 
	}
	// m >= n

	 = .make()
	for  := 0;  < ; ++ {
		[] = [] &^ []
	}
	copy([:], [:])

	return .norm()
}

func ( nat) (,  nat) nat {
	 := len()
	 := len()
	 := 
	if  <  {
		,  = , 
		 = 
	}
	// m >= n

	 = .make()
	for  := 0;  < ; ++ {
		[] = [] | []
	}
	copy([:], [:])

	return .norm()
}

func ( nat) (,  nat) nat {
	 := len()
	 := len()
	 := 
	if  <  {
		,  = , 
		 = 
	}
	// m >= n

	 = .make()
	for  := 0;  < ; ++ {
		[] = [] ^ []
	}
	copy([:], [:])

	return .norm()
}

// random creates a random integer in [0..limit), using the space in z if
// possible. n is the bit length of limit.
func ( nat) ( *rand.Rand,  nat,  int) nat {
	if alias(, ) {
		 = nil // z is an alias for limit - cannot reuse
	}
	 = .make(len())

	 := uint( % _W)
	if  == 0 {
		 = _W
	}
	 := Word((1 << ) - 1)

	for {
		switch _W {
		case 32:
			for  := range  {
				[] = Word(.Uint32())
			}
		case 64:
			for  := range  {
				[] = Word(.Uint32()) | Word(.Uint32())<<32
			}
		default:
			panic("unknown word size")
		}
		[len()-1] &= 
		if .cmp() < 0 {
			break
		}
	}

	return .norm()
}

// If m != 0 (i.e., len(m) != 0), expNN sets z to x**y mod m;
// otherwise it sets z to x**y. The result is the value of z.
// The caller may pass stk == nil to request that expNN obtain and release one itself.
func ( nat) ( *stack, , ,  nat,  bool) nat {
	if alias(, ) || alias(, ) {
		// We cannot allow in-place modification of x or y.
		 = nil
	}

	// x**y mod 1 == 0
	if len() == 1 && [0] == 1 {
		return .setWord(0)
	}
	// m == 0 || m > 1

	// x**0 == 1
	if len() == 0 {
		return .setWord(1)
	}
	// y > 0

	// 0**y = 0
	if len() == 0 {
		return .setWord(0)
	}
	// x > 0

	// 1**y = 1
	if len() == 1 && [0] == 1 {
		return .setWord(1)
	}
	// x > 1

	// x**1 == x
	if len() == 1 && [0] == 1 && len() == 0 {
		return .set()
	}
	if  == nil {
		 = getStack()
		defer .free()
	}
	if len() == 1 && [0] == 1 { // len(m) > 0
		return .rem(, , )
	}

	// y > 1

	if len() != 0 {
		// We likely end up being as long as the modulus.
		 = .make(len())

		// If the exponent is large, we use the Montgomery method for odd values,
		// and a 4-bit, windowed exponentiation for powers of two,
		// and a CRT-decomposed Montgomery method for the remaining values
		// (even values times non-trivial odd values, which decompose into one
		// instance of each of the first two cases).
		if len() > 1 && ! {
			if [0]&1 == 1 {
				return .expNNMontgomery(, , , )
			}
			if ,  := .isPow2();  {
				return .expNNWindowed(, , , )
			}
			return .expNNMontgomeryEven(, , , )
		}
	}

	 = .set()
	 := [len()-1] // v > 0 because y is normalized and y > 0
	 := nlz() + 1
	 <<= 
	var  nat

	const  = 1 << (_W - 1)

	// We walk through the bits of the exponent one by one. Each time we
	// see a bit, we square, thus doubling the power. If the bit is a one,
	// we also multiply by x, thus adding one to the power.

	 := _W - int()
	// zz and r are used to avoid allocating in mul and div as
	// otherwise the arguments would alias.
	var ,  nat
	for  := 0;  < ; ++ {
		 = .sqr(, )
		,  = , 

		if & != 0 {
			 = .mul(, , )
			,  = , 
		}

		if len() != 0 {
			,  = .div(, , , )
			, , ,  = , , , 
		}

		 <<= 1
	}

	for  := len() - 2;  >= 0; -- {
		 = []

		for  := 0;  < _W; ++ {
			 = .sqr(, )
			,  = , 

			if & != 0 {
				 = .mul(, , )
				,  = , 
			}

			if len() != 0 {
				,  = .div(, , , )
				, , ,  = , , , 
			}

			 <<= 1
		}
	}

	return .norm()
}

// expNNMontgomeryEven calculates x**y mod m where m = m1 × m2 for m1 = 2ⁿ and m2 odd.
// It uses two recursive calls to expNN for x**y mod m1 and x**y mod m2
// and then uses the Chinese Remainder Theorem to combine the results.
// The recursive call using m1 will use expNNWindowed,
// while the recursive call using m2 will use expNNMontgomery.
// For more details, see Ç. K. Koç, “Montgomery Reduction with Even Modulus”,
// IEE Proceedings: Computers and Digital Techniques, 141(5) 314-316, September 1994.
// http://www.people.vcu.edu/~jwang3/CMSC691/j34monex.pdf
func ( nat) ( *stack, , ,  nat) nat {
	// Split m = m₁ × m₂ where m₁ = 2ⁿ
	 := .trailingZeroBits()
	 := nat(nil).lsh(natOne, )
	 := nat(nil).rsh(, )

	// We want z = x**y mod m.
	// z₁ = x**y mod m1 = (x**y mod m) mod m1 = z mod m1
	// z₂ = x**y mod m2 = (x**y mod m) mod m2 = z mod m2
	// (We are using the math/big convention for names here,
	// where the computation is z = x**y mod m, so its parts are z1 and z2.
	// The paper is computing x = a**e mod n; it refers to these as x2 and z1.)
	 := nat(nil).expNN(, , , , false)
	 := nat(nil).expNN(, , , , false)

	// Reconstruct z from z₁, z₂ using CRT, using algorithm from paper,
	// which uses only a single modInverse (and an easy one at that).
	//	p = (z₁ - z₂) × m₂⁻¹ (mod m₁)
	//	z = z₂ + p × m₂
	// The final addition is in range because:
	//	z = z₂ + p × m₂
	//	  ≤ z₂ + (m₁-1) × m₂
	//	  < m₂ + (m₁-1) × m₂
	//	  = m₁ × m₂
	//	  = m.
	 = .set()

	// Compute (z₁ - z₂) mod m1 [m1 == 2**n] into z1.
	 = .subMod2N(, , )

	// Reuse z2 for p = (z₁ - z₂) [in z1] * m2⁻¹ (mod m₁ [= 2ⁿ]).
	 := nat(nil).modInverse(, )
	 = .mul(, , )
	 = .trunc(, )

	// Reuse z1 for p * m2.
	 = .add(, .mul(, , ))

	return 
}

// expNNWindowed calculates x**y mod m using a fixed, 4-bit window,
// where m = 2**logM.
func ( nat) ( *stack, ,  nat,  uint) nat {
	if len() <= 1 {
		panic("big: misuse of expNNWindowed")
	}
	if [0]&1 == 0 {
		// len(y) > 1, so y  > logM.
		// x is even, so x**y is a multiple of 2**y which is a multiple of 2**logM.
		return .setWord(0)
	}
	if  == 1 {
		return .setWord(1)
	}

	// zz is used to avoid allocating in mul as otherwise
	// the arguments would alias.
	defer .restore(.save())
	 := int(( + _W - 1) / _W)
	 := .nat()

	const  = 4
	// powers[i] contains x^i.
	var  [1 << ]nat
	for  := range  {
		[] = .nat()
	}
	[0] = [0].set(natOne)
	[1] = [1].trunc(, )
	for  := 2;  < 1<<;  += 2 {
		, ,  := &[/2], &[], &[+1]
		* = .sqr(, *)
		* = .trunc(*, )
		* = .mul(, *, )
		* = .trunc(*, )
	}

	// Because phi(2**logM) = 2**(logM-1), x**(2**(logM-1)) = 1,
	// so we can compute x**(y mod 2**(logM-1)) instead of x**y.
	// That is, we can throw away all but the bottom logM-1 bits of y.
	// Instead of allocating a new y, we start reading y at the right word
	// and truncate it appropriately at the start of the loop.
	 := len() - 1
	 := int(( - 2) / _W) // -2 because the top word of N bits is the (N-1)/W'th word.
	 := ^Word(0)
	if  := ( - 1) & (_W - 1);  != 0 {
		 = (1 << ) - 1
	}
	if  >  {
		 = 
	}
	 := false
	 = .setWord(1)
	for ;  >= 0; -- {
		 := []
		if  ==  {
			 &= 
		}
		for  := 0;  < _W;  +=  {
			if  {
				// Account for use of 4 bits in previous iteration.
				// Unrolled loop for significant performance
				// gain. Use go test -bench=".*" in crypto/rsa
				// to check performance before making changes.
				 = .sqr(, )
				,  = , 
				 = .trunc(, )

				 = .sqr(, )
				,  = , 
				 = .trunc(, )

				 = .sqr(, )
				,  = , 
				 = .trunc(, )

				 = .sqr(, )
				,  = , 
				 = .trunc(, )
			}

			 = .mul(, , [>>(_W-)])
			,  = , 
			 = .trunc(, )

			 <<= 
			 = true
		}
	}

	return .norm()
}

// expNNMontgomery calculates x**y mod m using a fixed, 4-bit window.
// Uses Montgomery representation.
func ( nat) ( *stack, , ,  nat) nat {
	 := len()

	// We want the lengths of x and m to be equal.
	// It is OK if x >= m as long as len(x) == len(m).
	if len() >  {
		_,  = nat(nil).div(, nil, , )
		// Note: now len(x) <= numWords, not guaranteed ==.
	}
	if len() <  {
		 := make(nat, )
		copy(, )
		 = 
	}

	// Ideally the precomputations would be performed outside, and reused
	// k0 = -m**-1 mod 2**_W. Algorithm from: Dumas, J.G. "On Newton–Raphson
	// Iteration for Multiplicative Inverses Modulo Prime Powers".
	 := 2 - [0]
	 := [0] - 1
	for  := 1;  < _W;  <<= 1 {
		 *= 
		 *= ( + 1)
	}
	 = -

	// RR = 2**(2*_W*len(m)) mod m
	 := nat(nil).setWord(1)
	 := nat(nil).lsh(, uint(2**_W))
	_,  = nat(nil).div(, , , )
	if len() <  {
		 = .make()
		copy(, )
		 = 
	}
	// one = 1, with equal length to that of m
	 := make(nat, )
	[0] = 1

	const  = 4
	// powers[i] contains x^i
	var  [1 << ]nat
	[0] = [0].montgomery(, , , , )
	[1] = [1].montgomery(, , , , )
	for  := 2;  < 1<<; ++ {
		[] = [].montgomery([-1], [1], , , )
	}

	// initialize z = 1 (Montgomery 1)
	 = .make()
	copy(, [0])

	 = .make()

	// same windowed exponent, but with Montgomery multiplications
	for  := len() - 1;  >= 0; -- {
		 := []
		for  := 0;  < _W;  +=  {
			if  != len()-1 ||  != 0 {
				 = .montgomery(, , , , )
				 = .montgomery(, , , , )
				 = .montgomery(, , , , )
				 = .montgomery(, , , , )
			}
			 = .montgomery(, [>>(_W-)], , , )
			,  = , 
			 <<= 
		}
	}
	// convert to regular number
	 = .montgomery(, , , , )

	// One last reduction, just in case.
	// See golang.org/issue/13907.
	if .cmp() >= 0 {
		// Common case is m has high bit set; in that case,
		// since zz is the same length as m, there can be just
		// one multiple of m to remove. Just subtract.
		// We think that the subtract should be sufficient in general,
		// so do that unconditionally, but double-check,
		// in case our beliefs are wrong.
		// The div is not expected to be reached.
		 = .sub(, )
		if .cmp() >= 0 {
			_,  = nat(nil).div(, nil, , )
		}
	}

	return .norm()
}

// bytes writes the value of z into buf using big-endian encoding.
// The value of z is encoded in the slice buf[i:]. If the value of z
// cannot be represented in buf, bytes panics. The number i of unused
// bytes at the beginning of buf is returned as result.
func ( nat) ( []byte) ( int) {
	// This function is used in cryptographic operations. It must not leak
	// anything but the Int's sign and bit size through side-channels. Any
	// changes must be reviewed by a security expert.
	 = len()
	for ,  := range  {
		for  := 0;  < _S; ++ {
			--
			if  >= 0 {
				[] = byte()
			} else if byte() != 0 {
				panic("math/big: buffer too small to fit value")
			}
			 >>= 8
		}
	}

	if  < 0 {
		 = 0
	}
	for  < len() && [] == 0 {
		++
	}

	return
}

// bigEndianWord returns the contents of buf interpreted as a big-endian encoded Word value.
func bigEndianWord( []byte) Word {
	if _W == 64 {
		return Word(byteorder.BEUint64())
	}
	return Word(byteorder.BEUint32())
}

// setBytes interprets buf as the bytes of a big-endian unsigned
// integer, sets z to that value, and returns z.
func ( nat) ( []byte) nat {
	 = .make((len() + _S - 1) / _S)

	 := len()
	for  := 0;  >= _S; ++ {
		[] = bigEndianWord([-_S : ])
		 -= _S
	}
	if  > 0 {
		var  Word
		for  := uint(0);  > 0;  += 8 {
			 |= Word([-1]) << 
			--
		}
		[len()-1] = 
	}

	return .norm()
}

// sqrt sets z = ⌊√x⌋
// The caller may pass stk == nil to request that sqrt obtain and release one itself.
func ( nat) ( *stack,  nat) nat {
	if .cmp(natOne) <= 0 {
		return .set()
	}
	if alias(, ) {
		 = nil
	}

	if  == nil {
		 = getStack()
		defer .free()
	}

	// Start with value known to be too large and repeat "z = ⌊(z + ⌊x/z⌋)/2⌋" until it stops getting smaller.
	// See Brent and Zimmermann, Modern Computer Arithmetic, Algorithm 1.13 (SqrtInt).
	// https://members.loria.fr/PZimmermann/mca/pub226.html
	// If x is one less than a perfect square, the sequence oscillates between the correct z and z+1;
	// otherwise it converges to the correct z and stays there.
	var ,  nat
	 = 
	 = .setUint64(1)
	 = .lsh(, uint(.bitLen()+1)/2) // must be ≥ √x
	for  := 0; ; ++ {
		, _ = .div(, nil, , )
		 = .add(, )
		 = .rsh(, 1)
		if .cmp() >= 0 {
			// z1 is answer.
			// Figure out whether z1 or z2 is currently aliased to z by looking at loop count.
			if &1 == 0 {
				return 
			}
			return .set()
		}
		,  = , 
	}
}

// subMod2N returns z = (x - y) mod 2ⁿ.
func ( nat) (,  nat,  uint) nat {
	if uint(.bitLen()) >  {
		if alias(, ) {
			// ok to overwrite x in place
			 = .trunc(, )
		} else {
			 = nat(nil).trunc(, )
		}
	}
	if uint(.bitLen()) >  {
		if alias(, ) {
			// ok to overwrite y in place
			 = .trunc(, )
		} else {
			 = nat(nil).trunc(, )
		}
	}
	if .cmp() >= 0 {
		return .sub(, )
	}
	// x - y < 0; x - y mod 2ⁿ = x - y + 2ⁿ = 2ⁿ - (y - x) = 1 + 2ⁿ-1 - (y - x) = 1 + ^(y - x).
	 = .sub(, )
	for uint(len())*_W <  {
		 = append(, 0)
	}
	for  := range  {
		[] = ^[]
	}
	 = .trunc(, )
	return .add(, natOne)
}