element.go

// Copyright 2020 ConsenSys Software Inc.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

// Code generated by consensys/gnark-crypto DO NOT EDIT

package fp

// /!\ WARNING /!\
// this code has not been audited and is provided as-is. In particular,
// there is no security guarantees such as constant time implementation
// or side-channel attack resistance
// /!\ WARNING /!\

import (
	"crypto/rand"
	"encoding/binary"
	"errors"
	"io"
	"math/big"
	"math/bits"
	"reflect"
	"strconv"
	"strings"
	"sync"
)

// Element represents a field element stored on 10 words (uint64)
// Element are assumed to be in Montgomery form in all methods
// field modulus q =
//
// 20494478644167774678813387386538961497669590920908778075528754551012016751717791778743535050360001387419576570244406805463255765034468441182772056330021723098661967429339971741066259394985997
type Element [10]uint64

// Limbs number of 64 bits words needed to represent Element
const Limbs = 10

// Bits number bits needed to represent Element
const Bits = 633

// Bytes number bytes needed to represent Element
const Bytes = Limbs * 8

// field modulus stored as big.Int
var _modulus big.Int

// Modulus returns q as a big.Int
// q =
//
// 20494478644167774678813387386538961497669590920908778075528754551012016751717791778743535050360001387419576570244406805463255765034468441182772056330021723098661967429339971741066259394985997
func Modulus() *big.Int {
	return new(big.Int).Set(&_modulus)
}

// q (modulus)
const qElementWord0 uint64 = 15512955586897510413
const qElementWord1 uint64 = 4410884215886313276
const qElementWord2 uint64 = 15543556715411259941
const qElementWord3 uint64 = 9083347379620258823
const qElementWord4 uint64 = 13320134076191308873
const qElementWord5 uint64 = 9318693926755804304
const qElementWord6 uint64 = 5645674015335635503
const qElementWord7 uint64 = 12176845843281334983
const qElementWord8 uint64 = 18165857675053050549
const qElementWord9 uint64 = 82862755739295587

var qElement = Element{
	qElementWord0,
	qElementWord1,
	qElementWord2,
	qElementWord3,
	qElementWord4,
	qElementWord5,
	qElementWord6,
	qElementWord7,
	qElementWord8,
	qElementWord9,
}

// Used for Montgomery reduction. (qInvNeg) q + r'.r = 1, i.e., qInvNeg = - q⁻¹ mod r
const qInvNegLsw uint64 = 13046692460116554043

// rSquare
var rSquare = Element{
	7358459907925294924,
	14414180951914241931,
	16619482658146888203,
	760736596725344926,
	12753071240931896792,
	13425190760400245818,
	12591714441439252728,
	15325516497554583360,
	5301152003049442834,
	35368377961363834,
}

var bigIntPool = sync.Pool{
	New: func() interface{} {
		return new(big.Int)
	},
}

func init() {
	_modulus.SetString("20494478644167774678813387386538961497669590920908778075528754551012016751717791778743535050360001387419576570244406805463255765034468441182772056330021723098661967429339971741066259394985997", 10)
}

// NewElement returns a new Element from a uint64 value
//
// it is equivalent to
// 		var v NewElement
// 		v.SetUint64(...)
func NewElement(v uint64) Element {
	z := Element{v}
	z.Mul(&z, &rSquare)
	return z
}

// SetUint64 sets z to v and returns z
func (z *Element) SetUint64(v uint64) *Element {
	//  sets z LSB to v (non-Montgomery form) and convert z to Montgomery form
	*z = Element{v}
	return z.Mul(z, &rSquare) // z.ToMont()
}

// SetInt64 sets z to v and returns z
func (z *Element) SetInt64(v int64) *Element {

	// absolute value of v
	m := v >> 63
	z.SetUint64(uint64((v ^ m) - m))

	if m != 0 {
		// v is negative
		z.Neg(z)
	}

	return z
}

// Set z = x
func (z *Element) Set(x *Element) *Element {
	z[0] = x[0]
	z[1] = x[1]
	z[2] = x[2]
	z[3] = x[3]
	z[4] = x[4]
	z[5] = x[5]
	z[6] = x[6]
	z[7] = x[7]
	z[8] = x[8]
	z[9] = x[9]
	return z
}

// SetInterface converts provided interface into Element
// returns an error if provided type is not supported
// supported types: Element, *Element, uint64, int, string (interpreted as base10 integer),
// *big.Int, big.Int, []byte
func (z *Element) SetInterface(i1 interface{}) (*Element, error) {
	switch c1 := i1.(type) {
	case Element:
		return z.Set(&c1), nil
	case *Element:
		if c1 == nil {
			return nil, errors.New("can't set fp.Element with <nil>")
		}
		return z.Set(c1), nil
	case uint8:
		return z.SetUint64(uint64(c1)), nil
	case uint16:
		return z.SetUint64(uint64(c1)), nil
	case uint32:
		return z.SetUint64(uint64(c1)), nil
	case uint:
		return z.SetUint64(uint64(c1)), nil
	case uint64:
		return z.SetUint64(c1), nil
	case int8:
		return z.SetInt64(int64(c1)), nil
	case int16:
		return z.SetInt64(int64(c1)), nil
	case int32:
		return z.SetInt64(int64(c1)), nil
	case int64:
		return z.SetInt64(c1), nil
	case int:
		return z.SetInt64(int64(c1)), nil
	case string:
		return z.SetString(c1), nil
	case *big.Int:
		if c1 == nil {
			return nil, errors.New("can't set fp.Element with <nil>")
		}
		return z.SetBigInt(c1), nil
	case big.Int:
		return z.SetBigInt(&c1), nil
	case []byte:
		return z.SetBytes(c1), nil
	default:
		return nil, errors.New("can't set fp.Element from type " + reflect.TypeOf(i1).String())
	}
}

// SetZero z = 0
func (z *Element) SetZero() *Element {
	z[0] = 0
	z[1] = 0
	z[2] = 0
	z[3] = 0
	z[4] = 0
	z[5] = 0
	z[6] = 0
	z[7] = 0
	z[8] = 0
	z[9] = 0
	return z
}

// SetOne z = 1 (in Montgomery form)
func (z *Element) SetOne() *Element {
	z[0] = 5665001492438840506
	z[1] = 16907884053554239805
	z[2] = 17318295036095996852
	z[3] = 12638729832353218866
	z[4] = 12856030952767240260
	z[5] = 15732028589390776959
	z[6] = 1038965607738428109
	z[7] = 8411601626847721258
	z[8] = 7016548280614581879
	z[9] = 51212299585931083
	return z
}

// Div z = x*y^-1 mod q
func (z *Element) Div(x, y *Element) *Element {
	var yInv Element
	yInv.Inverse(y)
	z.Mul(x, &yInv)
	return z
}

// Bit returns the i'th bit, with lsb == bit 0.
// It is the responsibility of the caller to convert from Montgomery to Regular form if needed
func (z *Element) Bit(i uint64) uint64 {
	j := i / 64
	if j >= 10 {
		return 0
	}
	return uint64(z[j] >> (i % 64) & 1)
}

// Equal returns z == x; constant-time
func (z *Element) Equal(x *Element) bool {
	return z.NotEqual(x) == 0
}

// NotEqual returns 0 if and only if z == x; constant-time
func (z *Element) NotEqual(x *Element) uint64 {
	return (z[9] ^ x[9]) | (z[8] ^ x[8]) | (z[7] ^ x[7]) | (z[6] ^ x[6]) | (z[5] ^ x[5]) | (z[4] ^ x[4]) | (z[3] ^ x[3]) | (z[2] ^ x[2]) | (z[1] ^ x[1]) | (z[0] ^ x[0])
}

// IsZero returns z == 0
func (z *Element) IsZero() bool {
	return (z[9] | z[8] | z[7] | z[6] | z[5] | z[4] | z[3] | z[2] | z[1] | z[0]) == 0
}

// IsUint64 reports whether z can be represented as an uint64.
func (z *Element) IsUint64() bool {
	return (z[9] | z[8] | z[7] | z[6] | z[5] | z[4] | z[3] | z[2] | z[1]) == 0
}

// Cmp compares (lexicographic order) z and x and returns:
//
//   -1 if z <  x
//    0 if z == x
//   +1 if z >  x
//
func (z *Element) Cmp(x *Element) int {
	_z := *z
	_x := *x
	_z.FromMont()
	_x.FromMont()
	if _z[9] > _x[9] {
		return 1
	} else if _z[9] < _x[9] {
		return -1
	}
	if _z[8] > _x[8] {
		return 1
	} else if _z[8] < _x[8] {
		return -1
	}
	if _z[7] > _x[7] {
		return 1
	} else if _z[7] < _x[7] {
		return -1
	}
	if _z[6] > _x[6] {
		return 1
	} else if _z[6] < _x[6] {
		return -1
	}
	if _z[5] > _x[5] {
		return 1
	} else if _z[5] < _x[5] {
		return -1
	}
	if _z[4] > _x[4] {
		return 1
	} else if _z[4] < _x[4] {
		return -1
	}
	if _z[3] > _x[3] {
		return 1
	} else if _z[3] < _x[3] {
		return -1
	}
	if _z[2] > _x[2] {
		return 1
	} else if _z[2] < _x[2] {
		return -1
	}
	if _z[1] > _x[1] {
		return 1
	} else if _z[1] < _x[1] {
		return -1
	}
	if _z[0] > _x[0] {
		return 1
	} else if _z[0] < _x[0] {
		return -1
	}
	return 0
}

// LexicographicallyLargest returns true if this element is strictly lexicographically
// larger than its negation, false otherwise
func (z *Element) LexicographicallyLargest() bool {
	// adapted from github.com/zkcrypto/bls12_381
	// we check if the element is larger than (q-1) / 2
	// if z - (((q -1) / 2) + 1) have no underflow, then z > (q-1) / 2

	_z := *z
	_z.FromMont()

	var b uint64
	_, b = bits.Sub64(_z[0], 7756477793448755207, 0)
	_, b = bits.Sub64(_z[1], 11428814144797932446, b)
	_, b = bits.Sub64(_z[2], 16995150394560405778, b)
	_, b = bits.Sub64(_z[3], 13765045726664905219, b)
	_, b = bits.Sub64(_z[4], 6660067038095654436, b)
	_, b = bits.Sub64(_z[5], 13882719000232677960, b)
	_, b = bits.Sub64(_z[6], 12046209044522593559, b)
	_, b = bits.Sub64(_z[7], 15311794958495443299, b)
	_, b = bits.Sub64(_z[8], 18306300874381301082, b)
	_, b = bits.Sub64(_z[9], 41431377869647793, b)

	return b == 0
}

// SetRandom sets z to a random element < q
func (z *Element) SetRandom() (*Element, error) {
	var bytes [80]byte
	if _, err := io.ReadFull(rand.Reader, bytes[:]); err != nil {
		return nil, err
	}
	z[0] = binary.BigEndian.Uint64(bytes[0:8])
	z[1] = binary.BigEndian.Uint64(bytes[8:16])
	z[2] = binary.BigEndian.Uint64(bytes[16:24])
	z[3] = binary.BigEndian.Uint64(bytes[24:32])
	z[4] = binary.BigEndian.Uint64(bytes[32:40])
	z[5] = binary.BigEndian.Uint64(bytes[40:48])
	z[6] = binary.BigEndian.Uint64(bytes[48:56])
	z[7] = binary.BigEndian.Uint64(bytes[56:64])
	z[8] = binary.BigEndian.Uint64(bytes[64:72])
	z[9] = binary.BigEndian.Uint64(bytes[72:80])
	z[9] %= 82862755739295587

	// if z > q → z -= q
	// note: this is NOT constant time
	if !(z[9] < 82862755739295587 || (z[9] == 82862755739295587 && (z[8] < 18165857675053050549 || (z[8] == 18165857675053050549 && (z[7] < 12176845843281334983 || (z[7] == 12176845843281334983 && (z[6] < 5645674015335635503 || (z[6] == 5645674015335635503 && (z[5] < 9318693926755804304 || (z[5] == 9318693926755804304 && (z[4] < 13320134076191308873 || (z[4] == 13320134076191308873 && (z[3] < 9083347379620258823 || (z[3] == 9083347379620258823 && (z[2] < 15543556715411259941 || (z[2] == 15543556715411259941 && (z[1] < 4410884215886313276 || (z[1] == 4410884215886313276 && (z[0] < 15512955586897510413))))))))))))))))))) {
		var b uint64
		z[0], b = bits.Sub64(z[0], 15512955586897510413, 0)
		z[1], b = bits.Sub64(z[1], 4410884215886313276, b)
		z[2], b = bits.Sub64(z[2], 15543556715411259941, b)
		z[3], b = bits.Sub64(z[3], 9083347379620258823, b)
		z[4], b = bits.Sub64(z[4], 13320134076191308873, b)
		z[5], b = bits.Sub64(z[5], 9318693926755804304, b)
		z[6], b = bits.Sub64(z[6], 5645674015335635503, b)
		z[7], b = bits.Sub64(z[7], 12176845843281334983, b)
		z[8], b = bits.Sub64(z[8], 18165857675053050549, b)
		z[9], _ = bits.Sub64(z[9], 82862755739295587, b)
	}

	return z, nil
}

// One returns 1 (in montgommery form)
func One() Element {
	var one Element
	one.SetOne()
	return one
}

