// Use of this source code is governed by a BSD-style

// 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) () {
for  := range  {
[] = 0
}
}

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)
if  >  {
}
[] =

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([0:], , )
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
}

func ( nat) ( nat, ,  Word) nat {
:= len()
if  == 0 ||  == 0 {
return .setWord() // result is r
}
// m > 0

= .make( + 1)
[] = mulAddVWW([0:], , , )

return .norm()
}

// basicMul multiplies x and y and leaves the result in z.
// The (non-normalized) result is placed in z[0 : len(x) + len(y)].
func basicMul(, ,  nat) {
[0 : len()+len()].clear() // initialize z
for ,  := range  {
if  != 0 {
}
}
}

// 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;  < ; ++ {
:= []
:= [] *
:=  +
:=  +
[+] =
if  <  ||  <  {
= 1
} else {
= 0
}
}
if  != 0 {
subVV([:], [:], )
} else {
copy([:], [:])
}
return [:]
}

// Fast version of z[0:n+n>>1].add(z[0:n+n>>1], x[0:n]) w/o bounds checks.
// Factored out for readability - do not use outside karatsuba.
if  := addVV([0:], , );  != 0 {
}
}

// Like karatsubaAdd, but does subtract.
func karatsubaSub(,  nat,  int) {
if  := subVV([0:], , );  != 0 {
subVW([:+>>1], [:], )
}
}

// Operands that are shorter than karatsubaThreshold are multiplied using
// "grade school" multiplication; for longer operands the Karatsuba algorithm
// is used.
var karatsubaThreshold = 40 // computed by calibrate_test.go

// karatsuba multiplies x and y and leaves the result in z.
// Both x and y must have the same length n and n must be a
// power of 2. The result vector z must have len(z) >= 6*n.
// The (non-normalized) result is placed in z[0 : 2*n].
func karatsuba(, ,  nat) {
:= len()

// Switch to basic multiplication if numbers are odd or small.
// (n is always even if karatsubaThreshold is even, but be
// conservative)
if &1 != 0 ||  < karatsubaThreshold ||  < 2 {
basicMul(, , )
return
}
// n&1 == 0 && n >= karatsubaThreshold && n >= 2

// Karatsuba multiplication is based on the observation that
// for two numbers x and y with:
//
//   x = x1*b + x0
//   y = y1*b + y0
//
// the product x*y can be obtained with 3 products z2, z1, z0
//
//   x*y = x1*y1*b*b + (x1*y0 + x0*y1)*b + x0*y0
//       =    z2*b*b +              z1*b +    z0
//
// with:
//
//   xd = x1 - x0
//   yd = y0 - y1
//
//   z1 =      xd*yd                    + z2 + z0
//      = (x1-x0)*(y0 - y1)             + z2 + z0
//      = x1*y0 - x1*y1 - x0*y0 + x0*y1 + z2 + z0
//      = x1*y0 -    z2 -    z0 + x0*y1 + z2 + z0
//      = x1*y0                 + x0*y1

// split x, y into "digits"
:=  >> 1              // n2 >= 1
,  := [:], [0:] // x = x1*b + y0
,  := [:], [0:] // y = y1*b + y0

// z is used for the result and temporary storage:
//
//   6*n     5*n     4*n     3*n     2*n     1*n     0*n
// z = [z2 copy|z0 copy| xd*yd | yd:xd | x1*y1 | x0*y0 ]
//
// For each recursive call of karatsuba, an unused slice of
// z is passed in that has (at least) half the length of the
// caller's z.

// compute z0 and z2 with the result "in place" in z
(, , )     // z0 = x0*y0
([:], , ) // z2 = x1*y1

// compute xd (or the negative value if underflow occurs)
:= 1 // sign of product xd*yd
:= [2* : 2*+]
if subVV(, , ) != 0 { // x1-x0
= -
subVV(, , ) // x0-x1
}

// compute yd (or the negative value if underflow occurs)
:= [2*+ : 3*]
if subVV(, , ) != 0 { // y0-y1
= -
subVV(, , ) // y1-y0
}

// p = (x1-x0)*(y0-y1) == x1*y0 - x1*y1 - x0*y0 + x0*y1 for s > 0
// p = (x0-x1)*(y0-y1) == x0*y0 - x0*y1 - x1*y0 + x1*y1 for s < 0
:= [*3:]
(, , )

// save original z2:z0
// (ok to use upper half of z since we're done recurring)
:= [*4:]
copy(, [:*2])

// add up all partial products
//
//   2*n     n     0
// z = [ z2  | z0  ]
//   +    [ z0  ]
//   +    [ z2  ]
//   +    [  p  ]
//
if  > 0 {
} else {
karatsubaSub([:], , )
}
}

