1 // Copyright 2009 The Go Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
5 // This file implements signed multi-precision integers.
16 // An Int represents a signed multi-precision integer.
17 // The zero value for an Int represents the value 0.
19 // Operations always take pointer arguments (*Int) rather
20 // than Int values, and each unique Int value requires
21 // its own unique *Int pointer. To "copy" an Int value,
22 // an existing (or newly allocated) Int must be set to
23 // a new value using the Int.Set method; shallow copies
24 // of Ints are not supported and may lead to errors.
27 abs nat
// absolute value of the integer
30 var intOne
= &Int
{false, natOne
}
38 func (x
*Int
) Sign() int {
48 // SetInt64 sets z to x and returns z.
49 func (z
*Int
) SetInt64(x
int64) *Int
{
55 z
.abs
= z
.abs
.setUint64(uint64(x
))
60 // SetUint64 sets z to x and returns z.
61 func (z
*Int
) SetUint64(x
uint64) *Int
{
62 z
.abs
= z
.abs
.setUint64(x
)
67 // NewInt allocates and returns a new Int set to x.
68 func NewInt(x
int64) *Int
{
69 return new(Int
).SetInt64(x
)
72 // Set sets z to x and returns z.
73 func (z
*Int
) Set(x
*Int
) *Int
{
75 z
.abs
= z
.abs
.set(x
.abs
)
81 // Bits provides raw (unchecked but fast) access to x by returning its
82 // absolute value as a little-endian Word slice. The result and x share
83 // the same underlying array.
84 // Bits is intended to support implementation of missing low-level Int
85 // functionality outside this package; it should be avoided otherwise.
86 func (x
*Int
) Bits() []Word
{
90 // SetBits provides raw (unchecked but fast) access to z by setting its
91 // value to abs, interpreted as a little-endian Word slice, and returning
92 // z. The result and abs share the same underlying array.
93 // SetBits is intended to support implementation of missing low-level Int
94 // functionality outside this package; it should be avoided otherwise.
95 func (z
*Int
) SetBits(abs
[]Word
) *Int
{
96 z
.abs
= nat(abs
).norm()
101 // Abs sets z to |x| (the absolute value of x) and returns z.
102 func (z
*Int
) Abs(x
*Int
) *Int
{
108 // Neg sets z to -x and returns z.
109 func (z
*Int
) Neg(x
*Int
) *Int
{
111 z
.neg
= len(z
.abs
) > 0 && !z
.neg
// 0 has no sign
115 // Add sets z to the sum x+y and returns z.
116 func (z
*Int
) Add(x
, y
*Int
) *Int
{
120 // (-x) + (-y) == -(x + y)
121 z
.abs
= z
.abs
.add(x
.abs
, y
.abs
)
123 // x + (-y) == x - y == -(y - x)
124 // (-x) + y == y - x == -(x - y)
125 if x
.abs
.cmp(y
.abs
) >= 0 {
126 z
.abs
= z
.abs
.sub(x
.abs
, y
.abs
)
129 z
.abs
= z
.abs
.sub(y
.abs
, x
.abs
)
132 z
.neg
= len(z
.abs
) > 0 && neg
// 0 has no sign
136 // Sub sets z to the difference x-y and returns z.
137 func (z
*Int
) Sub(x
, y
*Int
) *Int
{
141 // (-x) - y == -(x + y)
142 z
.abs
= z
.abs
.add(x
.abs
, y
.abs
)
144 // x - y == x - y == -(y - x)
145 // (-x) - (-y) == y - x == -(x - y)
146 if x
.abs
.cmp(y
.abs
) >= 0 {
147 z
.abs
= z
.abs
.sub(x
.abs
, y
.abs
)
150 z
.abs
= z
.abs
.sub(y
.abs
, x
.abs
)
153 z
.neg
= len(z
.abs
) > 0 && neg
// 0 has no sign
157 // Mul sets z to the product x*y and returns z.
158 func (z
*Int
) Mul(x
, y
*Int
) *Int
{
160 // x * (-y) == -(x * y)
161 // (-x) * y == -(x * y)
162 // (-x) * (-y) == x * y
164 z
.abs
= z
.abs
.sqr(x
.abs
)
168 z
.abs
= z
.abs
.mul(x
.abs
, y
.abs
)
169 z
.neg
= len(z
.abs
) > 0 && x
.neg
!= y
.neg
// 0 has no sign
173 // MulRange sets z to the product of all integers
174 // in the range [a, b] inclusively and returns z.
175 // If a > b (empty range), the result is 1.
176 func (z
*Int
) MulRange(a
, b
int64) *Int
{
179 return z
.SetInt64(1) // empty range
180 case a
<= 0 && b
>= 0:
181 return z
.SetInt64(0) // range includes 0
183 // a <= b && (b < 0 || a > 0)
191 z
.abs
= z
.abs
.mulRange(uint64(a
), uint64(b
))
196 // Binomial sets z to the binomial coefficient of (n, k) and returns z.
197 func (z
*Int
) Binomial(n
, k
int64) *Int
{
198 // reduce the number of multiplications by reducing k
199 if n
/2 < k
&& k
<= n
{
200 k
= n
- k
// Binomial(n, k) == Binomial(n, n-k)
208 // Quo sets z to the quotient x/y for y != 0 and returns z.
209 // If y == 0, a division-by-zero run-time panic occurs.
210 // Quo implements truncated division (like Go); see QuoRem for more details.
211 func (z
*Int
) Quo(x
, y
*Int
) *Int
{
212 z
.abs
, _
= z
.abs
.div(nil, x
.abs
, y
.abs
)
213 z
.neg
= len(z
.abs
) > 0 && x
.neg
!= y
.neg
// 0 has no sign
217 // Rem sets z to the remainder x%y for y != 0 and returns z.
218 // If y == 0, a division-by-zero run-time panic occurs.
219 // Rem implements truncated modulus (like Go); see QuoRem for more details.
220 func (z
*Int
) Rem(x
, y
*Int
) *Int
{
221 _
, z
.abs
= nat(nil).div(z
.abs
, x
.abs
, y
.abs
)
222 z
.neg
= len(z
.abs
) > 0 && x
.neg
// 0 has no sign
226 // QuoRem sets z to the quotient x/y and r to the remainder x%y
227 // and returns the pair (z, r) for y != 0.
228 // If y == 0, a division-by-zero run-time panic occurs.
230 // QuoRem implements T-division and modulus (like Go):
232 // q = x/y with the result truncated to zero
235 // (See Daan Leijen, ``Division and Modulus for Computer Scientists''.)
236 // See DivMod for Euclidean division and modulus (unlike Go).
238 func (z
*Int
) QuoRem(x
, y
, r
*Int
) (*Int
, *Int
) {
239 z
.abs
, r
.abs
= z
.abs
.div(r
.abs
, x
.abs
, y
.abs
)
240 z
.neg
, r
.neg
= len(z
.abs
) > 0 && x
.neg
!= y
.neg
, len(r
.abs
) > 0 && x
.neg
// 0 has no sign
244 // Div sets z to the quotient x/y for y != 0 and returns z.
245 // If y == 0, a division-by-zero run-time panic occurs.
246 // Div implements Euclidean division (unlike Go); see DivMod for more details.
247 func (z
*Int
) Div(x
, y
*Int
) *Int
{
248 y_neg
:= y
.neg
// z may be an alias for y
261 // Mod sets z to the modulus x%y for y != 0 and returns z.
262 // If y == 0, a division-by-zero run-time panic occurs.
263 // Mod implements Euclidean modulus (unlike Go); see DivMod for more details.
264 func (z
*Int
) Mod(x
, y
*Int
) *Int
{
266 if z
== y ||
alias(z
.abs
, y
.abs
) {
281 // DivMod sets z to the quotient x div y and m to the modulus x mod y
282 // and returns the pair (z, m) for y != 0.
283 // If y == 0, a division-by-zero run-time panic occurs.
285 // DivMod implements Euclidean division and modulus (unlike Go):
287 // q = x div y such that
288 // m = x - y*q with 0 <= m < |y|
290 // (See Raymond T. Boute, ``The Euclidean definition of the functions
291 // div and mod''. ACM Transactions on Programming Languages and
292 // Systems (TOPLAS), 14(2):127-144, New York, NY, USA, 4/1992.
294 // See QuoRem for T-division and modulus (like Go).
296 func (z
*Int
) DivMod(x
, y
, m
*Int
) (*Int
, *Int
) {
298 if z
== y ||
alias(z
.abs
, y
.abs
) {
314 // Cmp compares x and y and returns:
320 func (x
*Int
) Cmp(y
*Int
) (r
int) {
321 // x cmp y == x cmp y
324 // (-x) cmp (-y) == -(x cmp y)
341 // CmpAbs compares the absolute values of x and y and returns:
347 func (x
*Int
) CmpAbs(y
*Int
) int {
348 return x
.abs
.cmp(y
.abs
)
351 // low32 returns the least significant 32 bits of x.
352 func low32(x nat
) uint32 {
359 // low64 returns the least significant 64 bits of x.
360 func low64(x nat
) uint64 {
365 if _W
== 32 && len(x
) > 1 {
366 return uint64(x
[1])<<32 | v
371 // Int64 returns the int64 representation of x.
372 // If x cannot be represented in an int64, the result is undefined.
373 func (x
*Int
) Int64() int64 {
374 v
:= int64(low64(x
.abs
))
381 // Uint64 returns the uint64 representation of x.
382 // If x cannot be represented in a uint64, the result is undefined.
383 func (x
*Int
) Uint64() uint64 {
387 // IsInt64 reports whether x can be represented as an int64.
388 func (x
*Int
) IsInt64() bool {
389 if len(x
.abs
) <= 64/_W
{
390 w
:= int64(low64(x
.abs
))
391 return w
>= 0 || x
.neg
&& w
== -w
396 // IsUint64 reports whether x can be represented as a uint64.
397 func (x
*Int
) IsUint64() bool {
398 return !x
.neg
&& len(x
.abs
) <= 64/_W
401 // SetString sets z to the value of s, interpreted in the given base,
402 // and returns z and a boolean indicating success. The entire string
403 // (not just a prefix) must be valid for success. If SetString fails,
404 // the value of z is undefined but the returned value is nil.
406 // The base argument must be 0 or a value between 2 and MaxBase.
407 // For base 0, the number prefix determines the actual base: A prefix of
408 // ``0b'' or ``0B'' selects base 2, ``0'', ``0o'' or ``0O'' selects base 8,
409 // and ``0x'' or ``0X'' selects base 16. Otherwise, the selected base is 10
410 // and no prefix is accepted.
412 // For bases <= 36, lower and upper case letters are considered the same:
413 // The letters 'a' to 'z' and 'A' to 'Z' represent digit values 10 to 35.
414 // For bases > 36, the upper case letters 'A' to 'Z' represent the digit
417 // For base 0, an underscore character ``_'' may appear between a base
418 // prefix and an adjacent digit, and between successive digits; such
419 // underscores do not change the value of the number.
420 // Incorrect placement of underscores is reported as an error if there
421 // are no other errors. If base != 0, underscores are not recognized
422 // and act like any other character that is not a valid digit.
424 func (z
*Int
) SetString(s
string, base
int) (*Int
, bool) {
425 return z
.setFromScanner(strings
.NewReader(s
), base
)
428 // setFromScanner implements SetString given an io.ByteScanner.
429 // For documentation see comments of SetString.
430 func (z
*Int
) setFromScanner(r io
.ByteScanner
, base
int) (*Int
, bool) {
431 if _
, _
, err
:= z
.scan(r
, base
); err
!= nil {
434 // entire content must have been consumed
435 if _
, err
:= r
.ReadByte(); err
!= io
.EOF
{
438 return z
, true // err == io.EOF => scan consumed all content of r
441 // SetBytes interprets buf as the bytes of a big-endian unsigned
442 // integer, sets z to that value, and returns z.
443 func (z
*Int
) SetBytes(buf
[]byte) *Int
{
444 z
.abs
= z
.abs
.setBytes(buf
)
449 // Bytes returns the absolute value of x as a big-endian byte slice.
451 // To use a fixed length slice, or a preallocated one, use FillBytes.
452 func (x
*Int
) Bytes() []byte {
453 buf
:= make([]byte, len(x
.abs
)*_S
)
454 return buf
[x
.abs
.bytes(buf
):]
457 // FillBytes sets buf to the absolute value of x, storing it as a zero-extended
458 // big-endian byte slice, and returns buf.
460 // If the absolute value of x doesn't fit in buf, FillBytes will panic.
461 func (x
*Int
) FillBytes(buf
[]byte) []byte {
462 // Clear whole buffer. (This gets optimized into a memclr.)
470 // BitLen returns the length of the absolute value of x in bits.
471 // The bit length of 0 is 0.
472 func (x
*Int
) BitLen() int {
473 return x
.abs
.bitLen()
476 // TrailingZeroBits returns the number of consecutive least significant zero
478 func (x
*Int
) TrailingZeroBits() uint {
479 return x
.abs
.trailingZeroBits()
482 // Exp sets z = x**y mod |m| (i.e. the sign of m is ignored), and returns z.
483 // If m == nil or m == 0, z = x**y unless y <= 0 then z = 1. If m != 0, y < 0,
484 // and x and m are not relatively prime, z is unchanged and nil is returned.
486 // Modular exponentiation of inputs of a particular size is not a
487 // cryptographically constant-time operation.
488 func (z
*Int
) Exp(x
, y
, m
*Int
) *Int
{
489 // See Knuth, volume 2, section 4.6.3.
492 if m
== nil ||
len(m
.abs
) == 0 {
495 // for y < 0: x**y mod m == (x**(-1))**|y| mod m
496 inverse
:= new(Int
).ModInverse(x
, m
)
506 mWords
= m
.abs
// m.abs may be nil for m == 0
509 z
.abs
= z
.abs
.expNN(xWords
, yWords
, mWords
)
510 z
.neg
= len(z
.abs
) > 0 && x
.neg
&& len(yWords
) > 0 && yWords
[0]&1 == 1 // 0 has no sign
511 if z
.neg
&& len(mWords
) > 0 {
512 // make modulus result positive
513 z
.abs
= z
.abs
.sub(mWords
, z
.abs
) // z == x**y mod |m| && 0 <= z < |m|
520 // GCD sets z to the greatest common divisor of a and b and returns z.
521 // If x or y are not nil, GCD sets their value such that z = a*x + b*y.
523 // a and b may be positive, zero or negative. (Before Go 1.14 both had
524 // to be > 0.) Regardless of the signs of a and b, z is always >= 0.
526 // If a == b == 0, GCD sets z = x = y = 0.
528 // If a == 0 and b != 0, GCD sets z = |b|, x = 0, y = sign(b) * 1.
530 // If a != 0 and b == 0, GCD sets z = |a|, x = sign(a) * 1, y = 0.
531 func (z
*Int
) GCD(x
, y
, a
, b
*Int
) *Int
{
532 if len(a
.abs
) == 0 ||
len(b
.abs
) == 0 {
533 lenA
, lenB
, negA
, negB
:= len(a
.abs
), len(b
.abs
), a
.neg
, b
.neg
559 return z
.lehmerGCD(x
, y
, a
, b
)
562 // lehmerSimulate attempts to simulate several Euclidean update steps
563 // using the leading digits of A and B. It returns u0, u1, v0, v1
564 // such that A and B can be updated as:
567 // Requirements: A >= B and len(B.abs) >= 2
568 // Since we are calculating with full words to avoid overflow,
569 // we use 'even' to track the sign of the cosequences.
570 // For even iterations: u0, v1 >= 0 && u1, v0 <= 0
571 // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0
572 func lehmerSimulate(A
, B
*Int
) (u0
, u1
, v0
, v1 Word
, even
bool) {
573 // initialize the digits
574 var a1
, a2
, u2
, v2 Word
576 m
:= len(B
.abs
) // m >= 2
577 n
:= len(A
.abs
) // n >= m >= 2
579 // extract the top Word of bits from A and B
581 a1
= A
.abs
[n
-1]<<h | A
.abs
[n
-2]>>(_W
-h
)
582 // B may have implicit zero words in the high bits if the lengths differ
585 a2
= B
.abs
[n
-1]<<h | B
.abs
[n
-2]>>(_W
-h
)
587 a2
= B
.abs
[n
-2] >> (_W
- h
)
592 // Since we are calculating with full words to avoid overflow,
593 // we use 'even' to track the sign of the cosequences.
594 // For even iterations: u0, v1 >= 0 && u1, v0 <= 0
595 // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0
596 // The first iteration starts with k=1 (odd).
598 // variables to track the cosequences
602 // Calculate the quotient and cosequences using Collins' stopping condition.
603 // Note that overflow of a Word is not possible when computing the remainder
604 // sequence and cosequences since the cosequence size is bounded by the input size.
605 // See section 4.2 of Jebelean for details.
606 for a2
>= v2
&& a1
-a2
>= v1
+v2
{
609 u0
, u1
, u2
= u1
, u2
, u1
+q
*u2
610 v0
, v1
, v2
= v1
, v2
, v1
+q
*v2
616 // lehmerUpdate updates the inputs A and B such that:
619 // where the signs of u0, u1, v0, v1 are given by even
620 // For even == true: u0, v1 >= 0 && u1, v0 <= 0
621 // For even == false: u0, v1 <= 0 && u1, v0 >= 0
622 // q, r, s, t are temporary variables to avoid allocations in the multiplication
623 func lehmerUpdate(A
, B
, q
, r
, s
, t
*Int
, u0
, u1
, v0
, v1 Word
, even
bool) {
625 t
.abs
= t
.abs
.setWord(u0
)
626 s
.abs
= s
.abs
.setWord(v0
)
633 r
.abs
= r
.abs
.setWord(u1
)
634 q
.abs
= q
.abs
.setWord(v1
)
645 // euclidUpdate performs a single step of the Euclidean GCD algorithm
646 // if extended is true, it also updates the cosequence Ua, Ub
647 func euclidUpdate(A
, B
, Ua
, Ub
, q
, r
, s
, t
*Int
, extended
bool) {
648 q
, r
= q
.QuoRem(A
, B
, r
)
650 *A
, *B
, *r
= *B
, *r
, *A
653 // Ua, Ub = Ub, Ua - q*Ub
661 // lehmerGCD sets z to the greatest common divisor of a and b,
662 // which both must be != 0, and returns z.
663 // If x or y are not nil, their values are set such that z = a*x + b*y.
664 // See Knuth, The Art of Computer Programming, Vol. 2, Section 4.5.2, Algorithm L.
665 // This implementation uses the improved condition by Collins requiring only one
666 // quotient and avoiding the possibility of single Word overflow.
667 // See Jebelean, "Improving the multiprecision Euclidean algorithm",
668 // Design and Implementation of Symbolic Computation Systems, pp 45-58.
669 // The cosequences are updated according to Algorithm 10.45 from
670 // Cohen et al. "Handbook of Elliptic and Hyperelliptic Curve Cryptography" pp 192.
671 func (z
*Int
) lehmerGCD(x
, y
, a
, b
*Int
) *Int
{
672 var A
, B
, Ua
, Ub
*Int
677 extended
:= x
!= nil || y
!= nil
680 // Ua (Ub) tracks how many times input a has been accumulated into A (B).
681 Ua
= new(Int
).SetInt64(1)
685 // temp variables for multiprecision update
692 if A
.abs
.cmp(B
.abs
) < 0 {
697 // loop invariant A >= B
699 // Attempt to calculate in single-precision using leading words of A and B.
700 u0
, u1
, v0
, v1
, even
:= lehmerSimulate(A
, B
)
702 // multiprecision Step
704 // Simulate the effect of the single-precision steps using the cosequences.
707 lehmerUpdate(A
, B
, q
, r
, s
, t
, u0
, u1
, v0
, v1
, even
)
710 // Ua = u0*Ua + v0*Ub
711 // Ub = u1*Ua + v1*Ub
712 lehmerUpdate(Ua
, Ub
, q
, r
, s
, t
, u0
, u1
, v0
, v1
, even
)
716 // Single-digit calculations failed to simulate any quotients.
717 // Do a standard Euclidean step.
718 euclidUpdate(A
, B
, Ua
, Ub
, q
, r
, s
, t
, extended
)
723 // extended Euclidean algorithm base case if B is a single Word
725 // A is longer than a single Word, so one update is needed.
726 euclidUpdate(A
, B
, Ua
, Ub
, q
, r
, s
, t
, extended
)
729 // A and B are both a single Word.
730 aWord
, bWord
:= A
.abs
[0], B
.abs
[0]
732 var ua
, ub
, va
, vb Word
737 q
, r
:= aWord
/bWord
, aWord%bWord
738 aWord
, bWord
= bWord
, r
744 t
.abs
= t
.abs
.setWord(ua
)
745 s
.abs
= s
.abs
.setWord(va
)
755 aWord
, bWord
= bWord
, aWord%bWord
763 // avoid aliasing b needed in the division below
770 y
.Mul(a
, Ua
) // y can safely alias a
790 // Rand sets z to a pseudo-random number in [0, n) and returns z.
792 // As this uses the math/rand package, it must not be used for
793 // security-sensitive work. Use crypto/rand.Int instead.
794 func (z
*Int
) Rand(rnd
*rand
.Rand
, n
*Int
) *Int
{
796 if n
.neg ||
len(n
.abs
) == 0 {
800 z
.abs
= z
.abs
.random(rnd
, n
.abs
, n
.abs
.bitLen())
804 // ModInverse sets z to the multiplicative inverse of g in the ring ℤ/nℤ
805 // and returns z. If g and n are not relatively prime, g has no multiplicative
806 // inverse in the ring ℤ/nℤ. In this case, z is unchanged and the return value
808 func (z
*Int
) ModInverse(g
, n
*Int
) *Int
{
809 // GCD expects parameters a and b to be > 0.
821 // if and only if d==1, g and n are relatively prime
822 if d
.Cmp(intOne
) != 0 {
826 // x and y are such that g*x + n*y = 1, therefore x is the inverse element,
827 // but it may be negative, so convert to the range 0 <= z < |n|
836 // Jacobi returns the Jacobi symbol (x/y), either +1, -1, or 0.
837 // The y argument must be an odd integer.
838 func Jacobi(x
, y
*Int
) int {
839 if len(y
.abs
) == 0 || y
.abs
[0]&1 == 0 {
840 panic(fmt
.Sprintf("big: invalid 2nd argument to Int.Jacobi: need odd integer but got %s", y
))
843 // We use the formulation described in chapter 2, section 2.4,
844 // "The Yacas Book of Algorithms":
845 // http://yacas.sourceforge.net/Algo.book.pdf
860 if b
.Cmp(intOne
) == 0 {
872 // handle factors of 2 in 'a'
873 s
:= a
.abs
.trailingZeroBits()
875 bmod8
:= b
.abs
[0] & 7
876 if bmod8
== 3 || bmod8
== 5 {
880 c
.Rsh(&a
, s
) // a = 2^s*c
882 // swap numerator and denominator
883 if b
.abs
[0]&3 == 3 && c
.abs
[0]&3 == 3 {
891 // modSqrt3Mod4 uses the identity
892 // (a^((p+1)/4))^2 mod p
895 // to calculate the square root of any quadratic residue mod p quickly for 3
897 func (z
*Int
) modSqrt3Mod4Prime(x
, p
*Int
) *Int
{
898 e
:= new(Int
).Add(p
, intOne
) // e = p + 1
899 e
.Rsh(e
, 2) // e = (p + 1) / 4
900 z
.Exp(x
, e
, p
) // z = x^e mod p
904 // modSqrt5Mod8 uses Atkin's observation that 2 is not a square mod p
905 // alpha == (2*a)^((p-5)/8) mod p
906 // beta == 2*a*alpha^2 mod p is a square root of -1
907 // b == a*alpha*(beta-1) mod p is a square root of a
908 // to calculate the square root of any quadratic residue mod p quickly for 5
910 func (z
*Int
) modSqrt5Mod8Prime(x
, p
*Int
) *Int
{
911 // p == 5 mod 8 implies p = e*8 + 5
912 // e is the quotient and 5 the remainder on division by 8
913 e
:= new(Int
).Rsh(p
, 3) // e = (p - 5) / 8
914 tx
:= new(Int
).Lsh(x
, 1) // tx = 2*x
915 alpha
:= new(Int
).Exp(tx
, e
, p
)
916 beta
:= new(Int
).Mul(alpha
, alpha
)
920 beta
.Sub(beta
, intOne
)
923 beta
.Mul(beta
, alpha
)
928 // modSqrtTonelliShanks uses the Tonelli-Shanks algorithm to find the square
929 // root of a quadratic residue modulo any prime.
930 func (z
*Int
) modSqrtTonelliShanks(x
, p
*Int
) *Int
{
931 // Break p-1 into s*2^e such that s is odd.
934 e
:= s
.abs
.trailingZeroBits()
937 // find some non-square n
940 for Jacobi(&n
, p
) != -1 {
944 // Core of the Tonelli-Shanks algorithm. Follows the description in
945 // section 6 of "Square roots from 1; 24, 51, 10 to Dan Shanks" by Ezra
947 // https://www.maa.org/sites/default/files/pdf/upload_library/22/Polya/07468342.di020786.02p0470a.pdf
951 y
.Exp(x
, &y
, p
) // y = x^((s+1)/2)
952 b
.Exp(x
, &s
, p
) // b = x^s
953 g
.Exp(&n
, &s
, p
) // g = n^s
956 // find the least m such that ord_p(b) = 2^m
959 for t
.Cmp(intOne
) != 0 {
960 t
.Mul(&t
, &t
).Mod(&t
, p
)
968 t
.SetInt64(0).SetBit(&t
, int(r
-m
-1), 1).Exp(&g
, &t
, p
)
969 // t = g^(2^(r-m-1)) mod p
970 g
.Mul(&t
, &t
).Mod(&g
, p
) // g = g^(2^(r-m)) mod p
971 y
.Mul(&y
, &t
).Mod(&y
, p
)
972 b
.Mul(&b
, &g
).Mod(&b
, p
)
977 // ModSqrt sets z to a square root of x mod p if such a square root exists, and
978 // returns z. The modulus p must be an odd prime. If x is not a square mod p,
979 // ModSqrt leaves z unchanged and returns nil. This function panics if p is
980 // not an odd integer.
981 func (z
*Int
) ModSqrt(x
, p
*Int
) *Int
{
982 switch Jacobi(x
, p
) {
984 return nil // x is not a square mod p
986 return z
.SetInt64(0) // sqrt(0) mod p = 0
990 if x
.neg || x
.Cmp(p
) >= 0 { // ensure 0 <= x < p
991 x
= new(Int
).Mod(x
, p
)
995 case p
.abs
[0]%4
== 3:
996 // Check whether p is 3 mod 4, and if so, use the faster algorithm.
997 return z
.modSqrt3Mod4Prime(x
, p
)
998 case p
.abs
[0]%8
== 5:
999 // Check whether p is 5 mod 8, use Atkin's algorithm.
1000 return z
.modSqrt5Mod8Prime(x
, p
)
1002 // Otherwise, use Tonelli-Shanks.
1003 return z
.modSqrtTonelliShanks(x
, p
)
1007 // Lsh sets z = x << n and returns z.
1008 func (z
*Int
) Lsh(x
*Int
, n
uint) *Int
{
1009 z
.abs
= z
.abs
.shl(x
.abs
, n
)
1014 // Rsh sets z = x >> n and returns z.
1015 func (z
*Int
) Rsh(x
*Int
, n
uint) *Int
{
1017 // (-x) >> s == ^(x-1) >> s == ^((x-1) >> s) == -(((x-1) >> s) + 1)
1018 t
:= z
.abs
.sub(x
.abs
, natOne
) // no underflow because |x| > 0
1020 z
.abs
= t
.add(t
, natOne
)
1021 z
.neg
= true // z cannot be zero if x is negative
1025 z
.abs
= z
.abs
.shr(x
.abs
, n
)
1030 // Bit returns the value of the i'th bit of x. That is, it
1031 // returns (x>>i)&1. The bit index i must be >= 0.
1032 func (x
*Int
) Bit(i
int) uint {
1034 // optimization for common case: odd/even test of x
1036 return uint(x
.abs
[0] & 1) // bit 0 is same for -x
1041 panic("negative bit index")
1044 t
:= nat(nil).sub(x
.abs
, natOne
)
1045 return t
.bit(uint(i
)) ^ 1
1048 return x
.abs
.bit(uint(i
))
1051 // SetBit sets z to x, with x's i'th bit set to b (0 or 1).
1052 // That is, if b is 1 SetBit sets z = x | (1 << i);
1053 // if b is 0 SetBit sets z = x &^ (1 << i). If b is not 0 or 1,
1054 // SetBit will panic.
1055 func (z
*Int
) SetBit(x
*Int
, i
int, b
uint) *Int
{
1057 panic("negative bit index")
1060 t
:= z
.abs
.sub(x
.abs
, natOne
)
1061 t
= t
.setBit(t
, uint(i
), b
^1)
1062 z
.abs
= t
.add(t
, natOne
)
1063 z
.neg
= len(z
.abs
) > 0
1066 z
.abs
= z
.abs
.setBit(x
.abs
, uint(i
), b
)
1071 // And sets z = x & y and returns z.
1072 func (z
*Int
) And(x
, y
*Int
) *Int
{
1075 // (-x) & (-y) == ^(x-1) & ^(y-1) == ^((x-1) | (y-1)) == -(((x-1) | (y-1)) + 1)
1076 x1
:= nat(nil).sub(x
.abs
, natOne
)
1077 y1
:= nat(nil).sub(y
.abs
, natOne
)
1078 z
.abs
= z
.abs
.add(z
.abs
.or(x1
, y1
), natOne
)
1079 z
.neg
= true // z cannot be zero if x and y are negative
1084 z
.abs
= z
.abs
.and(x
.abs
, y
.abs
)
1091 x
, y
= y
, x
// & is symmetric
1094 // x & (-y) == x & ^(y-1) == x &^ (y-1)
1095 y1
:= nat(nil).sub(y
.abs
, natOne
)
1096 z
.abs
= z
.abs
.andNot(x
.abs
, y1
)
1101 // AndNot sets z = x &^ y and returns z.
1102 func (z
*Int
) AndNot(x
, y
*Int
) *Int
{
1105 // (-x) &^ (-y) == ^(x-1) &^ ^(y-1) == ^(x-1) & (y-1) == (y-1) &^ (x-1)
1106 x1
:= nat(nil).sub(x
.abs
, natOne
)
1107 y1
:= nat(nil).sub(y
.abs
, natOne
)
1108 z
.abs
= z
.abs
.andNot(y1
, x1
)
1114 z
.abs
= z
.abs
.andNot(x
.abs
, y
.abs
)
1120 // (-x) &^ y == ^(x-1) &^ y == ^(x-1) & ^y == ^((x-1) | y) == -(((x-1) | y) + 1)
1121 x1
:= nat(nil).sub(x
.abs
, natOne
)
1122 z
.abs
= z
.abs
.add(z
.abs
.or(x1
, y
.abs
), natOne
)
1123 z
.neg
= true // z cannot be zero if x is negative and y is positive
1127 // x &^ (-y) == x &^ ^(y-1) == x & (y-1)
1128 y1
:= nat(nil).sub(y
.abs
, natOne
)
1129 z
.abs
= z
.abs
.and(x
.abs
, y1
)
1134 // Or sets z = x | y and returns z.
1135 func (z
*Int
) Or(x
, y
*Int
) *Int
{
1138 // (-x) | (-y) == ^(x-1) | ^(y-1) == ^((x-1) & (y-1)) == -(((x-1) & (y-1)) + 1)
1139 x1
:= nat(nil).sub(x
.abs
, natOne
)
1140 y1
:= nat(nil).sub(y
.abs
, natOne
)
1141 z
.abs
= z
.abs
.add(z
.abs
.and(x1
, y1
), natOne
)
1142 z
.neg
= true // z cannot be zero if x and y are negative
1147 z
.abs
= z
.abs
.or(x
.abs
, y
.abs
)
1154 x
, y
= y
, x
// | is symmetric
1157 // x | (-y) == x | ^(y-1) == ^((y-1) &^ x) == -(^((y-1) &^ x) + 1)
1158 y1
:= nat(nil).sub(y
.abs
, natOne
)
1159 z
.abs
= z
.abs
.add(z
.abs
.andNot(y1
, x
.abs
), natOne
)
1160 z
.neg
= true // z cannot be zero if one of x or y is negative
1164 // Xor sets z = x ^ y and returns z.
1165 func (z
*Int
) Xor(x
, y
*Int
) *Int
{
1168 // (-x) ^ (-y) == ^(x-1) ^ ^(y-1) == (x-1) ^ (y-1)
1169 x1
:= nat(nil).sub(x
.abs
, natOne
)
1170 y1
:= nat(nil).sub(y
.abs
, natOne
)
1171 z
.abs
= z
.abs
.xor(x1
, y1
)
1177 z
.abs
= z
.abs
.xor(x
.abs
, y
.abs
)
1184 x
, y
= y
, x
// ^ is symmetric
1187 // x ^ (-y) == x ^ ^(y-1) == ^(x ^ (y-1)) == -((x ^ (y-1)) + 1)
1188 y1
:= nat(nil).sub(y
.abs
, natOne
)
1189 z
.abs
= z
.abs
.add(z
.abs
.xor(x
.abs
, y1
), natOne
)
1190 z
.neg
= true // z cannot be zero if only one of x or y is negative
1194 // Not sets z = ^x and returns z.
1195 func (z
*Int
) Not(x
*Int
) *Int
{
1197 // ^(-x) == ^(^(x-1)) == x-1
1198 z
.abs
= z
.abs
.sub(x
.abs
, natOne
)
1203 // ^x == -x-1 == -(x+1)
1204 z
.abs
= z
.abs
.add(x
.abs
, natOne
)
1205 z
.neg
= true // z cannot be zero if x is positive
1209 // Sqrt sets z to ⌊√x⌋, the largest integer such that z² ≤ x, and returns z.
1210 // It panics if x is negative.
1211 func (z
*Int
) Sqrt(x
*Int
) *Int
{
1213 panic("square root of negative number")
1216 z
.abs
= z
.abs
.sqrt(x
.abs
)