// Halve sets z to z / 2 (mod p)
func (z *Element) Halve() {
	if z[0]&1 == 1 {
		var carry uint64

		// z = z + q
		z[0], carry = bits.Add64(z[0], 15512955586897510413, 0)
		z[1], carry = bits.Add64(z[1], 4410884215886313276, carry)
		z[2], carry = bits.Add64(z[2], 15543556715411259941, carry)
		z[3], carry = bits.Add64(z[3], 9083347379620258823, carry)
		z[4], carry = bits.Add64(z[4], 13320134076191308873, carry)
		z[5], carry = bits.Add64(z[5], 9318693926755804304, carry)
		z[6], carry = bits.Add64(z[6], 5645674015335635503, carry)
		z[7], carry = bits.Add64(z[7], 12176845843281334983, carry)
		z[8], carry = bits.Add64(z[8], 18165857675053050549, carry)
		z[9], _ = bits.Add64(z[9], 82862755739295587, carry)

	}

	// z = z >> 1

	z[0] = z[0]>>1 | z[1]<<63
	z[1] = z[1]>>1 | z[2]<<63
	z[2] = z[2]>>1 | z[3]<<63
	z[3] = z[3]>>1 | z[4]<<63
	z[4] = z[4]>>1 | z[5]<<63
	z[5] = z[5]>>1 | z[6]<<63
	z[6] = z[6]>>1 | z[7]<<63
	z[7] = z[7]>>1 | z[8]<<63
	z[8] = z[8]>>1 | z[9]<<63
	z[9] >>= 1

}

// API with assembly impl

// Mul z = x * y mod q
// see https://hackmd.io/@gnark/modular_multiplication
func (z *Element) Mul(x, y *Element) *Element {
	mul(z, x, y)
	return z
}

// Square z = x * x mod q
// see https://hackmd.io/@gnark/modular_multiplication
func (z *Element) Square(x *Element) *Element {
	mul(z, x, x)
	return z
}

// FromMont converts z in place (i.e. mutates) from Montgomery to regular representation
// sets and returns z = z * 1
func (z *Element) FromMont() *Element {
	fromMont(z)
	return z
}

// Add z = x + y mod q
func (z *Element) Add(x, y *Element) *Element {
	add(z, x, y)
	return z
}

// Double z = x + x mod q, aka Lsh 1
func (z *Element) Double(x *Element) *Element {
	double(z, x)
	return z
}

// Sub  z = x - y mod q
func (z *Element) Sub(x, y *Element) *Element {
	sub(z, x, y)
	return z
}

// Neg z = q - x
func (z *Element) Neg(x *Element) *Element {
	neg(z, x)
	return z
}

// Select is a constant-time conditional move.
// If c=0, z = x0. Else z = x1
func (z *Element) Select(c int, x0 *Element, x1 *Element) *Element {
	cC := uint64((int64(c) | -int64(c)) >> 63) // "canonicized" into: 0 if c=0, -1 otherwise
	z[0] = x0[0] ^ cC&(x0[0]^x1[0])
	z[1] = x0[1] ^ cC&(x0[1]^x1[1])
	z[2] = x0[2] ^ cC&(x0[2]^x1[2])
	z[3] = x0[3] ^ cC&(x0[3]^x1[3])
	z[4] = x0[4] ^ cC&(x0[4]^x1[4])
	z[5] = x0[5] ^ cC&(x0[5]^x1[5])
	z[6] = x0[6] ^ cC&(x0[6]^x1[6])
	z[7] = x0[7] ^ cC&(x0[7]^x1[7])
	z[8] = x0[8] ^ cC&(x0[8]^x1[8])
	z[9] = x0[9] ^ cC&(x0[9]^x1[9])
	return z
}

// Generic (no ADX instructions, no AMD64) versions of multiplication and squaring algorithms

func _mulGeneric(z, x, y *Element) {

	var t [10]uint64
	var c [3]uint64
	{
		// round 0
		v := x[0]
		c[1], c[0] = bits.Mul64(v, y[0])
		m := c[0] * 13046692460116554043
		c[2] = madd0(m, 15512955586897510413, c[0])
		c[1], c[0] = madd1(v, y[1], c[1])
		c[2], t[0] = madd2(m, 4410884215886313276, c[2], c[0])
		c[1], c[0] = madd1(v, y[2], c[1])
		c[2], t[1] = madd2(m, 15543556715411259941, c[2], c[0])
		c[1], c[0] = madd1(v, y[3], c[1])
		c[2], t[2] = madd2(m, 9083347379620258823, c[2], c[0])
		c[1], c[0] = madd1(v, y[4], c[1])
		c[2], t[3] = madd2(m, 13320134076191308873, c[2], c[0])
		c[1], c[0] = madd1(v, y[5], c[1])
		c[2], t[4] = madd2(m, 9318693926755804304, c[2], c[0])
		c[1], c[0] = madd1(v, y[6], c[1])
		c[2], t[5] = madd2(m, 5645674015335635503, c[2], c[0])
		c[1], c[0] = madd1(v, y[7], c[1])
		c[2], t[6] = madd2(m, 12176845843281334983, c[2], c[0])
		c[1], c[0] = madd1(v, y[8], c[1])
		c[2], t[7] = madd2(m, 18165857675053050549, c[2], c[0])
		c[1], c[0] = madd1(v, y[9], c[1])
		t[9], t[8] = madd3(m, 82862755739295587, c[0], c[2], c[1])
	}
	{
		// round 1
		v := x[1]
		c[1], c[0] = madd1(v, y[0], t[0])
		m := c[0] * 13046692460116554043
		c[2] = madd0(m, 15512955586897510413, c[0])
		c[1], c[0] = madd2(v, y[1], c[1], t[1])
		c[2], t[0] = madd2(m, 4410884215886313276, c[2], c[0])
		c[1], c[0] = madd2(v, y[2], c[1], t[2])
		c[2], t[1] = madd2(m, 15543556715411259941, c[2], c[0])
		c[1], c[0] = madd2(v, y[3], c[1], t[3])
		c[2], t[2] = madd2(m, 9083347379620258823, c[2], c[0])
		c[1], c[0] = madd2(v, y[4], c[1], t[4])
		c[2], t[3] = madd2(m, 13320134076191308873, c[2], c[0])
		c[1], c[0] = madd2(v, y[5], c[1], t[5])
		c[2], t[4] = madd2(m, 9318693926755804304, c[2], c[0])
		c[1], c[0] = madd2(v, y[6], c[1], t[6])
		c[2], t[5] = madd2(m, 5645674015335635503, c[2], c[0])
		c[1], c[0] = madd2(v, y[7], c[1], t[7])
		c[2], t[6] = madd2(m, 12176845843281334983, c[2], c[0])
		c[1], c[0] = madd2(v, y[8], c[1], t[8])
		c[2], t[7] = madd2(m, 18165857675053050549, c[2], c[0])
		c[1], c[0] = madd2(v, y[9], c[1], t[9])
		t[9], t[8] = madd3(m, 82862755739295587, c[0], c[2], c[1])
	}
	{
		// round 2
		v := x[2]
		c[1], c[0] = madd1(v, y[0], t[0])
		m := c[0] * 13046692460116554043
		c[2] = madd0(m, 15512955586897510413, c[0])
		c[1], c[0] = madd2(v, y[1], c[1], t[1])
		c[2], t[0] = madd2(m, 4410884215886313276, c[2], c[0])
		c[1], c[0] = madd2(v, y[2], c[1], t[2])
		c[2], t[1] = madd2(m, 15543556715411259941, c[2], c[0])
		c[1], c[0] = madd2(v, y[3], c[1], t[3])
		c[2], t[2] = madd2(m, 9083347379620258823, c[2], c[0])
		c[1], c[0] = madd2(v, y[4], c[1], t[4])
		c[2], t[3] = madd2(m, 13320134076191308873, c[2], c[0])
		c[1], c[0] = madd2(v, y[5], c[1], t[5])
		c[2], t[4] = madd2(m, 9318693926755804304, c[2], c[0])
		c[1], c[0] = madd2(v, y[6], c[1], t[6])
		c[2], t[5] = madd2(m, 5645674015335635503, c[2], c[0])
		c[1], c[0] = madd2(v, y[7], c[1], t[7])
		c[2], t[6] = madd2(m, 12176845843281334983, c[2], c[0])
		c[1], c[0] = madd2(v, y[8], c[1], t[8])
		c[2], t[7] = madd2(m, 18165857675053050549, c[2], c[0])
		c[1], c[0] = madd2(v, y[9], c[1], t[9])
		t[9], t[8] = madd3(m, 82862755739295587, c[0], c[2], c[1])
	}
	{
		// round 3
		v := x[3]
		c[1], c[0] = madd1(v, y[0], t[0])
		m := c[0] * 13046692460116554043
		c[2] = madd0(m, 15512955586897510413, c[0])
		c[1], c[0] = madd2(v, y[1], c[1], t[1])
		c[2], t[0] = madd2(m, 4410884215886313276, c[2], c[0])
		c[1], c[0] = madd2(v, y[2], c[1], t[2])
		c[2], t[1] = madd2(m, 15543556715411259941, c[2], c[0])
		c[1], c[0] = madd2(v, y[3], c[1], t[3])
		c[2], t[2] = madd2(m, 9083347379620258823, c[2], c[0])
		c[1], c[0] = madd2(v, y[4], c[1], t[4])
		c[2], t[3] = madd2(m, 13320134076191308873, c[2], c[0])
		c[1], c[0] = madd2(v, y[5], c[1], t[5])
		c[2], t[4] = madd2(m, 9318693926755804304, c[2], c[0])
		c[1], c[0] = madd2(v, y[6], c[1], t[6])
		c[2], t[5] = madd2(m, 5645674015335635503, c[2], c[0])
		c[1], c[0] = madd2(v, y[7], c[1], t[7])
		c[2], t[6] = madd2(m, 12176845843281334983, c[2], c[0])
		c[1], c[0] = madd2(v, y[8], c[1], t[8])
		c[2], t[7] = madd2(m, 18165857675053050549, c[2], c[0])
		c[1], c[0] = madd2(v, y[9], c[1], t[9])
		t[9], t[8] = madd3(m, 82862755739295587, c[0], c[2], c[1])
	}
	{
		// round 4
		v := x[4]
		c[1], c[0] = madd1(v, y[0], t[0])
		m := c[0] * 13046692460116554043
		c[2] = madd0(m, 15512955586897510413, c[0])
		c[1], c[0] = madd2(v, y[1], c[1], t[1])
		c[2], t[0] = madd2(m, 4410884215886313276, c[2], c[0])
		c[1], c[0] = madd2(v, y[2], c[1], t[2])
		c[2], t[1] = madd2(m, 15543556715411259941, c[2], c[0])
		c[1], c[0] = madd2(v, y[3], c[1], t[3])
		c[2], t[2] = madd2(m, 9083347379620258823, c[2], c[0])
		c[1], c[0] = madd2(v, y[4], c[1], t[4])
		c[2], t[3] = madd2(m, 13320134076191308873, c[2], c[0])
		c[1], c[0] = madd2(v, y[5], c[1], t[5])
		c[2], t[4] = madd2(m, 9318693926755804304, c[2], c[0])
		c[1], c[0] = madd2(v, y[6], c[1], t[6])
		c[2], t[5] = madd2(m, 5645674015335635503, c[2], c[0])
		c[1], c[0] = madd2(v, y[7], c[1], t[7])
		c[2], t[6] = madd2(m, 12176845843281334983, c[2], c[0])
		c[1], c[0] = madd2(v, y[8], c[1], t[8])
		c[2], t[7] = madd2(m, 18165857675053050549, c[2], c[0])
		c[1], c[0] = madd2(v, y[9], c[1], t[9])
		t[9], t[8] = madd3(m, 82862755739295587, c[0], c[2], c[1])
	}
	{
		// round 5
		v := x[5]
		c[1], c[0] = madd1(v, y[0], t[0])
		m := c[0] * 13046692460116554043
		c[2] = madd0(m, 15512955586897510413, c[0])
		c[1], c[0] = madd2(v, y[1], c[1], t[1])
		c[2], t[0] = madd2(m, 4410884215886313276, c[2], c[0])
		c[1], c[0] = madd2(v, y[2], c[1], t[2])
		c[2], t[1] = madd2(m, 15543556715411259941, c[2], c[0])
		c[1], c[0] = madd2(v, y[3], c[1], t[3])
		c[2], t[2] = madd2(m, 9083347379620258823, c[2], c[0])
		c[1], c[0] = madd2(v, y[4], c[1], t[4])
		c[2], t[3] = madd2(m, 13320134076191308873, c[2], c[0])
		c[1], c[0] = madd2(v, y[5], c[1], t[5])
		c[2], t[4] = madd2(m, 9318693926755804304, c[2], c[0])
		c[1], c[0] = madd2(v, y[6], c[1], t[6])
		c[2], t[5] = madd2(m, 5645674015335635503, c[2], c[0])
		c[1], c[0] = madd2(v, y[7], c[1], t[7])
		c[2], t[6] = madd2(m, 12176845843281334983, c[2], c[0])
		c[1], c[0] = madd2(v, y[8], c[1], t[8])
		c[2], t[7] = madd2(m, 18165857675053050549, c[2], c[0])
		c[1], c[0] = madd2(v, y[9], c[1], t[9])
		t[9], t[8] = madd3(m, 82862755739295587, c[0], c[2], c[1])
	}
	{
		// round 6
		v := x[6]
		c[1], c[0] = madd1(v, y[0], t[0])
		m := c[0] * 13046692460116554043
		c[2] = madd0(m, 15512955586897510413, c[0])
		c[1], c[0] = madd2(v, y[1], c[1], t[1])
		c[2], t[0] = madd2(m, 4410884215886313276, c[2], c[0])
		c[1], c[0] = madd2(v, y[2], c[1], t[2])
		c[2], t[1] = madd2(m, 15543556715411259941, c[2], c[0])
		c[1], c[0] = madd2(v, y[3], c[1], t[3])
		c[2], t[2] = madd2(m, 9083347379620258823, c[2], c[0])
		c[1], c[0] = madd2(v, y[4], c[1], t[4])
		c[2], t[3] = madd2(m, 13320134076191308873, c[2], c[0])
		c[1], c[0] = madd2(v, y[5], c[1], t[5])
		c[2], t[4] = madd2(m, 9318693926755804304, c[2], c[0])
		c[1], c[0] = madd2(v, y[6], c[1], t[6])
		c[2], t[5] = madd2(m, 5645674015335635503, c[2], c[0])
		c[1], c[0] = madd2(v, y[7], c[1], t[7])
		c[2], t[6] = madd2(m, 12176845843281334983, c[2], c[0])
		c[1], c[0] = madd2(v, y[8], c[1], t[8])
		c[2], t[7] = madd2(m, 18165857675053050549, c[2], c[0])
		c[1], c[0] = madd2(v, y[9], c[1], t[9])
		t[9], t[8] = madd3(m, 82862755739295587, c[0], c[2], c[1])
	}
	{
		// round 7
		v := x[7]
		c[1], c[0] = madd1(v, y[0], t[0])
		m := c[0] * 13046692460116554043
		c[2] = madd0(m, 15512955586897510413, c[0])
		c[1], c[0] = madd2(v, y[1], c[1], t[1])
		c[2], t[0] = madd2(m, 4410884215886313276, c[2], c[0])
		c[1], c[0] = madd2(v, y[2], c[1], t[2])
		c[2], t[1] = madd2(m, 15543556715411259941, c[2], c[0])
		c[1], c[0] = madd2(v, y[3], c[1], t[3])
		c[2], t[2] = madd2(m, 9083347379620258823, c[2], c[0])
		c[1], c[0] = madd2(v, y[4], c[1], t[4])
		c[2], t[3] = madd2(m, 13320134076191308873, c[2], c[0])
		c[1], c[0] = madd2(v, y[5], c[1], t[5])
		c[2], t[4] = madd2(m, 9318693926755804304, c[2], c[0])
		c[1], c[0] = madd2(v, y[6], c[1], t[6])
		c[2], t[5] = madd2(m, 5645674015335635503, c[2], c[0])
		c[1], c[0] = madd2(v, y[7], c[1], t[7])
		c[2], t[6] = madd2(m, 12176845843281334983, c[2], c[0])
		c[1], c[0] = madd2(v, y[8], c[1], t[8])
		c[2], t[7] = madd2(m, 18165857675053050549, c[2], c[0])
		c[1], c[0] = madd2(v, y[9], c[1], t[9])
		t[9], t[8] = madd3(m, 82862755739295587, c[0], c[2], c[1])
	}
	{
		// round 8
		v := x[8]
		c[1], c[0] = madd1(v, y[0], t[0])
		m := c[0] * 13046692460116554043
		c[2] = madd0(m, 15512955586897510413, c[0])
		c[1], c[0] = madd2(v, y[1], c[1], t[1])
		c[2], t[0] = madd2(m, 4410884215886313276, c[2], c[0])
		c[1], c[0] = madd2(v, y[2], c[1], t[2])
		c[2], t[1] = madd2(m, 15543556715411259941, c[2], c[0])
		c[1], c[0] = madd2(v, y[3], c[1], t[3])
		c[2], t[2] = madd2(m, 9083347379620258823, c[2], c[0])
		c[1], c[0] = madd2(v, y[4], c[1], t[4])
		c[2], t[3] = madd2(m, 13320134076191308873, c[2], c[0])
		c[1], c[0] = madd2(v, y[5], c[1], t[5])
		c[2], t[4] = madd2(m, 9318693926755804304, c[2], c[0])
		c[1], c[0] = madd2(v, y[6], c[1], t[6])
		c[2], t[5] = madd2(m, 5645674015335635503, c[2], c[0])
		c[1], c[0] = madd2(v, y[7], c[1], t[7])
		c[2], t[6] = madd2(m, 12176845843281334983, c[2], c[0])
		c[1], c[0] = madd2(v, y[8], c[1], t[8])
		c[2], t[7] = madd2(m, 18165857675053050549, c[2], c[0])
		c[1], c[0] = madd2(v, y[9], c[1], t[9])
		t[9], t[8] = madd3(m, 82862755739295587, c[0], c[2], c[1])
	}
	{
		// round 9
		v := x[9]
		c[1], c[0] = madd1(v, y[0], t[0])
		m := c[0] * 13046692460116554043
		c[2] = madd0(m, 15512955586897510413, c[0])
		c[1], c[0] = madd2(v, y[1], c[1], t[1])
		c[2], z[0] = madd2(m, 4410884215886313276, c[2], c[0])
		c[1], c[0] = madd2(v, y[2], c[1], t[2])
		c[2], z[1] = madd2(m, 15543556715411259941, c[2], c[0])
		c[1], c[0] = madd2(v, y[3], c[1], t[3])
		c[2], z[2] = madd2(m, 9083347379620258823, c[2], c[0])
		c[1], c[0] = madd2(v, y[4], c[1], t[4])
		c[2], z[3] = madd2(m, 13320134076191308873, c[2], c[0])
		c[1], c[0] = madd2(v, y[5], c[1], t[5])
		c[2], z[4] = madd2(m, 9318693926755804304, c[2], c[0])
		c[1], c[0] = madd2(v, y[6], c[1], t[6])
		c[2], z[5] = madd2(m, 5645674015335635503, c[2], c[0])
		c[1], c[0] = madd2(v, y[7], c[1], t[7])
		c[2], z[6] = madd2(m, 12176845843281334983, c[2], c[0])
		c[1], c[0] = madd2(v, y[8], c[1], t[8])
		c[2], z[7] = madd2(m, 18165857675053050549, c[2], c[0])
		c[1], c[0] = madd2(v, y[9], c[1], t[9])
		z[9], z[8] = madd3(m, 82862755739295587, c[0], c[2], c[1])
	}

	// if z > q → z -= q
	// note: this is NOT constant time
	if !(z[9] < 82862755739295587 || (z[9] == 82862755739295587 && (z[8] < 18165857675053050549 || (z[8] == 18165857675053050549 && (z[7] < 12176845843281334983 || (z[7] == 12176845843281334983 && (z[6] < 5645674015335635503 || (z[6] == 5645674015335635503 && (z[5] < 9318693926755804304 || (z[5] == 9318693926755804304 && (z[4] < 13320134076191308873 || (z[4] == 13320134076191308873 && (z[3] < 9083347379620258823 || (z[3] == 9083347379620258823 && (z[2] < 15543556715411259941 || (z[2] == 15543556715411259941 && (z[1] < 4410884215886313276 || (z[1] == 4410884215886313276 && (z[0] < 15512955586897510413))))))))))))))))))) {
		var b uint64
		z[0], b = bits.Sub64(z[0], 15512955586897510413, 0)
		z[1], b = bits.Sub64(z[1], 4410884215886313276, b)
		z[2], b = bits.Sub64(z[2], 15543556715411259941, b)
		z[3], b = bits.Sub64(z[3], 9083347379620258823, b)
		z[4], b = bits.Sub64(z[4], 13320134076191308873, b)
		z[5], b = bits.Sub64(z[5], 9318693926755804304, b)
		z[6], b = bits.Sub64(z[6], 5645674015335635503, b)
		z[7], b = bits.Sub64(z[7], 12176845843281334983, b)
		z[8], b = bits.Sub64(z[8], 18165857675053050549, b)
		z[9], _ = bits.Sub64(z[9], 82862755739295587, b)
	}
}