// 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]
}

// addAt implements z += x<<(_W*i); 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)
if  := len();  > 0 {
if  := addVV([:+], [:], );  != 0 {
:=  +
if  < len() {
}
}
}
}

// karatsubaLen computes an approximation to the maximum k <= n such that
// k = p<<i for a number p <= threshold and an i >= 0. Thus, the
// result is the largest number that can be divided repeatedly by 2 before
// becoming about the value of threshold.
func karatsubaLen(,  int) int {
:= uint(0)
for  >  {
>>= 1
++
}
return  <<
}

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

switch {
case  < :
return .(, )
case  == 0 ||  == 0:
return [:0]
case  == 1:
}
// m >= n > 1

// determine if z can be reused
if alias(, ) || alias(, ) {
= nil // z is an alias for x or y - cannot reuse
}

// use basic multiplication if the numbers are small
if  < karatsubaThreshold {
= .make( + )
basicMul(, , )
return .norm()
}
// m >= n && n >= karatsubaThreshold && n >= 2

// determine Karatsuba length k such that
//
//   x = xh*b + x0  (0 <= x0 < b)
//   y = yh*b + y0  (0 <= y0 < b)
//   b = 1<<(_W*k)  ("base" of digits xi, yi)
//
:= karatsubaLen(, karatsubaThreshold)
// k <= n

// multiply x0 and y0 via Karatsuba
:= [0:]              // x0 is not normalized
:= [0:]              // y0 is not normalized
= .make(max(6*, +)) // enough space for karatsuba of x0*y0 and full result of x*y
karatsuba(, , )
= [0 : +]  // z has final length but may be incomplete
[2*:].clear() // upper portion of z is garbage (and 2*k <= m+n since k <= n <= m)

// If xh != 0 or yh != 0, add the missing terms to z. For
//
//   xh = xi*b^i + ... + x2*b^2 + x1*b (0 <= xi < b)
//   yh =                         y1*b (0 <= y1 < b)
//
// the missing terms are
//
//   x0*y1*b and xi*y0*b^i, xi*y1*b^(i+1) for i > 0
//
// since all the yi for i > 1 are 0 by choice of k: If any of them
// were > 0, then yh >= b^2 and thus y >= b^2. Then k' = k*2 would
//
if  <  ||  !=  {
:= getNat(3 * )
:= *

:= .norm()
:= [:]       // y1 is normalized because y is
= .(, ) // update t so we don't lose t's underlying array

:= .norm()
for  := ;  < len();  +=  {
:= [:]
if len() >  {
= [:]
}
= .norm()
= .(, )
= .(, )
}

putNat()
}

return .norm()
}

// basicSqr sets z = x*x and is asymptotically faster than basicMul
// by about a factor of 2, but slower for small arguments due to overhead.
// Requirements: len(x) > 0, len(z) == 2*len(x)
// The (non-normalized) result is placed in z.
func basicSqr(,  nat) {
:= len()
:= getNat(2 * )
:= * // temporary variable to hold the products
.clear()
[1], [0] = mulWW([0], [0]) // the initial square
for  := 1;  < ; ++ {
:= []
// z collects the squares x[i] * x[i]
[2*+1], [2*] = mulWW(, )
// t collects the products x[i] * x[j] where j < i
}
[2*-1] = shlVU([1:2*-1], [1:2*-1], 1) // double the j < i products
addVV(, , )                              // combine the result
putNat()
}

// karatsubaSqr squares x and leaves the result in z.
// len(x) must be a power of 2 and len(z) >= 6*len(x).
// The (non-normalized) result is placed in z[0 : 2*len(x)].
//
// The algorithm and the layout of z are the same as for karatsuba.
func karatsubaSqr(,  nat) {
:= len()

if &1 != 0 ||  < karatsubaSqrThreshold ||  < 2 {
basicSqr([:2*], )
return
}

:=  >> 1
,  := [:], [0:]

(, )
([:], )

// s = sign(xd*yd) == -1 for xd != 0; s == 1 for xd == 0
:= [2* : 2*+]
if subVV(, , ) != 0 {
subVV(, , )
}

:= [*3:]
(, )

:= [*4:]
copy(, [:*2])

karatsubaSub([:], , ) // s == -1 for p != 0; s == 1 for p == 0
}