func _mulWGeneric(z, x *Element, y uint64) {

	var t [10]uint64
	{
		// round 0
		c1, c0 := bits.Mul64(y, x[0])
		m := c0 * 13046692460116554043
		c2 := madd0(m, 15512955586897510413, c0)
		c1, c0 = madd1(y, x[1], c1)
		c2, t[0] = madd2(m, 4410884215886313276, c2, c0)
		c1, c0 = madd1(y, x[2], c1)
		c2, t[1] = madd2(m, 15543556715411259941, c2, c0)
		c1, c0 = madd1(y, x[3], c1)
		c2, t[2] = madd2(m, 9083347379620258823, c2, c0)
		c1, c0 = madd1(y, x[4], c1)
		c2, t[3] = madd2(m, 13320134076191308873, c2, c0)
		c1, c0 = madd1(y, x[5], c1)
		c2, t[4] = madd2(m, 9318693926755804304, c2, c0)
		c1, c0 = madd1(y, x[6], c1)
		c2, t[5] = madd2(m, 5645674015335635503, c2, c0)
		c1, c0 = madd1(y, x[7], c1)
		c2, t[6] = madd2(m, 12176845843281334983, c2, c0)
		c1, c0 = madd1(y, x[8], c1)
		c2, t[7] = madd2(m, 18165857675053050549, c2, c0)
		c1, c0 = madd1(y, x[9], c1)
		t[9], t[8] = madd3(m, 82862755739295587, c0, c2, c1)
	}
	{
		// round 1
		m := t[0] * 13046692460116554043
		c2 := madd0(m, 15512955586897510413, t[0])
		c2, t[0] = madd2(m, 4410884215886313276, c2, t[1])
		c2, t[1] = madd2(m, 15543556715411259941, c2, t[2])
		c2, t[2] = madd2(m, 9083347379620258823, c2, t[3])
		c2, t[3] = madd2(m, 13320134076191308873, c2, t[4])
		c2, t[4] = madd2(m, 9318693926755804304, c2, t[5])
		c2, t[5] = madd2(m, 5645674015335635503, c2, t[6])
		c2, t[6] = madd2(m, 12176845843281334983, c2, t[7])
		c2, t[7] = madd2(m, 18165857675053050549, c2, t[8])
		t[9], t[8] = madd2(m, 82862755739295587, t[9], c2)
	}
	{
		// round 2
		m := t[0] * 13046692460116554043
		c2 := madd0(m, 15512955586897510413, t[0])
		c2, t[0] = madd2(m, 4410884215886313276, c2, t[1])
		c2, t[1] = madd2(m, 15543556715411259941, c2, t[2])
		c2, t[2] = madd2(m, 9083347379620258823, c2, t[3])
		c2, t[3] = madd2(m, 13320134076191308873, c2, t[4])
		c2, t[4] = madd2(m, 9318693926755804304, c2, t[5])
		c2, t[5] = madd2(m, 5645674015335635503, c2, t[6])
		c2, t[6] = madd2(m, 12176845843281334983, c2, t[7])
		c2, t[7] = madd2(m, 18165857675053050549, c2, t[8])
		t[9], t[8] = madd2(m, 82862755739295587, t[9], c2)
	}
	{
		// round 3
		m := t[0] * 13046692460116554043
		c2 := madd0(m, 15512955586897510413, t[0])
		c2, t[0] = madd2(m, 4410884215886313276, c2, t[1])
		c2, t[1] = madd2(m, 15543556715411259941, c2, t[2])
		c2, t[2] = madd2(m, 9083347379620258823, c2, t[3])
		c2, t[3] = madd2(m, 13320134076191308873, c2, t[4])
		c2, t[4] = madd2(m, 9318693926755804304, c2, t[5])
		c2, t[5] = madd2(m, 5645674015335635503, c2, t[6])
		c2, t[6] = madd2(m, 12176845843281334983, c2, t[7])
		c2, t[7] = madd2(m, 18165857675053050549, c2, t[8])
		t[9], t[8] = madd2(m, 82862755739295587, t[9], c2)
	}
	{
		// round 4
		m := t[0] * 13046692460116554043
		c2 := madd0(m, 15512955586897510413, t[0])
		c2, t[0] = madd2(m, 4410884215886313276, c2, t[1])
		c2, t[1] = madd2(m, 15543556715411259941, c2, t[2])
		c2, t[2] = madd2(m, 9083347379620258823, c2, t[3])
		c2, t[3] = madd2(m, 13320134076191308873, c2, t[4])
		c2, t[4] = madd2(m, 9318693926755804304, c2, t[5])
		c2, t[5] = madd2(m, 5645674015335635503, c2, t[6])
		c2, t[6] = madd2(m, 12176845843281334983, c2, t[7])
		c2, t[7] = madd2(m, 18165857675053050549, c2, t[8])
		t[9], t[8] = madd2(m, 82862755739295587, t[9], c2)
	}
	{
		// round 5
		m := t[0] * 13046692460116554043
		c2 := madd0(m, 15512955586897510413, t[0])
		c2, t[0] = madd2(m, 4410884215886313276, c2, t[1])
		c2, t[1] = madd2(m, 15543556715411259941, c2, t[2])
		c2, t[2] = madd2(m, 9083347379620258823, c2, t[3])
		c2, t[3] = madd2(m, 13320134076191308873, c2, t[4])
		c2, t[4] = madd2(m, 9318693926755804304, c2, t[5])
		c2, t[5] = madd2(m, 5645674015335635503, c2, t[6])
		c2, t[6] = madd2(m, 12176845843281334983, c2, t[7])
		c2, t[7] = madd2(m, 18165857675053050549, c2, t[8])
		t[9], t[8] = madd2(m, 82862755739295587, t[9], c2)
	}
	{
		// round 6
		m := t[0] * 13046692460116554043
		c2 := madd0(m, 15512955586897510413, t[0])
		c2, t[0] = madd2(m, 4410884215886313276, c2, t[1])
		c2, t[1] = madd2(m, 15543556715411259941, c2, t[2])
		c2, t[2] = madd2(m, 9083347379620258823, c2, t[3])
		c2, t[3] = madd2(m, 13320134076191308873, c2, t[4])
		c2, t[4] = madd2(m, 9318693926755804304, c2, t[5])
		c2, t[5] = madd2(m, 5645674015335635503, c2, t[6])
		c2, t[6] = madd2(m, 12176845843281334983, c2, t[7])
		c2, t[7] = madd2(m, 18165857675053050549, c2, t[8])
		t[9], t[8] = madd2(m, 82862755739295587, t[9], c2)
	}
	{
		// round 7
		m := t[0] * 13046692460116554043
		c2 := madd0(m, 15512955586897510413, t[0])
		c2, t[0] = madd2(m, 4410884215886313276, c2, t[1])
		c2, t[1] = madd2(m, 15543556715411259941, c2, t[2])
		c2, t[2] = madd2(m, 9083347379620258823, c2, t[3])
		c2, t[3] = madd2(m, 13320134076191308873, c2, t[4])
		c2, t[4] = madd2(m, 9318693926755804304, c2, t[5])
		c2, t[5] = madd2(m, 5645674015335635503, c2, t[6])
		c2, t[6] = madd2(m, 12176845843281334983, c2, t[7])
		c2, t[7] = madd2(m, 18165857675053050549, c2, t[8])
		t[9], t[8] = madd2(m, 82862755739295587, t[9], c2)
	}
	{
		// round 8
		m := t[0] * 13046692460116554043
		c2 := madd0(m, 15512955586897510413, t[0])
		c2, t[0] = madd2(m, 4410884215886313276, c2, t[1])
		c2, t[1] = madd2(m, 15543556715411259941, c2, t[2])
		c2, t[2] = madd2(m, 9083347379620258823, c2, t[3])
		c2, t[3] = madd2(m, 13320134076191308873, c2, t[4])
		c2, t[4] = madd2(m, 9318693926755804304, c2, t[5])
		c2, t[5] = madd2(m, 5645674015335635503, c2, t[6])
		c2, t[6] = madd2(m, 12176845843281334983, c2, t[7])
		c2, t[7] = madd2(m, 18165857675053050549, c2, t[8])
		t[9], t[8] = madd2(m, 82862755739295587, t[9], c2)
	}
	{
		// round 9
		m := t[0] * 13046692460116554043
		c2 := madd0(m, 15512955586897510413, t[0])
		c2, z[0] = madd2(m, 4410884215886313276, c2, t[1])
		c2, z[1] = madd2(m, 15543556715411259941, c2, t[2])
		c2, z[2] = madd2(m, 9083347379620258823, c2, t[3])
		c2, z[3] = madd2(m, 13320134076191308873, c2, t[4])
		c2, z[4] = madd2(m, 9318693926755804304, c2, t[5])
		c2, z[5] = madd2(m, 5645674015335635503, c2, t[6])
		c2, z[6] = madd2(m, 12176845843281334983, c2, t[7])
		c2, z[7] = madd2(m, 18165857675053050549, c2, t[8])
		z[9], z[8] = madd2(m, 82862755739295587, t[9], c2)
	}

	// if z > q → z -= q
	// note: this is NOT constant time
	if !(z[9] < 82862755739295587 || (z[9] == 82862755739295587 && (z[8] < 18165857675053050549 || (z[8] == 18165857675053050549 && (z[7] < 12176845843281334983 || (z[7] == 12176845843281334983 && (z[6] < 5645674015335635503 || (z[6] == 5645674015335635503 && (z[5] < 9318693926755804304 || (z[5] == 9318693926755804304 && (z[4] < 13320134076191308873 || (z[4] == 13320134076191308873 && (z[3] < 9083347379620258823 || (z[3] == 9083347379620258823 && (z[2] < 15543556715411259941 || (z[2] == 15543556715411259941 && (z[1] < 4410884215886313276 || (z[1] == 4410884215886313276 && (z[0] < 15512955586897510413))))))))))))))))))) {
		var b uint64
		z[0], b = bits.Sub64(z[0], 15512955586897510413, 0)
		z[1], b = bits.Sub64(z[1], 4410884215886313276, b)
		z[2], b = bits.Sub64(z[2], 15543556715411259941, b)
		z[3], b = bits.Sub64(z[3], 9083347379620258823, b)
		z[4], b = bits.Sub64(z[4], 13320134076191308873, b)
		z[5], b = bits.Sub64(z[5], 9318693926755804304, b)
		z[6], b = bits.Sub64(z[6], 5645674015335635503, b)
		z[7], b = bits.Sub64(z[7], 12176845843281334983, b)
		z[8], b = bits.Sub64(z[8], 18165857675053050549, b)
		z[9], _ = bits.Sub64(z[9], 82862755739295587, b)
	}
}

func _fromMontGeneric(z *Element) {
	// the following lines implement z = z * 1
	// with a modified CIOS montgomery multiplication
	{
		// m = z[0]n'[0] mod W
		m := z[0] * 13046692460116554043
		C := madd0(m, 15512955586897510413, z[0])
		C, z[0] = madd2(m, 4410884215886313276, z[1], C)
		C, z[1] = madd2(m, 15543556715411259941, z[2], C)
		C, z[2] = madd2(m, 9083347379620258823, z[3], C)
		C, z[3] = madd2(m, 13320134076191308873, z[4], C)
		C, z[4] = madd2(m, 9318693926755804304, z[5], C)
		C, z[5] = madd2(m, 5645674015335635503, z[6], C)
		C, z[6] = madd2(m, 12176845843281334983, z[7], C)
		C, z[7] = madd2(m, 18165857675053050549, z[8], C)
		C, z[8] = madd2(m, 82862755739295587, z[9], C)
		z[9] = C
	}
	{
		// m = z[0]n'[0] mod W
		m := z[0] * 13046692460116554043
		C := madd0(m, 15512955586897510413, z[0])
		C, z[0] = madd2(m, 4410884215886313276, z[1], C)
		C, z[1] = madd2(m, 15543556715411259941, z[2], C)
		C, z[2] = madd2(m, 9083347379620258823, z[3], C)
		C, z[3] = madd2(m, 13320134076191308873, z[4], C)
		C, z[4] = madd2(m, 9318693926755804304, z[5], C)
		C, z[5] = madd2(m, 5645674015335635503, z[6], C)
		C, z[6] = madd2(m, 12176845843281334983, z[7], C)
		C, z[7] = madd2(m, 18165857675053050549, z[8], C)
		C, z[8] = madd2(m, 82862755739295587, z[9], C)
		z[9] = C
	}
	{
		// m = z[0]n'[0] mod W
		m := z[0] * 13046692460116554043
		C := madd0(m, 15512955586897510413, z[0])
		C, z[0] = madd2(m, 4410884215886313276, z[1], C)
		C, z[1] = madd2(m, 15543556715411259941, z[2], C)
		C, z[2] = madd2(m, 9083347379620258823, z[3], C)
		C, z[3] = madd2(m, 13320134076191308873, z[4], C)
		C, z[4] = madd2(m, 9318693926755804304, z[5], C)
		C, z[5] = madd2(m, 5645674015335635503, z[6], C)
		C, z[6] = madd2(m, 12176845843281334983, z[7], C)
		C, z[7] = madd2(m, 18165857675053050549, z[8], C)
		C, z[8] = madd2(m, 82862755739295587, z[9], C)
		z[9] = C
	}
	{
		// m = z[0]n'[0] mod W
		m := z[0] * 13046692460116554043
		C := madd0(m, 15512955586897510413, z[0])
		C, z[0] = madd2(m, 4410884215886313276, z[1], C)
		C, z[1] = madd2(m, 15543556715411259941, z[2], C)
		C, z[2] = madd2(m, 9083347379620258823, z[3], C)
		C, z[3] = madd2(m, 13320134076191308873, z[4], C)
		C, z[4] = madd2(m, 9318693926755804304, z[5], C)
		C, z[5] = madd2(m, 5645674015335635503, z[6], C)
		C, z[6] = madd2(m, 12176845843281334983, z[7], C)
		C, z[7] = madd2(m, 18165857675053050549, z[8], C)
		C, z[8] = madd2(m, 82862755739295587, z[9], C)
		z[9] = C
	}
	{
		// m = z[0]n'[0] mod W
		m := z[0] * 13046692460116554043
		C := madd0(m, 15512955586897510413, z[0])
		C, z[0] = madd2(m, 4410884215886313276, z[1], C)
		C, z[1] = madd2(m, 15543556715411259941, z[2], C)
		C, z[2] = madd2(m, 9083347379620258823, z[3], C)
		C, z[3] = madd2(m, 13320134076191308873, z[4], C)
		C, z[4] = madd2(m, 9318693926755804304, z[5], C)
		C, z[5] = madd2(m, 5645674015335635503, z[6], C)
		C, z[6] = madd2(m, 12176845843281334983, z[7], C)
		C, z[7] = madd2(m, 18165857675053050549, z[8], C)
		C, z[8] = madd2(m, 82862755739295587, z[9], C)
		z[9] = C
	}
	{
		// m = z[0]n'[0] mod W
		m := z[0] * 13046692460116554043
		C := madd0(m, 15512955586897510413, z[0])
		C, z[0] = madd2(m, 4410884215886313276, z[1], C)
		C, z[1] = madd2(m, 15543556715411259941, z[2], C)
		C, z[2] = madd2(m, 9083347379620258823, z[3], C)
		C, z[3] = madd2(m, 13320134076191308873, z[4], C)
		C, z[4] = madd2(m, 9318693926755804304, z[5], C)
		C, z[5] = madd2(m, 5645674015335635503, z[6], C)
		C, z[6] = madd2(m, 12176845843281334983, z[7], C)
		C, z[7] = madd2(m, 18165857675053050549, z[8], C)
		C, z[8] = madd2(m, 82862755739295587, z[9], C)
		z[9] = C
	}
	{
		// m = z[0]n'[0] mod W
		m := z[0] * 13046692460116554043
		C := madd0(m, 15512955586897510413, z[0])
		C, z[0] = madd2(m, 4410884215886313276, z[1], C)
		C, z[1] = madd2(m, 15543556715411259941, z[2], C)
		C, z[2] = madd2(m, 9083347379620258823, z[3], C)
		C, z[3] = madd2(m, 13320134076191308873, z[4], C)
		C, z[4] = madd2(m, 9318693926755804304, z[5], C)
		C, z[5] = madd2(m, 5645674015335635503, z[6], C)
		C, z[6] = madd2(m, 12176845843281334983, z[7], C)
		C, z[7] = madd2(m, 18165857675053050549, z[8], C)
		C, z[8] = madd2(m, 82862755739295587, z[9], C)
		z[9] = C
	}
	{
		// m = z[0]n'[0] mod W
		m := z[0] * 13046692460116554043
		C := madd0(m, 15512955586897510413, z[0])
		C, z[0] = madd2(m, 4410884215886313276, z[1], C)
		C, z[1] = madd2(m, 15543556715411259941, z[2], C)
		C, z[2] = madd2(m, 9083347379620258823, z[3], C)
		C, z[3] = madd2(m, 13320134076191308873, z[4], C)
		C, z[4] = madd2(m, 9318693926755804304, z[5], C)
		C, z[5] = madd2(m, 5645674015335635503, z[6], C)
		C, z[6] = madd2(m, 12176845843281334983, z[7], C)
		C, z[7] = madd2(m, 18165857675053050549, z[8], C)
		C, z[8] = madd2(m, 82862755739295587, z[9], C)
		z[9] = C
	}
	{
		// m = z[0]n'[0] mod W
		m := z[0] * 13046692460116554043
		C := madd0(m, 15512955586897510413, z[0])
		C, z[0] = madd2(m, 4410884215886313276, z[1], C)
		C, z[1] = madd2(m, 15543556715411259941, z[2], C)
		C, z[2] = madd2(m, 9083347379620258823, z[3], C)
		C, z[3] = madd2(m, 13320134076191308873, z[4], C)
		C, z[4] = madd2(m, 9318693926755804304, z[5], C)
		C, z[5] = madd2(m, 5645674015335635503, z[6], C)
		C, z[6] = madd2(m, 12176845843281334983, z[7], C)
		C, z[7] = madd2(m, 18165857675053050549, z[8], C)
		C, z[8] = madd2(m, 82862755739295587, z[9], C)
		z[9] = C
	}
	{
		// m = z[0]n'[0] mod W
		m := z[0] * 13046692460116554043
		C := madd0(m, 15512955586897510413, z[0])
		C, z[0] = madd2(m, 4410884215886313276, z[1], C)
		C, z[1] = madd2(m, 15543556715411259941, z[2], C)
		C, z[2] = madd2(m, 9083347379620258823, z[3], C)
		C, z[3] = madd2(m, 13320134076191308873, z[4], C)
		C, z[4] = madd2(m, 9318693926755804304, z[5], C)
		C, z[5] = madd2(m, 5645674015335635503, z[6], C)
		C, z[6] = madd2(m, 12176845843281334983, z[7], C)
		C, z[7] = madd2(m, 18165857675053050549, z[8], C)
		C, z[8] = madd2(m, 82862755739295587, z[9], C)
		z[9] = C
	}

	// if z > q → z -= q
	// note: this is NOT constant time
	if !(z[9] < 82862755739295587 || (z[9] == 82862755739295587 && (z[8] < 18165857675053050549 || (z[8] == 18165857675053050549 && (z[7] < 12176845843281334983 || (z[7] == 12176845843281334983 && (z[6] < 5645674015335635503 || (z[6] == 5645674015335635503 && (z[5] < 9318693926755804304 || (z[5] == 9318693926755804304 && (z[4] < 13320134076191308873 || (z[4] == 13320134076191308873 && (z[3] < 9083347379620258823 || (z[3] == 9083347379620258823 && (z[2] < 15543556715411259941 || (z[2] == 15543556715411259941 && (z[1] < 4410884215886313276 || (z[1] == 4410884215886313276 && (z[0] < 15512955586897510413))))))))))))))))))) {
		var b uint64
		z[0], b = bits.Sub64(z[0], 15512955586897510413, 0)
		z[1], b = bits.Sub64(z[1], 4410884215886313276, b)
		z[2], b = bits.Sub64(z[2], 15543556715411259941, b)
		z[3], b = bits.Sub64(z[3], 9083347379620258823, b)
		z[4], b = bits.Sub64(z[4], 13320134076191308873, b)
		z[5], b = bits.Sub64(z[5], 9318693926755804304, b)
		z[6], b = bits.Sub64(z[6], 5645674015335635503, b)
		z[7], b = bits.Sub64(z[7], 12176845843281334983, b)
		z[8], b = bits.Sub64(z[8], 18165857675053050549, b)
		z[9], _ = bits.Sub64(z[9], 82862755739295587, b)
	}
}