// Operands that are shorter than basicSqrThreshold are squared using
// "grade school" multiplication; for operands longer than karatsubaSqrThreshold
// we use the Karatsuba algorithm optimized for x == y.
var basicSqrThreshold = 20      // computed by calibrate_test.go
var karatsubaSqrThreshold = 260 // computed by calibrate_test.go

// z = x*x
func ( nat) ( nat) nat {
:= len()
switch {
case  == 0:
return [:0]
case  == 1:
:= [0]
= .make(2)
[1], [0] = mulWW(, )
return .norm()
}

if alias(, ) {
= nil // z is an alias for x - cannot reuse
}

if  < basicSqrThreshold {
= .make(2 * )
basicMul(, , )
return .norm()
}
if  < karatsubaSqrThreshold {
= .make(2 * )
basicSqr(, )
return .norm()
}

// Use Karatsuba multiplication optimized for x == y.
// The algorithm and layout of z are the same as for mul.

// z = (x1*b + x0)^2 = x1^2*b^2 + 2*x1*x0*b + x0^2

:= karatsubaLen(, karatsubaSqrThreshold)

:= [0:]
= .make(max(6*, 2*))
karatsubaSqr(, ) // z = x0^2
= [0 : 2*]
[2*:].clear()

if  <  {
:= getNat(2 * )
:= *
:= .norm()
:= [:]
= .mul(, )
addAt(, , ) // z = 2*x1*x0*b + x0^2
= .()
addAt(, , 2*) // z = x1^2*b^2 + 2*x1*x0*b + x0^2
putNat()
}

return .norm()
}

// mulRange computes the product of all the unsigned integers in the
// range [a, b] inclusively. If a > b (empty range), the result is 1.
func ( nat) (,  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())
}
:=  + (-)/2 // avoid overflow
return .mul(nat(nil).(, ), nat(nil).(+1, ))
}

// getNat returns a *nat of len n. The contents may not be zero.
// The pool holds *nat to avoid allocation when converting to interface{}.
func getNat( int) *nat {
var  *nat
if  := natPool.Get();  != nil {
= .(*nat)
}
if  == nil {
= new(nat)
}
* = .make()
if  > 0 {
(*)[0] = 0xfedcb // break code expecting zero
}
return
}

func putNat( *nat) {
natPool.Put()
}

var natPool sync.Pool

// 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)
[] = shlVU([-:], , %_W)
[0 : -].clear()

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()
shrVU(, [-:], %_W)

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.
func ( nat) (, ,  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 {
if len() != 0 {
return .rem(, )
}
return .set()
}
// 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) (, ,  nat) nat {
// Split m = m₁ × m₂ where m₁ = 2ⁿ
:= .trailingZeroBits()
:= nat(nil).shl(natOne, )
:= nat(nil).shr(, )

// 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.

return
}

// expNNWindowed calculates x**y mod m using a fixed, 4-bit window,
// where m = 2**logM.
func ( nat) (,  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.
:= int(( + _W - 1) / _W)
:= getNat()
:= *

const  = 4
// powers[i] contains x^i.
var  [1 << ]*nat
for  := range  {
[] = getNat()
}
*[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
}
}

* =
putNat()
for  := range  {
putNat([])
}

return .norm()
}

// expNNMontgomery calculates x**y mod m using a fixed, 4-bit window.
// Uses Montgomery representation.
func ( nat) (, ,  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).shl(, 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(binary.BigEndian.Uint64())
}
return Word(binary.BigEndian.Uint32())
}

// 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⌋
func ( nat) ( nat) nat {
if .cmp(natOne) <= 0 {
return .set()
}
if alias(, ) {
= nil
}

// 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)
= .shl(, uint(.bitLen()+1)/2) // must be ≥ √x
for  := 0; ; ++ {
, _ = .div(nil, , )
= .shr(, 1)
if .cmp() >= 0 {
// 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(, )