func _addGeneric(z, x, y *Element) {
	var carry uint64

	z[0], carry = bits.Add64(x[0], y[0], 0)
	z[1], carry = bits.Add64(x[1], y[1], carry)
	z[2], carry = bits.Add64(x[2], y[2], carry)
	z[3], carry = bits.Add64(x[3], y[3], carry)
	z[4], carry = bits.Add64(x[4], y[4], carry)
	z[5], carry = bits.Add64(x[5], y[5], carry)
	z[6], carry = bits.Add64(x[6], y[6], carry)
	z[7], carry = bits.Add64(x[7], y[7], carry)
	z[8], carry = bits.Add64(x[8], y[8], carry)
	z[9], _ = bits.Add64(x[9], y[9], carry)

	// if z > q → z -= q
	// note: this is NOT constant time
	if !(z[9] < 82862755739295587 || (z[9] == 82862755739295587 && (z[8] < 18165857675053050549 || (z[8] == 18165857675053050549 && (z[7] < 12176845843281334983 || (z[7] == 12176845843281334983 && (z[6] < 5645674015335635503 || (z[6] == 5645674015335635503 && (z[5] < 9318693926755804304 || (z[5] == 9318693926755804304 && (z[4] < 13320134076191308873 || (z[4] == 13320134076191308873 && (z[3] < 9083347379620258823 || (z[3] == 9083347379620258823 && (z[2] < 15543556715411259941 || (z[2] == 15543556715411259941 && (z[1] < 4410884215886313276 || (z[1] == 4410884215886313276 && (z[0] < 15512955586897510413))))))))))))))))))) {
		var b uint64
		z[0], b = bits.Sub64(z[0], 15512955586897510413, 0)
		z[1], b = bits.Sub64(z[1], 4410884215886313276, b)
		z[2], b = bits.Sub64(z[2], 15543556715411259941, b)
		z[3], b = bits.Sub64(z[3], 9083347379620258823, b)
		z[4], b = bits.Sub64(z[4], 13320134076191308873, b)
		z[5], b = bits.Sub64(z[5], 9318693926755804304, b)
		z[6], b = bits.Sub64(z[6], 5645674015335635503, b)
		z[7], b = bits.Sub64(z[7], 12176845843281334983, b)
		z[8], b = bits.Sub64(z[8], 18165857675053050549, b)
		z[9], _ = bits.Sub64(z[9], 82862755739295587, b)
	}
}

func _doubleGeneric(z, x *Element) {
	var carry uint64

	z[0], carry = bits.Add64(x[0], x[0], 0)
	z[1], carry = bits.Add64(x[1], x[1], carry)
	z[2], carry = bits.Add64(x[2], x[2], carry)
	z[3], carry = bits.Add64(x[3], x[3], carry)
	z[4], carry = bits.Add64(x[4], x[4], carry)
	z[5], carry = bits.Add64(x[5], x[5], carry)
	z[6], carry = bits.Add64(x[6], x[6], carry)
	z[7], carry = bits.Add64(x[7], x[7], carry)
	z[8], carry = bits.Add64(x[8], x[8], carry)
	z[9], _ = bits.Add64(x[9], x[9], carry)

	// if z > q → z -= q
	// note: this is NOT constant time
	if !(z[9] < 82862755739295587 || (z[9] == 82862755739295587 && (z[8] < 18165857675053050549 || (z[8] == 18165857675053050549 && (z[7] < 12176845843281334983 || (z[7] == 12176845843281334983 && (z[6] < 5645674015335635503 || (z[6] == 5645674015335635503 && (z[5] < 9318693926755804304 || (z[5] == 9318693926755804304 && (z[4] < 13320134076191308873 || (z[4] == 13320134076191308873 && (z[3] < 9083347379620258823 || (z[3] == 9083347379620258823 && (z[2] < 15543556715411259941 || (z[2] == 15543556715411259941 && (z[1] < 4410884215886313276 || (z[1] == 4410884215886313276 && (z[0] < 15512955586897510413))))))))))))))))))) {
		var b uint64
		z[0], b = bits.Sub64(z[0], 15512955586897510413, 0)
		z[1], b = bits.Sub64(z[1], 4410884215886313276, b)
		z[2], b = bits.Sub64(z[2], 15543556715411259941, b)
		z[3], b = bits.Sub64(z[3], 9083347379620258823, b)
		z[4], b = bits.Sub64(z[4], 13320134076191308873, b)
		z[5], b = bits.Sub64(z[5], 9318693926755804304, b)
		z[6], b = bits.Sub64(z[6], 5645674015335635503, b)
		z[7], b = bits.Sub64(z[7], 12176845843281334983, b)
		z[8], b = bits.Sub64(z[8], 18165857675053050549, b)
		z[9], _ = bits.Sub64(z[9], 82862755739295587, b)
	}
}

func _subGeneric(z, x, y *Element) {
	var b uint64
	z[0], b = bits.Sub64(x[0], y[0], 0)
	z[1], b = bits.Sub64(x[1], y[1], b)
	z[2], b = bits.Sub64(x[2], y[2], b)
	z[3], b = bits.Sub64(x[3], y[3], b)
	z[4], b = bits.Sub64(x[4], y[4], b)
	z[5], b = bits.Sub64(x[5], y[5], b)
	z[6], b = bits.Sub64(x[6], y[6], b)
	z[7], b = bits.Sub64(x[7], y[7], b)
	z[8], b = bits.Sub64(x[8], y[8], b)
	z[9], b = bits.Sub64(x[9], y[9], b)
	if b != 0 {
		var c uint64
		z[0], c = bits.Add64(z[0], 15512955586897510413, 0)
		z[1], c = bits.Add64(z[1], 4410884215886313276, c)
		z[2], c = bits.Add64(z[2], 15543556715411259941, c)
		z[3], c = bits.Add64(z[3], 9083347379620258823, c)
		z[4], c = bits.Add64(z[4], 13320134076191308873, c)
		z[5], c = bits.Add64(z[5], 9318693926755804304, c)
		z[6], c = bits.Add64(z[6], 5645674015335635503, c)
		z[7], c = bits.Add64(z[7], 12176845843281334983, c)
		z[8], c = bits.Add64(z[8], 18165857675053050549, c)
		z[9], _ = bits.Add64(z[9], 82862755739295587, c)
	}
}

func _negGeneric(z, x *Element) {
	if x.IsZero() {
		z.SetZero()
		return
	}
	var borrow uint64
	z[0], borrow = bits.Sub64(15512955586897510413, x[0], 0)
	z[1], borrow = bits.Sub64(4410884215886313276, x[1], borrow)
	z[2], borrow = bits.Sub64(15543556715411259941, x[2], borrow)
	z[3], borrow = bits.Sub64(9083347379620258823, x[3], borrow)
	z[4], borrow = bits.Sub64(13320134076191308873, x[4], borrow)
	z[5], borrow = bits.Sub64(9318693926755804304, x[5], borrow)
	z[6], borrow = bits.Sub64(5645674015335635503, x[6], borrow)
	z[7], borrow = bits.Sub64(12176845843281334983, x[7], borrow)
	z[8], borrow = bits.Sub64(18165857675053050549, x[8], borrow)
	z[9], _ = bits.Sub64(82862755739295587, x[9], borrow)
}

func _reduceGeneric(z *Element) {

	// if z > q → z -= q
	// note: this is NOT constant time
	if !(z[9] < 82862755739295587 || (z[9] == 82862755739295587 && (z[8] < 18165857675053050549 || (z[8] == 18165857675053050549 && (z[7] < 12176845843281334983 || (z[7] == 12176845843281334983 && (z[6] < 5645674015335635503 || (z[6] == 5645674015335635503 && (z[5] < 9318693926755804304 || (z[5] == 9318693926755804304 && (z[4] < 13320134076191308873 || (z[4] == 13320134076191308873 && (z[3] < 9083347379620258823 || (z[3] == 9083347379620258823 && (z[2] < 15543556715411259941 || (z[2] == 15543556715411259941 && (z[1] < 4410884215886313276 || (z[1] == 4410884215886313276 && (z[0] < 15512955586897510413))))))))))))))))))) {
		var b uint64
		z[0], b = bits.Sub64(z[0], 15512955586897510413, 0)
		z[1], b = bits.Sub64(z[1], 4410884215886313276, b)
		z[2], b = bits.Sub64(z[2], 15543556715411259941, b)
		z[3], b = bits.Sub64(z[3], 9083347379620258823, b)
		z[4], b = bits.Sub64(z[4], 13320134076191308873, b)
		z[5], b = bits.Sub64(z[5], 9318693926755804304, b)
		z[6], b = bits.Sub64(z[6], 5645674015335635503, b)
		z[7], b = bits.Sub64(z[7], 12176845843281334983, b)
		z[8], b = bits.Sub64(z[8], 18165857675053050549, b)
		z[9], _ = bits.Sub64(z[9], 82862755739295587, b)
	}
}

func mulByConstant(z *Element, c uint8) {
	switch c {
	case 0:
		z.SetZero()
		return
	case 1:
		return
	case 2:
		z.Double(z)
		return
	case 3:
		_z := *z
		z.Double(z).Add(z, &_z)
	case 5:
		_z := *z
		z.Double(z).Double(z).Add(z, &_z)
	case 11:
		_z := *z
		z.Double(z).Double(z).Add(z, &_z).Double(z).Add(z, &_z)
	default:
		var y Element
		y.SetUint64(uint64(c))
		z.Mul(z, &y)
	}
}

// BatchInvert returns a new slice with every element inverted.
// Uses Montgomery batch inversion trick
func BatchInvert(a []Element) []Element {
	res := make([]Element, len(a))
	if len(a) == 0 {
		return res
	}

	zeroes := make([]bool, len(a))
	accumulator := One()

	for i := 0; i < len(a); i++ {
		if a[i].IsZero() {
			zeroes[i] = true
			continue
		}
		res[i] = accumulator
		accumulator.Mul(&accumulator, &a[i])
	}

	accumulator.Inverse(&accumulator)

	for i := len(a) - 1; i >= 0; i-- {
		if zeroes[i] {
			continue
		}
		res[i].Mul(&res[i], &accumulator)
		accumulator.Mul(&accumulator, &a[i])
	}

	return res
}

func _butterflyGeneric(a, b *Element) {
	t := *a
	a.Add(a, b)
	b.Sub(&t, b)
}

// BitLen returns the minimum number of bits needed to represent z
// returns 0 if z == 0
func (z *Element) BitLen() int {
	if z[9] != 0 {
		return 576 + bits.Len64(z[9])
	}
	if z[8] != 0 {
		return 512 + bits.Len64(z[8])
	}
	if z[7] != 0 {
		return 448 + bits.Len64(z[7])
	}
	if z[6] != 0 {
		return 384 + bits.Len64(z[6])
	}
	if z[5] != 0 {
		return 320 + bits.Len64(z[5])
	}
	if z[4] != 0 {
		return 256 + bits.Len64(z[4])
	}
	if z[3] != 0 {
		return 192 + bits.Len64(z[3])
	}
	if z[2] != 0 {
		return 128 + bits.Len64(z[2])
	}
	if z[1] != 0 {
		return 64 + bits.Len64(z[1])
	}
	return bits.Len64(z[0])
}

// Exp z = x^exponent mod q
func (z *Element) Exp(x Element, exponent *big.Int) *Element {
	var bZero big.Int
	if exponent.Cmp(&bZero) == 0 {
		return z.SetOne()
	}

	z.Set(&x)

	for i := exponent.BitLen() - 2; i >= 0; i-- {
		z.Square(z)
		if exponent.Bit(i) == 1 {
			z.Mul(z, &x)
		}
	}

	return z
}

// ToMont converts z to Montgomery form
// sets and returns z = z * r²
func (z *Element) ToMont() *Element {
	return z.Mul(z, &rSquare)
}

// ToRegular returns z in regular form (doesn't mutate z)
func (z Element) ToRegular() Element {
	return *z.FromMont()
}

// String returns the decimal representation of z as generated by
// z.Text(10).
func (z *Element) String() string {
	return z.Text(10)
}

// Text returns the string representation of z in the given base.
// Base must be between 2 and 36, inclusive. The result uses the
// lower-case letters 'a' to 'z' for digit values 10 to 35.
// No prefix (such as "0x") is added to the string. If z is a nil
// pointer it returns "<nil>".
// If base == 10 and -z fits in a uint64 prefix "-" is added to the string.
func (z *Element) Text(base int) string {
	if base < 2 || base > 36 {
		panic("invalid base")
	}
	if z == nil {
		return "<nil>"
	}
	zz := *z
	zz.FromMont()
	if zz.IsUint64() {
		return strconv.FormatUint(zz[0], base)
	} else if base == 10 {
		var zzNeg Element
		zzNeg.Neg(z)
		zzNeg.FromMont()
		if zzNeg.IsUint64() {
			return "-" + strconv.FormatUint(zzNeg[0], base)
		}
	}
	vv := bigIntPool.Get().(*big.Int)
	r := zz.ToBigInt(vv).Text(base)
	bigIntPool.Put(vv)
	return r
}

// ToBigInt returns z as a big.Int in Montgomery form
func (z *Element) ToBigInt(res *big.Int) *big.Int {
	var b [Limbs * 8]byte
	binary.BigEndian.PutUint64(b[72:80], z[0])
	binary.BigEndian.PutUint64(b[64:72], z[1])
	binary.BigEndian.PutUint64(b[56:64], z[2])
	binary.BigEndian.PutUint64(b[48:56], z[3])
	binary.BigEndian.PutUint64(b[40:48], z[4])
	binary.BigEndian.PutUint64(b[32:40], z[5])
	binary.BigEndian.PutUint64(b[24:32], z[6])
	binary.BigEndian.PutUint64(b[16:24], z[7])
	binary.BigEndian.PutUint64(b[8:16], z[8])
	binary.BigEndian.PutUint64(b[0:8], z[9])

	return res.SetBytes(b[:])
}

// ToBigIntRegular returns z as a big.Int in regular form
func (z Element) ToBigIntRegular(res *big.Int) *big.Int {
	z.FromMont()
	return z.ToBigInt(res)
}

// Bytes returns the regular (non montgomery) value
// of z as a big-endian byte array.
func (z *Element) Bytes() (res [Limbs * 8]byte) {
	_z := z.ToRegular()
	binary.BigEndian.PutUint64(res[72:80], _z[0])
	binary.BigEndian.PutUint64(res[64:72], _z[1])
	binary.BigEndian.PutUint64(res[56:64], _z[2])
	binary.BigEndian.PutUint64(res[48:56], _z[3])
	binary.BigEndian.PutUint64(res[40:48], _z[4])
	binary.BigEndian.PutUint64(res[32:40], _z[5])
	binary.BigEndian.PutUint64(res[24:32], _z[6])
	binary.BigEndian.PutUint64(res[16:24], _z[7])
	binary.BigEndian.PutUint64(res[8:16], _z[8])
	binary.BigEndian.PutUint64(res[0:8], _z[9])

	return
}

// Marshal returns the regular (non montgomery) value
// of z as a big-endian byte slice.
func (z *Element) Marshal() []byte {
	b := z.Bytes()
	return b[:]
}

// SetBytes interprets e as the bytes of a big-endian unsigned integer,
// sets z to that value (in Montgomery form), and returns z.
func (z *Element) SetBytes(e []byte) *Element {
	// get a big int from our pool
	vv := bigIntPool.Get().(*big.Int)
	vv.SetBytes(e)

	// set big int
	z.SetBigInt(vv)

	// put temporary object back in pool
	bigIntPool.Put(vv)

	return z
}

// SetBigInt sets z to v (regular form) and returns z in Montgomery form
func (z *Element) SetBigInt(v *big.Int) *Element {
	z.SetZero()

	var zero big.Int

	// fast path
	c := v.Cmp(&_modulus)
	if c == 0 {
		// v == 0
		return z
	} else if c != 1 && v.Cmp(&zero) != -1 {
		// 0 < v < q
		return z.setBigInt(v)
	}

	// get temporary big int from the pool
	vv := bigIntPool.Get().(*big.Int)

	// copy input + modular reduction
	vv.Set(v)
	vv.Mod(v, &_modulus)

	// set big int byte value
	z.setBigInt(vv)

	// release object into pool
	bigIntPool.Put(vv)
	return z
}

// setBigInt assumes 0 ⩽ v < q
func (z *Element) setBigInt(v *big.Int) *Element {
	vBits := v.Bits()

	if bits.UintSize == 64 {
		for i := 0; i < len(vBits); i++ {
			z[i] = uint64(vBits[i])
		}
	} else {
		for i := 0; i < len(vBits); i++ {
			if i%2 == 0 {
				z[i/2] = uint64(vBits[i])
			} else {
				z[i/2] |= uint64(vBits[i]) << 32
			}
		}
	}

	return z.ToMont()
}

// SetString creates a big.Int with number and calls SetBigInt on z
//
// The number prefix determines the actual base: A prefix of
// ''0b'' or ''0B'' selects base 2, ''0'', ''0o'' or ''0O'' selects base 8,
// and ''0x'' or ''0X'' selects base 16. Otherwise, the selected base is 10
// and no prefix is accepted.
//
// For base 16, lower and upper case letters are considered the same:
// The letters 'a' to 'f' and 'A' to 'F' represent digit values 10 to 15.
//
// An underscore character ''_'' may appear between a base
// prefix and an adjacent digit, and between successive digits; such
// underscores do not change the value of the number.
// Incorrect placement of underscores is reported as a panic if there
// are no other errors.
//
func (z *Element) SetString(number string) *Element {
	// get temporary big int from the pool
	vv := bigIntPool.Get().(*big.Int)

	if _, ok := vv.SetString(number, 0); !ok {
		panic("Element.SetString failed -> can't parse number into a big.Int " + number)
	}

	z.SetBigInt(vv)

	// release object into pool
	bigIntPool.Put(vv)

	return z
}

// MarshalJSON returns json encoding of z (z.Text(10))
// If z == nil, returns null
func (z *Element) MarshalJSON() ([]byte, error) {
	if z == nil {
		return []byte("null"), nil
	}
	const maxSafeBound = 15 // we encode it as number if it's small
	s := z.Text(10)
	if len(s) <= maxSafeBound {
		return []byte(s), nil
	}
	var sbb strings.Builder
	sbb.WriteByte('"')
	sbb.WriteString(s)
	sbb.WriteByte('"')
	return []byte(sbb.String()), nil
}

// UnmarshalJSON accepts numbers and strings as input
// See Element.SetString for valid prefixes (0x, 0b, ...)
func (z *Element) UnmarshalJSON(data []byte) error {
	s := string(data)
	if len(s) > Bits*3 {
		return errors.New("value too large (max = Element.Bits * 3)")
	}

	// we accept numbers and strings, remove leading and trailing quotes if any
	if len(s) > 0 && s[0] == '"' {
		s = s[1:]
	}
	if len(s) > 0 && s[len(s)-1] == '"' {
		s = s[:len(s)-1]
	}

	// get temporary big int from the pool
	vv := bigIntPool.Get().(*big.Int)

	if _, ok := vv.SetString(s, 0); !ok {
		return errors.New("can't parse into a big.Int: " + s)
	}

	z.SetBigInt(vv)

	// release object into pool
	bigIntPool.Put(vv)
	return nil
}

// Legendre returns the Legendre symbol of z (either +1, -1, or 0.)
func (z *Element) Legendre() int {
	var l Element
	// z^((q-1)/2)
	l.expByLegendreExp(*z)

	if l.IsZero() {
		return 0
	}

	// if l == 1
	if (l[9] == 51212299585931083) && (l[8] == 7016548280614581879) && (l[7] == 8411601626847721258) && (l[6] == 1038965607738428109) && (l[5] == 15732028589390776959) && (l[4] == 12856030952767240260) && (l[3] == 12638729832353218866) && (l[2] == 17318295036095996852) && (l[1] == 16907884053554239805) && (l[0] == 5665001492438840506) {
		return 1
	}
	return -1
}

// Sqrt z = √x mod q
// if the square root doesn't exist (x is not a square mod q)
// Sqrt leaves z unchanged and returns nil
func (z *Element) Sqrt(x *Element) *Element {
	// q ≡ 5 (mod 8)
	// see modSqrt5Mod8Prime in math/big/int.go
	var one, alpha, beta, tx, square Element
	one.SetOne()
	tx.Double(x)
	alpha.expBySqrtExp(tx)

	beta.Square(&alpha).
		Mul(&beta, &tx).
		Sub(&beta, &one).
		Mul(&beta, x).
		Mul(&beta, &alpha)

	// as we didn't compute the legendre symbol, ensure we found beta such that beta * beta = x
	square.Square(&beta)
	if square.Equal(x) {
		return z.Set(&beta)
	}
	return nil
}

func max(a int, b int) int {
	if a > b {
		return a
	}
	return b
}

func min(a int, b int) int {
	if a < b {
		return a
	}
	return b
}

const updateFactorsConversionBias int64 = 0x7fffffff7fffffff // (2³¹ - 1)(2³² + 1)
const updateFactorIdentityMatrixRow0 = 1
const updateFactorIdentityMatrixRow1 = 1 << 32

func updateFactorsDecompose(c int64) (int64, int64) {
	c += updateFactorsConversionBias
	const low32BitsFilter int64 = 0xFFFFFFFF
	f := c&low32BitsFilter - 0x7FFFFFFF
	g := c>>32&low32BitsFilter - 0x7FFFFFFF
	return f, g
}

const k = 32 // word size / 2
const signBitSelector = uint64(1) << 63
const approxLowBitsN = k - 1
const approxHighBitsN = k + 1
const inversionCorrectionFactorWord0 = 17335095338408674528
const inversionCorrectionFactorWord1 = 1935156146725576072
const inversionCorrectionFactorWord2 = 12310223143035529855
const inversionCorrectionFactorWord3 = 14776388015283991997
const inversionCorrectionFactorWord4 = 13807356859349388480
const inversionCorrectionFactorWord5 = 10412247811534140886
const inversionCorrectionFactorWord6 = 1537112855455741892
const inversionCorrectionFactorWord7 = 5281081904757642912
const inversionCorrectionFactorWord8 = 14734303888675989218
const inversionCorrectionFactorWord9 = 64202171737444348

const invIterationsN = 42

// Inverse z = x⁻¹ mod q
// Implements "Optimized Binary GCD for Modular Inversion"
// https://github.com/pornin/bingcd/blob/main/doc/bingcd.pdf
func (z *Element) Inverse(x *Element) *Element {
	a := *x
	b := Element{
		qElementWord0,
		qElementWord1,
		qElementWord2,
		qElementWord3,
		qElementWord4,
		qElementWord5,
		qElementWord6,
		qElementWord7,
		qElementWord8,
		qElementWord9,
	} // b := q

	u := Element{1}

	// Update factors: we get [u; v]:= [f0 g0; f1 g1] [u; v]
	// c_i = f_i + 2³¹ - 1 + 2³² * (g_i + 2³¹ - 1)
	var c0, c1 int64

	// Saved update factors to reduce the number of field multiplications
	var pf0, pf1, pg0, pg1 int64

	var i uint

	var v, s Element

	// Since u,v are updated every other iteration, we must make sure we terminate after evenly many iterations
	// This also lets us get away with half as many updates to u,v
	// To make this constant-time-ish, replace the condition with i < invIterationsN
	for i = 0; i&1 == 1 || !a.IsZero(); i++ {
		n := max(a.BitLen(), b.BitLen())
		aApprox, bApprox := approximate(&a, n), approximate(&b, n)

		// After 0 iterations, we have f₀ ≤ 2⁰ and f₁ < 2⁰
		// f0, g0, f1, g1 = 1, 0, 0, 1
		c0, c1 = updateFactorIdentityMatrixRow0, updateFactorIdentityMatrixRow1

		for j := 0; j < approxLowBitsN; j++ {

			if aApprox&1 == 0 {
				aApprox /= 2
			} else {
				s, borrow := bits.Sub64(aApprox, bApprox, 0)
				if borrow == 1 {
					s = bApprox - aApprox
					bApprox = aApprox
					c0, c1 = c1, c0
				}

				aApprox = s / 2
				c0 = c0 - c1

				// Now |f₀| < 2ʲ + 2ʲ = 2ʲ⁺¹
				// |f₁| ≤ 2ʲ still
			}

			c1 *= 2
			// |f₁| ≤ 2ʲ⁺¹
		}

		s = a

		var g0 int64
		// from this point on c0 aliases for f0
		c0, g0 = updateFactorsDecompose(c0)
		aHi := a.linearCombNonModular(&s, c0, &b, g0)
		if aHi&signBitSelector != 0 {
			// if aHi < 0
			c0, g0 = -c0, -g0
			aHi = a.neg(&a, aHi)
		}
		// right-shift a by k-1 bits
		a[0] = (a[0] >> approxLowBitsN) | ((a[1]) << approxHighBitsN)
		a[1] = (a[1] >> approxLowBitsN) | ((a[2]) << approxHighBitsN)
		a[2] = (a[2] >> approxLowBitsN) | ((a[3]) << approxHighBitsN)
		a[3] = (a[3] >> approxLowBitsN) | ((a[4]) << approxHighBitsN)
		a[4] = (a[4] >> approxLowBitsN) | ((a[5]) << approxHighBitsN)
		a[5] = (a[5] >> approxLowBitsN) | ((a[6]) << approxHighBitsN)
		a[6] = (a[6] >> approxLowBitsN) | ((a[7]) << approxHighBitsN)
		a[7] = (a[7] >> approxLowBitsN) | ((a[8]) << approxHighBitsN)
		a[8] = (a[8] >> approxLowBitsN) | ((a[9]) << approxHighBitsN)
		a[9] = (a[9] >> approxLowBitsN) | (aHi << approxHighBitsN)

		var f1 int64
		// from this point on c1 aliases for g0
		f1, c1 = updateFactorsDecompose(c1)
		bHi := b.linearCombNonModular(&s, f1, &b, c1)
		if bHi&signBitSelector != 0 {
			// if bHi < 0
			f1, c1 = -f1, -c1
			bHi = b.neg(&b, bHi)
		}
		// right-shift b by k-1 bits
		b[0] = (b[0] >> approxLowBitsN) | ((b[1]) << approxHighBitsN)
		b[1] = (b[1] >> approxLowBitsN) | ((b[2]) << approxHighBitsN)
		b[2] = (b[2] >> approxLowBitsN) | ((b[3]) << approxHighBitsN)
		b[3] = (b[3] >> approxLowBitsN) | ((b[4]) << approxHighBitsN)
		b[4] = (b[4] >> approxLowBitsN) | ((b[5]) << approxHighBitsN)
		b[5] = (b[5] >> approxLowBitsN) | ((b[6]) << approxHighBitsN)
		b[6] = (b[6] >> approxLowBitsN) | ((b[7]) << approxHighBitsN)
		b[7] = (b[7] >> approxLowBitsN) | ((b[8]) << approxHighBitsN)
		b[8] = (b[8] >> approxLowBitsN) | ((b[9]) << approxHighBitsN)
		b[9] = (b[9] >> approxLowBitsN) | (bHi << approxHighBitsN)

		if i&1 == 1 {
			// Combine current update factors with previously stored ones
			// [f₀, g₀; f₁, g₁] ← [f₀, g₀; f₁, g₀] [pf₀, pg₀; pf₀, pg₀]
			// We have |f₀|, |g₀|, |pf₀|, |pf₁| ≤ 2ᵏ⁻¹, and that |pf_i| < 2ᵏ⁻¹ for i ∈ {0, 1}
			// Then for the new value we get |f₀| < 2ᵏ⁻¹ × 2ᵏ⁻¹ + 2ᵏ⁻¹ × 2ᵏ⁻¹ = 2²ᵏ⁻¹
			// Which leaves us with an extra bit for the sign

			// c0 aliases f0, c1 aliases g1
			c0, g0, f1, c1 = c0*pf0+g0*pf1,
				c0*pg0+g0*pg1,
				f1*pf0+c1*pf1,
				f1*pg0+c1*pg1

			s = u
			u.linearCombSosSigned(&u, c0, &v, g0)
			v.linearCombSosSigned(&s, f1, &v, c1)

		} else {
			// Save update factors
			pf0, pg0, pf1, pg1 = c0, g0, f1, c1
		}
	}

	// For every iteration that we miss, v is not being multiplied by 2²ᵏ⁻²
	const pSq int64 = 1 << (2 * (k - 1))
	// If the function is constant-time ish, this loop will not run (probably no need to take it out explicitly)
	for ; i < invIterationsN; i += 2 {
		v.mulWSigned(&v, pSq)
	}

	z.Mul(&v, &Element{
		inversionCorrectionFactorWord0,
		inversionCorrectionFactorWord1,
		inversionCorrectionFactorWord2,
		inversionCorrectionFactorWord3,
		inversionCorrectionFactorWord4,
		inversionCorrectionFactorWord5,
		inversionCorrectionFactorWord6,
		inversionCorrectionFactorWord7,
		inversionCorrectionFactorWord8,
		inversionCorrectionFactorWord9,
	})
	return z
}

// approximate a big number x into a single 64 bit word using its uppermost and lowermost bits
// if x fits in a word as is, no approximation necessary
func approximate(x *Element, nBits int) uint64 {

	if nBits <= 64 {
		return x[0]
	}

	const mask = (uint64(1) << (k - 1)) - 1 // k-1 ones
	lo := mask & x[0]

	hiWordIndex := (nBits - 1) / 64

	hiWordBitsAvailable := nBits - hiWordIndex*64
	hiWordBitsUsed := min(hiWordBitsAvailable, approxHighBitsN)

	mask_ := uint64(^((1 << (hiWordBitsAvailable - hiWordBitsUsed)) - 1))
	hi := (x[hiWordIndex] & mask_) << (64 - hiWordBitsAvailable)

	mask_ = ^(1<<(approxLowBitsN+hiWordBitsUsed) - 1)
	mid := (mask_ & x[hiWordIndex-1]) >> hiWordBitsUsed

	return lo | mid | hi
}

func (z *Element) linearCombSosSigned(x *Element, xC int64, y *Element, yC int64) {
	hi := z.linearCombNonModular(x, xC, y, yC)
	z.montReduceSigned(z, hi)
}

// montReduceSigned SOS algorithm; xHi must be at most 63 bits long. Last bit of xHi may be used as a sign bit
func (z *Element) montReduceSigned(x *Element, xHi uint64) {

	const signBitRemover = ^signBitSelector
	neg := xHi&signBitSelector != 0
	// the SOS implementation requires that most significant bit is 0
	// Let X be xHi*r + x
	// note that if X is negative we would have initially stored it as 2⁶⁴ r + X
	xHi &= signBitRemover
	// with this a negative X is now represented as 2⁶³ r + X

	var t [2*Limbs - 1]uint64
	var C uint64

	m := x[0] * qInvNegLsw

	C = madd0(m, qElementWord0, x[0])
	C, t[1] = madd2(m, qElementWord1, x[1], C)
	C, t[2] = madd2(m, qElementWord2, x[2], C)
	C, t[3] = madd2(m, qElementWord3, x[3], C)
	C, t[4] = madd2(m, qElementWord4, x[4], C)
	C, t[5] = madd2(m, qElementWord5, x[5], C)
	C, t[6] = madd2(m, qElementWord6, x[6], C)
	C, t[7] = madd2(m, qElementWord7, x[7], C)
	C, t[8] = madd2(m, qElementWord8, x[8], C)
	C, t[9] = madd2(m, qElementWord9, x[9], C)

	// the high word of m * qElement[9] is at most 62 bits
	// x[9] + C is at most 65 bits (high word at most 1 bit)
	// Thus the resulting C will be at most 63 bits
	t[10] = xHi + C
	// xHi and C are 63 bits, therefore no overflow

	{
		const i = 1
		m = t[i] * qInvNegLsw

		C = madd0(m, qElementWord0, t[i+0])
		C, t[i+1] = madd2(m, qElementWord1, t[i+1], C)
		C, t[i+2] = madd2(m, qElementWord2, t[i+2], C)
		C, t[i+3] = madd2(m, qElementWord3, t[i+3], C)
		C, t[i+4] = madd2(m, qElementWord4, t[i+4], C)
		C, t[i+5] = madd2(m, qElementWord5, t[i+5], C)
		C, t[i+6] = madd2(m, qElementWord6, t[i+6], C)
		C, t[i+7] = madd2(m, qElementWord7, t[i+7], C)
		C, t[i+8] = madd2(m, qElementWord8, t[i+8], C)
		C, t[i+9] = madd2(m, qElementWord9, t[i+9], C)

		t[i+Limbs] += C
	}
	{
		const i = 2
		m = t[i] * qInvNegLsw

		C = madd0(m, qElementWord0, t[i+0])
		C, t[i+1] = madd2(m, qElementWord1, t[i+1], C)
		C, t[i+2] = madd2(m, qElementWord2, t[i+2], C)
		C, t[i+3] = madd2(m, qElementWord3, t[i+3], C)
		C, t[i+4] = madd2(m, qElementWord4, t[i+4], C)
		C, t[i+5] = madd2(m, qElementWord5, t[i+5], C)
		C, t[i+6] = madd2(m, qElementWord6, t[i+6], C)
		C, t[i+7] = madd2(m, qElementWord7, t[i+7], C)
		C, t[i+8] = madd2(m, qElementWord8, t[i+8], C)
		C, t[i+9] = madd2(m, qElementWord9, t[i+9], C)

		t[i+Limbs] += C
	}
	{
		const i = 3
		m = t[i] * qInvNegLsw

		C = madd0(m, qElementWord0, t[i+0])
		C, t[i+1] = madd2(m, qElementWord1, t[i+1], C)
		C, t[i+2] = madd2(m, qElementWord2, t[i+2], C)
		C, t[i+3] = madd2(m, qElementWord3, t[i+3], C)
		C, t[i+4] = madd2(m, qElementWord4, t[i+4], C)
		C, t[i+5] = madd2(m, qElementWord5, t[i+5], C)
		C, t[i+6] = madd2(m, qElementWord6, t[i+6], C)
		C, t[i+7] = madd2(m, qElementWord7, t[i+7], C)
		C, t[i+8] = madd2(m, qElementWord8, t[i+8], C)
		C, t[i+9] = madd2(m, qElementWord9, t[i+9], C)

		t[i+Limbs] += C
	}
	{
		const i = 4
		m = t[i] * qInvNegLsw

		C = madd0(m, qElementWord0, t[i+0])
		C, t[i+1] = madd2(m, qElementWord1, t[i+1], C)
		C, t[i+2] = madd2(m, qElementWord2, t[i+2], C)
		C, t[i+3] = madd2(m, qElementWord3, t[i+3], C)
		C, t[i+4] = madd2(m, qElementWord4, t[i+4], C)
		C, t[i+5] = madd2(m, qElementWord5, t[i+5], C)
		C, t[i+6] = madd2(m, qElementWord6, t[i+6], C)
		C, t[i+7] = madd2(m, qElementWord7, t[i+7], C)
		C, t[i+8] = madd2(m, qElementWord8, t[i+8], C)
		C, t[i+9] = madd2(m, qElementWord9, t[i+9], C)

		t[i+Limbs] += C
	}
	{
		const i = 5
		m = t[i] * qInvNegLsw

		C = madd0(m, qElementWord0, t[i+0])
		C, t[i+1] = madd2(m, qElementWord1, t[i+1], C)
		C, t[i+2] = madd2(m, qElementWord2, t[i+2], C)
		C, t[i+3] = madd2(m, qElementWord3, t[i+3], C)
		C, t[i+4] = madd2(m, qElementWord4, t[i+4], C)
		C, t[i+5] = madd2(m, qElementWord5, t[i+5], C)
		C, t[i+6] = madd2(m, qElementWord6, t[i+6], C)
		C, t[i+7] = madd2(m, qElementWord7, t[i+7], C)
		C, t[i+8] = madd2(m, qElementWord8, t[i+8], C)
		C, t[i+9] = madd2(m, qElementWord9, t[i+9], C)

		t[i+Limbs] += C
	}
	{
		const i = 6
		m = t[i] * qInvNegLsw

		C = madd0(m, qElementWord0, t[i+0])
		C, t[i+1] = madd2(m, qElementWord1, t[i+1], C)
		C, t[i+2] = madd2(m, qElementWord2, t[i+2], C)
		C, t[i+3] = madd2(m, qElementWord3, t[i+3], C)
		C, t[i+4] = madd2(m, qElementWord4, t[i+4], C)
		C, t[i+5] = madd2(m, qElementWord5, t[i+5], C)
		C, t[i+6] = madd2(m, qElementWord6, t[i+6], C)
		C, t[i+7] = madd2(m, qElementWord7, t[i+7], C)
		C, t[i+8] = madd2(m, qElementWord8, t[i+8], C)
		C, t[i+9] = madd2(m, qElementWord9, t[i+9], C)

		t[i+Limbs] += C
	}
	{
		const i = 7
		m = t[i] * qInvNegLsw

		C = madd0(m, qElementWord0, t[i+0])
		C, t[i+1] = madd2(m, qElementWord1, t[i+1], C)
		C, t[i+2] = madd2(m, qElementWord2, t[i+2], C)
		C, t[i+3] = madd2(m, qElementWord3, t[i+3], C)
		C, t[i+4] = madd2(m, qElementWord4, t[i+4], C)
		C, t[i+5] = madd2(m, qElementWord5, t[i+5], C)
		C, t[i+6] = madd2(m, qElementWord6, t[i+6], C)
		C, t[i+7] = madd2(m, qElementWord7, t[i+7], C)
		C, t[i+8] = madd2(m, qElementWord8, t[i+8], C)
		C, t[i+9] = madd2(m, qElementWord9, t[i+9], C)

		t[i+Limbs] += C
	}
	{
		const i = 8
		m = t[i] * qInvNegLsw

		C = madd0(m, qElementWord0, t[i+0])
		C, t[i+1] = madd2(m, qElementWord1, t[i+1], C)
		C, t[i+2] = madd2(m, qElementWord2, t[i+2], C)
		C, t[i+3] = madd2(m, qElementWord3, t[i+3], C)
		C, t[i+4] = madd2(m, qElementWord4, t[i+4], C)
		C, t[i+5] = madd2(m, qElementWord5, t[i+5], C)
		C, t[i+6] = madd2(m, qElementWord6, t[i+6], C)
		C, t[i+7] = madd2(m, qElementWord7, t[i+7], C)
		C, t[i+8] = madd2(m, qElementWord8, t[i+8], C)
		C, t[i+9] = madd2(m, qElementWord9, t[i+9], C)

		t[i+Limbs] += C
	}
	{
		const i = 9
		m := t[i] * qInvNegLsw

		C = madd0(m, qElementWord0, t[i+0])
		C, z[0] = madd2(m, qElementWord1, t[i+1], C)
		C, z[1] = madd2(m, qElementWord2, t[i+2], C)
		C, z[2] = madd2(m, qElementWord3, t[i+3], C)
		C, z[3] = madd2(m, qElementWord4, t[i+4], C)
		C, z[4] = madd2(m, qElementWord5, t[i+5], C)
		C, z[5] = madd2(m, qElementWord6, t[i+6], C)
		C, z[6] = madd2(m, qElementWord7, t[i+7], C)
		C, z[7] = madd2(m, qElementWord8, t[i+8], C)
		z[9], z[8] = madd2(m, qElementWord9, t[i+9], C)
	}

	// if z > q → z -= q
	// note: this is NOT constant time
	if !(z[9] < 82862755739295587 || (z[9] == 82862755739295587 && (z[8] < 18165857675053050549 || (z[8] == 18165857675053050549 && (z[7] < 12176845843281334983 || (z[7] == 12176845843281334983 && (z[6] < 5645674015335635503 || (z[6] == 5645674015335635503 && (z[5] < 9318693926755804304 || (z[5] == 9318693926755804304 && (z[4] < 13320134076191308873 || (z[4] == 13320134076191308873 && (z[3] < 9083347379620258823 || (z[3] == 9083347379620258823 && (z[2] < 15543556715411259941 || (z[2] == 15543556715411259941 && (z[1] < 4410884215886313276 || (z[1] == 4410884215886313276 && (z[0] < 15512955586897510413))))))))))))))))))) {
		var b uint64
		z[0], b = bits.Sub64(z[0], 15512955586897510413, 0)
		z[1], b = bits.Sub64(z[1], 4410884215886313276, b)
		z[2], b = bits.Sub64(z[2], 15543556715411259941, b)
		z[3], b = bits.Sub64(z[3], 9083347379620258823, b)
		z[4], b = bits.Sub64(z[4], 13320134076191308873, b)
		z[5], b = bits.Sub64(z[5], 9318693926755804304, b)
		z[6], b = bits.Sub64(z[6], 5645674015335635503, b)
		z[7], b = bits.Sub64(z[7], 12176845843281334983, b)
		z[8], b = bits.Sub64(z[8], 18165857675053050549, b)
		z[9], _ = bits.Sub64(z[9], 82862755739295587, b)
	}
	if neg {
		// We have computed ( 2⁶³ r + X ) r⁻¹ = 2⁶³ + X r⁻¹ instead
		var b uint64
		z[0], b = bits.Sub64(z[0], signBitSelector, 0)
		z[1], b = bits.Sub64(z[1], 0, b)
		z[2], b = bits.Sub64(z[2], 0, b)
		z[3], b = bits.Sub64(z[3], 0, b)
		z[4], b = bits.Sub64(z[4], 0, b)
		z[5], b = bits.Sub64(z[5], 0, b)
		z[6], b = bits.Sub64(z[6], 0, b)
		z[7], b = bits.Sub64(z[7], 0, b)
		z[8], b = bits.Sub64(z[8], 0, b)
		z[9], b = bits.Sub64(z[9], 0, b)

		// Occurs iff x == 0 && xHi < 0, i.e. X = rX' for -2⁶³ ≤ X' < 0
		if b != 0 {
			// z[9] = -1
			// negative: add q
			const neg1 = 0xFFFFFFFFFFFFFFFF

			b = 0
			z[0], b = bits.Add64(z[0], qElementWord0, b)
			z[1], b = bits.Add64(z[1], qElementWord1, b)
			z[2], b = bits.Add64(z[2], qElementWord2, b)
			z[3], b = bits.Add64(z[3], qElementWord3, b)
			z[4], b = bits.Add64(z[4], qElementWord4, b)
			z[5], b = bits.Add64(z[5], qElementWord5, b)
			z[6], b = bits.Add64(z[6], qElementWord6, b)
			z[7], b = bits.Add64(z[7], qElementWord7, b)
			z[8], b = bits.Add64(z[8], qElementWord8, b)
			z[9], _ = bits.Add64(neg1, qElementWord9, b)
		}
	}
}

// mulWSigned mul word signed (w/ montgomery reduction)
func (z *Element) mulWSigned(x *Element, y int64) {
	m := y >> 63
	_mulWGeneric(z, x, uint64((y^m)-m))
	// multiply by abs(y)
	if y < 0 {
		z.Neg(z)
	}
}

func (z *Element) neg(x *Element, xHi uint64) uint64 {
	var b uint64

	z[0], b = bits.Sub64(0, x[0], 0)
	z[1], b = bits.Sub64(0, x[1], b)
	z[2], b = bits.Sub64(0, x[2], b)
	z[3], b = bits.Sub64(0, x[3], b)
	z[4], b = bits.Sub64(0, x[4], b)
	z[5], b = bits.Sub64(0, x[5], b)
	z[6], b = bits.Sub64(0, x[6], b)
	z[7], b = bits.Sub64(0, x[7], b)
	z[8], b = bits.Sub64(0, x[8], b)
	z[9], b = bits.Sub64(0, x[9], b)
	xHi, _ = bits.Sub64(0, xHi, b)

	return xHi
}

// regular multiplication by one word regular (non montgomery)
// Fewer additions than the branch-free for positive y. Could be faster on some architectures
func (z *Element) mulWRegular(x *Element, y int64) uint64 {

	// w := abs(y)
	m := y >> 63
	w := uint64((y ^ m) - m)

	var c uint64
	c, z[0] = bits.Mul64(x[0], w)
	c, z[1] = madd1(x[1], w, c)
	c, z[2] = madd1(x[2], w, c)
	c, z[3] = madd1(x[3], w, c)
	c, z[4] = madd1(x[4], w, c)
	c, z[5] = madd1(x[5], w, c)
	c, z[6] = madd1(x[6], w, c)
	c, z[7] = madd1(x[7], w, c)
	c, z[8] = madd1(x[8], w, c)
	c, z[9] = madd1(x[9], w, c)

	if y < 0 {
		c = z.neg(z, c)
	}

	return c
}

/*
Removed: seems slower
// mulWRegular branch-free regular multiplication by one word (non montgomery)
func (z *Element) mulWRegularBf(x *Element, y int64) uint64 {

	w := uint64(y)
	allNeg := uint64(y >> 63)	// -1 if y < 0, 0 o.w

	// s[0], s[1] so results are not stored immediately in z.
	// x[i] will be needed in the i+1 th iteration. We don't want to overwrite it in case x = z
	var s [2]uint64
	var h [2]uint64

	h[0], s[0] = bits.Mul64(x[0], w)

	c := uint64(0)
	b := uint64(0)

		{
			const curI = 1 % 2
			const prevI = 1 - curI
			const iMinusOne = 1 - 1

			h[curI], s[curI] = bits.Mul64(x[1], w)
			s[curI], c = bits.Add64(s[curI], h[prevI], c)
			s[curI], b = bits.Sub64(s[curI], allNeg & x[iMinusOne], b)
			z[iMinusOne] = s[prevI]
		}

		{
			const curI = 2 % 2
			const prevI = 1 - curI
			const iMinusOne = 2 - 1

			h[curI], s[curI] = bits.Mul64(x[2], w)
			s[curI], c = bits.Add64(s[curI], h[prevI], c)
			s[curI], b = bits.Sub64(s[curI], allNeg & x[iMinusOne], b)
			z[iMinusOne] = s[prevI]
		}

		{
			const curI = 3 % 2
			const prevI = 1 - curI
			const iMinusOne = 3 - 1

			h[curI], s[curI] = bits.Mul64(x[3], w)
			s[curI], c = bits.Add64(s[curI], h[prevI], c)
			s[curI], b = bits.Sub64(s[curI], allNeg & x[iMinusOne], b)
			z[iMinusOne] = s[prevI]
		}

		{
			const curI = 4 % 2
			const prevI = 1 - curI
			const iMinusOne = 4 - 1

			h[curI], s[curI] = bits.Mul64(x[4], w)
			s[curI], c = bits.Add64(s[curI], h[prevI], c)
			s[curI], b = bits.Sub64(s[curI], allNeg & x[iMinusOne], b)
			z[iMinusOne] = s[prevI]
		}

		{
			const curI = 5 % 2
			const prevI = 1 - curI
			const iMinusOne = 5 - 1

			h[curI], s[curI] = bits.Mul64(x[5], w)
			s[curI], c = bits.Add64(s[curI], h[prevI], c)
			s[curI], b = bits.Sub64(s[curI], allNeg & x[iMinusOne], b)
			z[iMinusOne] = s[prevI]
		}

		{
			const curI = 6 % 2
			const prevI = 1 - curI
			const iMinusOne = 6 - 1

			h[curI], s[curI] = bits.Mul64(x[6], w)
			s[curI], c = bits.Add64(s[curI], h[prevI], c)
			s[curI], b = bits.Sub64(s[curI], allNeg & x[iMinusOne], b)
			z[iMinusOne] = s[prevI]
		}

		{
			const curI = 7 % 2
			const prevI = 1 - curI
			const iMinusOne = 7 - 1

			h[curI], s[curI] = bits.Mul64(x[7], w)
			s[curI], c = bits.Add64(s[curI], h[prevI], c)
			s[curI], b = bits.Sub64(s[curI], allNeg & x[iMinusOne], b)
			z[iMinusOne] = s[prevI]
		}

		{
			const curI = 8 % 2
			const prevI = 1 - curI
			const iMinusOne = 8 - 1

			h[curI], s[curI] = bits.Mul64(x[8], w)
			s[curI], c = bits.Add64(s[curI], h[prevI], c)
			s[curI], b = bits.Sub64(s[curI], allNeg & x[iMinusOne], b)
			z[iMinusOne] = s[prevI]
		}

		{
			const curI = 9 % 2
			const prevI = 1 - curI
			const iMinusOne = 9 - 1

			h[curI], s[curI] = bits.Mul64(x[9], w)
			s[curI], c = bits.Add64(s[curI], h[prevI], c)
			s[curI], b = bits.Sub64(s[curI], allNeg & x[iMinusOne], b)
			z[iMinusOne] = s[prevI]
		}
	{
		const curI = 10 % 2
		const prevI = 1 - curI
		const iMinusOne = 9

		s[curI], _ = bits.Sub64(h[prevI], allNeg & x[iMinusOne], b)
		z[iMinusOne] = s[prevI]

		return s[curI] + c
	}
}*/

// Requires NoCarry
func (z *Element) linearCombNonModular(x *Element, xC int64, y *Element, yC int64) uint64 {
	var yTimes Element

	yHi := yTimes.mulWRegular(y, yC)
	xHi := z.mulWRegular(x, xC)

	carry := uint64(0)
	z[0], carry = bits.Add64(z[0], yTimes[0], carry)
	z[1], carry = bits.Add64(z[1], yTimes[1], carry)
	z[2], carry = bits.Add64(z[2], yTimes[2], carry)
	z[3], carry = bits.Add64(z[3], yTimes[3], carry)
	z[4], carry = bits.Add64(z[4], yTimes[4], carry)
	z[5], carry = bits.Add64(z[5], yTimes[5], carry)
	z[6], carry = bits.Add64(z[6], yTimes[6], carry)
	z[7], carry = bits.Add64(z[7], yTimes[7], carry)
	z[8], carry = bits.Add64(z[8], yTimes[8], carry)
	z[9], carry = bits.Add64(z[9], yTimes[9], carry)

	yHi, _ = bits.Add64(xHi, yHi, carry)

	return yHi
}

func (z *Element) EvalPolynomial(monic bool, coefficients []Element, x *Element) {
	dst := coefficients[len(coefficients)-1]

	if monic {
		dst.Add(&dst, x)
	}

	for i := len(coefficients) - 2; i >= 0; i-- {
		dst.Mul(&dst, x)
		dst.Add(&dst, &coefficients[i])
	}

	*z = dst
}