From 2473e19f735f88a69609e7b5c96bbd0434e3a14a Mon Sep 17 00:00:00 2001 From: Jbop1626 <34898270+jbop1626@users.noreply.github.com> Date: Tue, 5 Nov 2019 14:22:15 -0500 Subject: [PATCH] Rewrite everything in C99 - Replace BigInteger with GNU MP - Update README --- README.md | 12 +- src/bigint/include/BigInteger.hpp | 215 - src/bigint/include/BigIntegerAlgorithms.hpp | 25 - src/bigint/include/BigIntegerLibrary.hpp | 8 - src/bigint/include/BigIntegerUtils.hpp | 72 - src/bigint/include/BigUnsigned.hpp | 418 -- src/bigint/include/BigUnsignedInABase.hpp | 122 - src/bigint/include/NumberlikeArray.hpp | 177 - src/bigint/src/BigInteger.cpp | 405 -- src/bigint/src/BigIntegerAlgorithms.cpp | 71 - src/bigint/src/BigIntegerUtils.cpp | 50 - src/bigint/src/BigUnsigned.cpp | 697 --- src/bigint/src/BigUnsignedInABase.cpp | 125 - src/ecc/ecc.c | 389 ++ src/ecc/ecc.cpp | 367 -- src/ecc/ecc.h | 103 + src/ecc/ecc.hpp | 95 - src/example.c | 113 + src/mini-gmp/mini-gmp.c | 4412 +++++++++++++++++++ src/mini-gmp/mini-gmp.h | 298 ++ src/ninty-233.c | 248 ++ src/ninty-233.cpp | 204 - src/ninty-233.h | 59 + src/ninty-233.hpp | 56 - src/sha1/sha1.h | 12 +- src/sha1/sha1.hpp | 32 - 26 files changed, 5642 insertions(+), 3143 deletions(-) delete mode 100644 src/bigint/include/BigInteger.hpp delete mode 100644 src/bigint/include/BigIntegerAlgorithms.hpp delete mode 100644 src/bigint/include/BigIntegerLibrary.hpp delete mode 100644 src/bigint/include/BigIntegerUtils.hpp delete mode 100644 src/bigint/include/BigUnsigned.hpp delete mode 100644 src/bigint/include/BigUnsignedInABase.hpp delete mode 100644 src/bigint/include/NumberlikeArray.hpp delete mode 100644 src/bigint/src/BigInteger.cpp delete mode 100644 src/bigint/src/BigIntegerAlgorithms.cpp delete mode 100644 src/bigint/src/BigIntegerUtils.cpp delete mode 100644 src/bigint/src/BigUnsigned.cpp delete mode 100644 src/bigint/src/BigUnsignedInABase.cpp create mode 100644 src/ecc/ecc.c delete mode 100644 src/ecc/ecc.cpp create mode 100644 src/ecc/ecc.h delete mode 100644 src/ecc/ecc.hpp create mode 100644 src/example.c create mode 100644 src/mini-gmp/mini-gmp.c create mode 100644 src/mini-gmp/mini-gmp.h create mode 100644 src/ninty-233.c delete mode 100644 src/ninty-233.cpp create mode 100644 src/ninty-233.h delete mode 100644 src/ninty-233.hpp delete mode 100644 src/sha1/sha1.hpp diff --git a/README.md b/README.md index 3eb667f..18827b3 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,15 @@ ### ninty-233 -ninty-233 is a library for [ECC](https://en.wikipedia.org/wiki/Elliptic-curve_cryptography) operations using sect233r1 / NIST B-233, the curve used in the [iQue Player](https://en.wikipedia.org/wiki/IQue_Player) and [Nintendo Wii](https://en.wikipedia.org/wiki/Wii). +ninty-233 is a C99 library for [ECC](https://en.wikipedia.org/wiki/Elliptic-curve_cryptography) operations using sect233r1 / NIST B-233, the curve used in the [iQue Player](https://en.wikipedia.org/wiki/IQue_Player) and [Nintendo Wii](https://en.wikipedia.org/wiki/Wii). It can be used for [ECDH](https://en.wikipedia.org/wiki/Elliptic-curve_Diffie%E2%80%93Hellman) (used to create some encryption keys on the iQue Player) and [ECDSA](https://en.wikipedia.org/wiki/Elliptic_Curve_Digital_Signature_Algorithm) signing/verification (used to sign game saves on both consoles and to sign recrypt.sys on the iQue Player). -Arbitrary-precision arithmetic is done using the public domain [C++ Big Integer Library](https://mattmccutchen.net/bigint/) by Matt McCutchen. +ninty-233 should **NOT** be expected to be secure; it is intended to be used as a tool for working with keys and/or data that are **already known** (made obvious by the fact that there is no function provided for generating keys). In its current state, I expect it to be vulnerable to various attacks (e.g. the GF(2^m) element and elliptic curve point arithmetic operations are not resistant to timing analysis). The ``generate_k()`` function in particular is not cryptographically secure, and is simply a trivial implementation that will allow, for example, signing homebrew apps. -SHA1 implementation is from the public domain [WjCryptLib](https://github.com/WaterJuice/WjCryptLib) by [WaterJuice](https://github.com/WaterJuice). +ninty-233 currently requires an architecture with unsigned integer types of exactly 8 and 32 bits - that is, both ``uint8_t`` and ``uint32_t`` must be defined. This likely isn't a problem for literally anyone, but it's probably good to know before compiling for a calculator or something. + +Arbitrary-precision arithmetic is done using the [GNU MP library](https://gmplib.org/) in the form of mini-gmp, a standalone subset of GMP. + +The SHA1 implementation is a slightly modified version of the public domain [WjCryptLib](https://github.com/WaterJuice/WjCryptLib) implementation from [WaterJuice](https://github.com/WaterJuice) et al. + +ninty-233 is licensed under the GPLv3 or (at your option) any later version (and I've opted to use GMP under the same license). diff --git a/src/bigint/include/BigInteger.hpp b/src/bigint/include/BigInteger.hpp deleted file mode 100644 index a5526e5..0000000 --- a/src/bigint/include/BigInteger.hpp +++ /dev/null @@ -1,215 +0,0 @@ -#ifndef BIGINTEGER_H -#define BIGINTEGER_H - -#include "BigUnsigned.hpp" - -/* A BigInteger object represents a signed integer of size limited only by - * available memory. BigUnsigneds support most mathematical operators and can - * be converted to and from most primitive integer types. - * - * A BigInteger is just an aggregate of a BigUnsigned and a sign. (It is no - * longer derived from BigUnsigned because that led to harmful implicit - * conversions.) */ -class BigInteger { - -public: - typedef BigUnsigned::Blk Blk; - typedef BigUnsigned::Index Index; - typedef BigUnsigned::CmpRes CmpRes; - static const CmpRes - less = BigUnsigned::less , - equal = BigUnsigned::equal , - greater = BigUnsigned::greater; - // Enumeration for the sign of a BigInteger. - enum Sign { negative = -1, zero = 0, positive = 1 }; - -protected: - Sign sign; - BigUnsigned mag; - -public: - // Constructs zero. - BigInteger() : sign(zero), mag() {} - - // Copy constructor - BigInteger(const BigInteger &x) : sign(x.sign), mag(x.mag) {}; - - // Assignment operator - void operator=(const BigInteger &x); - - // Constructor that copies from a given array of blocks with a sign. - BigInteger(const Blk *b, Index blen, Sign s); - - // Nonnegative constructor that copies from a given array of blocks. - BigInteger(const Blk *b, Index blen) : mag(b, blen) { - sign = mag.isZero() ? zero : positive; - } - - // Constructor from a BigUnsigned and a sign - BigInteger(const BigUnsigned &x, Sign s); - - // Nonnegative constructor from a BigUnsigned - BigInteger(const BigUnsigned &x) : mag(x) { - sign = mag.isZero() ? zero : positive; - } - - // Constructors from primitive integer types - BigInteger(unsigned long x); - BigInteger( long x); - BigInteger(unsigned int x); - BigInteger( int x); - BigInteger(unsigned short x); - BigInteger( short x); - - /* Converters to primitive integer types - * The implicit conversion operators caused trouble, so these are now - * named. */ - unsigned long toUnsignedLong () const; - long toLong () const; - unsigned int toUnsignedInt () const; - int toInt () const; - unsigned short toUnsignedShort() const; - short toShort () const; -protected: - // Helper - template X convertToUnsignedPrimitive() const; - template X convertToSignedPrimitive() const; -public: - - // ACCESSORS - Sign getSign() const { return sign; } - /* The client can't do any harm by holding a read-only reference to the - * magnitude. */ - const BigUnsigned &getMagnitude() const { return mag; } - - // Some accessors that go through to the magnitude - Index getLength() const { return mag.getLength(); } - Index getCapacity() const { return mag.getCapacity(); } - Blk getBlock(Index i) const { return mag.getBlock(i); } - bool isZero() const { return sign == zero; } // A bit special - - // COMPARISONS - - // Compares this to x like Perl's <=> - CmpRes compareTo(const BigInteger &x) const; - - // Ordinary comparison operators - bool operator ==(const BigInteger &x) const { - return sign == x.sign && mag == x.mag; - } - bool operator !=(const BigInteger &x) const { return !operator ==(x); }; - bool operator < (const BigInteger &x) const { return compareTo(x) == less ; } - bool operator <=(const BigInteger &x) const { return compareTo(x) != greater; } - bool operator >=(const BigInteger &x) const { return compareTo(x) != less ; } - bool operator > (const BigInteger &x) const { return compareTo(x) == greater; } - - // OPERATORS -- See the discussion in BigUnsigned.hh. - void add (const BigInteger &a, const BigInteger &b); - void subtract(const BigInteger &a, const BigInteger &b); - void multiply(const BigInteger &a, const BigInteger &b); - /* See the comment on BigUnsigned::divideWithRemainder. Semantics - * differ from those of primitive integers when negatives and/or zeros - * are involved. */ - void divideWithRemainder(const BigInteger &b, BigInteger &q); - void negate(const BigInteger &a); - - /* Bitwise operators are not provided for BigIntegers. Use - * getMagnitude to get the magnitude and operate on that instead. */ - - BigInteger operator +(const BigInteger &x) const; - BigInteger operator -(const BigInteger &x) const; - BigInteger operator *(const BigInteger &x) const; - BigInteger operator /(const BigInteger &x) const; - BigInteger operator %(const BigInteger &x) const; - BigInteger operator -() const; - - void operator +=(const BigInteger &x); - void operator -=(const BigInteger &x); - void operator *=(const BigInteger &x); - void operator /=(const BigInteger &x); - void operator %=(const BigInteger &x); - void flipSign(); - - // INCREMENT/DECREMENT OPERATORS - void operator ++( ); - void operator ++(int); - void operator --( ); - void operator --(int); -}; - -// NORMAL OPERATORS -/* These create an object to hold the result and invoke - * the appropriate put-here operation on it, passing - * this and x. The new object is then returned. */ -inline BigInteger BigInteger::operator +(const BigInteger &x) const { - BigInteger ans; - ans.add(*this, x); - return ans; -} -inline BigInteger BigInteger::operator -(const BigInteger &x) const { - BigInteger ans; - ans.subtract(*this, x); - return ans; -} -inline BigInteger BigInteger::operator *(const BigInteger &x) const { - BigInteger ans; - ans.multiply(*this, x); - return ans; -} -inline BigInteger BigInteger::operator /(const BigInteger &x) const { - if (x.isZero()) throw "BigInteger::operator /: division by zero"; - BigInteger q, r; - r = *this; - r.divideWithRemainder(x, q); - return q; -} -inline BigInteger BigInteger::operator %(const BigInteger &x) const { - if (x.isZero()) throw "BigInteger::operator %: division by zero"; - BigInteger q, r; - r = *this; - r.divideWithRemainder(x, q); - return r; -} -inline BigInteger BigInteger::operator -() const { - BigInteger ans; - ans.negate(*this); - return ans; -} - -/* - * ASSIGNMENT OPERATORS - * - * Now the responsibility for making a temporary copy if necessary - * belongs to the put-here operations. See Assignment Operators in - * BigUnsigned.hh. - */ -inline void BigInteger::operator +=(const BigInteger &x) { - add(*this, x); -} -inline void BigInteger::operator -=(const BigInteger &x) { - subtract(*this, x); -} -inline void BigInteger::operator *=(const BigInteger &x) { - multiply(*this, x); -} -inline void BigInteger::operator /=(const BigInteger &x) { - if (x.isZero()) throw "BigInteger::operator /=: division by zero"; - /* The following technique is slightly faster than copying *this first - * when x is large. */ - BigInteger q; - divideWithRemainder(x, q); - // *this contains the remainder, but we overwrite it with the quotient. - *this = q; -} -inline void BigInteger::operator %=(const BigInteger &x) { - if (x.isZero()) throw "BigInteger::operator %=: division by zero"; - BigInteger q; - // Mods *this by x. Don't care about quotient left in q. - divideWithRemainder(x, q); -} -// This one is trivial -inline void BigInteger::flipSign() { - sign = Sign(-sign); -} - -#endif diff --git a/src/bigint/include/BigIntegerAlgorithms.hpp b/src/bigint/include/BigIntegerAlgorithms.hpp deleted file mode 100644 index 0b31c41..0000000 --- a/src/bigint/include/BigIntegerAlgorithms.hpp +++ /dev/null @@ -1,25 +0,0 @@ -#ifndef BIGINTEGERALGORITHMS_H -#define BIGINTEGERALGORITHMS_H - -#include "BigInteger.hpp" - -/* Some mathematical algorithms for big integers. - * This code is new and, as such, experimental. */ - -// Returns the greatest common divisor of a and b. -BigUnsigned gcd(BigUnsigned a, BigUnsigned b); - -/* Extended Euclidean algorithm. - * Given m and n, finds gcd g and numbers r, s such that r*m + s*n == g. */ -void extendedEuclidean(BigInteger m, BigInteger n, - BigInteger &g, BigInteger &r, BigInteger &s); - -/* Returns the multiplicative inverse of x modulo n, or throws an exception if - * they have a common factor. */ -BigUnsigned modinv(const BigInteger &x, const BigUnsigned &n); - -// Returns (base ^ exponent) % modulus. -BigUnsigned modexp(const BigInteger &base, const BigUnsigned &exponent, - const BigUnsigned &modulus); - -#endif diff --git a/src/bigint/include/BigIntegerLibrary.hpp b/src/bigint/include/BigIntegerLibrary.hpp deleted file mode 100644 index 07b61ad..0000000 --- a/src/bigint/include/BigIntegerLibrary.hpp +++ /dev/null @@ -1,8 +0,0 @@ -// This header file includes all of the library header files. - -#include "NumberlikeArray.hpp" -#include "BigUnsigned.hpp" -#include "BigInteger.hpp" -#include "BigIntegerAlgorithms.hpp" -#include "BigUnsignedInABase.hpp" -#include "BigIntegerUtils.hpp" diff --git a/src/bigint/include/BigIntegerUtils.hpp b/src/bigint/include/BigIntegerUtils.hpp deleted file mode 100644 index 87cd3d0..0000000 --- a/src/bigint/include/BigIntegerUtils.hpp +++ /dev/null @@ -1,72 +0,0 @@ -#ifndef BIGINTEGERUTILS_H -#define BIGINTEGERUTILS_H - -#include "BigInteger.hpp" -#include -#include - -/* This file provides: - * - Convenient std::string <-> BigUnsigned/BigInteger conversion routines - * - std::ostream << operators for BigUnsigned/BigInteger */ - -// std::string conversion routines. Base 10 only. -std::string bigUnsignedToString(const BigUnsigned &x); -std::string bigIntegerToString(const BigInteger &x); -BigUnsigned stringToBigUnsigned(const std::string &s); -BigInteger stringToBigInteger(const std::string &s); - -// Creates a BigInteger from data such as `char's; read below for details. -template -BigInteger dataToBigInteger(const T* data, BigInteger::Index length, BigInteger::Sign sign); - -// Outputs x to os, obeying the flags `dec', `hex', `bin', and `showbase'. -std::ostream &operator <<(std::ostream &os, const BigUnsigned &x); - -// Outputs x to os, obeying the flags `dec', `hex', `bin', and `showbase'. -// My somewhat arbitrary policy: a negative sign comes before a base indicator (like -0xFF). -std::ostream &operator <<(std::ostream &os, const BigInteger &x); - -// BEGIN TEMPLATE DEFINITIONS. - -/* - * Converts binary data to a BigInteger. - * Pass an array `data', its length, and the desired sign. - * - * Elements of `data' may be of any type `T' that has the following - * two properties (this includes almost all integral types): - * - * (1) `sizeof(T)' correctly gives the amount of binary data in one - * value of `T' and is a factor of `sizeof(Blk)'. - * - * (2) When a value of `T' is casted to a `Blk', the low bytes of - * the result contain the desired binary data. - */ -template -BigInteger dataToBigInteger(const T* data, BigInteger::Index length, BigInteger::Sign sign) { - // really ceiling(numBytes / sizeof(BigInteger::Blk)) - unsigned int pieceSizeInBits = 8 * sizeof(T); - unsigned int piecesPerBlock = sizeof(BigInteger::Blk) / sizeof(T); - unsigned int numBlocks = (length + piecesPerBlock - 1) / piecesPerBlock; - - // Allocate our block array - BigInteger::Blk *blocks = new BigInteger::Blk[numBlocks]; - - BigInteger::Index blockNum, pieceNum, pieceNumHere; - - // Convert - for (blockNum = 0, pieceNum = 0; blockNum < numBlocks; blockNum++) { - BigInteger::Blk curBlock = 0; - for (pieceNumHere = 0; pieceNumHere < piecesPerBlock && pieceNum < length; - pieceNumHere++, pieceNum++) - curBlock |= (BigInteger::Blk(data[pieceNum]) << (pieceSizeInBits * pieceNumHere)); - blocks[blockNum] = curBlock; - } - - // Create the BigInteger. - BigInteger x(blocks, numBlocks, sign); - - delete [] blocks; - return x; -} - -#endif diff --git a/src/bigint/include/BigUnsigned.hpp b/src/bigint/include/BigUnsigned.hpp deleted file mode 100644 index c89bf67..0000000 --- a/src/bigint/include/BigUnsigned.hpp +++ /dev/null @@ -1,418 +0,0 @@ -#ifndef BIGUNSIGNED_H -#define BIGUNSIGNED_H - -#include "NumberlikeArray.hpp" - -/* A BigUnsigned object represents a nonnegative integer of size limited only by - * available memory. BigUnsigneds support most mathematical operators and can - * be converted to and from most primitive integer types. - * - * The number is stored as a NumberlikeArray of unsigned longs as if it were - * written in base 256^sizeof(unsigned long). The least significant block is - * first, and the length is such that the most significant block is nonzero. */ -class BigUnsigned : protected NumberlikeArray { - -public: - // Enumeration for the result of a comparison. - enum CmpRes { less = -1, equal = 0, greater = 1 }; - - // BigUnsigneds are built with a Blk type of unsigned long. - typedef unsigned long Blk; - - typedef NumberlikeArray::Index Index; - NumberlikeArray::N; - -protected: - // Creates a BigUnsigned with a capacity; for internal use. - BigUnsigned(int, Index c) : NumberlikeArray(0, c) {} - - // Decreases len to eliminate any leading zero blocks. - void zapLeadingZeros() { - while (len > 0 && blk[len - 1] == 0) - len--; - } - -public: - // Constructs zero. - BigUnsigned() : NumberlikeArray() {} - - // Copy constructor - BigUnsigned(const BigUnsigned &x) : NumberlikeArray(x) {} - - // Assignment operator - void operator=(const BigUnsigned &x) { - NumberlikeArray::operator =(x); - } - - // Constructor that copies from a given array of blocks. - BigUnsigned(const Blk *b, Index blen) : NumberlikeArray(b, blen) { - // Eliminate any leading zeros we may have been passed. - zapLeadingZeros(); - } - - // Destructor. NumberlikeArray does the delete for us. - ~BigUnsigned() {} - - // Constructors from primitive integer types - BigUnsigned(unsigned long x); - BigUnsigned( long x); - BigUnsigned(unsigned int x); - BigUnsigned( int x); - BigUnsigned(unsigned short x); - BigUnsigned( short x); -protected: - // Helpers - template void initFromPrimitive (X x); - template void initFromSignedPrimitive(X x); -public: - - /* Converters to primitive integer types - * The implicit conversion operators caused trouble, so these are now - * named. */ - unsigned long toUnsignedLong () const; - long toLong () const; - unsigned int toUnsignedInt () const; - int toInt () const; - unsigned short toUnsignedShort() const; - short toShort () const; -protected: - // Helpers - template X convertToSignedPrimitive() const; - template X convertToPrimitive () const; -public: - - // BIT/BLOCK ACCESSORS - - // Expose these from NumberlikeArray directly. - NumberlikeArray::getCapacity; - NumberlikeArray::getLength; - - /* Returns the requested block, or 0 if it is beyond the length (as if - * the number had 0s infinitely to the left). */ - Blk getBlock(Index i) const { return i >= len ? 0 : blk[i]; } - /* Sets the requested block. The number grows or shrinks as necessary. */ - void setBlock(Index i, Blk newBlock); - - // The number is zero if and only if the canonical length is zero. - bool isZero() const { return NumberlikeArray::isEmpty(); } - - /* Returns the length of the number in bits, i.e., zero if the number - * is zero and otherwise one more than the largest value of bi for - * which getBit(bi) returns true. */ - Index bitLength() const; - /* Get the state of bit bi, which has value 2^bi. Bits beyond the - * number's length are considered to be 0. */ - bool getBit(Index bi) const { - return (getBlock(bi / N) & (Blk(1) << (bi % N))) != 0; - } - /* Sets the state of bit bi to newBit. The number grows or shrinks as - * necessary. */ - void setBit(Index bi, bool newBit); - - // COMPARISONS - - // Compares this to x like Perl's <=> - CmpRes compareTo(const BigUnsigned &x) const; - - // Ordinary comparison operators - bool operator ==(const BigUnsigned &x) const { - return NumberlikeArray::operator ==(x); - } - bool operator !=(const BigUnsigned &x) const { - return NumberlikeArray::operator !=(x); - } - bool operator < (const BigUnsigned &x) const { return compareTo(x) == less ; } - bool operator <=(const BigUnsigned &x) const { return compareTo(x) != greater; } - bool operator >=(const BigUnsigned &x) const { return compareTo(x) != less ; } - bool operator > (const BigUnsigned &x) const { return compareTo(x) == greater; } - - /* - * BigUnsigned and BigInteger both provide three kinds of operators. - * Here ``big-integer'' refers to BigInteger or BigUnsigned. - * - * (1) Overloaded ``return-by-value'' operators: - * +, -, *, /, %, unary -, &, |, ^, <<, >>. - * Big-integer code using these operators looks identical to code using - * the primitive integer types. These operators take one or two - * big-integer inputs and return a big-integer result, which can then - * be assigned to a BigInteger variable or used in an expression. - * Example: - * BigInteger a(1), b = 1; - * BigInteger c = a + b; - * - * (2) Overloaded assignment operators: - * +=, -=, *=, /=, %=, flipSign, &=, |=, ^=, <<=, >>=, ++, --. - * Again, these are used on big integers just like on ints. They take - * one writable big integer that both provides an operand and receives a - * result. Most also take a second read-only operand. - * Example: - * BigInteger a(1), b(1); - * a += b; - * - * (3) Copy-less operations: `add', `subtract', etc. - * These named methods take operands as arguments and store the result - * in the receiver (*this), avoiding unnecessary copies and allocations. - * `divideWithRemainder' is special: it both takes the dividend from and - * stores the remainder into the receiver, and it takes a separate - * object in which to store the quotient. NOTE: If you are wondering - * why these don't return a value, you probably mean to use the - * overloaded return-by-value operators instead. - * - * Examples: - * BigInteger a(43), b(7), c, d; - * - * c = a + b; // Now c == 50. - * c.add(a, b); // Same effect but without the two copies. - * - * c.divideWithRemainder(b, d); - * // 50 / 7; now d == 7 (quotient) and c == 1 (remainder). - * - * // ``Aliased'' calls now do the right thing using a temporary - * // copy, but see note on `divideWithRemainder'. - * a.add(a, b); - */ - - // COPY-LESS OPERATIONS - - // These 8: Arguments are read-only operands, result is saved in *this. - void add(const BigUnsigned &a, const BigUnsigned &b); - void subtract(const BigUnsigned &a, const BigUnsigned &b); - void multiply(const BigUnsigned &a, const BigUnsigned &b); - void bitAnd(const BigUnsigned &a, const BigUnsigned &b); - void bitOr(const BigUnsigned &a, const BigUnsigned &b); - void bitXor(const BigUnsigned &a, const BigUnsigned &b); - /* Negative shift amounts translate to opposite-direction shifts, - * except for -2^(8*sizeof(int)-1) which is unimplemented. */ - void bitShiftLeft(const BigUnsigned &a, int b); - void bitShiftRight(const BigUnsigned &a, int b); - - /* `a.divideWithRemainder(b, q)' is like `q = a / b, a %= b'. - * / and % use semantics similar to Knuth's, which differ from the - * primitive integer semantics under division by zero. See the - * implementation in BigUnsigned.cc for details. - * `a.divideWithRemainder(b, a)' throws an exception: it doesn't make - * sense to write quotient and remainder into the same variable. */ - void divideWithRemainder(const BigUnsigned &b, BigUnsigned &q); - - /* `divide' and `modulo' are no longer offered. Use - * `divideWithRemainder' instead. */ - - // OVERLOADED RETURN-BY-VALUE OPERATORS - BigUnsigned operator +(const BigUnsigned &x) const; - BigUnsigned operator -(const BigUnsigned &x) const; - BigUnsigned operator *(const BigUnsigned &x) const; - BigUnsigned operator /(const BigUnsigned &x) const; - BigUnsigned operator %(const BigUnsigned &x) const; - /* OK, maybe unary minus could succeed in one case, but it really - * shouldn't be used, so it isn't provided. */ - BigUnsigned operator &(const BigUnsigned &x) const; - BigUnsigned operator |(const BigUnsigned &x) const; - BigUnsigned operator ^(const BigUnsigned &x) const; - BigUnsigned operator <<(int b) const; - BigUnsigned operator >>(int b) const; - - // OVERLOADED ASSIGNMENT OPERATORS - void operator +=(const BigUnsigned &x); - void operator -=(const BigUnsigned &x); - void operator *=(const BigUnsigned &x); - void operator /=(const BigUnsigned &x); - void operator %=(const BigUnsigned &x); - void operator &=(const BigUnsigned &x); - void operator |=(const BigUnsigned &x); - void operator ^=(const BigUnsigned &x); - void operator <<=(int b); - void operator >>=(int b); - - /* INCREMENT/DECREMENT OPERATORS - * To discourage messy coding, these do not return *this, so prefix - * and postfix behave the same. */ - void operator ++( ); - void operator ++(int); - void operator --( ); - void operator --(int); - - // Helper function that needs access to BigUnsigned internals - friend Blk getShiftedBlock(const BigUnsigned &num, Index x, - unsigned int y); - - // See BigInteger.cc. - template - friend X convertBigUnsignedToPrimitiveAccess(const BigUnsigned &a); -}; - -/* Implementing the return-by-value and assignment operators in terms of the - * copy-less operations. The copy-less operations are responsible for making - * any necessary temporary copies to work around aliasing. */ - -inline BigUnsigned BigUnsigned::operator +(const BigUnsigned &x) const { - BigUnsigned ans; - ans.add(*this, x); - return ans; -} -inline BigUnsigned BigUnsigned::operator -(const BigUnsigned &x) const { - BigUnsigned ans; - ans.subtract(*this, x); - return ans; -} -inline BigUnsigned BigUnsigned::operator *(const BigUnsigned &x) const { - BigUnsigned ans; - ans.multiply(*this, x); - return ans; -} -inline BigUnsigned BigUnsigned::operator /(const BigUnsigned &x) const { - if (x.isZero()) throw "BigUnsigned::operator /: division by zero"; - BigUnsigned q, r; - r = *this; - r.divideWithRemainder(x, q); - return q; -} -inline BigUnsigned BigUnsigned::operator %(const BigUnsigned &x) const { - if (x.isZero()) throw "BigUnsigned::operator %: division by zero"; - BigUnsigned q, r; - r = *this; - r.divideWithRemainder(x, q); - return r; -} -inline BigUnsigned BigUnsigned::operator &(const BigUnsigned &x) const { - BigUnsigned ans; - ans.bitAnd(*this, x); - return ans; -} -inline BigUnsigned BigUnsigned::operator |(const BigUnsigned &x) const { - BigUnsigned ans; - ans.bitOr(*this, x); - return ans; -} -inline BigUnsigned BigUnsigned::operator ^(const BigUnsigned &x) const { - BigUnsigned ans; - ans.bitXor(*this, x); - return ans; -} -inline BigUnsigned BigUnsigned::operator <<(int b) const { - BigUnsigned ans; - ans.bitShiftLeft(*this, b); - return ans; -} -inline BigUnsigned BigUnsigned::operator >>(int b) const { - BigUnsigned ans; - ans.bitShiftRight(*this, b); - return ans; -} - -inline void BigUnsigned::operator +=(const BigUnsigned &x) { - add(*this, x); -} -inline void BigUnsigned::operator -=(const BigUnsigned &x) { - subtract(*this, x); -} -inline void BigUnsigned::operator *=(const BigUnsigned &x) { - multiply(*this, x); -} -inline void BigUnsigned::operator /=(const BigUnsigned &x) { - if (x.isZero()) throw "BigUnsigned::operator /=: division by zero"; - /* The following technique is slightly faster than copying *this first - * when x is large. */ - BigUnsigned q; - divideWithRemainder(x, q); - // *this contains the remainder, but we overwrite it with the quotient. - *this = q; -} -inline void BigUnsigned::operator %=(const BigUnsigned &x) { - if (x.isZero()) throw "BigUnsigned::operator %=: division by zero"; - BigUnsigned q; - // Mods *this by x. Don't care about quotient left in q. - divideWithRemainder(x, q); -} -inline void BigUnsigned::operator &=(const BigUnsigned &x) { - bitAnd(*this, x); -} -inline void BigUnsigned::operator |=(const BigUnsigned &x) { - bitOr(*this, x); -} -inline void BigUnsigned::operator ^=(const BigUnsigned &x) { - bitXor(*this, x); -} -inline void BigUnsigned::operator <<=(int b) { - bitShiftLeft(*this, b); -} -inline void BigUnsigned::operator >>=(int b) { - bitShiftRight(*this, b); -} - -/* Templates for conversions of BigUnsigned to and from primitive integers. - * BigInteger.cc needs to instantiate convertToPrimitive, and the uses in - * BigUnsigned.cc didn't do the trick; I think g++ inlined convertToPrimitive - * instead of generating linkable instantiations. So for consistency, I put - * all the templates here. */ - -// CONSTRUCTION FROM PRIMITIVE INTEGERS - -/* Initialize this BigUnsigned from the given primitive integer. The same - * pattern works for all primitive integer types, so I put it into a template to - * reduce code duplication. (Don't worry: this is protected and we instantiate - * it only with primitive integer types.) Type X could be signed, but x is - * known to be nonnegative. */ -template -void BigUnsigned::initFromPrimitive(X x) { - if (x == 0) - ; // NumberlikeArray already initialized us to zero. - else { - // Create a single block. blk is NULL; no need to delete it. - cap = 1; - blk = new Blk[1]; - len = 1; - blk[0] = Blk(x); - } -} - -/* Ditto, but first check that x is nonnegative. I could have put the check in - * initFromPrimitive and let the compiler optimize it out for unsigned-type - * instantiations, but I wanted to avoid the warning stupidly issued by g++ for - * a condition that is constant in *any* instantiation, even if not in all. */ -template -void BigUnsigned::initFromSignedPrimitive(X x) { - if (x < 0) - throw "BigUnsigned constructor: " - "Cannot construct a BigUnsigned from a negative number"; - else - initFromPrimitive(x); -} - -// CONVERSION TO PRIMITIVE INTEGERS - -/* Template with the same idea as initFromPrimitive. This might be slightly - * slower than the previous version with the masks, but it's much shorter and - * clearer, which is the library's stated goal. */ -template -X BigUnsigned::convertToPrimitive() const { - if (len == 0) - // The number is zero; return zero. - return 0; - else if (len == 1) { - // The single block might fit in an X. Try the conversion. - X x = X(blk[0]); - // Make sure the result accurately represents the block. - if (Blk(x) == blk[0]) - // Successful conversion. - return x; - // Otherwise fall through. - } - throw "BigUnsigned::to: " - "Value is too big to fit in the requested type"; -} - -/* Wrap the above in an x >= 0 test to make sure we got a nonnegative result, - * not a negative one that happened to convert back into the correct nonnegative - * one. (E.g., catch incorrect conversion of 2^31 to the long -2^31.) Again, - * separated to avoid a g++ warning. */ -template -X BigUnsigned::convertToSignedPrimitive() const { - X x = convertToPrimitive(); - if (x >= 0) - return x; - else - throw "BigUnsigned::to(Primitive): " - "Value is too big to fit in the requested type"; -} - -#endif diff --git a/src/bigint/include/BigUnsignedInABase.hpp b/src/bigint/include/BigUnsignedInABase.hpp deleted file mode 100644 index 8893ef9..0000000 --- a/src/bigint/include/BigUnsignedInABase.hpp +++ /dev/null @@ -1,122 +0,0 @@ -#ifndef BIGUNSIGNEDINABASE_H -#define BIGUNSIGNEDINABASE_H - -#include "NumberlikeArray.hpp" -#include "BigUnsigned.hpp" -#include - -/* - * A BigUnsignedInABase object represents a nonnegative integer of size limited - * only by available memory, represented in a user-specified base that can fit - * in an `unsigned short' (most can, and this saves memory). - * - * BigUnsignedInABase is intended as an intermediary class with little - * functionality of its own. BigUnsignedInABase objects can be constructed - * from, and converted to, BigUnsigneds (requiring multiplication, mods, etc.) - * and `std::string's (by switching digit values for appropriate characters). - * - * BigUnsignedInABase is similar to BigUnsigned. Note the following: - * - * (1) They represent the number in exactly the same way, except that - * BigUnsignedInABase uses ``digits'' (or Digit) where BigUnsigned uses - * ``blocks'' (or Blk). - * - * (2) Both use the management features of NumberlikeArray. (In fact, my desire - * to add a BigUnsignedInABase class without duplicating a lot of code led me to - * introduce NumberlikeArray.) - * - * (3) The only arithmetic operation supported by BigUnsignedInABase is an - * equality test. Use BigUnsigned for arithmetic. - */ - -class BigUnsignedInABase : protected NumberlikeArray { - -public: - // The digits of a BigUnsignedInABase are unsigned shorts. - typedef unsigned short Digit; - // That's also the type of a base. - typedef Digit Base; - -protected: - // The base in which this BigUnsignedInABase is expressed - Base base; - - // Creates a BigUnsignedInABase with a capacity; for internal use. - BigUnsignedInABase(int, Index c) : NumberlikeArray(0, c) {} - - // Decreases len to eliminate any leading zero digits. - void zapLeadingZeros() { - while (len > 0 && blk[len - 1] == 0) - len--; - } - -public: - // Constructs zero in base 2. - BigUnsignedInABase() : NumberlikeArray(), base(2) {} - - // Copy constructor - BigUnsignedInABase(const BigUnsignedInABase &x) : NumberlikeArray(x), base(x.base) {} - - // Assignment operator - void operator =(const BigUnsignedInABase &x) { - NumberlikeArray::operator =(x); - base = x.base; - } - - // Constructor that copies from a given array of digits. - BigUnsignedInABase(const Digit *d, Index l, Base base); - - // Destructor. NumberlikeArray does the delete for us. - ~BigUnsignedInABase() {} - - // LINKS TO BIGUNSIGNED - BigUnsignedInABase(const BigUnsigned &x, Base base); - operator BigUnsigned() const; - - /* LINKS TO STRINGS - * - * These use the symbols ``0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'' to - * represent digits of 0 through 35. When parsing strings, lowercase is - * also accepted. - * - * All string representations are big-endian (big-place-value digits - * first). (Computer scientists have adopted zero-based counting; why - * can't they tolerate little-endian numbers?) - * - * No string representation has a ``base indicator'' like ``0x''. - * - * An exception is made for zero: it is converted to ``0'' and not the - * empty string. - * - * If you want different conventions, write your own routines to go - * between BigUnsignedInABase and strings. It's not hard. - */ - operator std::string() const; - BigUnsignedInABase(const std::string &s, Base base); - -public: - - // ACCESSORS - Base getBase() const { return base; } - - // Expose these from NumberlikeArray directly. - NumberlikeArray::getCapacity; - NumberlikeArray::getLength; - - /* Returns the requested digit, or 0 if it is beyond the length (as if - * the number had 0s infinitely to the left). */ - Digit getDigit(Index i) const { return i >= len ? 0 : blk[i]; } - - // The number is zero if and only if the canonical length is zero. - bool isZero() const { return NumberlikeArray::isEmpty(); } - - /* Equality test. For the purposes of this test, two BigUnsignedInABase - * values must have the same base to be equal. */ - bool operator ==(const BigUnsignedInABase &x) const { - return base == x.base && NumberlikeArray::operator ==(x); - } - bool operator !=(const BigUnsignedInABase &x) const { return !operator ==(x); } - -}; - -#endif diff --git a/src/bigint/include/NumberlikeArray.hpp b/src/bigint/include/NumberlikeArray.hpp deleted file mode 100644 index 53c8e5b..0000000 --- a/src/bigint/include/NumberlikeArray.hpp +++ /dev/null @@ -1,177 +0,0 @@ -#ifndef NUMBERLIKEARRAY_H -#define NUMBERLIKEARRAY_H - -// Make sure we have NULL. -#ifndef NULL -#define NULL 0 -#endif - -/* A NumberlikeArray object holds a heap-allocated array of Blk with a - * length and a capacity and provides basic memory management features. - * BigUnsigned and BigUnsignedInABase both subclass it. - * - * NumberlikeArray provides no information hiding. Subclasses should use - * nonpublic inheritance and manually expose members as desired using - * declarations like this: - * - * public: - * NumberlikeArray< the-type-argument >::getLength; - */ -template -class NumberlikeArray { -public: - - // Type for the index of a block in the array - typedef unsigned int Index; - // The number of bits in a block, defined below. - static const unsigned int N; - - // The current allocated capacity of this NumberlikeArray (in blocks) - Index cap; - // The actual length of the value stored in this NumberlikeArray (in blocks) - Index len; - // Heap-allocated array of the blocks (can be NULL if len == 0) - Blk *blk; - - // Constructs a ``zero'' NumberlikeArray with the given capacity. - NumberlikeArray(Index c) : cap(c), len(0) { - blk = (cap > 0) ? (new Blk[cap]) : NULL; - } - - /* Constructs a zero NumberlikeArray without allocating a backing array. - * A subclass that doesn't know the needed capacity at initialization - * time can use this constructor and then overwrite blk without first - * deleting it. */ - NumberlikeArray() : cap(0), len(0) { - blk = NULL; - } - - // Destructor. Note that `delete NULL' is a no-op. - ~NumberlikeArray() { - delete [] blk; - } - - /* Ensures that the array has at least the requested capacity; may - * destroy the contents. */ - void allocate(Index c); - - /* Ensures that the array has at least the requested capacity; does not - * destroy the contents. */ - void allocateAndCopy(Index c); - - // Copy constructor - NumberlikeArray(const NumberlikeArray &x); - - // Assignment operator - void operator=(const NumberlikeArray &x); - - // Constructor that copies from a given array of blocks - NumberlikeArray(const Blk *b, Index blen); - - // ACCESSORS - Index getCapacity() const { return cap; } - Index getLength() const { return len; } - Blk getBlock(Index i) const { return blk[i]; } - bool isEmpty() const { return len == 0; } - - /* Equality comparison: checks if both objects have the same length and - * equal (==) array elements to that length. Subclasses may wish to - * override. */ - bool operator ==(const NumberlikeArray &x) const; - - bool operator !=(const NumberlikeArray &x) const { - return !operator ==(x); - } -}; - -/* BEGIN TEMPLATE DEFINITIONS. They are present here so that source files that - * include this header file can generate the necessary real definitions. */ - -template -const unsigned int NumberlikeArray::N = 8 * sizeof(Blk); - -template -void NumberlikeArray::allocate(Index c) { - // If the requested capacity is more than the current capacity... - if (c > cap) { - // Delete the old number array - delete [] blk; - // Allocate the new array - cap = c; - blk = new Blk[cap]; - } -} - -template -void NumberlikeArray::allocateAndCopy(Index c) { - // If the requested capacity is more than the current capacity... - if (c > cap) { - Blk *oldBlk = blk; - // Allocate the new number array - cap = c; - blk = new Blk[cap]; - // Copy number blocks - Index i; - for (i = 0; i < len; i++) - blk[i] = oldBlk[i]; - // Delete the old array - delete [] oldBlk; - } -} - -template -NumberlikeArray::NumberlikeArray(const NumberlikeArray &x) - : len(x.len) { - // Create array - cap = len; - blk = new Blk[cap]; - // Copy blocks - Index i; - for (i = 0; i < len; i++) - blk[i] = x.blk[i]; -} - -template -void NumberlikeArray::operator=(const NumberlikeArray &x) { - /* Calls like a = a have no effect; catch them before the aliasing - * causes a problem */ - if (this == &x) - return; - // Copy length - len = x.len; - // Expand array if necessary - allocate(len); - // Copy number blocks - Index i; - for (i = 0; i < len; i++) - blk[i] = x.blk[i]; -} - -template -NumberlikeArray::NumberlikeArray(const Blk *b, Index blen) - : cap(blen), len(blen) { - // Create array - blk = new Blk[cap]; - // Copy blocks - Index i; - for (i = 0; i < len; i++) - blk[i] = b[i]; -} - -template -bool NumberlikeArray::operator ==(const NumberlikeArray &x) const { - if (len != x.len) - // Definitely unequal. - return false; - else { - // Compare corresponding blocks one by one. - Index i; - for (i = 0; i < len; i++) - if (blk[i] != x.blk[i]) - return false; - // No blocks differed, so the objects are equal. - return true; - } -} - -#endif diff --git a/src/bigint/src/BigInteger.cpp b/src/bigint/src/BigInteger.cpp deleted file mode 100644 index d12d4d3..0000000 --- a/src/bigint/src/BigInteger.cpp +++ /dev/null @@ -1,405 +0,0 @@ -#include "../include/BigInteger.hpp" - -void BigInteger::operator =(const BigInteger &x) { - // Calls like a = a have no effect - if (this == &x) - return; - // Copy sign - sign = x.sign; - // Copy the rest - mag = x.mag; -} - -BigInteger::BigInteger(const Blk *b, Index blen, Sign s) : mag(b, blen) { - switch (s) { - case zero: - if (!mag.isZero()) - throw "BigInteger::BigInteger(const Blk *, Index, Sign): Cannot use a sign of zero with a nonzero magnitude"; - sign = zero; - break; - case positive: - case negative: - // If the magnitude is zero, force the sign to zero. - sign = mag.isZero() ? zero : s; - break; - default: - /* g++ seems to be optimizing out this case on the assumption - * that the sign is a valid member of the enumeration. Oh well. */ - throw "BigInteger::BigInteger(const Blk *, Index, Sign): Invalid sign"; - } -} - -BigInteger::BigInteger(const BigUnsigned &x, Sign s) : mag(x) { - switch (s) { - case zero: - if (!mag.isZero()) - throw "BigInteger::BigInteger(const BigUnsigned &, Sign): Cannot use a sign of zero with a nonzero magnitude"; - sign = zero; - break; - case positive: - case negative: - // If the magnitude is zero, force the sign to zero. - sign = mag.isZero() ? zero : s; - break; - default: - /* g++ seems to be optimizing out this case on the assumption - * that the sign is a valid member of the enumeration. Oh well. */ - throw "BigInteger::BigInteger(const BigUnsigned &, Sign): Invalid sign"; - } -} - -/* CONSTRUCTION FROM PRIMITIVE INTEGERS - * Same idea as in BigUnsigned.cc, except that negative input results in a - * negative BigInteger instead of an exception. */ - -// Done longhand to let us use initialization. -BigInteger::BigInteger(unsigned long x) : mag(x) { sign = mag.isZero() ? zero : positive; } -BigInteger::BigInteger(unsigned int x) : mag(x) { sign = mag.isZero() ? zero : positive; } -BigInteger::BigInteger(unsigned short x) : mag(x) { sign = mag.isZero() ? zero : positive; } - -// For signed input, determine the desired magnitude and sign separately. - -namespace { - template - BigInteger::Blk magOf(X x) { - /* UX(...) cast needed to stop short(-2^15), which negates to - * itself, from sign-extending in the conversion to Blk. */ - return BigInteger::Blk(x < 0 ? UX(-x) : x); - } - template - BigInteger::Sign signOf(X x) { - return (x == 0) ? BigInteger::zero - : (x > 0) ? BigInteger::positive - : BigInteger::negative; - } -} - -BigInteger::BigInteger(long x) : sign(signOf(x)), mag(magOf(x)) {} -BigInteger::BigInteger(int x) : sign(signOf(x)), mag(magOf(x)) {} -BigInteger::BigInteger(short x) : sign(signOf(x)), mag(magOf(x)) {} - -// CONVERSION TO PRIMITIVE INTEGERS - -/* Reuse BigUnsigned's conversion to an unsigned primitive integer. - * The friend is a separate function rather than - * BigInteger::convertToUnsignedPrimitive to avoid requiring BigUnsigned to - * declare BigInteger. */ -template -inline X convertBigUnsignedToPrimitiveAccess(const BigUnsigned &a) { - return a.convertToPrimitive(); -} - -template -X BigInteger::convertToUnsignedPrimitive() const { - if (sign == negative) - throw "BigInteger::to: " - "Cannot convert a negative integer to an unsigned type"; - else - return convertBigUnsignedToPrimitiveAccess(mag); -} - -/* Similar to BigUnsigned::convertToPrimitive, but split into two cases for - * nonnegative and negative numbers. */ -template -X BigInteger::convertToSignedPrimitive() const { - if (sign == zero) - return 0; - else if (mag.getLength() == 1) { - // The single block might fit in an X. Try the conversion. - Blk b = mag.getBlock(0); - if (sign == positive) { - X x = X(b); - if (x >= 0 && Blk(x) == b) - return x; - } else { - X x = -X(b); - /* UX(...) needed to avoid rejecting conversion of - * -2^15 to a short. */ - if (x < 0 && Blk(UX(-x)) == b) - return x; - } - // Otherwise fall through. - } - throw "BigInteger::to: " - "Value is too big to fit in the requested type"; -} - -unsigned long BigInteger::toUnsignedLong () const { return convertToUnsignedPrimitive (); } -unsigned int BigInteger::toUnsignedInt () const { return convertToUnsignedPrimitive (); } -unsigned short BigInteger::toUnsignedShort() const { return convertToUnsignedPrimitive (); } -long BigInteger::toLong () const { return convertToSignedPrimitive (); } -int BigInteger::toInt () const { return convertToSignedPrimitive (); } -short BigInteger::toShort () const { return convertToSignedPrimitive (); } - -// COMPARISON -BigInteger::CmpRes BigInteger::compareTo(const BigInteger &x) const { - // A greater sign implies a greater number - if (sign < x.sign) - return less; - else if (sign > x.sign) - return greater; - else switch (sign) { - // If the signs are the same... - case zero: - return equal; // Two zeros are equal - case positive: - // Compare the magnitudes - return mag.compareTo(x.mag); - case negative: - // Compare the magnitudes, but return the opposite result - return CmpRes(-mag.compareTo(x.mag)); - default: - throw "BigInteger internal error"; - } -} - -/* COPY-LESS OPERATIONS - * These do some messing around to determine the sign of the result, - * then call one of BigUnsigned's copy-less operations. */ - -// See remarks about aliased calls in BigUnsigned.cc . -#define DTRT_ALIASED(cond, op) \ - if (cond) { \ - BigInteger tmpThis; \ - tmpThis.op; \ - *this = tmpThis; \ - return; \ - } - -void BigInteger::add(const BigInteger &a, const BigInteger &b) { - DTRT_ALIASED(this == &a || this == &b, add(a, b)); - // If one argument is zero, copy the other. - if (a.sign == zero) - operator =(b); - else if (b.sign == zero) - operator =(a); - // If the arguments have the same sign, take the - // common sign and add their magnitudes. - else if (a.sign == b.sign) { - sign = a.sign; - mag.add(a.mag, b.mag); - } else { - // Otherwise, their magnitudes must be compared. - switch (a.mag.compareTo(b.mag)) { - case equal: - // If their magnitudes are the same, copy zero. - mag = 0; - sign = zero; - break; - // Otherwise, take the sign of the greater, and subtract - // the lesser magnitude from the greater magnitude. - case greater: - sign = a.sign; - mag.subtract(a.mag, b.mag); - break; - case less: - sign = b.sign; - mag.subtract(b.mag, a.mag); - break; - } - } -} - -void BigInteger::subtract(const BigInteger &a, const BigInteger &b) { - // Notice that this routine is identical to BigInteger::add, - // if one replaces b.sign by its opposite. - DTRT_ALIASED(this == &a || this == &b, subtract(a, b)); - // If a is zero, copy b and flip its sign. If b is zero, copy a. - if (a.sign == zero) { - mag = b.mag; - // Take the negative of _b_'s, sign, not ours. - // Bug pointed out by Sam Larkin on 2005.03.30. - sign = Sign(-b.sign); - } else if (b.sign == zero) - operator =(a); - // If their signs differ, take a.sign and add the magnitudes. - else if (a.sign != b.sign) { - sign = a.sign; - mag.add(a.mag, b.mag); - } else { - // Otherwise, their magnitudes must be compared. - switch (a.mag.compareTo(b.mag)) { - // If their magnitudes are the same, copy zero. - case equal: - mag = 0; - sign = zero; - break; - // If a's magnitude is greater, take a.sign and - // subtract a from b. - case greater: - sign = a.sign; - mag.subtract(a.mag, b.mag); - break; - // If b's magnitude is greater, take the opposite - // of b.sign and subtract b from a. - case less: - sign = Sign(-b.sign); - mag.subtract(b.mag, a.mag); - break; - } - } -} - -void BigInteger::multiply(const BigInteger &a, const BigInteger &b) { - DTRT_ALIASED(this == &a || this == &b, multiply(a, b)); - // If one object is zero, copy zero and return. - if (a.sign == zero || b.sign == zero) { - sign = zero; - mag = 0; - return; - } - // If the signs of the arguments are the same, the result - // is positive, otherwise it is negative. - sign = (a.sign == b.sign) ? positive : negative; - // Multiply the magnitudes. - mag.multiply(a.mag, b.mag); -} - -/* - * DIVISION WITH REMAINDER - * Please read the comments before the definition of - * `BigUnsigned::divideWithRemainder' in `BigUnsigned.cc' for lots of - * information you should know before reading this function. - * - * Following Knuth, I decree that x / y is to be - * 0 if y==0 and floor(real-number x / y) if y!=0. - * Then x % y shall be x - y*(integer x / y). - * - * Note that x = y * (x / y) + (x % y) always holds. - * In addition, (x % y) is from 0 to y - 1 if y > 0, - * and from -(|y| - 1) to 0 if y < 0. (x % y) = x if y = 0. - * - * Examples: (q = a / b, r = a % b) - * a b q r - * === === === === - * 4 3 1 1 - * -4 3 -2 2 - * 4 -3 -2 -2 - * -4 -3 1 -1 - */ -void BigInteger::divideWithRemainder(const BigInteger &b, BigInteger &q) { - // Defend against aliased calls; - // same idea as in BigUnsigned::divideWithRemainder . - if (this == &q) - throw "BigInteger::divideWithRemainder: Cannot write quotient and remainder into the same variable"; - if (this == &b || &q == &b) { - BigInteger tmpB(b); - divideWithRemainder(tmpB, q); - return; - } - - // Division by zero gives quotient 0 and remainder *this - if (b.sign == zero) { - q.mag = 0; - q.sign = zero; - return; - } - // 0 / b gives quotient 0 and remainder 0 - if (sign == zero) { - q.mag = 0; - q.sign = zero; - return; - } - - // Here *this != 0, b != 0. - - // Do the operands have the same sign? - if (sign == b.sign) { - // Yes: easy case. Quotient is zero or positive. - q.sign = positive; - } else { - // No: harder case. Quotient is negative. - q.sign = negative; - // Decrease the magnitude of the dividend by one. - mag--; - /* - * We tinker with the dividend before and with the - * quotient and remainder after so that the result - * comes out right. To see why it works, consider the following - * list of examples, where A is the magnitude-decreased - * a, Q and R are the results of BigUnsigned division - * with remainder on A and |b|, and q and r are the - * final results we want: - * - * a A b Q R q r - * -3 -2 3 0 2 -1 0 - * -4 -3 3 1 0 -2 2 - * -5 -4 3 1 1 -2 1 - * -6 -5 3 1 2 -2 0 - * - * It appears that we need a total of 3 corrections: - * Decrease the magnitude of a to get A. Increase the - * magnitude of Q to get q (and make it negative). - * Find r = (b - 1) - R and give it the desired sign. - */ - } - - // Divide the magnitudes. - mag.divideWithRemainder(b.mag, q.mag); - - if (sign != b.sign) { - // More for the harder case (as described): - // Increase the magnitude of the quotient by one. - q.mag++; - // Modify the remainder. - mag.subtract(b.mag, mag); - mag--; - } - - // Sign of the remainder is always the sign of the divisor b. - sign = b.sign; - - // Set signs to zero as necessary. (Thanks David Allen!) - if (mag.isZero()) - sign = zero; - if (q.mag.isZero()) - q.sign = zero; - - // WHEW!!! -} - -// Negation -void BigInteger::negate(const BigInteger &a) { - DTRT_ALIASED(this == &a, negate(a)); - // Copy a's magnitude - mag = a.mag; - // Copy the opposite of a.sign - sign = Sign(-a.sign); -} - -// INCREMENT/DECREMENT OPERATORS - -// Prefix increment -void BigInteger::operator ++() { - if (sign == negative) { - mag--; - if (mag == 0) - sign = zero; - } else { - mag++; - sign = positive; // if not already - } -} - -// Postfix increment: same as prefix -void BigInteger::operator ++(int) { - operator ++(); -} - -// Prefix decrement -void BigInteger::operator --() { - if (sign == positive) { - mag--; - if (mag == 0) - sign = zero; - } else { - mag++; - sign = negative; - } -} - -// Postfix decrement: same as prefix -void BigInteger::operator --(int) { - operator --(); -} - diff --git a/src/bigint/src/BigIntegerAlgorithms.cpp b/src/bigint/src/BigIntegerAlgorithms.cpp deleted file mode 100644 index c84799d..0000000 --- a/src/bigint/src/BigIntegerAlgorithms.cpp +++ /dev/null @@ -1,71 +0,0 @@ -#include "../include/BigIntegerAlgorithms.hpp" - -BigUnsigned gcd(BigUnsigned a, BigUnsigned b) { - BigUnsigned trash; - // Neat in-place alternating technique. - for (;;) { - if (b.isZero()) - return a; - a.divideWithRemainder(b, trash); - if (a.isZero()) - return b; - b.divideWithRemainder(a, trash); - } -} - -void extendedEuclidean(BigInteger m, BigInteger n, - BigInteger &g, BigInteger &r, BigInteger &s) { - if (&g == &r || &g == &s || &r == &s) { - throw "BigInteger extendedEuclidean: Outputs are aliased"; - } - BigInteger r1(1), s1(0), r2(0), s2(1), q; - /* Invariants: - * r1*m(orig) + s1*n(orig) == m(current) - * r2*m(orig) + s2*n(orig) == n(current) */ - for (;;) { - if (n.isZero()) { - r = r1; s = s1; g = m; - return; - } - // Subtract q times the second invariant from the first invariant. - m.divideWithRemainder(n, q); - r1 -= q*r2; s1 -= q*s2; - - if (m.isZero()) { - r = r2; s = s2; g = n; - return; - } - // Subtract q times the first invariant from the second invariant. - n.divideWithRemainder(m, q); - r2 -= q*r1; s2 -= q*s1; - } -} - -BigUnsigned modinv(const BigInteger &x, const BigUnsigned &n) { - BigInteger g, r, s; - extendedEuclidean(x, n, g, r, s); - if (g == 1) - // r*x + s*n == 1, so r*x === 1 (mod n), so r is the answer. - return (r % n).getMagnitude(); // (r % n) will be nonnegative - else - throw "BigInteger modinv: x and n have a common factor"; -} - -BigUnsigned modexp(const BigInteger &base, const BigUnsigned &exponent, - const BigUnsigned &modulus) { - BigUnsigned ans = 1, base2 = (base % modulus).getMagnitude(); - BigUnsigned::Index i = exponent.bitLength(); - // For each bit of the exponent, most to least significant... - while (i > 0) { - i--; - // Square. - ans *= ans; - ans %= modulus; - // And multiply if the bit is a 1. - if (exponent.getBit(i)) { - ans *= base2; - ans %= modulus; - } - } - return ans; -} diff --git a/src/bigint/src/BigIntegerUtils.cpp b/src/bigint/src/BigIntegerUtils.cpp deleted file mode 100644 index 7157b8a..0000000 --- a/src/bigint/src/BigIntegerUtils.cpp +++ /dev/null @@ -1,50 +0,0 @@ -#include "../include/BigIntegerUtils.hpp" -#include "../include/BigUnsignedInABase.hpp" - -std::string bigUnsignedToString(const BigUnsigned &x) { - return std::string(BigUnsignedInABase(x, 10)); -} - -std::string bigIntegerToString(const BigInteger &x) { - return (x.getSign() == BigInteger::negative) - ? (std::string("-") + bigUnsignedToString(x.getMagnitude())) - : (bigUnsignedToString(x.getMagnitude())); -} - -BigUnsigned stringToBigUnsigned(const std::string &s) { - return BigUnsigned(BigUnsignedInABase(s, 10)); -} - -BigInteger stringToBigInteger(const std::string &s) { - // Recognize a sign followed by a BigUnsigned. - return (s[0] == '-') ? BigInteger(stringToBigUnsigned(s.substr(1, s.length() - 1)), BigInteger::negative) - : (s[0] == '+') ? BigInteger(stringToBigUnsigned(s.substr(1, s.length() - 1))) - : BigInteger(stringToBigUnsigned(s)); -} - -std::ostream &operator <<(std::ostream &os, const BigUnsigned &x) { - BigUnsignedInABase::Base base; - long osFlags = os.flags(); - if (osFlags & os.dec) - base = 10; - else if (osFlags & os.hex) { - base = 16; - if (osFlags & os.showbase) - os << "0x"; - } else if (osFlags & os.oct) { - base = 8; - if (osFlags & os.showbase) - os << '0'; - } else - throw "std::ostream << BigUnsigned: Could not determine the desired base from output-stream flags"; - std::string s = std::string(BigUnsignedInABase(x, base)); - os << s; - return os; -} - -std::ostream &operator <<(std::ostream &os, const BigInteger &x) { - if (x.getSign() == BigInteger::negative) - os << '-'; - os << x.getMagnitude(); - return os; -} diff --git a/src/bigint/src/BigUnsigned.cpp b/src/bigint/src/BigUnsigned.cpp deleted file mode 100644 index f5b05ed..0000000 --- a/src/bigint/src/BigUnsigned.cpp +++ /dev/null @@ -1,697 +0,0 @@ -#include "../include/BigUnsigned.hpp" - -// Memory management definitions have moved to the bottom of NumberlikeArray.hh. - -// The templates used by these constructors and converters are at the bottom of -// BigUnsigned.hh. - -BigUnsigned::BigUnsigned(unsigned long x) { initFromPrimitive (x); } -BigUnsigned::BigUnsigned(unsigned int x) { initFromPrimitive (x); } -BigUnsigned::BigUnsigned(unsigned short x) { initFromPrimitive (x); } -BigUnsigned::BigUnsigned( long x) { initFromSignedPrimitive(x); } -BigUnsigned::BigUnsigned( int x) { initFromSignedPrimitive(x); } -BigUnsigned::BigUnsigned( short x) { initFromSignedPrimitive(x); } - -unsigned long BigUnsigned::toUnsignedLong () const { return convertToPrimitive (); } -unsigned int BigUnsigned::toUnsignedInt () const { return convertToPrimitive (); } -unsigned short BigUnsigned::toUnsignedShort() const { return convertToPrimitive (); } -long BigUnsigned::toLong () const { return convertToSignedPrimitive< long >(); } -int BigUnsigned::toInt () const { return convertToSignedPrimitive< int >(); } -short BigUnsigned::toShort () const { return convertToSignedPrimitive< short>(); } - -// BIT/BLOCK ACCESSORS - -void BigUnsigned::setBlock(Index i, Blk newBlock) { - if (newBlock == 0) { - if (i < len) { - blk[i] = 0; - zapLeadingZeros(); - } - // If i >= len, no effect. - } else { - if (i >= len) { - // The nonzero block extends the number. - allocateAndCopy(i+1); - // Zero any added blocks that we aren't setting. - for (Index j = len; j < i; j++) - blk[j] = 0; - len = i+1; - } - blk[i] = newBlock; - } -} - -/* Evidently the compiler wants BigUnsigned:: on the return type because, at - * that point, it hasn't yet parsed the BigUnsigned:: on the name to get the - * proper scope. */ -BigUnsigned::Index BigUnsigned::bitLength() const { - if (isZero()) - return 0; - else { - Blk leftmostBlock = getBlock(len - 1); - Index leftmostBlockLen = 0; - while (leftmostBlock != 0) { - leftmostBlock >>= 1; - leftmostBlockLen++; - } - return leftmostBlockLen + (len - 1) * N; - } -} - -void BigUnsigned::setBit(Index bi, bool newBit) { - Index blockI = bi / N; - Blk block = getBlock(blockI), mask = Blk(1) << (bi % N); - block = newBit ? (block | mask) : (block & ~mask); - setBlock(blockI, block); -} - -// COMPARISON -BigUnsigned::CmpRes BigUnsigned::compareTo(const BigUnsigned &x) const { - // A bigger length implies a bigger number. - if (len < x.len) - return less; - else if (len > x.len) - return greater; - else { - // Compare blocks one by one from left to right. - Index i = len; - while (i > 0) { - i--; - if (blk[i] == x.blk[i]) - continue; - else if (blk[i] > x.blk[i]) - return greater; - else - return less; - } - // If no blocks differed, the numbers are equal. - return equal; - } -} - -// COPY-LESS OPERATIONS - -/* - * On most calls to copy-less operations, it's safe to read the inputs little by - * little and write the outputs little by little. However, if one of the - * inputs is coming from the same variable into which the output is to be - * stored (an "aliased" call), we risk overwriting the input before we read it. - * In this case, we first compute the result into a temporary BigUnsigned - * variable and then copy it into the requested output variable *this. - * Each put-here operation uses the DTRT_ALIASED macro (Do The Right Thing on - * aliased calls) to generate code for this check. - * - * I adopted this approach on 2007.02.13 (see Assignment Operators in - * BigUnsigned.hh). Before then, put-here operations rejected aliased calls - * with an exception. I think doing the right thing is better. - * - * Some of the put-here operations can probably handle aliased calls safely - * without the extra copy because (for example) they process blocks strictly - * right-to-left. At some point I might determine which ones don't need the - * copy, but my reasoning would need to be verified very carefully. For now - * I'll leave in the copy. - */ -#define DTRT_ALIASED(cond, op) \ - if (cond) { \ - BigUnsigned tmpThis; \ - tmpThis.op; \ - *this = tmpThis; \ - return; \ - } - - - -void BigUnsigned::add(const BigUnsigned &a, const BigUnsigned &b) { - DTRT_ALIASED(this == &a || this == &b, add(a, b)); - // If one argument is zero, copy the other. - if (a.len == 0) { - operator =(b); - return; - } else if (b.len == 0) { - operator =(a); - return; - } - // Some variables... - // Carries in and out of an addition stage - bool carryIn, carryOut; - Blk temp; - Index i; - // a2 points to the longer input, b2 points to the shorter - const BigUnsigned *a2, *b2; - if (a.len >= b.len) { - a2 = &a; - b2 = &b; - } else { - a2 = &b; - b2 = &a; - } - // Set prelimiary length and make room in this BigUnsigned - len = a2->len + 1; - allocate(len); - // For each block index that is present in both inputs... - for (i = 0, carryIn = false; i < b2->len; i++) { - // Add input blocks - temp = a2->blk[i] + b2->blk[i]; - // If a rollover occurred, the result is less than either input. - // This test is used many times in the BigUnsigned code. - carryOut = (temp < a2->blk[i]); - // If a carry was input, handle it - if (carryIn) { - temp++; - carryOut |= (temp == 0); - } - blk[i] = temp; // Save the addition result - carryIn = carryOut; // Pass the carry along - } - // If there is a carry left over, increase blocks until - // one does not roll over. - for (; i < a2->len && carryIn; i++) { - temp = a2->blk[i] + 1; - carryIn = (temp == 0); - blk[i] = temp; - } - // If the carry was resolved but the larger number - // still has blocks, copy them over. - for (; i < a2->len; i++) - blk[i] = a2->blk[i]; - // Set the extra block if there's still a carry, decrease length otherwise - if (carryIn) - blk[i] = 1; - else - len--; -} - -void BigUnsigned::subtract(const BigUnsigned &a, const BigUnsigned &b) { - DTRT_ALIASED(this == &a || this == &b, subtract(a, b)); - if (b.len == 0) { - // If b is zero, copy a. - operator =(a); - return; - } else if (a.len < b.len) - // If a is shorter than b, the result is negative. - throw "BigUnsigned::subtract: " - "Negative result in unsigned calculation"; - // Some variables... - bool borrowIn, borrowOut; - Blk temp; - Index i; - // Set preliminary length and make room - len = a.len; - allocate(len); - // For each block index that is present in both inputs... - for (i = 0, borrowIn = false; i < b.len; i++) { - temp = a.blk[i] - b.blk[i]; - // If a reverse rollover occurred, - // the result is greater than the block from a. - borrowOut = (temp > a.blk[i]); - // Handle an incoming borrow - if (borrowIn) { - borrowOut |= (temp == 0); - temp--; - } - blk[i] = temp; // Save the subtraction result - borrowIn = borrowOut; // Pass the borrow along - } - // If there is a borrow left over, decrease blocks until - // one does not reverse rollover. - for (; i < a.len && borrowIn; i++) { - borrowIn = (a.blk[i] == 0); - blk[i] = a.blk[i] - 1; - } - /* If there's still a borrow, the result is negative. - * Throw an exception, but zero out this object so as to leave it in a - * predictable state. */ - if (borrowIn) { - len = 0; - throw "BigUnsigned::subtract: Negative result in unsigned calculation"; - } else - // Copy over the rest of the blocks - for (; i < a.len; i++) - blk[i] = a.blk[i]; - // Zap leading zeros - zapLeadingZeros(); -} - -/* - * About the multiplication and division algorithms: - * - * I searched unsucessfully for fast C++ built-in operations like the `b_0' - * and `c_0' Knuth describes in Section 4.3.1 of ``The Art of Computer - * Programming'' (replace `place' by `Blk'): - * - * ``b_0[:] multiplication of a one-place integer by another one-place - * integer, giving a two-place answer; - * - * ``c_0[:] division of a two-place integer by a one-place integer, - * provided that the quotient is a one-place integer, and yielding - * also a one-place remainder.'' - * - * I also missed his note that ``[b]y adjusting the word size, if - * necessary, nearly all computers will have these three operations - * available'', so I gave up on trying to use algorithms similar to his. - * A future version of the library might include such algorithms; I - * would welcome contributions from others for this. - * - * I eventually decided to use bit-shifting algorithms. To multiply `a' - * and `b', we zero out the result. Then, for each `1' bit in `a', we - * shift `b' left the appropriate amount and add it to the result. - * Similarly, to divide `a' by `b', we shift `b' left varying amounts, - * repeatedly trying to subtract it from `a'. When we succeed, we note - * the fact by setting a bit in the quotient. While these algorithms - * have the same O(n^2) time complexity as Knuth's, the ``constant factor'' - * is likely to be larger. - * - * Because I used these algorithms, which require single-block addition - * and subtraction rather than single-block multiplication and division, - * the innermost loops of all four routines are very similar. Study one - * of them and all will become clear. - */ - -/* - * This is a little inline function used by both the multiplication - * routine and the division routine. - * - * `getShiftedBlock' returns the `x'th block of `num << y'. - * `y' may be anything from 0 to N - 1, and `x' may be anything from - * 0 to `num.len'. - * - * Two things contribute to this block: - * - * (1) The `N - y' low bits of `num.blk[x]', shifted `y' bits left. - * - * (2) The `y' high bits of `num.blk[x-1]', shifted `N - y' bits right. - * - * But we must be careful if `x == 0' or `x == num.len', in - * which case we should use 0 instead of (2) or (1), respectively. - * - * If `y == 0', then (2) contributes 0, as it should. However, - * in some computer environments, for a reason I cannot understand, - * `a >> b' means `a >> (b % N)'. This means `num.blk[x-1] >> (N - y)' - * will return `num.blk[x-1]' instead of the desired 0 when `y == 0'; - * the test `y == 0' handles this case specially. - */ -inline BigUnsigned::Blk getShiftedBlock(const BigUnsigned &num, - BigUnsigned::Index x, unsigned int y) { - BigUnsigned::Blk part1 = (x == 0 || y == 0) ? 0 : (num.blk[x - 1] >> (BigUnsigned::N - y)); - BigUnsigned::Blk part2 = (x == num.len) ? 0 : (num.blk[x] << y); - return part1 | part2; -} - -void BigUnsigned::multiply(const BigUnsigned &a, const BigUnsigned &b) { - DTRT_ALIASED(this == &a || this == &b, multiply(a, b)); - // If either a or b is zero, set to zero. - if (a.len == 0 || b.len == 0) { - len = 0; - return; - } - /* - * Overall method: - * - * Set this = 0. - * For each 1-bit of `a' (say the `i2'th bit of block `i'): - * Add `b << (i blocks and i2 bits)' to *this. - */ - // Variables for the calculation - Index i, j, k; - unsigned int i2; - Blk temp; - bool carryIn, carryOut; - // Set preliminary length and make room - len = a.len + b.len; - allocate(len); - // Zero out this object - for (i = 0; i < len; i++) - blk[i] = 0; - // For each block of the first number... - for (i = 0; i < a.len; i++) { - // For each 1-bit of that block... - for (i2 = 0; i2 < N; i2++) { - if ((a.blk[i] & (Blk(1) << i2)) == 0) - continue; - /* - * Add b to this, shifted left i blocks and i2 bits. - * j is the index in b, and k = i + j is the index in this. - * - * `getShiftedBlock', a short inline function defined above, - * is now used for the bit handling. It replaces the more - * complex `bHigh' code, in which each run of the loop dealt - * immediately with the low bits and saved the high bits to - * be picked up next time. The last run of the loop used to - * leave leftover high bits, which were handled separately. - * Instead, this loop runs an additional time with j == b.len. - * These changes were made on 2005.01.11. - */ - for (j = 0, k = i, carryIn = false; j <= b.len; j++, k++) { - /* - * The body of this loop is very similar to the body of the first loop - * in `add', except that this loop does a `+=' instead of a `+'. - */ - temp = blk[k] + getShiftedBlock(b, j, i2); - carryOut = (temp < blk[k]); - if (carryIn) { - temp++; - carryOut |= (temp == 0); - } - blk[k] = temp; - carryIn = carryOut; - } - // No more extra iteration to deal with `bHigh'. - // Roll-over a carry as necessary. - for (; carryIn; k++) { - blk[k]++; - carryIn = (blk[k] == 0); - } - } - } - // Zap possible leading zero - if (blk[len - 1] == 0) - len--; -} - -/* - * DIVISION WITH REMAINDER - * This monstrous function mods *this by the given divisor b while storing the - * quotient in the given object q; at the end, *this contains the remainder. - * The seemingly bizarre pattern of inputs and outputs was chosen so that the - * function copies as little as possible (since it is implemented by repeated - * subtraction of multiples of b from *this). - * - * "modWithQuotient" might be a better name for this function, but I would - * rather not change the name now. - */ -void BigUnsigned::divideWithRemainder(const BigUnsigned &b, BigUnsigned &q) { - /* Defending against aliased calls is more complex than usual because we - * are writing to both *this and q. - * - * It would be silly to try to write quotient and remainder to the - * same variable. Rule that out right away. */ - if (this == &q) - throw "BigUnsigned::divideWithRemainder: Cannot write quotient and remainder into the same variable"; - /* Now *this and q are separate, so the only concern is that b might be - * aliased to one of them. If so, use a temporary copy of b. */ - if (this == &b || &q == &b) { - BigUnsigned tmpB(b); - divideWithRemainder(tmpB, q); - return; - } - - /* - * Knuth's definition of mod (which this function uses) is somewhat - * different from the C++ definition of % in case of division by 0. - * - * We let a / 0 == 0 (it doesn't matter much) and a % 0 == a, no - * exceptions thrown. This allows us to preserve both Knuth's demand - * that a mod 0 == a and the useful property that - * (a / b) * b + (a % b) == a. - */ - if (b.len == 0) { - q.len = 0; - return; - } - - /* - * If *this.len < b.len, then *this < b, and we can be sure that b doesn't go into - * *this at all. The quotient is 0 and *this is already the remainder (so leave it alone). - */ - if (len < b.len) { - q.len = 0; - return; - } - - // At this point we know (*this).len >= b.len > 0. (Whew!) - - /* - * Overall method: - * - * For each appropriate i and i2, decreasing: - * Subtract (b << (i blocks and i2 bits)) from *this, storing the - * result in subtractBuf. - * If the subtraction succeeds with a nonnegative result: - * Turn on bit i2 of block i of the quotient q. - * Copy subtractBuf back into *this. - * Otherwise bit i2 of block i remains off, and *this is unchanged. - * - * Eventually q will contain the entire quotient, and *this will - * be left with the remainder. - * - * subtractBuf[x] corresponds to blk[x], not blk[x+i], since 2005.01.11. - * But on a single iteration, we don't touch the i lowest blocks of blk - * (and don't use those of subtractBuf) because these blocks are - * unaffected by the subtraction: we are subtracting - * (b << (i blocks and i2 bits)), which ends in at least `i' zero - * blocks. */ - // Variables for the calculation - Index i, j, k; - unsigned int i2; - Blk temp; - bool borrowIn, borrowOut; - - /* - * Make sure we have an extra zero block just past the value. - * - * When we attempt a subtraction, we might shift `b' so - * its first block begins a few bits left of the dividend, - * and then we'll try to compare these extra bits with - * a nonexistent block to the left of the dividend. The - * extra zero block ensures sensible behavior; we need - * an extra block in `subtractBuf' for exactly the same reason. - */ - Index origLen = len; // Save real length. - /* To avoid an out-of-bounds access in case of reallocation, allocate - * first and then increment the logical length. */ - allocateAndCopy(len + 1); - len++; - blk[origLen] = 0; // Zero the added block. - - // subtractBuf holds part of the result of a subtraction; see above. - Blk *subtractBuf = new Blk[len]; - - // Set preliminary length for quotient and make room - q.len = origLen - b.len + 1; - q.allocate(q.len); - // Zero out the quotient - for (i = 0; i < q.len; i++) - q.blk[i] = 0; - - // For each possible left-shift of b in blocks... - i = q.len; - while (i > 0) { - i--; - // For each possible left-shift of b in bits... - // (Remember, N is the number of bits in a Blk.) - q.blk[i] = 0; - i2 = N; - while (i2 > 0) { - i2--; - /* - * Subtract b, shifted left i blocks and i2 bits, from *this, - * and store the answer in subtractBuf. In the for loop, `k == i + j'. - * - * Compare this to the middle section of `multiply'. They - * are in many ways analogous. See especially the discussion - * of `getShiftedBlock'. - */ - for (j = 0, k = i, borrowIn = false; j <= b.len; j++, k++) { - temp = blk[k] - getShiftedBlock(b, j, i2); - borrowOut = (temp > blk[k]); - if (borrowIn) { - borrowOut |= (temp == 0); - temp--; - } - // Since 2005.01.11, indices of `subtractBuf' directly match those of `blk', so use `k'. - subtractBuf[k] = temp; - borrowIn = borrowOut; - } - // No more extra iteration to deal with `bHigh'. - // Roll-over a borrow as necessary. - for (; k < origLen && borrowIn; k++) { - borrowIn = (blk[k] == 0); - subtractBuf[k] = blk[k] - 1; - } - /* - * If the subtraction was performed successfully (!borrowIn), - * set bit i2 in block i of the quotient. - * - * Then, copy the portion of subtractBuf filled by the subtraction - * back to *this. This portion starts with block i and ends-- - * where? Not necessarily at block `i + b.len'! Well, we - * increased k every time we saved a block into subtractBuf, so - * the region of subtractBuf we copy is just [i, k). - */ - if (!borrowIn) { - q.blk[i] |= (Blk(1) << i2); - while (k > i) { - k--; - blk[k] = subtractBuf[k]; - } - } - } - } - // Zap possible leading zero in quotient - if (q.blk[q.len - 1] == 0) - q.len--; - // Zap any/all leading zeros in remainder - zapLeadingZeros(); - // Deallocate subtractBuf. - // (Thanks to Brad Spencer for noticing my accidental omission of this!) - delete [] subtractBuf; -} - -/* BITWISE OPERATORS - * These are straightforward blockwise operations except that they differ in - * the output length and the necessity of zapLeadingZeros. */ - -void BigUnsigned::bitAnd(const BigUnsigned &a, const BigUnsigned &b) { - DTRT_ALIASED(this == &a || this == &b, bitAnd(a, b)); - // The bitwise & can't be longer than either operand. - len = (a.len >= b.len) ? b.len : a.len; - allocate(len); - Index i; - for (i = 0; i < len; i++) - blk[i] = a.blk[i] & b.blk[i]; - zapLeadingZeros(); -} - -void BigUnsigned::bitOr(const BigUnsigned &a, const BigUnsigned &b) { - DTRT_ALIASED(this == &a || this == &b, bitOr(a, b)); - Index i; - const BigUnsigned *a2, *b2; - if (a.len >= b.len) { - a2 = &a; - b2 = &b; - } else { - a2 = &b; - b2 = &a; - } - allocate(a2->len); - for (i = 0; i < b2->len; i++) - blk[i] = a2->blk[i] | b2->blk[i]; - for (; i < a2->len; i++) - blk[i] = a2->blk[i]; - len = a2->len; - // Doesn't need zapLeadingZeros. -} - -void BigUnsigned::bitXor(const BigUnsigned &a, const BigUnsigned &b) { - DTRT_ALIASED(this == &a || this == &b, bitXor(a, b)); - Index i; - const BigUnsigned *a2, *b2; - if (a.len >= b.len) { - a2 = &a; - b2 = &b; - } else { - a2 = &b; - b2 = &a; - } - allocate(a2->len); - for (i = 0; i < b2->len; i++) - blk[i] = a2->blk[i] ^ b2->blk[i]; - for (; i < a2->len; i++) - blk[i] = a2->blk[i]; - len = a2->len; - zapLeadingZeros(); -} - -void BigUnsigned::bitShiftLeft(const BigUnsigned &a, int b) { - DTRT_ALIASED(this == &a, bitShiftLeft(a, b)); - if (b < 0) { - if (b << 1 == 0) - throw "BigUnsigned::bitShiftLeft: " - "Pathological shift amount not implemented"; - else { - bitShiftRight(a, -b); - return; - } - } - Index shiftBlocks = b / N; - unsigned int shiftBits = b % N; - // + 1: room for high bits nudged left into another block - len = a.len + shiftBlocks + 1; - allocate(len); - Index i, j; - for (i = 0; i < shiftBlocks; i++) - blk[i] = 0; - for (j = 0, i = shiftBlocks; j <= a.len; j++, i++) - blk[i] = getShiftedBlock(a, j, shiftBits); - // Zap possible leading zero - if (blk[len - 1] == 0) - len--; -} - -void BigUnsigned::bitShiftRight(const BigUnsigned &a, int b) { - DTRT_ALIASED(this == &a, bitShiftRight(a, b)); - if (b < 0) { - if (b << 1 == 0) - throw "BigUnsigned::bitShiftRight: " - "Pathological shift amount not implemented"; - else { - bitShiftLeft(a, -b); - return; - } - } - // This calculation is wacky, but expressing the shift as a left bit shift - // within each block lets us use getShiftedBlock. - Index rightShiftBlocks = (b + N - 1) / N; - unsigned int leftShiftBits = N * rightShiftBlocks - b; - // Now (N * rightShiftBlocks - leftShiftBits) == b - // and 0 <= leftShiftBits < N. - if (rightShiftBlocks >= a.len + 1) { - // All of a is guaranteed to be shifted off, even considering the left - // bit shift. - len = 0; - return; - } - // Now we're allocating a positive amount. - // + 1: room for high bits nudged left into another block - len = a.len + 1 - rightShiftBlocks; - allocate(len); - Index i, j; - for (j = rightShiftBlocks, i = 0; j <= a.len; j++, i++) - blk[i] = getShiftedBlock(a, j, leftShiftBits); - // Zap possible leading zero - if (blk[len - 1] == 0) - len--; -} - -// INCREMENT/DECREMENT OPERATORS - -// Prefix increment -void BigUnsigned::operator ++() { - Index i; - bool carry = true; - for (i = 0; i < len && carry; i++) { - blk[i]++; - carry = (blk[i] == 0); - } - if (carry) { - // Allocate and then increase length, as in divideWithRemainder - allocateAndCopy(len + 1); - len++; - blk[i] = 1; - } -} - -// Postfix increment: same as prefix -void BigUnsigned::operator ++(int) { - operator ++(); -} - -// Prefix decrement -void BigUnsigned::operator --() { - if (len == 0) - throw "BigUnsigned::operator --(): Cannot decrement an unsigned zero"; - Index i; - bool borrow = true; - for (i = 0; borrow; i++) { - borrow = (blk[i] == 0); - blk[i]--; - } - // Zap possible leading zero (there can only be one) - if (blk[len - 1] == 0) - len--; -} - -// Postfix decrement: same as prefix -void BigUnsigned::operator --(int) { - operator --(); -} diff --git a/src/bigint/src/BigUnsignedInABase.cpp b/src/bigint/src/BigUnsignedInABase.cpp deleted file mode 100644 index 0289e94..0000000 --- a/src/bigint/src/BigUnsignedInABase.cpp +++ /dev/null @@ -1,125 +0,0 @@ -#include "../include/BigUnsignedInABase.hpp" - -BigUnsignedInABase::BigUnsignedInABase(const Digit *d, Index l, Base base) - : NumberlikeArray(d, l), base(base) { - // Check the base - if (base < 2) - throw "BigUnsignedInABase::BigUnsignedInABase(const Digit *, Index, Base): The base must be at least 2"; - - // Validate the digits. - for (Index i = 0; i < l; i++) - if (blk[i] >= base) - throw "BigUnsignedInABase::BigUnsignedInABase(const Digit *, Index, Base): A digit is too large for the specified base"; - - // Eliminate any leading zeros we may have been passed. - zapLeadingZeros(); -} - -namespace { - unsigned int bitLen(unsigned int x) { - unsigned int len = 0; - while (x > 0) { - x >>= 1; - len++; - } - return len; - } - unsigned int ceilingDiv(unsigned int a, unsigned int b) { - return (a + b - 1) / b; - } -} - -BigUnsignedInABase::BigUnsignedInABase(const BigUnsigned &x, Base base) { - // Check the base - if (base < 2) - throw "BigUnsignedInABase(BigUnsigned, Base): The base must be at least 2"; - this->base = base; - - // Get an upper bound on how much space we need - int maxBitLenOfX = x.getLength() * BigUnsigned::N; - int minBitsPerDigit = bitLen(base) - 1; - int maxDigitLenOfX = ceilingDiv(maxBitLenOfX, minBitsPerDigit); - len = maxDigitLenOfX; // Another change to comply with `staying in bounds'. - allocate(len); // Get the space - - BigUnsigned x2(x), buBase(base); - Index digitNum = 0; - - while (!x2.isZero()) { - // Get last digit. This is like `lastDigit = x2 % buBase, x2 /= buBase'. - BigUnsigned lastDigit(x2); - lastDigit.divideWithRemainder(buBase, x2); - // Save the digit. - blk[digitNum] = lastDigit.toUnsignedShort(); - // Move on. We can't run out of room: we figured it out above. - digitNum++; - } - - // Save the actual length. - len = digitNum; -} - -BigUnsignedInABase::operator BigUnsigned() const { - BigUnsigned ans(0), buBase(base), temp; - Index digitNum = len; - while (digitNum > 0) { - digitNum--; - temp.multiply(ans, buBase); - ans.add(temp, BigUnsigned(blk[digitNum])); - } - return ans; -} - -BigUnsignedInABase::BigUnsignedInABase(const std::string &s, Base base) { - // Check the base. - if (base > 36) - throw "BigUnsignedInABase(std::string, Base): The default string conversion routines use the symbol set 0-9, A-Z and therefore support only up to base 36. You tried a conversion with a base over 36; write your own string conversion routine."; - // Save the base. - // This pattern is seldom seen in C++, but the analogous ``this.'' is common in Java. - this->base = base; - - // `s.length()' is a `size_t', while `len' is a `NumberlikeArray::Index', - // also known as an `unsigned int'. Some compilers warn without this cast. - len = Index(s.length()); - allocate(len); - - Index digitNum, symbolNumInString; - for (digitNum = 0; digitNum < len; digitNum++) { - symbolNumInString = len - 1 - digitNum; - char theSymbol = s[symbolNumInString]; - if (theSymbol >= '0' && theSymbol <= '9') - blk[digitNum] = theSymbol - '0'; - else if (theSymbol >= 'A' && theSymbol <= 'Z') - blk[digitNum] = theSymbol - 'A' + 10; - else if (theSymbol >= 'a' && theSymbol <= 'z') - blk[digitNum] = theSymbol - 'a' + 10; - else - throw "BigUnsignedInABase(std::string, Base): Bad symbol in input. Only 0-9, A-Z, a-z are accepted."; - - if (blk[digitNum] >= base) - throw "BigUnsignedInABase::BigUnsignedInABase(const Digit *, Index, Base): A digit is too large for the specified base"; - } - zapLeadingZeros(); -} - -BigUnsignedInABase::operator std::string() const { - if (base > 36) - throw "BigUnsignedInABase ==> std::string: The default string conversion routines use the symbol set 0-9, A-Z and therefore support only up to base 36. You tried a conversion with a base over 36; write your own string conversion routine."; - if (len == 0) - return std::string("0"); - // Some compilers don't have push_back, so use a char * buffer instead. - char *s = new char[len + 1]; - s[len] = '\0'; - Index digitNum, symbolNumInString; - for (symbolNumInString = 0; symbolNumInString < len; symbolNumInString++) { - digitNum = len - 1 - symbolNumInString; - Digit theDigit = blk[digitNum]; - if (theDigit < 10) - s[symbolNumInString] = char('0' + theDigit); - else - s[symbolNumInString] = char('A' + theDigit - 10); - } - std::string s2(s); - delete [] s; - return s2; -} diff --git a/src/ecc/ecc.c b/src/ecc/ecc.c new file mode 100644 index 0000000..94ef88e --- /dev/null +++ b/src/ecc/ecc.c @@ -0,0 +1,389 @@ +/* + ecc.c - inplementations of ECC operations using keys + defined with sect233r1 / NIST B-233 + + This is NOT intended to be used in an actual cryptographic + scheme; as written, it is vulnerable to several attacks. + This might or might not change in the future. It is intended + to be used for doing operations on keys which are already known. + + Copyright © 2018, 2019 Jbop (https://github.com/jbop1626) + + This file is a part of ninty-233. + + ninty-233 is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ninty-233 is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#include +#include +#include + +#include "ecc.h" + +/* + sect233r1 - domain parameters over GF(2^m). + Defined in "Standards for Efficient Cryptography 2 (SEC 2)" v2.0, pp. 19-20 + Not all are currently used. +*/ +const element POLY_F = {0x0200, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000400, 0x00000000, 0x00000001}; +const element POLY_R = {0x0000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000400, 0x00000000, 0x00000001}; +const element A_COEFF = {0x0000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000001}; +const element B_COEFF = {0x0066, 0x647EDE6C, 0x332C7F8C, 0x0923BB58, 0x213B333B, 0x20E9CE42, 0x81FE115F, 0x7D8F90AD}; +const element G_X = {0x00FA, 0xC9DFCBAC, 0x8313BB21, 0x39F1BB75, 0x5FEF65BC, 0x391F8B36, 0xF8F8EB73, 0x71FD558B}; +const element G_Y = {0x0100, 0x6A08A419, 0x03350678, 0xE58528BE, 0xBF8A0BEF, 0xF867A7CA, 0x36716F7E, 0x01F81052}; +const element G_ORDER = {0x0100, 0x00000000, 0x00000000, 0x00000000, 0x0013E974, 0xE72F8A69, 0x22031D26, 0x03CFE0D7}; /* +const uint32_t COFACTOR = 0x02; */ + + +/* + Printing +*/ +void print_element(const element a) { + for (int i = 0; i < 8; ++i) { + printf("%08"PRIX32" ", a[i]); + } + printf("\n"); +} + +void print_point(const ec_point * a) { + printf("x: "); + print_element(a->x); + printf("y: "); + print_element(a->y); + printf("\n"); +} + +/* + Helper functions for working with elements in GF(2^m) +*/ +int gf2m_is_equal(const element a, const element b) { + for (int i = 0; i < 8; ++i) { + if (a[i] != b[i]) { + return 0; + } + } + return 1; +} + +void gf2m_set_zero(element a) { + for (int i = 0; i < 8; ++i) { + a[i] = 0; + } +} + +void gf2m_copy(const element src, element dst) { + for (int i = 0; i < 8; ++i) { + dst[i] = src[i]; + } +} + +int gf2m_get_bit(const element a, int index) { + if (index >= 233 || index < 0) return -1; + int word_index = ((index / 32) - 7) * -1; + int shift = index - (32 * (7 - word_index)); + return (a[word_index] >> shift) & 1; +} + +void gf2m_left_shift(element a, int shift) { + if (shift <= 0) { + a[0] &= 0x1FF; + return; + } + for (int i = 0; i < 7; ++i) { + a[i] <<= 1; + if (a[i + 1] >= 0x80000000) { + a[i] |= 1; + } + } + a[7] <<= 1; + gf2m_left_shift(a, shift - 1); +} + +int gf2m_is_one(const element a) { + if (a[7] != 1) { + return 0; + } + else { + for (int i = 0; i < 7; ++i) { + if (a[i] != 0) { + return 0; + } + } + } + return 1; +} + +int gf2m_degree(const element a) { + int degree = 0; + int i = 0; + while ((i < 7) && (a[i] == 0)) { + i++; + } + degree = (7 - i) * 32; + uint32_t most_significant_word = a[i]; + while (most_significant_word != 0) { + most_significant_word >>= 1; + degree += 1; + } + return degree - 1; +} + +void gf2m_swap(element a, element b) { + element temp; + gf2m_copy(a, temp); + gf2m_copy(b, a); + gf2m_copy(temp, b); +} + +/* + Arithmetic operations on elements in GF(2^m) +*/ +void gf2m_add(const element a, const element b, element result) { + for (int i = 0; i < 8; ++i) { + result[i] = a[i] ^ b[i]; + } +} + +void gf2m_inv(const element a, element result) { + element u, v, g_1, g_2, temp; + gf2m_copy(a, u); + gf2m_copy(POLY_F, v); + gf2m_set_zero(g_1); + g_1[7] |= 1; + gf2m_set_zero(g_2); + int j = gf2m_degree(u) - 233; + while (!gf2m_is_one(u)) { + if (j < 0) { + gf2m_swap(u, v); + gf2m_swap(g_1, g_2); + j = -j; + } + gf2m_copy(v, temp); + gf2m_left_shift(temp, j); + gf2m_add(u, temp, u); + gf2m_copy(g_2, temp); + gf2m_left_shift(temp, j); + gf2m_add(g_1, temp, g_1); + + u[0] &= 0x1FF; + g_1[0] &= 0x1FF; + + j = gf2m_degree(u) - gf2m_degree(v); + } + gf2m_copy(g_1, result); +} + +// basic implementation +void gf2m_mul(const element a, const element b, element result) { + element t1, t2, t3; + gf2m_copy(a, t1); + gf2m_copy(b, t2); + gf2m_set_zero(t3); + for (int i = 0; i < 233; ++i) { + if (gf2m_get_bit(t2, i)) { + gf2m_add(t3, t1, t3); + } + int carry = gf2m_get_bit(t1, 232); + gf2m_left_shift(t1, 1); + if (carry == 1) { + gf2m_add(POLY_R, t1, t1); + } + } + gf2m_copy(t3, result); +} + +void gf2m_div(const element a, const element b, element result) { + element temp; + gf2m_inv(b, temp); + gf2m_mul(a, temp, result); +} + +/* + Operations on points on the elliptic curve + y^2 + xy = x^3 + ax^2 + b over GF(2^m) +*/ +void ec_point_copy(const ec_point * src, ec_point * dst) { + gf2m_copy(src->x, dst->x); + gf2m_copy(src->y, dst->y); +} + +int ec_point_is_equal(const ec_point * a, const ec_point * b) { + return gf2m_is_equal(a->x, b->x) && gf2m_is_equal(a->y, b->y); +} + +void ec_point_neg(const ec_point * a, ec_point * result) { + gf2m_copy(a->x, result->x); + gf2m_add(a->x, a->y, result->y); +} + +void ec_point_double(const ec_point * a, ec_point * result) { + ec_point temp, zero; + gf2m_set_zero(zero.x); + gf2m_set_zero(zero.y); + + ec_point_neg(a, &temp); + if (ec_point_is_equal(a, &temp) || ec_point_is_equal(a, &zero)) { + ec_point_copy(&zero, result); + return; + } + + element lambda, x, y, t, t2; + // Compute lambda (a.x + (a.y / a.x)) + gf2m_div(a->y, a->x, t); + gf2m_add(a->x, t, lambda); + // Compute X (lambda^2 + lambda + A_COEFF) + gf2m_mul(lambda, lambda, t); + gf2m_add(t, lambda, t); + gf2m_add(t, A_COEFF, x); + // Compute Y (a.x^2 + (lambda * X) + X) + gf2m_mul(a->x, a->x, t); + gf2m_mul(lambda, x, t2); + gf2m_add(t, t2, t); + gf2m_add(t, x, y); + // Copy X,Y to output point result + gf2m_copy(x, result->x); + gf2m_copy(y, result->y); +} + +void ec_point_add(const ec_point * a, const ec_point * b, ec_point * result) { + if (!ec_point_is_equal(a, b)) { + ec_point temp, zero; + gf2m_set_zero(zero.x); + gf2m_set_zero(zero.y); + ec_point_neg(b, &temp); + if (ec_point_is_equal(a, &temp)) { + ec_point_copy(&zero, result); + return; + } + else if (ec_point_is_equal(a, &zero)) { + ec_point_copy(b, result); + return; + } + else if (ec_point_is_equal(b, &zero)) { + ec_point_copy(a, result); + return; + } + else { + element lambda, x, y, t, t2; + // Compute lambda ((b.y + a.y) / (b.x + a.x)) + gf2m_add(b->y, a->y, t); + gf2m_add(b->x, a->x, t2); + gf2m_div(t, t2, lambda); + // Compute X (lambda^2 + lambda + a.x + b.x + A_COEFF) + gf2m_mul(lambda, lambda, t); + gf2m_add(t, lambda, t2); + gf2m_add(t2, a->x, t); + gf2m_add(t, b->x, t2); + gf2m_add(t2, A_COEFF, x); + // Compute Y ((lambda * (a.x + X)) + X + a.y) + gf2m_add(a->x, x, t); + gf2m_mul(lambda, t, t2); + gf2m_add(t2, x, t); + gf2m_add(t, a->y, y); + // Copy X,Y to output point result + gf2m_copy(x, result->x); + gf2m_copy(y, result->y); + return; + } + } + else { + ec_point_double(a, result); + } +} + +void ec_point_mul(const element a, const ec_point * b, ec_point * result) { + element k; + ec_point P, Q; + + gf2m_copy(a, k); + ec_point_copy(b, &P); + gf2m_set_zero(Q.x); + gf2m_set_zero(Q.y); + for (int i = 0; i < 233; ++i) { + if (gf2m_get_bit(k, i)) { + ec_point_add(&Q, &P, &Q); + } + ec_point_double(&P, &P); + } + ec_point_copy(&Q, result); +} + +int ec_point_on_curve(const ec_point * P) { + // y^2 + xy = x^3 + ax^2 + b + element xx, yy, xy, lhs, rhs; + + // lhs = y^2 + xy + gf2m_mul(P->y, P->y, yy); + gf2m_mul(P->x, P->y, xy); + gf2m_add(yy, xy, lhs); + + // rhs = x^3 + ax^2 + b = (x^2)(x + a) + b + gf2m_mul(P->x, P->x, xx); + gf2m_add(P->x, A_COEFF, rhs); + gf2m_mul(xx, rhs, rhs); + gf2m_add(rhs, B_COEFF, rhs); + + return gf2m_is_equal(lhs, rhs); +} + +/* + I/O Helpers + Private keys are expected to be 32 bytes; Public keys + are expected to be 64 bytes and in uncompressed form. + + Wii keys will need to be padded - two 0 bytes at the + start of the private key, and two 0 bytes before each + coordinate in the public key. Keys on the iQue Player + are already stored with padding and need no changes. + + These functions are mainly intended for reading/writing + *keys* as byte arrays or octet streams, but they will + work fine for any input with the correct length. +*/ +// (32-byte) octet stream to GF(2^m) element +void os_to_elem(const uint8_t * src_os, element dst_elem) { + int j = 0; + for (int i = 0; i < 8; ++i) { + uint32_t result = src_os[j]; + result = (result << 8) | src_os[j + 1]; + result = (result << 8) | src_os[j + 2]; + result = (result << 8) | src_os[j + 3]; + dst_elem[i] = result; + j += 4; + } +} + +// (64-byte) octet stream to elliptic curve point +void os_to_point(const uint8_t * src_os, ec_point * dst_point) { + os_to_elem(src_os, dst_point->x); + os_to_elem(src_os + 32, dst_point->y); +} + +// GF(2^m) element to (32-byte) octet stream +void elem_to_os(const element src_elem, uint8_t * dst_os) { + int j = 0; + for (int i = 0; i < 8; ++i) { + dst_os[j] = (src_elem[i] & 0xFF000000) >> 24; + dst_os[j + 1] = (src_elem[i] & 0x00FF0000) >> 16; + dst_os[j + 2] = (src_elem[i] & 0x0000FF00) >> 8; + dst_os[j + 3] = src_elem[i] & 0x000000FF; + j += 4; + } +} + +// Elliptic curve point to (64-byte) octet stream +void point_to_os(const ec_point * src_point, uint8_t * dst_os) { + elem_to_os(src_point->x, dst_os); + elem_to_os(src_point->y, dst_os + 32); +} diff --git a/src/ecc/ecc.cpp b/src/ecc/ecc.cpp deleted file mode 100644 index 01e033a..0000000 --- a/src/ecc/ecc.cpp +++ /dev/null @@ -1,367 +0,0 @@ -/* - ecc.cpp - inplementations of ECC operations using keys - defined with sect233r1 / NIST B-233 - This is NOT intended to be used in an actual cryptographic - scheme; as written, it is vulnerable to several attacks. - This might or might not change in the future. It is intended - to be used for doing operations on keys which are already known. - - Copyright © 2018 Jbop (https://github.com/jbop1626); - Modification of a part of iQueCrypt - (https://github.com/jbop1626/iquecrypt) - - This file is a part of ninty-233. - - ninty-233 is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - ninty-233 is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . -*/ -#include -#include -#include -#include "ecc.hpp" - -/* - Printing -*/ -void print_element(const element a) { - for (int i = 0; i < 8; ++i) { - std::cout << std::setw(8) << std::setfill('0') << std::hex << a[i] << " "; - } - std::cout << std::endl << std::endl;; -} - -void print_point(const ec_point & a) { - std::cout << "x: "; - print_element(a.x); - std::cout << "y: "; - print_element(a.y); - std::cout << std::endl; -} - -/* - Helper functions for working with elements in GF(2^m) -*/ -bool gf2m_is_equal(const element a, const element b) { - for (int i = 0; i < 7; ++i) { - if (a[i] != b[i]) return false; - } - return true; -} - -void gf2m_set_zero(element a) { - for (int i = 0; i < 8; ++i) { - a[i] = 0; - } -} - -void gf2m_copy(const element src, element dst) { - std::memcpy(dst, src, 32); -} - -int gf2m_get_bit(const element a, int index) { - int word_index = ((index / 32) - 7) * -1; - int shift = index - (32 * (7 - word_index)); - return (a[word_index] >> shift) & 1; -} - -void gf2m_left_shift(element a, int shift) { - if (!shift) { - a[0] &= 0x1FF; - return; - } - for (int i = 0; i < 7; ++i) { - a[i] <<= 1; - if (a[i + 1] >= 0x80000000) a[i] |= 1; - } - a[7] <<= 1; - gf2m_left_shift(a, shift - 1); -} - -bool gf2m_is_one(const element a) { - if (a[7] != 1) return false; - else { - for (int i = 0; i < 7; ++i) { - if (a[i] != 0) return false; - } - } - return true; -} - -int gf2m_degree(const element a) { - int degree = 0; - int i = 0; - while (a[i] == 0) { - i++; - } - degree = (7 - i) * 32; - uint32_t MSW = a[i]; - while (MSW != 0) { - MSW >>= 1; - degree += 1; - } - return degree - 1; -} - -void gf2m_swap(element a, element b) { - element temp; - gf2m_copy(a, temp); - gf2m_copy(b, a); - gf2m_copy(temp, b); -} - -/* - Arithmetic operations on elements in GF(2^m) -*/ -void gf2m_add(const element a, const element b, element c) { - for (int i = 0; i < 8; ++i) { - c[i] = a[i] ^ b[i]; - } -} - -void gf2m_inv(const element a, element c) { - element u, v, g_1, g_2, temp; - gf2m_copy(a, u); - gf2m_copy(poly_f, v); - gf2m_set_zero(g_1); - g_1[7] |= 1; - gf2m_set_zero(g_2); - int j = gf2m_degree(u) - 233; - while (!gf2m_is_one(u)) { - if (j < 0) { - gf2m_swap(u, v); - gf2m_swap(g_1, g_2); - j = -j; - } - gf2m_copy(v, temp); - gf2m_left_shift(temp, j); - gf2m_add(u, temp, u); - gf2m_copy(g_2, temp); - gf2m_left_shift(temp, j); - gf2m_add(g_1, temp, g_1); - - u[0] &= 0x1FF; - g_1[0] &= 0x1FF; - - j = gf2m_degree(u) - gf2m_degree(v); - } - gf2m_copy(g_1, c); -} - -// basic implementation -void gf2m_mul(const element a, const element b, element c) { - element t1, t2, t3; - gf2m_copy(a, t1); - gf2m_copy(b, t2); - gf2m_set_zero(t3); - for (int i = 0; i < 233; ++i) { - if (gf2m_get_bit(t2, i)) { - gf2m_add(t3, t1, t3); - } - int carry = gf2m_get_bit(t1, 232); - gf2m_left_shift(t1, 1); - if (carry == 1) { - gf2m_add(poly_r, t1, t1); - } - } - gf2m_copy(t3, c); -} - -void gf2m_div(const element a, const element b, element c) { - element temp; - gf2m_inv(b, temp); - gf2m_mul(a, temp, c); -} -// void gf2m_reduce(element c) -// void gf2m_square(const element a, element c) - -/* - Operations on points on the elliptic curve - y^2 + xy = x^3 + ax^2 + b over GF(2^m) -*/ -void ec_point_copy(const ec_point & src, ec_point & dst) { - gf2m_copy(src.x, dst.x); - gf2m_copy(src.y, dst.y); -} - -bool ec_point_is_equal(const ec_point & a, const ec_point & c) { - return gf2m_is_equal(a.x, c.x) && gf2m_is_equal(a.y, c.y); -} - -void ec_point_neg(const ec_point & a, ec_point & c) { - element temp; - gf2m_copy(a.x, c.x); - gf2m_add(a.x, a.y, temp); - gf2m_copy(temp, c.y); -} - -void ec_point_double(const ec_point & a, ec_point & c) { - ec_point temp; - ec_point zero; - gf2m_set_zero(zero.x); - gf2m_set_zero(zero.y); - - ec_point_neg(a, temp); - if (ec_point_is_equal(a, temp) || ec_point_is_equal(a, zero)) { - ec_point_copy(zero, c); - return; - } - - element lambda, x, y, t, t2; - // Compute lambda (a.x + (a.y / a.x)) - gf2m_div(a.y, a.x, t); - gf2m_add(a.x, t, lambda); - // Compute X (lambda^2 + lambda + a_coeff) - gf2m_mul(lambda, lambda, t); - gf2m_add(t, lambda, t); - gf2m_add(t, a_coeff, x); - // Compute Y (a.x^2 + (lambda * X) + X) - gf2m_mul(a.x, a.x, t); - gf2m_mul(lambda, x, t2); - gf2m_add(t, t2, t); - gf2m_add(t, x, y); - // Copy X,Y to output point c - gf2m_copy(x, c.x); - gf2m_copy(y, c.y); -} - -void ec_point_add(const ec_point & a, const ec_point & b, ec_point & c) { - if (!ec_point_is_equal(a, b)) { - ec_point temp; - ec_point zero; - gf2m_set_zero(zero.x); - gf2m_set_zero(zero.y); - ec_point_neg(b, temp); - if (ec_point_is_equal(a, temp)) { - ec_point_copy(zero, c); - return; - } - else if (ec_point_is_equal(a, zero)) { - ec_point_copy(b, c); - return; - } - else if (ec_point_is_equal(b, zero)) { - ec_point_copy(a, c); - return; - } - else { - element lambda, x, y, t, t2; - // Compute lambda ((b.y + a.y) / (b.x + a.x)) - gf2m_add(b.y, a.y, t); - gf2m_add(b.x, a.x, t2); - gf2m_div(t, t2, lambda); - // Compute X (lambda^2 + lambda + a.x + b.x + a_coeff) - gf2m_mul(lambda, lambda, t); - gf2m_add(t, lambda, t2); - gf2m_add(t2, a.x, t); - gf2m_add(t, b.x, t2); - gf2m_add(t2, a_coeff, x); - // Compute Y ((lambda * (a.x + X)) + X + a.y) - gf2m_add(a.x, x, t); - gf2m_mul(lambda, t, t2); - gf2m_add(t2, x, t); - gf2m_add(t, a.y, y); - // Copy X,Y to output point c - gf2m_copy(x, c.x); - gf2m_copy(y, c.y); - return; - } - } - else { - ec_point_double(a, c); - } -} - -void ec_point_mul(const element a, const ec_point & b, ec_point & c) { - element k; - ec_point P; - ec_point Q; - - gf2m_copy(a, k); - ec_point_copy(b, P); - gf2m_set_zero(Q.x); - gf2m_set_zero(Q.y); - for (int i = 0; i < 233; ++i) { - if (gf2m_get_bit(k, i)) { - ec_point_add(Q, P, Q); - } - ec_point_double(P, P); - } - ec_point_copy(Q, c); -} - -bool ec_point_on_curve(const ec_point & P) { - // y^2 + xy = x^3 + ax^2 + b - element xx, yy, xy, lhs, rhs; - - gf2m_mul(P.y, P.y, yy); - gf2m_mul(P.x, P.y, xy); - gf2m_add(yy, xy, lhs); - - gf2m_mul(P.x, P.x, xx); - gf2m_add(P.x, a_coeff, rhs); - gf2m_mul(xx, rhs, rhs); - gf2m_add(rhs, b_coeff, rhs); - - return gf2m_is_equal(lhs, rhs); -} - -/* - I/O Helpers - Private keys are expected to be 32 bytes; Public keys - are expected to be 64 bytes and in uncompressed form. - - Wii keys will need to be padded - two 0 bytes at the - start of the private key, and two 0 bytes before each - coordinate in the public key. - - These functions are mainly intended for reading/writing - *keys* as byte arrays or octet streams, but they will - work fine for any input with the correct length. -*/ -// (32-byte) octet stream to GF(2^m) element -void os_to_elem(const uint8_t * os, element elem) { - int j = 0; - for (int i = 0; i < 8; ++i) { - uint32_t temp = 0; - temp |= (os[j] << 24); - temp |= (os[j + 1] << 16); - temp |= (os[j + 2] << 8); - temp |= os[j + 3]; - elem[i] = temp; - j += 4; - } -} - -// (64-byte) octet stream to elliptic curve point -void os_to_point(const uint8_t * os, ec_point & point) { - os_to_elem(os, point.x); - os_to_elem(os + 32, point.y); -} - -// GF(2^m) element to (32-byte) octet stream -void elem_to_os(const element src, uint8_t * output_os) { - int j = 0; - for (int i = 0; i < 7; ++i) { - output_os[j] = ((src[i] & 0xFF000000) >> 24); - output_os[j + 1] = ((src[i] & 0x00FF0000) >> 16); - output_os[j + 2] = ((src[i] & 0x0000FF00) >> 8); - output_os[j + 3] = src[i] & 0x000000FF; - j += 4; - } -} - -// Elliptic curve point to (64-byte) octet stream -void point_to_os(const ec_point & src, uint8_t * output_os) { - elem_to_os(src.x, output_os); - elem_to_os(src.y, output_os + 32); -} diff --git a/src/ecc/ecc.h b/src/ecc/ecc.h new file mode 100644 index 0000000..92aff54 --- /dev/null +++ b/src/ecc/ecc.h @@ -0,0 +1,103 @@ +/* + ecc.h - definitions required for ECC operations using keys + defined with sect233r1 / NIST B-233 + + Copyright © 2018, 2019 Jbop (https://github.com/jbop1626) + + This file is a part of ninty-233. + + ninty-233 is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ninty-233 is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ +#ifndef NINTY_233_ECC_H +#define NINTY_233_ECC_H + +#include + +#if defined (__cplusplus) +extern "C" { +#endif + +typedef uint32_t element[8]; + +typedef struct { + element x; + element y; +} ec_point; + +/* + sect233r1 - domain parameters over GF(2^m). + Defined in "Standards for Efficient Cryptography 2 (SEC 2)" v2.0, pp. 19-20 + Not all are currently used. + (Actual definitions are in ecc.c) +*/ +extern const element POLY_F; +extern const element POLY_R; +extern const element A_COEFF; +extern const element B_COEFF; +extern const element G_X; +extern const element G_Y; +extern const element G_ORDER; /* +extern const uint32_t COFACTOR; */ + +/* + Printing +*/ +void print_element(const element a); +void print_point(const ec_point * a); + +/* + Helper functions for working with elements in GF(2^m) +*/ +int gf2m_is_equal(const element a, const element b); +void gf2m_set_zero(element a); +void gf2m_copy(const element src, element dst); +int gf2m_get_bit(const element a, int index); +void gf2m_left_shift(element a, int shift); +int gf2m_is_one(const element a); +int gf2m_degree(const element a); +void gf2m_swap(element a, element b); + +/* + Arithmetic operations on elements in GF(2^m) +*/ +void gf2m_add(const element a, const element b, element result); +void gf2m_inv(const element a, element result); +void gf2m_mul(const element a, const element b, element result); +void gf2m_div(const element a, const element b, element result); + +/* + Operations on points on the elliptic curve + y^2 + xy = x^3 + ax^2 + b over GF(2^m) +*/ +void ec_point_copy(const ec_point * src, ec_point * dst); +int ec_point_is_equal(const ec_point * a, const ec_point * b); +void ec_point_neg(const ec_point * a, ec_point * result); +void ec_point_double(const ec_point * a, ec_point * result); +void ec_point_add(const ec_point * a, const ec_point * b, ec_point * result); +void ec_point_mul(const element a, const ec_point * b, ec_point * result); +int ec_point_on_curve(const ec_point * a); + +/* + I/O Helpers +*/ +void os_to_elem(const uint8_t * src_os, element dst_elem); +void os_to_point(const uint8_t * src_os, ec_point * dst_point); +void elem_to_os(const element src_elem, uint8_t * dst_os); +void point_to_os(const ec_point * src_point, uint8_t * dst_os); + +#if defined (__cplusplus) +} +#endif + +#endif diff --git a/src/ecc/ecc.hpp b/src/ecc/ecc.hpp deleted file mode 100644 index 86ecfb6..0000000 --- a/src/ecc/ecc.hpp +++ /dev/null @@ -1,95 +0,0 @@ -/* - ecc.hpp - definitions required for ECC operations using keys - defined with sect233r1 / NIST B-233 - - Copyright © 2018 Jbop (https://github.com/jbop1626); - Modification of a part of iQueCrypt - (https://github.com/jbop1626/iquecrypt) - - This file is a part of ninty-233. - - ninty-233 is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - ninty-233 is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . -*/ -#ifndef NINTY_233_ECC_HPP -#define NINTY_233_ECC_HPP - -typedef uint32_t element[8]; - -typedef struct { - element x; - element y; -} ec_point; - -/* - sect233r1 - domain parameters over GF(2^m). Defined in SEC 2 v2.0, pp. 19-20 - Not all are currently used. -*/ -const element poly_f = {0x0200, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000400, 0x00000000, 0x00000001}; -const element poly_r = {0x0000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000400, 0x00000000, 0x00000001}; -const element a_coeff = {0x0000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000001}; -const element b_coeff = {0x0066, 0x647EDE6C, 0x332C7F8C, 0x0923BB58, 0x213B333B, 0x20E9CE42, 0x81FE115F, 0x7D8F90AD}; -const element G_x = {0x00FA, 0xC9DFCBAC, 0x8313BB21, 0x39F1BB75, 0x5FEF65BC, 0x391F8B36, 0xF8F8EB73, 0x71FD558B}; -const element G_y = {0x0100, 0x6A08A419, 0x03350678, 0xE58528BE, 0xBF8A0BEF, 0xF867A7CA, 0x36716F7E, 0x01F81052}; -const element G_order = {0x0100, 0x00000000, 0x00000000, 0x00000000, 0x0013E974, 0xE72F8A69, 0x22031D26, 0x03CFE0D7}; /* -const uint32_t cofactor = 0x02; */ - -/* - Printing -*/ -void print_element(const element a); -void print_point(const ec_point & a); - -/* - Helper functions for working with elements in GF(2^m) -*/ -bool gf2m_is_equal(const element a, const element b); -void gf2m_set_zero(element a); -void gf2m_copy(const element src, element dst); -int gf2m_get_bit(const element a, int index); -void gf2m_left_shift(element a, int shift); -bool gf2m_is_one(const element a); -int gf2m_degree(const element a); -void gf2m_swap(element a, element b); - -/* - Arithmetic operations on elements in GF(2^m) -*/ -void gf2m_add(const element a, const element b, element c); -void gf2m_inv(const element a, element c); -void gf2m_mul(const element a, const element b, element c); -void gf2m_div(const element a, const element b, element c); -// void gf2m_reduce(element c); -// void gf2m_square(const element a, element c); - -/* - Operations on points on the elliptic curve - y^2 + xy = x^3 + ax^2 + b over GF(2^m) -*/ -void ec_point_copy(const ec_point & src, ec_point & dst); -bool ec_point_is_equal(const ec_point & a, const ec_point & c); -void ec_point_neg(const ec_point & a, ec_point & c); -void ec_point_double(const ec_point & a, ec_point & c); -void ec_point_add(const ec_point & a, const ec_point & b, ec_point & c); -void ec_point_mul(const element a, const ec_point & b, ec_point & c); -bool ec_point_on_curve(const ec_point & a); - -/* - I/O Helpers -*/ -void os_to_elem(const uint8_t * os, element elem); -void os_to_point(const uint8_t * os, ec_point & point); -void elem_to_os(const element src, uint8_t * output_os); -void point_to_os(const ec_point & src, uint8_t * output_os); - -#endif diff --git a/src/example.c b/src/example.c new file mode 100644 index 0000000..05aa37e --- /dev/null +++ b/src/example.c @@ -0,0 +1,113 @@ +/* + example.c + Simple example of how to use ninty-233 + + Copyright © 2019 Jbop (https://github.com/jbop1626) + + This file is a part of ninty-233. + + ninty-233 is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ninty-233 is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#include +#include + +#include "ninty-233.h" + +int main(int argc, char * argv[]) { + uint8_t private_key_A[32] = { 0x00, 0x00, 0x01, 0x5F, 0xC2, 0xC6, 0x53, 0xCD, + 0xBA, 0xDC, 0xD8, 0x24, 0x23, 0xFE, 0xA2, 0xE8, + 0xCC, 0x75, 0x7C, 0xC3, 0x5A, 0x27, 0x46, 0x73, + 0xF5, 0x70, 0xAE, 0x7A, 0xD2, 0xBB, 0x59, 0xB0 }; + + uint8_t public_key_A[64] = { 0x00, 0x00, 0x00, 0x27, 0x64, 0xB9, 0x37, 0x21, + 0xDA, 0xA4, 0x33, 0xC1, 0xC4, 0x45, 0x0B, 0xFA, + 0x8E, 0x4D, 0x36, 0x9C, 0x41, 0x16, 0xD6, 0xED, + 0xBE, 0x03, 0x0D, 0x9F, 0x12, 0x3B, 0xB1, 0x18, + 0x00, 0x00, 0x01, 0x87, 0xBE, 0xC4, 0xF5, 0x1A, + 0x9F, 0x4B, 0x0B, 0x3F, 0xD3, 0x08, 0x01, 0xBC, + 0x6F, 0x11, 0xAC, 0x4C, 0x26, 0xDB, 0x3B, 0x6C, + 0x70, 0x5D, 0xC1, 0x7E, 0x9C, 0x00, 0x84, 0x27 }; + + uint8_t private_key_B[32] = { 0x00, 0x00, 0x00, 0xE4, 0xC0, 0x89, 0x52, 0x2C, + 0xA3, 0x38, 0x6C, 0xC8, 0xF6, 0x29, 0x80, 0x1E, + 0x0F, 0xE9, 0xD0, 0x92, 0xA5, 0x61, 0x27, 0x48, + 0xC1, 0xE9, 0x51, 0x1D, 0x82, 0xDB, 0x93, 0xE0 }; + + uint8_t public_key_B[64] = { 0x00, 0x00, 0x00, 0xDD, 0x56, 0x98, 0x2B, 0xED, + 0x1F, 0x91, 0x0C, 0x20, 0x2D, 0x91, 0x38, 0xE8, + 0x6B, 0xFC, 0x60, 0x77, 0x3F, 0x38, 0xF5, 0x4A, + 0x08, 0xEC, 0xB3, 0xD6, 0xEB, 0x40, 0x10, 0xD1, + 0x00, 0x00, 0x00, 0xBB, 0xFD, 0x3C, 0xA6, 0x76, + 0x6F, 0xB1, 0x19, 0xCE, 0xC4, 0xEB, 0x65, 0x74, + 0x8D, 0x54, 0x9B, 0xD6, 0x94, 0x0F, 0x70, 0x44, + 0x00, 0x0F, 0x8E, 0xA1, 0xD5, 0x1B, 0x47, 0x1A}; + + uint8_t data[48] = { 0x6B, 0xFC, 0x60, 0x77, 0x3F, 0x38, 0xF5, 0x4A, + 0x08, 0xEC, 0xB3, 0xD6, 0xEB, 0x40, 0x10, 0xD1, + 0x00, 0x00, 0x00, 0xBB, 0xFD, 0x3C, 0xA6, 0x76, + 0xA3, 0x38, 0x6C, 0xC8, 0xF6, 0x29, 0x80, 0x1E, + 0x8E, 0x4D, 0x36, 0x9C, 0x41, 0x16, 0xD6, 0xED, + 0xBE, 0x03, 0x0D, 0x9F, 0x12, 0x3B, 0xB1, 0x18 }; + + // Convert bytes to GF(2^m) elements and elliptic curve points + element priv_key_A, priv_key_B; + ec_point pub_key_A, pub_key_B; + os_to_elem(private_key_A, priv_key_A); + os_to_elem(private_key_B, priv_key_B); + os_to_point(public_key_A, &pub_key_A); + os_to_point(public_key_B, &pub_key_B); + + /* ECDH */ + printf("ECDH with private key A and public key B:\n"); + ec_point shared_secret1; + ecdh(priv_key_A, &pub_key_B, &shared_secret1); + print_point(&shared_secret1); + + printf("ECDH with private key B and public key A:\n"); + ec_point shared_secret2; + ecdh(priv_key_B, &pub_key_A, &shared_secret2); + print_point(&shared_secret2); + + if (ec_point_is_equal(&shared_secret1, &shared_secret2)) { + printf("Success! Shared secret outputs are equal.\n"); + } + else { + printf("Failure! Shared secret outputs are not equal.\n"); + } + + + /* HASHING */ + printf("\n\nSign data with private key A:\n"); + mpz_t hash; + mpz_init(hash); + + // If we wanted to hash in the way the iQue Player does it (by adding + // certain magic data to the SHA1 state), we would pass in the aptly-named + // IQUE_HASH flag instead of NOT_IQUE_HASH as the 3rd argument. + sha1(data, 48, NOT_IQUE_HASH, hash); + + /* SIGNING */ + element r, s; + ecdsa_sign(hash, priv_key_A, r, s); + printf("Complete!\n"); + + /* SIGNATURE VERIFICATION */ + int result1 = ecdsa_verify(hash, &pub_key_A, r, s); + printf("Verify signature with public key A: %d\n", result1); + + int result2 = ecdsa_verify(hash, &pub_key_B, r, s); + printf("Verify signature with public key B: %d\n", result2); + mpz_clear(hash); +} diff --git a/src/mini-gmp/mini-gmp.c b/src/mini-gmp/mini-gmp.c new file mode 100644 index 0000000..38f4f19 --- /dev/null +++ b/src/mini-gmp/mini-gmp.c @@ -0,0 +1,4412 @@ +/* mini-gmp, a minimalistic implementation of a GNU GMP subset. + + Contributed to the GNU project by Niels Möller + +Copyright 1991-1997, 1999-2016 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +/* NOTE: All functions in this file which are not declared in + mini-gmp.h are internal, and are not intended to be compatible + neither with GMP nor with future versions of mini-gmp. */ + +/* Much of the material copied from GMP files, including: gmp-impl.h, + longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c, + mpn/generic/lshift.c, mpn/generic/mul_1.c, + mpn/generic/mul_basecase.c, mpn/generic/rshift.c, + mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c, + mpn/generic/submul_1.c. */ + +#include +#include +#include +#include +#include +#include + +#include "mini-gmp.h" + + +/* Macros */ +#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT) + +#define GMP_LIMB_MAX (~ (mp_limb_t) 0) +#define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1)) + +#define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2)) +#define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1) + +#define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT) +#define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1)) + +#define GMP_ABS(x) ((x) >= 0 ? (x) : -(x)) +#define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1)) + +#define GMP_MIN(a, b) ((a) < (b) ? (a) : (b)) +#define GMP_MAX(a, b) ((a) > (b) ? (a) : (b)) + +#define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b))) + +#define gmp_assert_nocarry(x) do { \ + mp_limb_t cy_ = (x); \ + assert (cy_ == 0); \ + } while (0) + +#define gmp_clz(count, x) do { \ + mp_limb_t clz_x_ = (x); \ + unsigned clz_c_; \ + for (clz_c_ = 0; \ + (clz_x_ & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \ + clz_c_ += 8) \ + clz_x_ <<= 8; \ + for (; (clz_x_ & GMP_LIMB_HIGHBIT) == 0; clz_c_++) \ + clz_x_ <<= 1; \ + (count) = clz_c_; \ + } while (0) + +#define gmp_ctz(count, x) do { \ + mp_limb_t ctz_x_ = (x); \ + unsigned ctz_c_ = 0; \ + gmp_clz (ctz_c_, ctz_x_ & - ctz_x_); \ + (count) = GMP_LIMB_BITS - 1 - ctz_c_; \ + } while (0) + +#define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \ + do { \ + mp_limb_t x_; \ + x_ = (al) + (bl); \ + (sh) = (ah) + (bh) + (x_ < (al)); \ + (sl) = x_; \ + } while (0) + +#define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \ + do { \ + mp_limb_t x_; \ + x_ = (al) - (bl); \ + (sh) = (ah) - (bh) - ((al) < (bl)); \ + (sl) = x_; \ + } while (0) + +#define gmp_umul_ppmm(w1, w0, u, v) \ + do { \ + mp_limb_t x0_, x1_, x2_, x3_; \ + unsigned ul_, vl_, uh_, vh_; \ + mp_limb_t u_ = (u), v_ = (v); \ + \ + ul_ = u_ & GMP_LLIMB_MASK; \ + uh_ = u_ >> (GMP_LIMB_BITS / 2); \ + vl_ = v_ & GMP_LLIMB_MASK; \ + vh_ = v_ >> (GMP_LIMB_BITS / 2); \ + \ + x0_ = (mp_limb_t) ul_ * vl_; \ + x1_ = (mp_limb_t) ul_ * vh_; \ + x2_ = (mp_limb_t) uh_ * vl_; \ + x3_ = (mp_limb_t) uh_ * vh_; \ + \ + x1_ += x0_ >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \ + x1_ += x2_; /* but this indeed can */ \ + if (x1_ < x2_) /* did we get it? */ \ + x3_ += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \ + \ + (w1) = x3_ + (x1_ >> (GMP_LIMB_BITS / 2)); \ + (w0) = (x1_ << (GMP_LIMB_BITS / 2)) + (x0_ & GMP_LLIMB_MASK); \ + } while (0) + +#define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \ + do { \ + mp_limb_t qh_, ql_, r_, mask_; \ + gmp_umul_ppmm (qh_, ql_, (nh), (di)); \ + gmp_add_ssaaaa (qh_, ql_, qh_, ql_, (nh) + 1, (nl)); \ + r_ = (nl) - qh_ * (d); \ + mask_ = -(mp_limb_t) (r_ > ql_); /* both > and >= are OK */ \ + qh_ += mask_; \ + r_ += mask_ & (d); \ + if (r_ >= (d)) \ + { \ + r_ -= (d); \ + qh_++; \ + } \ + \ + (r) = r_; \ + (q) = qh_; \ + } while (0) + +#define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \ + do { \ + mp_limb_t q0_, t1_, t0_, mask_; \ + gmp_umul_ppmm ((q), q0_, (n2), (dinv)); \ + gmp_add_ssaaaa ((q), q0_, (q), q0_, (n2), (n1)); \ + \ + /* Compute the two most significant limbs of n - q'd */ \ + (r1) = (n1) - (d1) * (q); \ + gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0)); \ + gmp_umul_ppmm (t1_, t0_, (d0), (q)); \ + gmp_sub_ddmmss ((r1), (r0), (r1), (r0), t1_, t0_); \ + (q)++; \ + \ + /* Conditionally adjust q and the remainders */ \ + mask_ = - (mp_limb_t) ((r1) >= q0_); \ + (q) += mask_; \ + gmp_add_ssaaaa ((r1), (r0), (r1), (r0), mask_ & (d1), mask_ & (d0)); \ + if ((r1) >= (d1)) \ + { \ + if ((r1) > (d1) || (r0) >= (d0)) \ + { \ + (q)++; \ + gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \ + } \ + } \ + } while (0) + +/* Swap macros. */ +#define MP_LIMB_T_SWAP(x, y) \ + do { \ + mp_limb_t mp_limb_t_swap_tmp_ = (x); \ + (x) = (y); \ + (y) = mp_limb_t_swap_tmp_; \ + } while (0) +#define MP_SIZE_T_SWAP(x, y) \ + do { \ + mp_size_t mp_size_t_swap_tmp_ = (x); \ + (x) = (y); \ + (y) = mp_size_t_swap_tmp_; \ + } while (0) +#define MP_BITCNT_T_SWAP(x,y) \ + do { \ + mp_bitcnt_t mp_bitcnt_t_swap_tmp_ = (x); \ + (x) = (y); \ + (y) = mp_bitcnt_t_swap_tmp_; \ + } while (0) +#define MP_PTR_SWAP(x, y) \ + do { \ + mp_ptr mp_ptr_swap_tmp_ = (x); \ + (x) = (y); \ + (y) = mp_ptr_swap_tmp_; \ + } while (0) +#define MP_SRCPTR_SWAP(x, y) \ + do { \ + mp_srcptr mp_srcptr_swap_tmp_ = (x); \ + (x) = (y); \ + (y) = mp_srcptr_swap_tmp_; \ + } while (0) + +#define MPN_PTR_SWAP(xp,xs, yp,ys) \ + do { \ + MP_PTR_SWAP (xp, yp); \ + MP_SIZE_T_SWAP (xs, ys); \ + } while(0) +#define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \ + do { \ + MP_SRCPTR_SWAP (xp, yp); \ + MP_SIZE_T_SWAP (xs, ys); \ + } while(0) + +#define MPZ_PTR_SWAP(x, y) \ + do { \ + mpz_ptr mpz_ptr_swap_tmp_ = (x); \ + (x) = (y); \ + (y) = mpz_ptr_swap_tmp_; \ + } while (0) +#define MPZ_SRCPTR_SWAP(x, y) \ + do { \ + mpz_srcptr mpz_srcptr_swap_tmp_ = (x); \ + (x) = (y); \ + (y) = mpz_srcptr_swap_tmp_; \ + } while (0) + +const int mp_bits_per_limb = GMP_LIMB_BITS; + + +/* Memory allocation and other helper functions. */ +static void +gmp_die (const char *msg) +{ + fprintf (stderr, "%s\n", msg); + abort(); +} + +static void * +gmp_default_alloc (size_t size) +{ + void *p; + + assert (size > 0); + + p = malloc (size); + if (!p) + gmp_die("gmp_default_alloc: Virtual memory exhausted."); + + return p; +} + +static void * +gmp_default_realloc (void *old, size_t old_size, size_t new_size) +{ + void * p; + + p = realloc (old, new_size); + + if (!p) + gmp_die("gmp_default_realloc: Virtual memory exhausted."); + + return p; +} + +static void +gmp_default_free (void *p, size_t size) +{ + free (p); +} + +static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc; +static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc; +static void (*gmp_free_func) (void *, size_t) = gmp_default_free; + +void +mp_get_memory_functions (void *(**alloc_func) (size_t), + void *(**realloc_func) (void *, size_t, size_t), + void (**free_func) (void *, size_t)) +{ + if (alloc_func) + *alloc_func = gmp_allocate_func; + + if (realloc_func) + *realloc_func = gmp_reallocate_func; + + if (free_func) + *free_func = gmp_free_func; +} + +void +mp_set_memory_functions (void *(*alloc_func) (size_t), + void *(*realloc_func) (void *, size_t, size_t), + void (*free_func) (void *, size_t)) +{ + if (!alloc_func) + alloc_func = gmp_default_alloc; + if (!realloc_func) + realloc_func = gmp_default_realloc; + if (!free_func) + free_func = gmp_default_free; + + gmp_allocate_func = alloc_func; + gmp_reallocate_func = realloc_func; + gmp_free_func = free_func; +} + +#define gmp_xalloc(size) ((*gmp_allocate_func)((size))) +#define gmp_free(p) ((*gmp_free_func) ((p), 0)) + +static mp_ptr +gmp_xalloc_limbs (mp_size_t size) +{ + return (mp_ptr) gmp_xalloc (size * sizeof (mp_limb_t)); +} + +static mp_ptr +gmp_xrealloc_limbs (mp_ptr old, mp_size_t size) +{ + assert (size > 0); + return (mp_ptr) (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t)); +} + + +/* MPN interface */ + +void +mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n) +{ + mp_size_t i; + for (i = 0; i < n; i++) + d[i] = s[i]; +} + +void +mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n) +{ + while (--n >= 0) + d[n] = s[n]; +} + +int +mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + while (--n >= 0) + { + if (ap[n] != bp[n]) + return ap[n] > bp[n] ? 1 : -1; + } + return 0; +} + +static int +mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) +{ + if (an != bn) + return an < bn ? -1 : 1; + else + return mpn_cmp (ap, bp, an); +} + +static mp_size_t +mpn_normalized_size (mp_srcptr xp, mp_size_t n) +{ + while (n > 0 && xp[n-1] == 0) + --n; + return n; +} + +int +mpn_zero_p(mp_srcptr rp, mp_size_t n) +{ + return mpn_normalized_size (rp, n) == 0; +} + +void +mpn_zero (mp_ptr rp, mp_size_t n) +{ + while (--n >= 0) + rp[n] = 0; +} + +mp_limb_t +mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b) +{ + mp_size_t i; + + assert (n > 0); + i = 0; + do + { + mp_limb_t r = ap[i] + b; + /* Carry out */ + b = (r < b); + rp[i] = r; + } + while (++i < n); + + return b; +} + +mp_limb_t +mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + mp_size_t i; + mp_limb_t cy; + + for (i = 0, cy = 0; i < n; i++) + { + mp_limb_t a, b, r; + a = ap[i]; b = bp[i]; + r = a + cy; + cy = (r < cy); + r += b; + cy += (r < b); + rp[i] = r; + } + return cy; +} + +mp_limb_t +mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) +{ + mp_limb_t cy; + + assert (an >= bn); + + cy = mpn_add_n (rp, ap, bp, bn); + if (an > bn) + cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy); + return cy; +} + +mp_limb_t +mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b) +{ + mp_size_t i; + + assert (n > 0); + + i = 0; + do + { + mp_limb_t a = ap[i]; + /* Carry out */ + mp_limb_t cy = a < b; + rp[i] = a - b; + b = cy; + } + while (++i < n); + + return b; +} + +mp_limb_t +mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + mp_size_t i; + mp_limb_t cy; + + for (i = 0, cy = 0; i < n; i++) + { + mp_limb_t a, b; + a = ap[i]; b = bp[i]; + b += cy; + cy = (b < cy); + cy += (a < b); + rp[i] = a - b; + } + return cy; +} + +mp_limb_t +mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn) +{ + mp_limb_t cy; + + assert (an >= bn); + + cy = mpn_sub_n (rp, ap, bp, bn); + if (an > bn) + cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy); + return cy; +} + +mp_limb_t +mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) +{ + mp_limb_t ul, cl, hpl, lpl; + + assert (n >= 1); + + cl = 0; + do + { + ul = *up++; + gmp_umul_ppmm (hpl, lpl, ul, vl); + + lpl += cl; + cl = (lpl < cl) + hpl; + + *rp++ = lpl; + } + while (--n != 0); + + return cl; +} + +mp_limb_t +mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) +{ + mp_limb_t ul, cl, hpl, lpl, rl; + + assert (n >= 1); + + cl = 0; + do + { + ul = *up++; + gmp_umul_ppmm (hpl, lpl, ul, vl); + + lpl += cl; + cl = (lpl < cl) + hpl; + + rl = *rp; + lpl = rl + lpl; + cl += lpl < rl; + *rp++ = lpl; + } + while (--n != 0); + + return cl; +} + +mp_limb_t +mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl) +{ + mp_limb_t ul, cl, hpl, lpl, rl; + + assert (n >= 1); + + cl = 0; + do + { + ul = *up++; + gmp_umul_ppmm (hpl, lpl, ul, vl); + + lpl += cl; + cl = (lpl < cl) + hpl; + + rl = *rp; + lpl = rl - lpl; + cl += lpl > rl; + *rp++ = lpl; + } + while (--n != 0); + + return cl; +} + +mp_limb_t +mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn) +{ + assert (un >= vn); + assert (vn >= 1); + + /* We first multiply by the low order limb. This result can be + stored, not added, to rp. We also avoid a loop for zeroing this + way. */ + + rp[un] = mpn_mul_1 (rp, up, un, vp[0]); + + /* Now accumulate the product of up[] and the next higher limb from + vp[]. */ + + while (--vn >= 1) + { + rp += 1, vp += 1; + rp[un] = mpn_addmul_1 (rp, up, un, vp[0]); + } + return rp[un]; +} + +void +mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + mpn_mul (rp, ap, n, bp, n); +} + +void +mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n) +{ + mpn_mul (rp, ap, n, ap, n); +} + +mp_limb_t +mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) +{ + mp_limb_t high_limb, low_limb; + unsigned int tnc; + mp_limb_t retval; + + assert (n >= 1); + assert (cnt >= 1); + assert (cnt < GMP_LIMB_BITS); + + up += n; + rp += n; + + tnc = GMP_LIMB_BITS - cnt; + low_limb = *--up; + retval = low_limb >> tnc; + high_limb = (low_limb << cnt); + + while (--n != 0) + { + low_limb = *--up; + *--rp = high_limb | (low_limb >> tnc); + high_limb = (low_limb << cnt); + } + *--rp = high_limb; + + return retval; +} + +mp_limb_t +mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt) +{ + mp_limb_t high_limb, low_limb; + unsigned int tnc; + mp_limb_t retval; + + assert (n >= 1); + assert (cnt >= 1); + assert (cnt < GMP_LIMB_BITS); + + tnc = GMP_LIMB_BITS - cnt; + high_limb = *up++; + retval = (high_limb << tnc); + low_limb = high_limb >> cnt; + + while (--n != 0) + { + high_limb = *up++; + *rp++ = low_limb | (high_limb << tnc); + low_limb = high_limb >> cnt; + } + *rp = low_limb; + + return retval; +} + +static mp_bitcnt_t +mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un, + mp_limb_t ux) +{ + unsigned cnt; + + assert (ux == 0 || ux == GMP_LIMB_MAX); + assert (0 <= i && i <= un ); + + while (limb == 0) + { + i++; + if (i == un) + return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS); + limb = ux ^ up[i]; + } + gmp_ctz (cnt, limb); + return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt; +} + +mp_bitcnt_t +mpn_scan1 (mp_srcptr ptr, mp_bitcnt_t bit) +{ + mp_size_t i; + i = bit / GMP_LIMB_BITS; + + return mpn_common_scan ( ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)), + i, ptr, i, 0); +} + +mp_bitcnt_t +mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit) +{ + mp_size_t i; + i = bit / GMP_LIMB_BITS; + + return mpn_common_scan (~ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)), + i, ptr, i, GMP_LIMB_MAX); +} + +void +mpn_com (mp_ptr rp, mp_srcptr up, mp_size_t n) +{ + while (--n >= 0) + *rp++ = ~ *up++; +} + +mp_limb_t +mpn_neg (mp_ptr rp, mp_srcptr up, mp_size_t n) +{ + while (*up == 0) + { + *rp = 0; + if (!--n) + return 0; + ++up; ++rp; + } + *rp = - *up; + mpn_com (++rp, ++up, --n); + return 1; +} + + +/* MPN division interface. */ + +/* The 3/2 inverse is defined as + + m = floor( (B^3-1) / (B u1 + u0)) - B +*/ +mp_limb_t +mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0) +{ + mp_limb_t r, p, m, ql; + unsigned ul, uh, qh; + + assert (u1 >= GMP_LIMB_HIGHBIT); + + /* For notation, let b denote the half-limb base, so that B = b^2. + Split u1 = b uh + ul. */ + ul = u1 & GMP_LLIMB_MASK; + uh = u1 >> (GMP_LIMB_BITS / 2); + + /* Approximation of the high half of quotient. Differs from the 2/1 + inverse of the half limb uh, since we have already subtracted + u0. */ + qh = ~u1 / uh; + + /* Adjust to get a half-limb 3/2 inverse, i.e., we want + + qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u + = floor( (b (~u) + b-1) / u), + + and the remainder + + r = b (~u) + b-1 - qh (b uh + ul) + = b (~u - qh uh) + b-1 - qh ul + + Subtraction of qh ul may underflow, which implies adjustments. + But by normalization, 2 u >= B > qh ul, so we need to adjust by + at most 2. + */ + + r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK; + + p = (mp_limb_t) qh * ul; + /* Adjustment steps taken from udiv_qrnnd_c */ + if (r < p) + { + qh--; + r += u1; + if (r >= u1) /* i.e. we didn't get carry when adding to r */ + if (r < p) + { + qh--; + r += u1; + } + } + r -= p; + + /* Low half of the quotient is + + ql = floor ( (b r + b-1) / u1). + + This is a 3/2 division (on half-limbs), for which qh is a + suitable inverse. */ + + p = (r >> (GMP_LIMB_BITS / 2)) * qh + r; + /* Unlike full-limb 3/2, we can add 1 without overflow. For this to + work, it is essential that ql is a full mp_limb_t. */ + ql = (p >> (GMP_LIMB_BITS / 2)) + 1; + + /* By the 3/2 trick, we don't need the high half limb. */ + r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1; + + if (r >= (p << (GMP_LIMB_BITS / 2))) + { + ql--; + r += u1; + } + m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql; + if (r >= u1) + { + m++; + r -= u1; + } + + /* Now m is the 2/1 invers of u1. If u0 > 0, adjust it to become a + 3/2 inverse. */ + if (u0 > 0) + { + mp_limb_t th, tl; + r = ~r; + r += u0; + if (r < u0) + { + m--; + if (r >= u1) + { + m--; + r -= u1; + } + r -= u1; + } + gmp_umul_ppmm (th, tl, u0, m); + r += th; + if (r < th) + { + m--; + m -= ((r > u1) | ((r == u1) & (tl > u0))); + } + } + + return m; +} + +struct gmp_div_inverse +{ + /* Normalization shift count. */ + unsigned shift; + /* Normalized divisor (d0 unused for mpn_div_qr_1) */ + mp_limb_t d1, d0; + /* Inverse, for 2/1 or 3/2. */ + mp_limb_t di; +}; + +static void +mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d) +{ + unsigned shift; + + assert (d > 0); + gmp_clz (shift, d); + inv->shift = shift; + inv->d1 = d << shift; + inv->di = mpn_invert_limb (inv->d1); +} + +static void +mpn_div_qr_2_invert (struct gmp_div_inverse *inv, + mp_limb_t d1, mp_limb_t d0) +{ + unsigned shift; + + assert (d1 > 0); + gmp_clz (shift, d1); + inv->shift = shift; + if (shift > 0) + { + d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift)); + d0 <<= shift; + } + inv->d1 = d1; + inv->d0 = d0; + inv->di = mpn_invert_3by2 (d1, d0); +} + +static void +mpn_div_qr_invert (struct gmp_div_inverse *inv, + mp_srcptr dp, mp_size_t dn) +{ + assert (dn > 0); + + if (dn == 1) + mpn_div_qr_1_invert (inv, dp[0]); + else if (dn == 2) + mpn_div_qr_2_invert (inv, dp[1], dp[0]); + else + { + unsigned shift; + mp_limb_t d1, d0; + + d1 = dp[dn-1]; + d0 = dp[dn-2]; + assert (d1 > 0); + gmp_clz (shift, d1); + inv->shift = shift; + if (shift > 0) + { + d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift)); + d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift)); + } + inv->d1 = d1; + inv->d0 = d0; + inv->di = mpn_invert_3by2 (d1, d0); + } +} + +/* Not matching current public gmp interface, rather corresponding to + the sbpi1_div_* functions. */ +static mp_limb_t +mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn, + const struct gmp_div_inverse *inv) +{ + mp_limb_t d, di; + mp_limb_t r; + mp_ptr tp = NULL; + + if (inv->shift > 0) + { + tp = gmp_xalloc_limbs (nn); + r = mpn_lshift (tp, np, nn, inv->shift); + np = tp; + } + else + r = 0; + + d = inv->d1; + di = inv->di; + while (--nn >= 0) + { + mp_limb_t q; + + gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di); + if (qp) + qp[nn] = q; + } + if (inv->shift > 0) + gmp_free (tp); + + return r >> inv->shift; +} + +static mp_limb_t +mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d) +{ + assert (d > 0); + + /* Special case for powers of two. */ + if ((d & (d-1)) == 0) + { + mp_limb_t r = np[0] & (d-1); + if (qp) + { + if (d <= 1) + mpn_copyi (qp, np, nn); + else + { + unsigned shift; + gmp_ctz (shift, d); + mpn_rshift (qp, np, nn, shift); + } + } + return r; + } + else + { + struct gmp_div_inverse inv; + mpn_div_qr_1_invert (&inv, d); + return mpn_div_qr_1_preinv (qp, np, nn, &inv); + } +} + +static void +mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn, + const struct gmp_div_inverse *inv) +{ + unsigned shift; + mp_size_t i; + mp_limb_t d1, d0, di, r1, r0; + mp_ptr tp; + + assert (nn >= 2); + shift = inv->shift; + d1 = inv->d1; + d0 = inv->d0; + di = inv->di; + + if (shift > 0) + { + tp = gmp_xalloc_limbs (nn); + r1 = mpn_lshift (tp, np, nn, shift); + np = tp; + } + else + r1 = 0; + + r0 = np[nn - 1]; + + i = nn - 2; + do + { + mp_limb_t n0, q; + n0 = np[i]; + gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di); + + if (qp) + qp[i] = q; + } + while (--i >= 0); + + if (shift > 0) + { + assert ((r0 << (GMP_LIMB_BITS - shift)) == 0); + r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift)); + r1 >>= shift; + + gmp_free (tp); + } + + rp[1] = r1; + rp[0] = r0; +} + +#if 0 +static void +mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn, + mp_limb_t d1, mp_limb_t d0) +{ + struct gmp_div_inverse inv; + assert (nn >= 2); + + mpn_div_qr_2_invert (&inv, d1, d0); + mpn_div_qr_2_preinv (qp, rp, np, nn, &inv); +} +#endif + +static void +mpn_div_qr_pi1 (mp_ptr qp, + mp_ptr np, mp_size_t nn, mp_limb_t n1, + mp_srcptr dp, mp_size_t dn, + mp_limb_t dinv) +{ + mp_size_t i; + + mp_limb_t d1, d0; + mp_limb_t cy, cy1; + mp_limb_t q; + + assert (dn > 2); + assert (nn >= dn); + + d1 = dp[dn - 1]; + d0 = dp[dn - 2]; + + assert ((d1 & GMP_LIMB_HIGHBIT) != 0); + /* Iteration variable is the index of the q limb. + * + * We divide + * by + */ + + i = nn - dn; + do + { + mp_limb_t n0 = np[dn-1+i]; + + if (n1 == d1 && n0 == d0) + { + q = GMP_LIMB_MAX; + mpn_submul_1 (np+i, dp, dn, q); + n1 = np[dn-1+i]; /* update n1, last loop's value will now be invalid */ + } + else + { + gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv); + + cy = mpn_submul_1 (np + i, dp, dn-2, q); + + cy1 = n0 < cy; + n0 = n0 - cy; + cy = n1 < cy1; + n1 = n1 - cy1; + np[dn-2+i] = n0; + + if (cy != 0) + { + n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1); + q--; + } + } + + if (qp) + qp[i] = q; + } + while (--i >= 0); + + np[dn - 1] = n1; +} + +static void +mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn, + mp_srcptr dp, mp_size_t dn, + const struct gmp_div_inverse *inv) +{ + assert (dn > 0); + assert (nn >= dn); + + if (dn == 1) + np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv); + else if (dn == 2) + mpn_div_qr_2_preinv (qp, np, np, nn, inv); + else + { + mp_limb_t nh; + unsigned shift; + + assert (inv->d1 == dp[dn-1]); + assert (inv->d0 == dp[dn-2]); + assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0); + + shift = inv->shift; + if (shift > 0) + nh = mpn_lshift (np, np, nn, shift); + else + nh = 0; + + mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di); + + if (shift > 0) + gmp_assert_nocarry (mpn_rshift (np, np, dn, shift)); + } +} + +static void +mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn) +{ + struct gmp_div_inverse inv; + mp_ptr tp = NULL; + + assert (dn > 0); + assert (nn >= dn); + + mpn_div_qr_invert (&inv, dp, dn); + if (dn > 2 && inv.shift > 0) + { + tp = gmp_xalloc_limbs (dn); + gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift)); + dp = tp; + } + mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv); + if (tp) + gmp_free (tp); +} + + +/* MPN base conversion. */ +static unsigned +mpn_base_power_of_two_p (unsigned b) +{ + switch (b) + { + case 2: return 1; + case 4: return 2; + case 8: return 3; + case 16: return 4; + case 32: return 5; + case 64: return 6; + case 128: return 7; + case 256: return 8; + default: return 0; + } +} + +struct mpn_base_info +{ + /* bb is the largest power of the base which fits in one limb, and + exp is the corresponding exponent. */ + unsigned exp; + mp_limb_t bb; +}; + +static void +mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b) +{ + mp_limb_t m; + mp_limb_t p; + unsigned exp; + + m = GMP_LIMB_MAX / b; + for (exp = 1, p = b; p <= m; exp++) + p *= b; + + info->exp = exp; + info->bb = p; +} + +static mp_bitcnt_t +mpn_limb_size_in_base_2 (mp_limb_t u) +{ + unsigned shift; + + assert (u > 0); + gmp_clz (shift, u); + return GMP_LIMB_BITS - shift; +} + +static size_t +mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un) +{ + unsigned char mask; + size_t sn, j; + mp_size_t i; + unsigned shift; + + sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]) + + bits - 1) / bits; + + mask = (1U << bits) - 1; + + for (i = 0, j = sn, shift = 0; j-- > 0;) + { + unsigned char digit = up[i] >> shift; + + shift += bits; + + if (shift >= GMP_LIMB_BITS && ++i < un) + { + shift -= GMP_LIMB_BITS; + digit |= up[i] << (bits - shift); + } + sp[j] = digit & mask; + } + return sn; +} + +/* We generate digits from the least significant end, and reverse at + the end. */ +static size_t +mpn_limb_get_str (unsigned char *sp, mp_limb_t w, + const struct gmp_div_inverse *binv) +{ + mp_size_t i; + for (i = 0; w > 0; i++) + { + mp_limb_t h, l, r; + + h = w >> (GMP_LIMB_BITS - binv->shift); + l = w << binv->shift; + + gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di); + assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0); + r >>= binv->shift; + + sp[i] = r; + } + return i; +} + +static size_t +mpn_get_str_other (unsigned char *sp, + int base, const struct mpn_base_info *info, + mp_ptr up, mp_size_t un) +{ + struct gmp_div_inverse binv; + size_t sn; + size_t i; + + mpn_div_qr_1_invert (&binv, base); + + sn = 0; + + if (un > 1) + { + struct gmp_div_inverse bbinv; + mpn_div_qr_1_invert (&bbinv, info->bb); + + do + { + mp_limb_t w; + size_t done; + w = mpn_div_qr_1_preinv (up, up, un, &bbinv); + un -= (up[un-1] == 0); + done = mpn_limb_get_str (sp + sn, w, &binv); + + for (sn += done; done < info->exp; done++) + sp[sn++] = 0; + } + while (un > 1); + } + sn += mpn_limb_get_str (sp + sn, up[0], &binv); + + /* Reverse order */ + for (i = 0; 2*i + 1 < sn; i++) + { + unsigned char t = sp[i]; + sp[i] = sp[sn - i - 1]; + sp[sn - i - 1] = t; + } + + return sn; +} + +size_t +mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un) +{ + unsigned bits; + + assert (un > 0); + assert (up[un-1] > 0); + + bits = mpn_base_power_of_two_p (base); + if (bits) + return mpn_get_str_bits (sp, bits, up, un); + else + { + struct mpn_base_info info; + + mpn_get_base_info (&info, base); + return mpn_get_str_other (sp, base, &info, up, un); + } +} + +static mp_size_t +mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn, + unsigned bits) +{ + mp_size_t rn; + size_t j; + unsigned shift; + + for (j = sn, rn = 0, shift = 0; j-- > 0; ) + { + if (shift == 0) + { + rp[rn++] = sp[j]; + shift += bits; + } + else + { + rp[rn-1] |= (mp_limb_t) sp[j] << shift; + shift += bits; + if (shift >= GMP_LIMB_BITS) + { + shift -= GMP_LIMB_BITS; + if (shift > 0) + rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift); + } + } + } + rn = mpn_normalized_size (rp, rn); + return rn; +} + +/* Result is usually normalized, except for all-zero input, in which + case a single zero limb is written at *RP, and 1 is returned. */ +static mp_size_t +mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn, + mp_limb_t b, const struct mpn_base_info *info) +{ + mp_size_t rn; + mp_limb_t w; + unsigned k; + size_t j; + + assert (sn > 0); + + k = 1 + (sn - 1) % info->exp; + + j = 0; + w = sp[j++]; + while (--k != 0) + w = w * b + sp[j++]; + + rp[0] = w; + + for (rn = 1; j < sn;) + { + mp_limb_t cy; + + w = sp[j++]; + for (k = 1; k < info->exp; k++) + w = w * b + sp[j++]; + + cy = mpn_mul_1 (rp, rp, rn, info->bb); + cy += mpn_add_1 (rp, rp, rn, w); + if (cy > 0) + rp[rn++] = cy; + } + assert (j == sn); + + return rn; +} + +mp_size_t +mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base) +{ + unsigned bits; + + if (sn == 0) + return 0; + + bits = mpn_base_power_of_two_p (base); + if (bits) + return mpn_set_str_bits (rp, sp, sn, bits); + else + { + struct mpn_base_info info; + + mpn_get_base_info (&info, base); + return mpn_set_str_other (rp, sp, sn, base, &info); + } +} + + +/* MPZ interface */ +void +mpz_init (mpz_t r) +{ + static const mp_limb_t dummy_limb = 0xc1a0; + + r->_mp_alloc = 0; + r->_mp_size = 0; + r->_mp_d = (mp_ptr) &dummy_limb; +} + +/* The utility of this function is a bit limited, since many functions + assigns the result variable using mpz_swap. */ +void +mpz_init2 (mpz_t r, mp_bitcnt_t bits) +{ + mp_size_t rn; + + bits -= (bits != 0); /* Round down, except if 0 */ + rn = 1 + bits / GMP_LIMB_BITS; + + r->_mp_alloc = rn; + r->_mp_size = 0; + r->_mp_d = gmp_xalloc_limbs (rn); +} + +void +mpz_clear (mpz_t r) +{ + if (r->_mp_alloc) + gmp_free (r->_mp_d); +} + +static mp_ptr +mpz_realloc (mpz_t r, mp_size_t size) +{ + size = GMP_MAX (size, 1); + + if (r->_mp_alloc) + r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size); + else + r->_mp_d = gmp_xalloc_limbs (size); + r->_mp_alloc = size; + + if (GMP_ABS (r->_mp_size) > size) + r->_mp_size = 0; + + return r->_mp_d; +} + +/* Realloc for an mpz_t WHAT if it has less than NEEDED limbs. */ +#define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc \ + ? mpz_realloc(z,n) \ + : (z)->_mp_d) + +/* MPZ assignment and basic conversions. */ +void +mpz_set_si (mpz_t r, signed long int x) +{ + if (x >= 0) + mpz_set_ui (r, x); + else /* (x < 0) */ + { + r->_mp_size = -1; + MPZ_REALLOC (r, 1)[0] = GMP_NEG_CAST (unsigned long int, x); + } +} + +void +mpz_set_ui (mpz_t r, unsigned long int x) +{ + if (x > 0) + { + r->_mp_size = 1; + MPZ_REALLOC (r, 1)[0] = x; + } + else + r->_mp_size = 0; +} + +void +mpz_set (mpz_t r, const mpz_t x) +{ + /* Allow the NOP r == x */ + if (r != x) + { + mp_size_t n; + mp_ptr rp; + + n = GMP_ABS (x->_mp_size); + rp = MPZ_REALLOC (r, n); + + mpn_copyi (rp, x->_mp_d, n); + r->_mp_size = x->_mp_size; + } +} + +void +mpz_init_set_si (mpz_t r, signed long int x) +{ + mpz_init (r); + mpz_set_si (r, x); +} + +void +mpz_init_set_ui (mpz_t r, unsigned long int x) +{ + mpz_init (r); + mpz_set_ui (r, x); +} + +void +mpz_init_set (mpz_t r, const mpz_t x) +{ + mpz_init (r); + mpz_set (r, x); +} + +int +mpz_fits_slong_p (const mpz_t u) +{ + mp_size_t us = u->_mp_size; + + if (us == 1) + return u->_mp_d[0] < GMP_LIMB_HIGHBIT; + else if (us == -1) + return u->_mp_d[0] <= GMP_LIMB_HIGHBIT; + else + return (us == 0); +} + +int +mpz_fits_ulong_p (const mpz_t u) +{ + mp_size_t us = u->_mp_size; + + return (us == (us > 0)); +} + +long int +mpz_get_si (const mpz_t u) +{ + if (u->_mp_size < 0) + /* This expression is necessary to properly handle 0x80000000 */ + return -1 - (long) ((u->_mp_d[0] - 1) & ~GMP_LIMB_HIGHBIT); + else + return (long) (mpz_get_ui (u) & ~GMP_LIMB_HIGHBIT); +} + +unsigned long int +mpz_get_ui (const mpz_t u) +{ + return u->_mp_size == 0 ? 0 : u->_mp_d[0]; +} + +size_t +mpz_size (const mpz_t u) +{ + return GMP_ABS (u->_mp_size); +} + +mp_limb_t +mpz_getlimbn (const mpz_t u, mp_size_t n) +{ + if (n >= 0 && n < GMP_ABS (u->_mp_size)) + return u->_mp_d[n]; + else + return 0; +} + +void +mpz_realloc2 (mpz_t x, mp_bitcnt_t n) +{ + mpz_realloc (x, 1 + (n - (n != 0)) / GMP_LIMB_BITS); +} + +mp_srcptr +mpz_limbs_read (mpz_srcptr x) +{ + return x->_mp_d; +} + +mp_ptr +mpz_limbs_modify (mpz_t x, mp_size_t n) +{ + assert (n > 0); + return MPZ_REALLOC (x, n); +} + +mp_ptr +mpz_limbs_write (mpz_t x, mp_size_t n) +{ + return mpz_limbs_modify (x, n); +} + +void +mpz_limbs_finish (mpz_t x, mp_size_t xs) +{ + mp_size_t xn; + xn = mpn_normalized_size (x->_mp_d, GMP_ABS (xs)); + x->_mp_size = xs < 0 ? -xn : xn; +} + +mpz_srcptr +mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs) +{ + x->_mp_alloc = 0; + x->_mp_d = (mp_ptr) xp; + mpz_limbs_finish (x, xs); + return x; +} + + +/* Conversions and comparison to double. */ +void +mpz_set_d (mpz_t r, double x) +{ + int sign; + mp_ptr rp; + mp_size_t rn, i; + double B; + double Bi; + mp_limb_t f; + + /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is + zero or infinity. */ + if (x != x || x == x * 0.5) + { + r->_mp_size = 0; + return; + } + + sign = x < 0.0 ; + if (sign) + x = - x; + + if (x < 1.0) + { + r->_mp_size = 0; + return; + } + B = 2.0 * (double) GMP_LIMB_HIGHBIT; + Bi = 1.0 / B; + for (rn = 1; x >= B; rn++) + x *= Bi; + + rp = MPZ_REALLOC (r, rn); + + f = (mp_limb_t) x; + x -= f; + assert (x < 1.0); + i = rn-1; + rp[i] = f; + while (--i >= 0) + { + x = B * x; + f = (mp_limb_t) x; + x -= f; + assert (x < 1.0); + rp[i] = f; + } + + r->_mp_size = sign ? - rn : rn; +} + +void +mpz_init_set_d (mpz_t r, double x) +{ + mpz_init (r); + mpz_set_d (r, x); +} + +double +mpz_get_d (const mpz_t u) +{ + mp_size_t un; + double x; + double B = 2.0 * (double) GMP_LIMB_HIGHBIT; + + un = GMP_ABS (u->_mp_size); + + if (un == 0) + return 0.0; + + x = u->_mp_d[--un]; + while (un > 0) + x = B*x + u->_mp_d[--un]; + + if (u->_mp_size < 0) + x = -x; + + return x; +} + +int +mpz_cmpabs_d (const mpz_t x, double d) +{ + mp_size_t xn; + double B, Bi; + mp_size_t i; + + xn = x->_mp_size; + d = GMP_ABS (d); + + if (xn != 0) + { + xn = GMP_ABS (xn); + + B = 2.0 * (double) GMP_LIMB_HIGHBIT; + Bi = 1.0 / B; + + /* Scale d so it can be compared with the top limb. */ + for (i = 1; i < xn; i++) + d *= Bi; + + if (d >= B) + return -1; + + /* Compare floor(d) to top limb, subtract and cancel when equal. */ + for (i = xn; i-- > 0;) + { + mp_limb_t f, xl; + + f = (mp_limb_t) d; + xl = x->_mp_d[i]; + if (xl > f) + return 1; + else if (xl < f) + return -1; + d = B * (d - f); + } + } + return - (d > 0.0); +} + +int +mpz_cmp_d (const mpz_t x, double d) +{ + if (x->_mp_size < 0) + { + if (d >= 0.0) + return -1; + else + return -mpz_cmpabs_d (x, d); + } + else + { + if (d < 0.0) + return 1; + else + return mpz_cmpabs_d (x, d); + } +} + + +/* MPZ comparisons and the like. */ +int +mpz_sgn (const mpz_t u) +{ + return GMP_CMP (u->_mp_size, 0); +} + +int +mpz_cmp_si (const mpz_t u, long v) +{ + mp_size_t usize = u->_mp_size; + + if (usize < -1) + return -1; + else if (v >= 0) + return mpz_cmp_ui (u, v); + else if (usize >= 0) + return 1; + else /* usize == -1 */ + return GMP_CMP (GMP_NEG_CAST (mp_limb_t, v), u->_mp_d[0]); +} + +int +mpz_cmp_ui (const mpz_t u, unsigned long v) +{ + mp_size_t usize = u->_mp_size; + + if (usize > 1) + return 1; + else if (usize < 0) + return -1; + else + return GMP_CMP (mpz_get_ui (u), v); +} + +int +mpz_cmp (const mpz_t a, const mpz_t b) +{ + mp_size_t asize = a->_mp_size; + mp_size_t bsize = b->_mp_size; + + if (asize != bsize) + return (asize < bsize) ? -1 : 1; + else if (asize >= 0) + return mpn_cmp (a->_mp_d, b->_mp_d, asize); + else + return mpn_cmp (b->_mp_d, a->_mp_d, -asize); +} + +int +mpz_cmpabs_ui (const mpz_t u, unsigned long v) +{ + if (GMP_ABS (u->_mp_size) > 1) + return 1; + else + return GMP_CMP (mpz_get_ui (u), v); +} + +int +mpz_cmpabs (const mpz_t u, const mpz_t v) +{ + return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size), + v->_mp_d, GMP_ABS (v->_mp_size)); +} + +void +mpz_abs (mpz_t r, const mpz_t u) +{ + mpz_set (r, u); + r->_mp_size = GMP_ABS (r->_mp_size); +} + +void +mpz_neg (mpz_t r, const mpz_t u) +{ + mpz_set (r, u); + r->_mp_size = -r->_mp_size; +} + +void +mpz_swap (mpz_t u, mpz_t v) +{ + MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size); + MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc); + MP_PTR_SWAP (u->_mp_d, v->_mp_d); +} + + +/* MPZ addition and subtraction */ + +/* Adds to the absolute value. Returns new size, but doesn't store it. */ +static mp_size_t +mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b) +{ + mp_size_t an; + mp_ptr rp; + mp_limb_t cy; + + an = GMP_ABS (a->_mp_size); + if (an == 0) + { + MPZ_REALLOC (r, 1)[0] = b; + return b > 0; + } + + rp = MPZ_REALLOC (r, an + 1); + + cy = mpn_add_1 (rp, a->_mp_d, an, b); + rp[an] = cy; + an += cy; + + return an; +} + +/* Subtract from the absolute value. Returns new size, (or -1 on underflow), + but doesn't store it. */ +static mp_size_t +mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b) +{ + mp_size_t an = GMP_ABS (a->_mp_size); + mp_ptr rp; + + if (an == 0) + { + MPZ_REALLOC (r, 1)[0] = b; + return -(b > 0); + } + rp = MPZ_REALLOC (r, an); + if (an == 1 && a->_mp_d[0] < b) + { + rp[0] = b - a->_mp_d[0]; + return -1; + } + else + { + gmp_assert_nocarry (mpn_sub_1 (rp, a->_mp_d, an, b)); + return mpn_normalized_size (rp, an); + } +} + +void +mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b) +{ + if (a->_mp_size >= 0) + r->_mp_size = mpz_abs_add_ui (r, a, b); + else + r->_mp_size = -mpz_abs_sub_ui (r, a, b); +} + +void +mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b) +{ + if (a->_mp_size < 0) + r->_mp_size = -mpz_abs_add_ui (r, a, b); + else + r->_mp_size = mpz_abs_sub_ui (r, a, b); +} + +void +mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b) +{ + if (b->_mp_size < 0) + r->_mp_size = mpz_abs_add_ui (r, b, a); + else + r->_mp_size = -mpz_abs_sub_ui (r, b, a); +} + +static mp_size_t +mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t an = GMP_ABS (a->_mp_size); + mp_size_t bn = GMP_ABS (b->_mp_size); + mp_ptr rp; + mp_limb_t cy; + + if (an < bn) + { + MPZ_SRCPTR_SWAP (a, b); + MP_SIZE_T_SWAP (an, bn); + } + + rp = MPZ_REALLOC (r, an + 1); + cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn); + + rp[an] = cy; + + return an + cy; +} + +static mp_size_t +mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t an = GMP_ABS (a->_mp_size); + mp_size_t bn = GMP_ABS (b->_mp_size); + int cmp; + mp_ptr rp; + + cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn); + if (cmp > 0) + { + rp = MPZ_REALLOC (r, an); + gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn)); + return mpn_normalized_size (rp, an); + } + else if (cmp < 0) + { + rp = MPZ_REALLOC (r, bn); + gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an)); + return -mpn_normalized_size (rp, bn); + } + else + return 0; +} + +void +mpz_add (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t rn; + + if ( (a->_mp_size ^ b->_mp_size) >= 0) + rn = mpz_abs_add (r, a, b); + else + rn = mpz_abs_sub (r, a, b); + + r->_mp_size = a->_mp_size >= 0 ? rn : - rn; +} + +void +mpz_sub (mpz_t r, const mpz_t a, const mpz_t b) +{ + mp_size_t rn; + + if ( (a->_mp_size ^ b->_mp_size) >= 0) + rn = mpz_abs_sub (r, a, b); + else + rn = mpz_abs_add (r, a, b); + + r->_mp_size = a->_mp_size >= 0 ? rn : - rn; +} + + +/* MPZ multiplication */ +void +mpz_mul_si (mpz_t r, const mpz_t u, long int v) +{ + if (v < 0) + { + mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v)); + mpz_neg (r, r); + } + else + mpz_mul_ui (r, u, (unsigned long int) v); +} + +void +mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v) +{ + mp_size_t un, us; + mp_ptr tp; + mp_limb_t cy; + + us = u->_mp_size; + + if (us == 0 || v == 0) + { + r->_mp_size = 0; + return; + } + + un = GMP_ABS (us); + + tp = MPZ_REALLOC (r, un + 1); + cy = mpn_mul_1 (tp, u->_mp_d, un, v); + tp[un] = cy; + + un += (cy > 0); + r->_mp_size = (us < 0) ? - un : un; +} + +void +mpz_mul (mpz_t r, const mpz_t u, const mpz_t v) +{ + int sign; + mp_size_t un, vn, rn; + mpz_t t; + mp_ptr tp; + + un = u->_mp_size; + vn = v->_mp_size; + + if (un == 0 || vn == 0) + { + r->_mp_size = 0; + return; + } + + sign = (un ^ vn) < 0; + + un = GMP_ABS (un); + vn = GMP_ABS (vn); + + mpz_init2 (t, (un + vn) * GMP_LIMB_BITS); + + tp = t->_mp_d; + if (un >= vn) + mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn); + else + mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un); + + rn = un + vn; + rn -= tp[rn-1] == 0; + + t->_mp_size = sign ? - rn : rn; + mpz_swap (r, t); + mpz_clear (t); +} + +void +mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits) +{ + mp_size_t un, rn; + mp_size_t limbs; + unsigned shift; + mp_ptr rp; + + un = GMP_ABS (u->_mp_size); + if (un == 0) + { + r->_mp_size = 0; + return; + } + + limbs = bits / GMP_LIMB_BITS; + shift = bits % GMP_LIMB_BITS; + + rn = un + limbs + (shift > 0); + rp = MPZ_REALLOC (r, rn); + if (shift > 0) + { + mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift); + rp[rn-1] = cy; + rn -= (cy == 0); + } + else + mpn_copyd (rp + limbs, u->_mp_d, un); + + mpn_zero (rp, limbs); + + r->_mp_size = (u->_mp_size < 0) ? - rn : rn; +} + +void +mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v) +{ + mpz_t t; + mpz_init (t); + mpz_mul_ui (t, u, v); + mpz_add (r, r, t); + mpz_clear (t); +} + +void +mpz_submul_ui (mpz_t r, const mpz_t u, unsigned long int v) +{ + mpz_t t; + mpz_init (t); + mpz_mul_ui (t, u, v); + mpz_sub (r, r, t); + mpz_clear (t); +} + +void +mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v) +{ + mpz_t t; + mpz_init (t); + mpz_mul (t, u, v); + mpz_add (r, r, t); + mpz_clear (t); +} + +void +mpz_submul (mpz_t r, const mpz_t u, const mpz_t v) +{ + mpz_t t; + mpz_init (t); + mpz_mul (t, u, v); + mpz_sub (r, r, t); + mpz_clear (t); +} + + +/* MPZ division */ +enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC }; + +/* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */ +static int +mpz_div_qr (mpz_t q, mpz_t r, + const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode) +{ + mp_size_t ns, ds, nn, dn, qs; + ns = n->_mp_size; + ds = d->_mp_size; + + if (ds == 0) + gmp_die("mpz_div_qr: Divide by zero."); + + if (ns == 0) + { + if (q) + q->_mp_size = 0; + if (r) + r->_mp_size = 0; + return 0; + } + + nn = GMP_ABS (ns); + dn = GMP_ABS (ds); + + qs = ds ^ ns; + + if (nn < dn) + { + if (mode == GMP_DIV_CEIL && qs >= 0) + { + /* q = 1, r = n - d */ + if (r) + mpz_sub (r, n, d); + if (q) + mpz_set_ui (q, 1); + } + else if (mode == GMP_DIV_FLOOR && qs < 0) + { + /* q = -1, r = n + d */ + if (r) + mpz_add (r, n, d); + if (q) + mpz_set_si (q, -1); + } + else + { + /* q = 0, r = d */ + if (r) + mpz_set (r, n); + if (q) + q->_mp_size = 0; + } + return 1; + } + else + { + mp_ptr np, qp; + mp_size_t qn, rn; + mpz_t tq, tr; + + mpz_init_set (tr, n); + np = tr->_mp_d; + + qn = nn - dn + 1; + + if (q) + { + mpz_init2 (tq, qn * GMP_LIMB_BITS); + qp = tq->_mp_d; + } + else + qp = NULL; + + mpn_div_qr (qp, np, nn, d->_mp_d, dn); + + if (qp) + { + qn -= (qp[qn-1] == 0); + + tq->_mp_size = qs < 0 ? -qn : qn; + } + rn = mpn_normalized_size (np, dn); + tr->_mp_size = ns < 0 ? - rn : rn; + + if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0) + { + if (q) + mpz_sub_ui (tq, tq, 1); + if (r) + mpz_add (tr, tr, d); + } + else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0) + { + if (q) + mpz_add_ui (tq, tq, 1); + if (r) + mpz_sub (tr, tr, d); + } + + if (q) + { + mpz_swap (tq, q); + mpz_clear (tq); + } + if (r) + mpz_swap (tr, r); + + mpz_clear (tr); + + return rn != 0; + } +} + +void +mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, r, n, d, GMP_DIV_CEIL); +} + +void +mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR); +} + +void +mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC); +} + +void +mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL); +} + +void +mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR); +} + +void +mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC); +} + +void +mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL); +} + +void +mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR); +} + +void +mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC); +} + +void +mpz_mod (mpz_t r, const mpz_t n, const mpz_t d) +{ + mpz_div_qr (NULL, r, n, d, d->_mp_size >= 0 ? GMP_DIV_FLOOR : GMP_DIV_CEIL); +} + +static void +mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index, + enum mpz_div_round_mode mode) +{ + mp_size_t un, qn; + mp_size_t limb_cnt; + mp_ptr qp; + int adjust; + + un = u->_mp_size; + if (un == 0) + { + q->_mp_size = 0; + return; + } + limb_cnt = bit_index / GMP_LIMB_BITS; + qn = GMP_ABS (un) - limb_cnt; + bit_index %= GMP_LIMB_BITS; + + if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */ + /* Note: Below, the final indexing at limb_cnt is valid because at + that point we have qn > 0. */ + adjust = (qn <= 0 + || !mpn_zero_p (u->_mp_d, limb_cnt) + || (u->_mp_d[limb_cnt] + & (((mp_limb_t) 1 << bit_index) - 1))); + else + adjust = 0; + + if (qn <= 0) + qn = 0; + else + { + qp = MPZ_REALLOC (q, qn); + + if (bit_index != 0) + { + mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index); + qn -= qp[qn - 1] == 0; + } + else + { + mpn_copyi (qp, u->_mp_d + limb_cnt, qn); + } + } + + q->_mp_size = qn; + + if (adjust) + mpz_add_ui (q, q, 1); + if (un < 0) + mpz_neg (q, q); +} + +static void +mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index, + enum mpz_div_round_mode mode) +{ + mp_size_t us, un, rn; + mp_ptr rp; + mp_limb_t mask; + + us = u->_mp_size; + if (us == 0 || bit_index == 0) + { + r->_mp_size = 0; + return; + } + rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; + assert (rn > 0); + + rp = MPZ_REALLOC (r, rn); + un = GMP_ABS (us); + + mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index); + + if (rn > un) + { + /* Quotient (with truncation) is zero, and remainder is + non-zero */ + if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */ + { + /* Have to negate and sign extend. */ + mp_size_t i; + + gmp_assert_nocarry (! mpn_neg (rp, u->_mp_d, un)); + for (i = un; i < rn - 1; i++) + rp[i] = GMP_LIMB_MAX; + + rp[rn-1] = mask; + us = -us; + } + else + { + /* Just copy */ + if (r != u) + mpn_copyi (rp, u->_mp_d, un); + + rn = un; + } + } + else + { + if (r != u) + mpn_copyi (rp, u->_mp_d, rn - 1); + + rp[rn-1] = u->_mp_d[rn-1] & mask; + + if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */ + { + /* If r != 0, compute 2^{bit_count} - r. */ + mpn_neg (rp, rp, rn); + + rp[rn-1] &= mask; + + /* us is not used for anything else, so we can modify it + here to indicate flipped sign. */ + us = -us; + } + } + rn = mpn_normalized_size (rp, rn); + r->_mp_size = us < 0 ? -rn : rn; +} + +void +mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL); +} + +void +mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR); +} + +void +mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC); +} + +void +mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL); +} + +void +mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR); +} + +void +mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt) +{ + mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC); +} + +void +mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d) +{ + gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC)); +} + +int +mpz_divisible_p (const mpz_t n, const mpz_t d) +{ + return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0; +} + +int +mpz_congruent_p (const mpz_t a, const mpz_t b, const mpz_t m) +{ + mpz_t t; + int res; + + /* a == b (mod 0) iff a == b */ + if (mpz_sgn (m) == 0) + return (mpz_cmp (a, b) == 0); + + mpz_init (t); + mpz_sub (t, a, b); + res = mpz_divisible_p (t, m); + mpz_clear (t); + + return res; +} + +static unsigned long +mpz_div_qr_ui (mpz_t q, mpz_t r, + const mpz_t n, unsigned long d, enum mpz_div_round_mode mode) +{ + mp_size_t ns, qn; + mp_ptr qp; + mp_limb_t rl; + mp_size_t rs; + + ns = n->_mp_size; + if (ns == 0) + { + if (q) + q->_mp_size = 0; + if (r) + r->_mp_size = 0; + return 0; + } + + qn = GMP_ABS (ns); + if (q) + qp = MPZ_REALLOC (q, qn); + else + qp = NULL; + + rl = mpn_div_qr_1 (qp, n->_mp_d, qn, d); + assert (rl < d); + + rs = rl > 0; + rs = (ns < 0) ? -rs : rs; + + if (rl > 0 && ( (mode == GMP_DIV_FLOOR && ns < 0) + || (mode == GMP_DIV_CEIL && ns >= 0))) + { + if (q) + gmp_assert_nocarry (mpn_add_1 (qp, qp, qn, 1)); + rl = d - rl; + rs = -rs; + } + + if (r) + { + MPZ_REALLOC (r, 1)[0] = rl; + r->_mp_size = rs; + } + if (q) + { + qn -= (qp[qn-1] == 0); + assert (qn == 0 || qp[qn-1] > 0); + + q->_mp_size = (ns < 0) ? - qn : qn; + } + + return rl; +} + +unsigned long +mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL); +} + +unsigned long +mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR); +} + +unsigned long +mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC); +} + +unsigned long +mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL); +} + +unsigned long +mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR); +} + +unsigned long +mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC); +} + +unsigned long +mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL); +} +unsigned long +mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR); +} +unsigned long +mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC); +} + +unsigned long +mpz_cdiv_ui (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL); +} + +unsigned long +mpz_fdiv_ui (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR); +} + +unsigned long +mpz_tdiv_ui (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC); +} + +unsigned long +mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR); +} + +void +mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d) +{ + gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC)); +} + +int +mpz_divisible_ui_p (const mpz_t n, unsigned long d) +{ + return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0; +} + + +/* GCD */ +static mp_limb_t +mpn_gcd_11 (mp_limb_t u, mp_limb_t v) +{ + unsigned shift; + + assert ( (u | v) > 0); + + if (u == 0) + return v; + else if (v == 0) + return u; + + gmp_ctz (shift, u | v); + + u >>= shift; + v >>= shift; + + if ( (u & 1) == 0) + MP_LIMB_T_SWAP (u, v); + + while ( (v & 1) == 0) + v >>= 1; + + while (u != v) + { + if (u > v) + { + u -= v; + do + u >>= 1; + while ( (u & 1) == 0); + } + else + { + v -= u; + do + v >>= 1; + while ( (v & 1) == 0); + } + } + return u << shift; +} + +unsigned long +mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v) +{ + mp_size_t un; + + if (v == 0) + { + if (g) + mpz_abs (g, u); + } + else + { + un = GMP_ABS (u->_mp_size); + if (un != 0) + v = mpn_gcd_11 (mpn_div_qr_1 (NULL, u->_mp_d, un, v), v); + + if (g) + mpz_set_ui (g, v); + } + + return v; +} + +static mp_bitcnt_t +mpz_make_odd (mpz_t r) +{ + mp_bitcnt_t shift; + + assert (r->_mp_size > 0); + /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */ + shift = mpn_common_scan (r->_mp_d[0], 0, r->_mp_d, 0, 0); + mpz_tdiv_q_2exp (r, r, shift); + + return shift; +} + +void +mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v) +{ + mpz_t tu, tv; + mp_bitcnt_t uz, vz, gz; + + if (u->_mp_size == 0) + { + mpz_abs (g, v); + return; + } + if (v->_mp_size == 0) + { + mpz_abs (g, u); + return; + } + + mpz_init (tu); + mpz_init (tv); + + mpz_abs (tu, u); + uz = mpz_make_odd (tu); + mpz_abs (tv, v); + vz = mpz_make_odd (tv); + gz = GMP_MIN (uz, vz); + + if (tu->_mp_size < tv->_mp_size) + mpz_swap (tu, tv); + + mpz_tdiv_r (tu, tu, tv); + if (tu->_mp_size == 0) + { + mpz_swap (g, tv); + } + else + for (;;) + { + int c; + + mpz_make_odd (tu); + c = mpz_cmp (tu, tv); + if (c == 0) + { + mpz_swap (g, tu); + break; + } + if (c < 0) + mpz_swap (tu, tv); + + if (tv->_mp_size == 1) + { + mp_limb_t vl = tv->_mp_d[0]; + mp_limb_t ul = mpz_tdiv_ui (tu, vl); + mpz_set_ui (g, mpn_gcd_11 (ul, vl)); + break; + } + mpz_sub (tu, tu, tv); + } + mpz_clear (tu); + mpz_clear (tv); + mpz_mul_2exp (g, g, gz); +} + +void +mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v) +{ + mpz_t tu, tv, s0, s1, t0, t1; + mp_bitcnt_t uz, vz, gz; + mp_bitcnt_t power; + + if (u->_mp_size == 0) + { + /* g = 0 u + sgn(v) v */ + signed long sign = mpz_sgn (v); + mpz_abs (g, v); + if (s) + mpz_set_ui (s, 0); + if (t) + mpz_set_si (t, sign); + return; + } + + if (v->_mp_size == 0) + { + /* g = sgn(u) u + 0 v */ + signed long sign = mpz_sgn (u); + mpz_abs (g, u); + if (s) + mpz_set_si (s, sign); + if (t) + mpz_set_ui (t, 0); + return; + } + + mpz_init (tu); + mpz_init (tv); + mpz_init (s0); + mpz_init (s1); + mpz_init (t0); + mpz_init (t1); + + mpz_abs (tu, u); + uz = mpz_make_odd (tu); + mpz_abs (tv, v); + vz = mpz_make_odd (tv); + gz = GMP_MIN (uz, vz); + + uz -= gz; + vz -= gz; + + /* Cofactors corresponding to odd gcd. gz handled later. */ + if (tu->_mp_size < tv->_mp_size) + { + mpz_swap (tu, tv); + MPZ_SRCPTR_SWAP (u, v); + MPZ_PTR_SWAP (s, t); + MP_BITCNT_T_SWAP (uz, vz); + } + + /* Maintain + * + * u = t0 tu + t1 tv + * v = s0 tu + s1 tv + * + * where u and v denote the inputs with common factors of two + * eliminated, and det (s0, t0; s1, t1) = 2^p. Then + * + * 2^p tu = s1 u - t1 v + * 2^p tv = -s0 u + t0 v + */ + + /* After initial division, tu = q tv + tu', we have + * + * u = 2^uz (tu' + q tv) + * v = 2^vz tv + * + * or + * + * t0 = 2^uz, t1 = 2^uz q + * s0 = 0, s1 = 2^vz + */ + + mpz_setbit (t0, uz); + mpz_tdiv_qr (t1, tu, tu, tv); + mpz_mul_2exp (t1, t1, uz); + + mpz_setbit (s1, vz); + power = uz + vz; + + if (tu->_mp_size > 0) + { + mp_bitcnt_t shift; + shift = mpz_make_odd (tu); + mpz_mul_2exp (t0, t0, shift); + mpz_mul_2exp (s0, s0, shift); + power += shift; + + for (;;) + { + int c; + c = mpz_cmp (tu, tv); + if (c == 0) + break; + + if (c < 0) + { + /* tv = tv' + tu + * + * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv' + * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */ + + mpz_sub (tv, tv, tu); + mpz_add (t0, t0, t1); + mpz_add (s0, s0, s1); + + shift = mpz_make_odd (tv); + mpz_mul_2exp (t1, t1, shift); + mpz_mul_2exp (s1, s1, shift); + } + else + { + mpz_sub (tu, tu, tv); + mpz_add (t1, t0, t1); + mpz_add (s1, s0, s1); + + shift = mpz_make_odd (tu); + mpz_mul_2exp (t0, t0, shift); + mpz_mul_2exp (s0, s0, shift); + } + power += shift; + } + } + + /* Now tv = odd part of gcd, and -s0 and t0 are corresponding + cofactors. */ + + mpz_mul_2exp (tv, tv, gz); + mpz_neg (s0, s0); + + /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To + adjust cofactors, we need u / g and v / g */ + + mpz_divexact (s1, v, tv); + mpz_abs (s1, s1); + mpz_divexact (t1, u, tv); + mpz_abs (t1, t1); + + while (power-- > 0) + { + /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */ + if (mpz_odd_p (s0) || mpz_odd_p (t0)) + { + mpz_sub (s0, s0, s1); + mpz_add (t0, t0, t1); + } + mpz_divexact_ui (s0, s0, 2); + mpz_divexact_ui (t0, t0, 2); + } + + /* Arrange so that |s| < |u| / 2g */ + mpz_add (s1, s0, s1); + if (mpz_cmpabs (s0, s1) > 0) + { + mpz_swap (s0, s1); + mpz_sub (t0, t0, t1); + } + if (u->_mp_size < 0) + mpz_neg (s0, s0); + if (v->_mp_size < 0) + mpz_neg (t0, t0); + + mpz_swap (g, tv); + if (s) + mpz_swap (s, s0); + if (t) + mpz_swap (t, t0); + + mpz_clear (tu); + mpz_clear (tv); + mpz_clear (s0); + mpz_clear (s1); + mpz_clear (t0); + mpz_clear (t1); +} + +void +mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v) +{ + mpz_t g; + + if (u->_mp_size == 0 || v->_mp_size == 0) + { + r->_mp_size = 0; + return; + } + + mpz_init (g); + + mpz_gcd (g, u, v); + mpz_divexact (g, u, g); + mpz_mul (r, g, v); + + mpz_clear (g); + mpz_abs (r, r); +} + +void +mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v) +{ + if (v == 0 || u->_mp_size == 0) + { + r->_mp_size = 0; + return; + } + + v /= mpz_gcd_ui (NULL, u, v); + mpz_mul_ui (r, u, v); + + mpz_abs (r, r); +} + +int +mpz_invert (mpz_t r, const mpz_t u, const mpz_t m) +{ + mpz_t g, tr; + int invertible; + + if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0) + return 0; + + mpz_init (g); + mpz_init (tr); + + mpz_gcdext (g, tr, NULL, u, m); + invertible = (mpz_cmp_ui (g, 1) == 0); + + if (invertible) + { + if (tr->_mp_size < 0) + { + if (m->_mp_size >= 0) + mpz_add (tr, tr, m); + else + mpz_sub (tr, tr, m); + } + mpz_swap (r, tr); + } + + mpz_clear (g); + mpz_clear (tr); + return invertible; +} + + +/* Higher level operations (sqrt, pow and root) */ + +void +mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e) +{ + unsigned long bit; + mpz_t tr; + mpz_init_set_ui (tr, 1); + + bit = GMP_ULONG_HIGHBIT; + do + { + mpz_mul (tr, tr, tr); + if (e & bit) + mpz_mul (tr, tr, b); + bit >>= 1; + } + while (bit > 0); + + mpz_swap (r, tr); + mpz_clear (tr); +} + +void +mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e) +{ + mpz_t b; + mpz_pow_ui (r, mpz_roinit_n (b, &blimb, 1), e); +} + +void +mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m) +{ + mpz_t tr; + mpz_t base; + mp_size_t en, mn; + mp_srcptr mp; + struct gmp_div_inverse minv; + unsigned shift; + mp_ptr tp = NULL; + + en = GMP_ABS (e->_mp_size); + mn = GMP_ABS (m->_mp_size); + if (mn == 0) + gmp_die ("mpz_powm: Zero modulo."); + + if (en == 0) + { + mpz_set_ui (r, 1); + return; + } + + mp = m->_mp_d; + mpn_div_qr_invert (&minv, mp, mn); + shift = minv.shift; + + if (shift > 0) + { + /* To avoid shifts, we do all our reductions, except the final + one, using a *normalized* m. */ + minv.shift = 0; + + tp = gmp_xalloc_limbs (mn); + gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift)); + mp = tp; + } + + mpz_init (base); + + if (e->_mp_size < 0) + { + if (!mpz_invert (base, b, m)) + gmp_die ("mpz_powm: Negative exponent and non-invertible base."); + } + else + { + mp_size_t bn; + mpz_abs (base, b); + + bn = base->_mp_size; + if (bn >= mn) + { + mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv); + bn = mn; + } + + /* We have reduced the absolute value. Now take care of the + sign. Note that we get zero represented non-canonically as + m. */ + if (b->_mp_size < 0) + { + mp_ptr bp = MPZ_REALLOC (base, mn); + gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn)); + bn = mn; + } + base->_mp_size = mpn_normalized_size (base->_mp_d, bn); + } + mpz_init_set_ui (tr, 1); + + while (--en >= 0) + { + mp_limb_t w = e->_mp_d[en]; + mp_limb_t bit; + + bit = GMP_LIMB_HIGHBIT; + do + { + mpz_mul (tr, tr, tr); + if (w & bit) + mpz_mul (tr, tr, base); + if (tr->_mp_size > mn) + { + mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv); + tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn); + } + bit >>= 1; + } + while (bit > 0); + } + + /* Final reduction */ + if (tr->_mp_size >= mn) + { + minv.shift = shift; + mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv); + tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn); + } + if (tp) + gmp_free (tp); + + mpz_swap (r, tr); + mpz_clear (tr); + mpz_clear (base); +} + +void +mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m) +{ + mpz_t e; + mpz_powm (r, b, mpz_roinit_n (e, &elimb, 1), m); +} + +/* x=trunc(y^(1/z)), r=y-x^z */ +void +mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z) +{ + int sgn; + mpz_t t, u; + + sgn = y->_mp_size < 0; + if ((~z & sgn) != 0) + gmp_die ("mpz_rootrem: Negative argument, with even root."); + if (z == 0) + gmp_die ("mpz_rootrem: Zeroth root."); + + if (mpz_cmpabs_ui (y, 1) <= 0) { + if (x) + mpz_set (x, y); + if (r) + r->_mp_size = 0; + return; + } + + mpz_init (u); + mpz_init (t); + mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1); + + if (z == 2) /* simplify sqrt loop: z-1 == 1 */ + do { + mpz_swap (u, t); /* u = x */ + mpz_tdiv_q (t, y, u); /* t = y/x */ + mpz_add (t, t, u); /* t = y/x + x */ + mpz_tdiv_q_2exp (t, t, 1); /* x'= (y/x + x)/2 */ + } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */ + else /* z != 2 */ { + mpz_t v; + + mpz_init (v); + if (sgn) + mpz_neg (t, t); + + do { + mpz_swap (u, t); /* u = x */ + mpz_pow_ui (t, u, z - 1); /* t = x^(z-1) */ + mpz_tdiv_q (t, y, t); /* t = y/x^(z-1) */ + mpz_mul_ui (v, u, z - 1); /* v = x*(z-1) */ + mpz_add (t, t, v); /* t = y/x^(z-1) + x*(z-1) */ + mpz_tdiv_q_ui (t, t, z); /* x'=(y/x^(z-1) + x*(z-1))/z */ + } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */ + + mpz_clear (v); + } + + if (r) { + mpz_pow_ui (t, u, z); + mpz_sub (r, y, t); + } + if (x) + mpz_swap (x, u); + mpz_clear (u); + mpz_clear (t); +} + +int +mpz_root (mpz_t x, const mpz_t y, unsigned long z) +{ + int res; + mpz_t r; + + mpz_init (r); + mpz_rootrem (x, r, y, z); + res = r->_mp_size == 0; + mpz_clear (r); + + return res; +} + +/* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */ +void +mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u) +{ + mpz_rootrem (s, r, u, 2); +} + +void +mpz_sqrt (mpz_t s, const mpz_t u) +{ + mpz_rootrem (s, NULL, u, 2); +} + +int +mpz_perfect_square_p (const mpz_t u) +{ + if (u->_mp_size <= 0) + return (u->_mp_size == 0); + else + return mpz_root (NULL, u, 2); +} + +int +mpn_perfect_square_p (mp_srcptr p, mp_size_t n) +{ + mpz_t t; + + assert (n > 0); + assert (p [n-1] != 0); + return mpz_root (NULL, mpz_roinit_n (t, p, n), 2); +} + +mp_size_t +mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n) +{ + mpz_t s, r, u; + mp_size_t res; + + assert (n > 0); + assert (p [n-1] != 0); + + mpz_init (r); + mpz_init (s); + mpz_rootrem (s, r, mpz_roinit_n (u, p, n), 2); + + assert (s->_mp_size == (n+1)/2); + mpn_copyd (sp, s->_mp_d, s->_mp_size); + mpz_clear (s); + res = r->_mp_size; + if (rp) + mpn_copyd (rp, r->_mp_d, res); + mpz_clear (r); + return res; +} + +/* Combinatorics */ + +void +mpz_fac_ui (mpz_t x, unsigned long n) +{ + mpz_set_ui (x, n + (n == 0)); + while (n > 2) + mpz_mul_ui (x, x, --n); +} + +void +mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k) +{ + mpz_t t; + + mpz_set_ui (r, k <= n); + + if (k > (n >> 1)) + k = (k <= n) ? n - k : 0; + + mpz_init (t); + mpz_fac_ui (t, k); + + for (; k > 0; k--) + mpz_mul_ui (r, r, n--); + + mpz_divexact (r, r, t); + mpz_clear (t); +} + + +/* Primality testing */ +static int +gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y, + const mpz_t q, mp_bitcnt_t k) +{ + assert (k > 0); + + /* Caller must initialize y to the base. */ + mpz_powm (y, y, q, n); + + if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0) + return 1; + + while (--k > 0) + { + mpz_powm_ui (y, y, 2, n); + if (mpz_cmp (y, nm1) == 0) + return 1; + /* y == 1 means that the previous y was a non-trivial square root + of 1 (mod n). y == 0 means that n is a power of the base. + In either case, n is not prime. */ + if (mpz_cmp_ui (y, 1) <= 0) + return 0; + } + return 0; +} + +/* This product is 0xc0cfd797, and fits in 32 bits. */ +#define GMP_PRIME_PRODUCT \ + (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL) + +/* Bit (p+1)/2 is set, for each odd prime <= 61 */ +#define GMP_PRIME_MASK 0xc96996dcUL + +int +mpz_probab_prime_p (const mpz_t n, int reps) +{ + mpz_t nm1; + mpz_t q; + mpz_t y; + mp_bitcnt_t k; + int is_prime; + int j; + + /* Note that we use the absolute value of n only, for compatibility + with the real GMP. */ + if (mpz_even_p (n)) + return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0; + + /* Above test excludes n == 0 */ + assert (n->_mp_size != 0); + + if (mpz_cmpabs_ui (n, 64) < 0) + return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2; + + if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1) + return 0; + + /* All prime factors are >= 31. */ + if (mpz_cmpabs_ui (n, 31*31) < 0) + return 2; + + /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] = + j^2 + j + 41 using Euler's polynomial. We potentially stop early, + if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps > + 30 (a[30] == 971 > 31*31 == 961). */ + + mpz_init (nm1); + mpz_init (q); + mpz_init (y); + + /* Find q and k, where q is odd and n = 1 + 2**k * q. */ + nm1->_mp_size = mpz_abs_sub_ui (nm1, n, 1); + k = mpz_scan1 (nm1, 0); + mpz_tdiv_q_2exp (q, nm1, k); + + for (j = 0, is_prime = 1; is_prime & (j < reps); j++) + { + mpz_set_ui (y, (unsigned long) j*j+j+41); + if (mpz_cmp (y, nm1) >= 0) + { + /* Don't try any further bases. This "early" break does not affect + the result for any reasonable reps value (<=5000 was tested) */ + assert (j >= 30); + break; + } + is_prime = gmp_millerrabin (n, nm1, y, q, k); + } + mpz_clear (nm1); + mpz_clear (q); + mpz_clear (y); + + return is_prime; +} + + +/* Logical operations and bit manipulation. */ + +/* Numbers are treated as if represented in two's complement (and + infinitely sign extended). For a negative values we get the two's + complement from -x = ~x + 1, where ~ is bitwise complement. + Negation transforms + + xxxx10...0 + + into + + yyyy10...0 + + where yyyy is the bitwise complement of xxxx. So least significant + bits, up to and including the first one bit, are unchanged, and + the more significant bits are all complemented. + + To change a bit from zero to one in a negative number, subtract the + corresponding power of two from the absolute value. This can never + underflow. To change a bit from one to zero, add the corresponding + power of two, and this might overflow. E.g., if x = -001111, the + two's complement is 110001. Clearing the least significant bit, we + get two's complement 110000, and -010000. */ + +int +mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index) +{ + mp_size_t limb_index; + unsigned shift; + mp_size_t ds; + mp_size_t dn; + mp_limb_t w; + int bit; + + ds = d->_mp_size; + dn = GMP_ABS (ds); + limb_index = bit_index / GMP_LIMB_BITS; + if (limb_index >= dn) + return ds < 0; + + shift = bit_index % GMP_LIMB_BITS; + w = d->_mp_d[limb_index]; + bit = (w >> shift) & 1; + + if (ds < 0) + { + /* d < 0. Check if any of the bits below is set: If so, our bit + must be complemented. */ + if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0) + return bit ^ 1; + while (--limb_index >= 0) + if (d->_mp_d[limb_index] > 0) + return bit ^ 1; + } + return bit; +} + +static void +mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index) +{ + mp_size_t dn, limb_index; + mp_limb_t bit; + mp_ptr dp; + + dn = GMP_ABS (d->_mp_size); + + limb_index = bit_index / GMP_LIMB_BITS; + bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS); + + if (limb_index >= dn) + { + mp_size_t i; + /* The bit should be set outside of the end of the number. + We have to increase the size of the number. */ + dp = MPZ_REALLOC (d, limb_index + 1); + + dp[limb_index] = bit; + for (i = dn; i < limb_index; i++) + dp[i] = 0; + dn = limb_index + 1; + } + else + { + mp_limb_t cy; + + dp = d->_mp_d; + + cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit); + if (cy > 0) + { + dp = MPZ_REALLOC (d, dn + 1); + dp[dn++] = cy; + } + } + + d->_mp_size = (d->_mp_size < 0) ? - dn : dn; +} + +static void +mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index) +{ + mp_size_t dn, limb_index; + mp_ptr dp; + mp_limb_t bit; + + dn = GMP_ABS (d->_mp_size); + dp = d->_mp_d; + + limb_index = bit_index / GMP_LIMB_BITS; + bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS); + + assert (limb_index < dn); + + gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index, + dn - limb_index, bit)); + dn = mpn_normalized_size (dp, dn); + d->_mp_size = (d->_mp_size < 0) ? - dn : dn; +} + +void +mpz_setbit (mpz_t d, mp_bitcnt_t bit_index) +{ + if (!mpz_tstbit (d, bit_index)) + { + if (d->_mp_size >= 0) + mpz_abs_add_bit (d, bit_index); + else + mpz_abs_sub_bit (d, bit_index); + } +} + +void +mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index) +{ + if (mpz_tstbit (d, bit_index)) + { + if (d->_mp_size >= 0) + mpz_abs_sub_bit (d, bit_index); + else + mpz_abs_add_bit (d, bit_index); + } +} + +void +mpz_combit (mpz_t d, mp_bitcnt_t bit_index) +{ + if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0)) + mpz_abs_sub_bit (d, bit_index); + else + mpz_abs_add_bit (d, bit_index); +} + +void +mpz_com (mpz_t r, const mpz_t u) +{ + mpz_neg (r, u); + mpz_sub_ui (r, r, 1); +} + +void +mpz_and (mpz_t r, const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, rn, i; + mp_ptr up, vp, rp; + + mp_limb_t ux, vx, rx; + mp_limb_t uc, vc, rc; + mp_limb_t ul, vl, rl; + + un = GMP_ABS (u->_mp_size); + vn = GMP_ABS (v->_mp_size); + if (un < vn) + { + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (un, vn); + } + if (vn == 0) + { + r->_mp_size = 0; + return; + } + + uc = u->_mp_size < 0; + vc = v->_mp_size < 0; + rc = uc & vc; + + ux = -uc; + vx = -vc; + rx = -rc; + + /* If the smaller input is positive, higher limbs don't matter. */ + rn = vx ? un : vn; + + rp = MPZ_REALLOC (r, rn + (mp_size_t) rc); + + up = u->_mp_d; + vp = v->_mp_d; + + i = 0; + do + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + vl = (vp[i] ^ vx) + vc; + vc = vl < vc; + + rl = ( (ul & vl) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + while (++i < vn); + assert (vc == 0); + + for (; i < rn; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + rl = ( (ul & vx) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + if (rc) + rp[rn++] = rc; + else + rn = mpn_normalized_size (rp, rn); + + r->_mp_size = rx ? -rn : rn; +} + +void +mpz_ior (mpz_t r, const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, rn, i; + mp_ptr up, vp, rp; + + mp_limb_t ux, vx, rx; + mp_limb_t uc, vc, rc; + mp_limb_t ul, vl, rl; + + un = GMP_ABS (u->_mp_size); + vn = GMP_ABS (v->_mp_size); + if (un < vn) + { + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (un, vn); + } + if (vn == 0) + { + mpz_set (r, u); + return; + } + + uc = u->_mp_size < 0; + vc = v->_mp_size < 0; + rc = uc | vc; + + ux = -uc; + vx = -vc; + rx = -rc; + + /* If the smaller input is negative, by sign extension higher limbs + don't matter. */ + rn = vx ? vn : un; + + rp = MPZ_REALLOC (r, rn + (mp_size_t) rc); + + up = u->_mp_d; + vp = v->_mp_d; + + i = 0; + do + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + vl = (vp[i] ^ vx) + vc; + vc = vl < vc; + + rl = ( (ul | vl) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + while (++i < vn); + assert (vc == 0); + + for (; i < rn; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + rl = ( (ul | vx) ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + if (rc) + rp[rn++] = rc; + else + rn = mpn_normalized_size (rp, rn); + + r->_mp_size = rx ? -rn : rn; +} + +void +mpz_xor (mpz_t r, const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, i; + mp_ptr up, vp, rp; + + mp_limb_t ux, vx, rx; + mp_limb_t uc, vc, rc; + mp_limb_t ul, vl, rl; + + un = GMP_ABS (u->_mp_size); + vn = GMP_ABS (v->_mp_size); + if (un < vn) + { + MPZ_SRCPTR_SWAP (u, v); + MP_SIZE_T_SWAP (un, vn); + } + if (vn == 0) + { + mpz_set (r, u); + return; + } + + uc = u->_mp_size < 0; + vc = v->_mp_size < 0; + rc = uc ^ vc; + + ux = -uc; + vx = -vc; + rx = -rc; + + rp = MPZ_REALLOC (r, un + (mp_size_t) rc); + + up = u->_mp_d; + vp = v->_mp_d; + + i = 0; + do + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + vl = (vp[i] ^ vx) + vc; + vc = vl < vc; + + rl = (ul ^ vl ^ rx) + rc; + rc = rl < rc; + rp[i] = rl; + } + while (++i < vn); + assert (vc == 0); + + for (; i < un; i++) + { + ul = (up[i] ^ ux) + uc; + uc = ul < uc; + + rl = (ul ^ ux) + rc; + rc = rl < rc; + rp[i] = rl; + } + if (rc) + rp[un++] = rc; + else + un = mpn_normalized_size (rp, un); + + r->_mp_size = rx ? -un : un; +} + +static unsigned +gmp_popcount_limb (mp_limb_t x) +{ + unsigned c; + + /* Do 16 bits at a time, to avoid limb-sized constants. */ + for (c = 0; x > 0; x >>= 16) + { + unsigned w = ((x >> 1) & 0x5555) + (x & 0x5555); + w = ((w >> 2) & 0x3333) + (w & 0x3333); + w = ((w >> 4) & 0x0f0f) + (w & 0x0f0f); + w = (w >> 8) + (w & 0x00ff); + c += w; + } + return c; +} + +mp_bitcnt_t +mpn_popcount (mp_srcptr p, mp_size_t n) +{ + mp_size_t i; + mp_bitcnt_t c; + + for (c = 0, i = 0; i < n; i++) + c += gmp_popcount_limb (p[i]); + + return c; +} + +mp_bitcnt_t +mpz_popcount (const mpz_t u) +{ + mp_size_t un; + + un = u->_mp_size; + + if (un < 0) + return ~(mp_bitcnt_t) 0; + + return mpn_popcount (u->_mp_d, un); +} + +mp_bitcnt_t +mpz_hamdist (const mpz_t u, const mpz_t v) +{ + mp_size_t un, vn, i; + mp_limb_t uc, vc, ul, vl, comp; + mp_srcptr up, vp; + mp_bitcnt_t c; + + un = u->_mp_size; + vn = v->_mp_size; + + if ( (un ^ vn) < 0) + return ~(mp_bitcnt_t) 0; + + comp = - (uc = vc = (un < 0)); + if (uc) + { + assert (vn < 0); + un = -un; + vn = -vn; + } + + up = u->_mp_d; + vp = v->_mp_d; + + if (un < vn) + MPN_SRCPTR_SWAP (up, un, vp, vn); + + for (i = 0, c = 0; i < vn; i++) + { + ul = (up[i] ^ comp) + uc; + uc = ul < uc; + + vl = (vp[i] ^ comp) + vc; + vc = vl < vc; + + c += gmp_popcount_limb (ul ^ vl); + } + assert (vc == 0); + + for (; i < un; i++) + { + ul = (up[i] ^ comp) + uc; + uc = ul < uc; + + c += gmp_popcount_limb (ul ^ comp); + } + + return c; +} + +mp_bitcnt_t +mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit) +{ + mp_ptr up; + mp_size_t us, un, i; + mp_limb_t limb, ux; + + us = u->_mp_size; + un = GMP_ABS (us); + i = starting_bit / GMP_LIMB_BITS; + + /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit + for u<0. Notice this test picks up any u==0 too. */ + if (i >= un) + return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit); + + up = u->_mp_d; + ux = 0; + limb = up[i]; + + if (starting_bit != 0) + { + if (us < 0) + { + ux = mpn_zero_p (up, i); + limb = ~ limb + ux; + ux = - (mp_limb_t) (limb >= ux); + } + + /* Mask to 0 all bits before starting_bit, thus ignoring them. */ + limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS)); + } + + return mpn_common_scan (limb, i, up, un, ux); +} + +mp_bitcnt_t +mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit) +{ + mp_ptr up; + mp_size_t us, un, i; + mp_limb_t limb, ux; + + us = u->_mp_size; + ux = - (mp_limb_t) (us >= 0); + un = GMP_ABS (us); + i = starting_bit / GMP_LIMB_BITS; + + /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for + u<0. Notice this test picks up all cases of u==0 too. */ + if (i >= un) + return (ux ? starting_bit : ~(mp_bitcnt_t) 0); + + up = u->_mp_d; + limb = up[i] ^ ux; + + if (ux == 0) + limb -= mpn_zero_p (up, i); /* limb = ~(~limb + zero_p) */ + + /* Mask all bits before starting_bit, thus ignoring them. */ + limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS)); + + return mpn_common_scan (limb, i, up, un, ux); +} + + +/* MPZ base conversion. */ + +size_t +mpz_sizeinbase (const mpz_t u, int base) +{ + mp_size_t un; + mp_srcptr up; + mp_ptr tp; + mp_bitcnt_t bits; + struct gmp_div_inverse bi; + size_t ndigits; + + assert (base >= 2); + assert (base <= 36); + + un = GMP_ABS (u->_mp_size); + if (un == 0) + return 1; + + up = u->_mp_d; + + bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]); + switch (base) + { + case 2: + return bits; + case 4: + return (bits + 1) / 2; + case 8: + return (bits + 2) / 3; + case 16: + return (bits + 3) / 4; + case 32: + return (bits + 4) / 5; + /* FIXME: Do something more clever for the common case of base + 10. */ + } + + tp = gmp_xalloc_limbs (un); + mpn_copyi (tp, up, un); + mpn_div_qr_1_invert (&bi, base); + + ndigits = 0; + do + { + ndigits++; + mpn_div_qr_1_preinv (tp, tp, un, &bi); + un -= (tp[un-1] == 0); + } + while (un > 0); + + gmp_free (tp); + return ndigits; +} + +char * +mpz_get_str (char *sp, int base, const mpz_t u) +{ + unsigned bits; + const char *digits; + mp_size_t un; + size_t i, sn; + + if (base >= 0) + { + digits = "0123456789abcdefghijklmnopqrstuvwxyz"; + } + else + { + base = -base; + digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + } + if (base <= 1) + base = 10; + if (base > 36) + return NULL; + + sn = 1 + mpz_sizeinbase (u, base); + if (!sp) + sp = (char *) gmp_xalloc (1 + sn); + + un = GMP_ABS (u->_mp_size); + + if (un == 0) + { + sp[0] = '0'; + sp[1] = '\0'; + return sp; + } + + i = 0; + + if (u->_mp_size < 0) + sp[i++] = '-'; + + bits = mpn_base_power_of_two_p (base); + + if (bits) + /* Not modified in this case. */ + sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un); + else + { + struct mpn_base_info info; + mp_ptr tp; + + mpn_get_base_info (&info, base); + tp = gmp_xalloc_limbs (un); + mpn_copyi (tp, u->_mp_d, un); + + sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un); + gmp_free (tp); + } + + for (; i < sn; i++) + sp[i] = digits[(unsigned char) sp[i]]; + + sp[sn] = '\0'; + return sp; +} + +int +mpz_set_str (mpz_t r, const char *sp, int base) +{ + unsigned bits; + mp_size_t rn, alloc; + mp_ptr rp; + size_t dn; + int sign; + unsigned char *dp; + + assert (base == 0 || (base >= 2 && base <= 36)); + + while (isspace( (unsigned char) *sp)) + sp++; + + sign = (*sp == '-'); + sp += sign; + + if (base == 0) + { + if (sp[0] == '0') + { + if (sp[1] == 'x' || sp[1] == 'X') + { + base = 16; + sp += 2; + } + else if (sp[1] == 'b' || sp[1] == 'B') + { + base = 2; + sp += 2; + } + else + base = 8; + } + else + base = 10; + } + + if (!*sp) + { + r->_mp_size = 0; + return -1; + } + dp = (unsigned char *) gmp_xalloc (strlen (sp)); + + for (dn = 0; *sp; sp++) + { + unsigned digit; + + if (isspace ((unsigned char) *sp)) + continue; + else if (*sp >= '0' && *sp <= '9') + digit = *sp - '0'; + else if (*sp >= 'a' && *sp <= 'z') + digit = *sp - 'a' + 10; + else if (*sp >= 'A' && *sp <= 'Z') + digit = *sp - 'A' + 10; + else + digit = base; /* fail */ + + if (digit >= (unsigned) base) + { + gmp_free (dp); + r->_mp_size = 0; + return -1; + } + + dp[dn++] = digit; + } + + if (!dn) + { + gmp_free (dp); + r->_mp_size = 0; + return -1; + } + bits = mpn_base_power_of_two_p (base); + + if (bits > 0) + { + alloc = (dn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; + rp = MPZ_REALLOC (r, alloc); + rn = mpn_set_str_bits (rp, dp, dn, bits); + } + else + { + struct mpn_base_info info; + mpn_get_base_info (&info, base); + alloc = (dn + info.exp - 1) / info.exp; + rp = MPZ_REALLOC (r, alloc); + rn = mpn_set_str_other (rp, dp, dn, base, &info); + /* Normalization, needed for all-zero input. */ + assert (rn > 0); + rn -= rp[rn-1] == 0; + } + assert (rn <= alloc); + gmp_free (dp); + + r->_mp_size = sign ? - rn : rn; + + return 0; +} + +int +mpz_init_set_str (mpz_t r, const char *sp, int base) +{ + mpz_init (r); + return mpz_set_str (r, sp, base); +} + +size_t +mpz_out_str (FILE *stream, int base, const mpz_t x) +{ + char *str; + size_t len; + + str = mpz_get_str (NULL, base, x); + len = strlen (str); + len = fwrite (str, 1, len, stream); + gmp_free (str); + return len; +} + + +static int +gmp_detect_endian (void) +{ + static const int i = 2; + const unsigned char *p = (const unsigned char *) &i; + return 1 - *p; +} + +/* Import and export. Does not support nails. */ +void +mpz_import (mpz_t r, size_t count, int order, size_t size, int endian, + size_t nails, const void *src) +{ + const unsigned char *p; + ptrdiff_t word_step; + mp_ptr rp; + mp_size_t rn; + + /* The current (partial) limb. */ + mp_limb_t limb; + /* The number of bytes already copied to this limb (starting from + the low end). */ + size_t bytes; + /* The index where the limb should be stored, when completed. */ + mp_size_t i; + + if (nails != 0) + gmp_die ("mpz_import: Nails not supported."); + + assert (order == 1 || order == -1); + assert (endian >= -1 && endian <= 1); + + if (endian == 0) + endian = gmp_detect_endian (); + + p = (unsigned char *) src; + + word_step = (order != endian) ? 2 * size : 0; + + /* Process bytes from the least significant end, so point p at the + least significant word. */ + if (order == 1) + { + p += size * (count - 1); + word_step = - word_step; + } + + /* And at least significant byte of that word. */ + if (endian == 1) + p += (size - 1); + + rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t); + rp = MPZ_REALLOC (r, rn); + + for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step) + { + size_t j; + for (j = 0; j < size; j++, p -= (ptrdiff_t) endian) + { + limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT); + if (bytes == sizeof(mp_limb_t)) + { + rp[i++] = limb; + bytes = 0; + limb = 0; + } + } + } + assert (i + (bytes > 0) == rn); + if (limb != 0) + rp[i++] = limb; + else + i = mpn_normalized_size (rp, i); + + r->_mp_size = i; +} + +void * +mpz_export (void *r, size_t *countp, int order, size_t size, int endian, + size_t nails, const mpz_t u) +{ + size_t count; + mp_size_t un; + + if (nails != 0) + gmp_die ("mpz_import: Nails not supported."); + + assert (order == 1 || order == -1); + assert (endian >= -1 && endian <= 1); + assert (size > 0 || u->_mp_size == 0); + + un = u->_mp_size; + count = 0; + if (un != 0) + { + size_t k; + unsigned char *p; + ptrdiff_t word_step; + /* The current (partial) limb. */ + mp_limb_t limb; + /* The number of bytes left to to in this limb. */ + size_t bytes; + /* The index where the limb was read. */ + mp_size_t i; + + un = GMP_ABS (un); + + /* Count bytes in top limb. */ + limb = u->_mp_d[un-1]; + assert (limb != 0); + + k = 0; + do { + k++; limb >>= CHAR_BIT; + } while (limb != 0); + + count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size; + + if (!r) + r = gmp_xalloc (count * size); + + if (endian == 0) + endian = gmp_detect_endian (); + + p = (unsigned char *) r; + + word_step = (order != endian) ? 2 * size : 0; + + /* Process bytes from the least significant end, so point p at the + least significant word. */ + if (order == 1) + { + p += size * (count - 1); + word_step = - word_step; + } + + /* And at least significant byte of that word. */ + if (endian == 1) + p += (size - 1); + + for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step) + { + size_t j; + for (j = 0; j < size; j++, p -= (ptrdiff_t) endian) + { + if (bytes == 0) + { + if (i < un) + limb = u->_mp_d[i++]; + bytes = sizeof (mp_limb_t); + } + *p = limb; + limb >>= CHAR_BIT; + bytes--; + } + } + assert (i == un); + assert (k == count); + } + + if (countp) + *countp = count; + + return r; +} diff --git a/src/mini-gmp/mini-gmp.h b/src/mini-gmp/mini-gmp.h new file mode 100644 index 0000000..18d94ed --- /dev/null +++ b/src/mini-gmp/mini-gmp.h @@ -0,0 +1,298 @@ +/* mini-gmp, a minimalistic implementation of a GNU GMP subset. + +Copyright 2011-2015 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +/* About mini-gmp: This is a minimal implementation of a subset of the + GMP interface. It is intended for inclusion into applications which + have modest bignums needs, as a fallback when the real GMP library + is not installed. + + This file defines the public interface. */ + +#ifndef MINI_GMP_H +#define MINI_GMP_H + +/* For size_t */ +#include + +#if defined (__cplusplus) +extern "C" { +#endif + +void mp_set_memory_functions (void *(*) (size_t), + void *(*) (void *, size_t, size_t), + void (*) (void *, size_t)); + +void mp_get_memory_functions (void *(**) (size_t), + void *(**) (void *, size_t, size_t), + void (**) (void *, size_t)); + +typedef unsigned long mp_limb_t; +typedef long mp_size_t; +typedef unsigned long mp_bitcnt_t; + +typedef mp_limb_t *mp_ptr; +typedef const mp_limb_t *mp_srcptr; + +typedef struct +{ + int _mp_alloc; /* Number of *limbs* allocated and pointed + to by the _mp_d field. */ + int _mp_size; /* abs(_mp_size) is the number of limbs the + last field points to. If _mp_size is + negative this is a negative number. */ + mp_limb_t *_mp_d; /* Pointer to the limbs. */ +} mpz_struct; + +typedef mpz_struct mpz_t[1]; + +typedef mpz_struct *mpz_ptr; +typedef const mpz_struct *mpz_srcptr; + +extern const int mp_bits_per_limb; + +void mpn_copyi (mp_ptr, mp_srcptr, mp_size_t); +void mpn_copyd (mp_ptr, mp_srcptr, mp_size_t); +void mpn_zero (mp_ptr, mp_size_t); + +int mpn_cmp (mp_srcptr, mp_srcptr, mp_size_t); +int mpn_zero_p (mp_srcptr, mp_size_t); + +mp_limb_t mpn_add_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); +mp_limb_t mpn_add_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t); +mp_limb_t mpn_add (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t); + +mp_limb_t mpn_sub_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); +mp_limb_t mpn_sub_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t); +mp_limb_t mpn_sub (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t); + +mp_limb_t mpn_mul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); +mp_limb_t mpn_addmul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); +mp_limb_t mpn_submul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t); + +mp_limb_t mpn_mul (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t); +void mpn_mul_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t); +void mpn_sqr (mp_ptr, mp_srcptr, mp_size_t); +int mpn_perfect_square_p (mp_srcptr, mp_size_t); +mp_size_t mpn_sqrtrem (mp_ptr, mp_ptr, mp_srcptr, mp_size_t); + +mp_limb_t mpn_lshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int); +mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int); + +mp_bitcnt_t mpn_scan0 (mp_srcptr, mp_bitcnt_t); +mp_bitcnt_t mpn_scan1 (mp_srcptr, mp_bitcnt_t); + +void mpn_com (mp_ptr, mp_srcptr, mp_size_t); +mp_limb_t mpn_neg (mp_ptr, mp_srcptr, mp_size_t); + +mp_bitcnt_t mpn_popcount (mp_srcptr, mp_size_t); + +mp_limb_t mpn_invert_3by2 (mp_limb_t, mp_limb_t); +#define mpn_invert_limb(x) mpn_invert_3by2 ((x), 0) + +size_t mpn_get_str (unsigned char *, int, mp_ptr, mp_size_t); +mp_size_t mpn_set_str (mp_ptr, const unsigned char *, size_t, int); + +void mpz_init (mpz_t); +void mpz_init2 (mpz_t, mp_bitcnt_t); +void mpz_clear (mpz_t); + +#define mpz_odd_p(z) (((z)->_mp_size != 0) & (int) (z)->_mp_d[0]) +#define mpz_even_p(z) (! mpz_odd_p (z)) + +int mpz_sgn (const mpz_t); +int mpz_cmp_si (const mpz_t, long); +int mpz_cmp_ui (const mpz_t, unsigned long); +int mpz_cmp (const mpz_t, const mpz_t); +int mpz_cmpabs_ui (const mpz_t, unsigned long); +int mpz_cmpabs (const mpz_t, const mpz_t); +int mpz_cmp_d (const mpz_t, double); +int mpz_cmpabs_d (const mpz_t, double); + +void mpz_abs (mpz_t, const mpz_t); +void mpz_neg (mpz_t, const mpz_t); +void mpz_swap (mpz_t, mpz_t); + +void mpz_add_ui (mpz_t, const mpz_t, unsigned long); +void mpz_add (mpz_t, const mpz_t, const mpz_t); +void mpz_sub_ui (mpz_t, const mpz_t, unsigned long); +void mpz_ui_sub (mpz_t, unsigned long, const mpz_t); +void mpz_sub (mpz_t, const mpz_t, const mpz_t); + +void mpz_mul_si (mpz_t, const mpz_t, long int); +void mpz_mul_ui (mpz_t, const mpz_t, unsigned long int); +void mpz_mul (mpz_t, const mpz_t, const mpz_t); +void mpz_mul_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_addmul_ui (mpz_t, const mpz_t, unsigned long int); +void mpz_addmul (mpz_t, const mpz_t, const mpz_t); +void mpz_submul_ui (mpz_t, const mpz_t, unsigned long int); +void mpz_submul (mpz_t, const mpz_t, const mpz_t); + +void mpz_cdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_fdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_tdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_cdiv_q (mpz_t, const mpz_t, const mpz_t); +void mpz_fdiv_q (mpz_t, const mpz_t, const mpz_t); +void mpz_tdiv_q (mpz_t, const mpz_t, const mpz_t); +void mpz_cdiv_r (mpz_t, const mpz_t, const mpz_t); +void mpz_fdiv_r (mpz_t, const mpz_t, const mpz_t); +void mpz_tdiv_r (mpz_t, const mpz_t, const mpz_t); + +void mpz_cdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_fdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_tdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_cdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_fdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t); +void mpz_tdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t); + +void mpz_mod (mpz_t, const mpz_t, const mpz_t); + +void mpz_divexact (mpz_t, const mpz_t, const mpz_t); + +int mpz_divisible_p (const mpz_t, const mpz_t); +int mpz_congruent_p (const mpz_t, const mpz_t, const mpz_t); + +unsigned long mpz_cdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long); +unsigned long mpz_fdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long); +unsigned long mpz_tdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long); +unsigned long mpz_cdiv_q_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_fdiv_q_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_tdiv_q_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_cdiv_r_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_fdiv_r_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_tdiv_r_ui (mpz_t, const mpz_t, unsigned long); +unsigned long mpz_cdiv_ui (const mpz_t, unsigned long); +unsigned long mpz_fdiv_ui (const mpz_t, unsigned long); +unsigned long mpz_tdiv_ui (const mpz_t, unsigned long); + +unsigned long mpz_mod_ui (mpz_t, const mpz_t, unsigned long); + +void mpz_divexact_ui (mpz_t, const mpz_t, unsigned long); + +int mpz_divisible_ui_p (const mpz_t, unsigned long); + +unsigned long mpz_gcd_ui (mpz_t, const mpz_t, unsigned long); +void mpz_gcd (mpz_t, const mpz_t, const mpz_t); +void mpz_gcdext (mpz_t, mpz_t, mpz_t, const mpz_t, const mpz_t); +void mpz_lcm_ui (mpz_t, const mpz_t, unsigned long); +void mpz_lcm (mpz_t, const mpz_t, const mpz_t); +int mpz_invert (mpz_t, const mpz_t, const mpz_t); + +void mpz_sqrtrem (mpz_t, mpz_t, const mpz_t); +void mpz_sqrt (mpz_t, const mpz_t); +int mpz_perfect_square_p (const mpz_t); + +void mpz_pow_ui (mpz_t, const mpz_t, unsigned long); +void mpz_ui_pow_ui (mpz_t, unsigned long, unsigned long); +void mpz_powm (mpz_t, const mpz_t, const mpz_t, const mpz_t); +void mpz_powm_ui (mpz_t, const mpz_t, unsigned long, const mpz_t); + +void mpz_rootrem (mpz_t, mpz_t, const mpz_t, unsigned long); +int mpz_root (mpz_t, const mpz_t, unsigned long); + +void mpz_fac_ui (mpz_t, unsigned long); +void mpz_bin_uiui (mpz_t, unsigned long, unsigned long); + +int mpz_probab_prime_p (const mpz_t, int); + +int mpz_tstbit (const mpz_t, mp_bitcnt_t); +void mpz_setbit (mpz_t, mp_bitcnt_t); +void mpz_clrbit (mpz_t, mp_bitcnt_t); +void mpz_combit (mpz_t, mp_bitcnt_t); + +void mpz_com (mpz_t, const mpz_t); +void mpz_and (mpz_t, const mpz_t, const mpz_t); +void mpz_ior (mpz_t, const mpz_t, const mpz_t); +void mpz_xor (mpz_t, const mpz_t, const mpz_t); + +mp_bitcnt_t mpz_popcount (const mpz_t); +mp_bitcnt_t mpz_hamdist (const mpz_t, const mpz_t); +mp_bitcnt_t mpz_scan0 (const mpz_t, mp_bitcnt_t); +mp_bitcnt_t mpz_scan1 (const mpz_t, mp_bitcnt_t); + +int mpz_fits_slong_p (const mpz_t); +int mpz_fits_ulong_p (const mpz_t); +long int mpz_get_si (const mpz_t); +unsigned long int mpz_get_ui (const mpz_t); +double mpz_get_d (const mpz_t); +size_t mpz_size (const mpz_t); +mp_limb_t mpz_getlimbn (const mpz_t, mp_size_t); + +void mpz_realloc2 (mpz_t, mp_bitcnt_t); +mp_srcptr mpz_limbs_read (mpz_srcptr); +mp_ptr mpz_limbs_modify (mpz_t, mp_size_t); +mp_ptr mpz_limbs_write (mpz_t, mp_size_t); +void mpz_limbs_finish (mpz_t, mp_size_t); +mpz_srcptr mpz_roinit_n (mpz_t, mp_srcptr, mp_size_t); + +#define MPZ_ROINIT_N(xp, xs) {{0, (xs),(xp) }} + +void mpz_set_si (mpz_t, signed long int); +void mpz_set_ui (mpz_t, unsigned long int); +void mpz_set (mpz_t, const mpz_t); +void mpz_set_d (mpz_t, double); + +void mpz_init_set_si (mpz_t, signed long int); +void mpz_init_set_ui (mpz_t, unsigned long int); +void mpz_init_set (mpz_t, const mpz_t); +void mpz_init_set_d (mpz_t, double); + +size_t mpz_sizeinbase (const mpz_t, int); +char *mpz_get_str (char *, int, const mpz_t); +int mpz_set_str (mpz_t, const char *, int); +int mpz_init_set_str (mpz_t, const char *, int); + +/* This long list taken from gmp.h. */ +/* For reference, "defined(EOF)" cannot be used here. In g++ 2.95.4, + defines EOF but not FILE. */ +#if defined (FILE) \ + || defined (H_STDIO) \ + || defined (_H_STDIO) /* AIX */ \ + || defined (_STDIO_H) /* glibc, Sun, SCO */ \ + || defined (_STDIO_H_) /* BSD, OSF */ \ + || defined (__STDIO_H) /* Borland */ \ + || defined (__STDIO_H__) /* IRIX */ \ + || defined (_STDIO_INCLUDED) /* HPUX */ \ + || defined (__dj_include_stdio_h_) /* DJGPP */ \ + || defined (_FILE_DEFINED) /* Microsoft */ \ + || defined (__STDIO__) /* Apple MPW MrC */ \ + || defined (_MSL_STDIO_H) /* Metrowerks */ \ + || defined (_STDIO_H_INCLUDED) /* QNX4 */ \ + || defined (_ISO_STDIO_ISO_H) /* Sun C++ */ \ + || defined (__STDIO_LOADED) /* VMS */ +size_t mpz_out_str (FILE *, int, const mpz_t); +#endif + +void mpz_import (mpz_t, size_t, int, size_t, int, size_t, const void *); +void *mpz_export (void *, size_t *, int, size_t, int, size_t, const mpz_t); + +#if defined (__cplusplus) +} +#endif +#endif /* MINI_GMP_H */ diff --git a/src/ninty-233.c b/src/ninty-233.c new file mode 100644 index 0000000..ca6377d --- /dev/null +++ b/src/ninty-233.c @@ -0,0 +1,248 @@ +/* + ninty-233.c + + Copyright © 2018, 2019 Jbop (https://github.com/jbop1626) + + This file is a part of ninty-233. + + ninty-233 is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ninty-233 is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#include +#include +#include +#include +#include +#include + +#include "ninty-233.h" +#include "ecc/ecc.h" +#include "sha1/sha1.h" +#include "mini-gmp/mini-gmp.h" + + +static void init_mpz_list(size_t count, mpz_ptr x, ...) { + va_list mpz_list; + va_start(mpz_list, x); + + size_t i = 0; + while (i < count) { + mpz_init(x); + x = va_arg(mpz_list, mpz_ptr); + i++; + } + + va_end(mpz_list); +} + +static void clear_mpz_list(size_t count, mpz_ptr x, ...) { + va_list mpz_list; + va_start(mpz_list, x); + + size_t i = 0; + while (i < count) { + mpz_clear(x); + x = va_arg(mpz_list, mpz_ptr); + i++; + } + + va_end(mpz_list); +} + +static void generate_k(const mpz_t n, const mpz_t hash, mpz_t k_out) { + // Do NOT use this implementation for generation of k + // when creating a signature that must be secure! + srand(time(NULL)); + + mpz_t random_mpz; + mpz_init(random_mpz); + + uint32_t buffer[8] = { 0 }; + for(int i = 0; i < 8; ++i) { + buffer[i] = rand() % UINT32_MAX; + } + + mpz_import(random_mpz, 8, 1, sizeof(buffer[0]), 0, 0, buffer); + mpz_mul(k_out, random_mpz, hash); + while (mpz_cmp(k_out, n) >= 0) { + mpz_tdiv_q_ui(k_out, k_out, 7); + } + + mpz_clear(random_mpz); +} + + + +void mpz_to_gf2m(const mpz_t src, element dst) { + uint32_t buffer[32] = { 0 }; + gf2m_set_zero(dst); + + size_t count = 0; + mpz_export((void *)buffer, &count, 1, sizeof(dst[0]), 0, 0, src); + if (count == 0 || count > INT_MAX) { + fprintf(stderr, "mpz_to_gf2m error! Element argument is now zero.\n"); + return; + } + + int i = 7; + int j = count - 1; + while(i >= 0 && j >= 0) { + dst[i] = buffer[j]; + i--; + j--; + } +} + +void gf2m_to_mpz(const element src, mpz_t dst) { + mpz_import(dst, 8, 1, sizeof(src[0]), 0, 0, src); +} + +void sha1(const uint8_t * input, uint32_t input_length, unsigned ique_flag, mpz_t hash_out) { + SHA1_HASH hash; + Sha1Context context; + + Sha1Initialise(&context); + Sha1Update(&context, input, input_length); + if (ique_flag) { + // When performing certain hashes, the iQue Player updates the + // SHA1 state with the following magic data. + uint8_t ique_magic[4] = { 0x06, 0x09, 0x19, 0x68 }; + Sha1Update(&context, &ique_magic, 4); + } + Sha1Finalise(&context, &hash); + + mpz_import(hash_out, 20, 1, sizeof(hash.bytes[0]), 0, 0, (void *)hash.bytes); +} + +void ecdh(const element private_key, const ec_point * public_key, ec_point * shared_secret_output) { + ec_point_mul(private_key, public_key, shared_secret_output); +} + +void ecdsa_sign(const mpz_t z, const element private_key, element r_out, element s_out) { + mpz_t r, s, n, D, zero, k, x_p, k_inv, med; + init_mpz_list(9, r, s, n, D, zero, k, x_p, k_inv, med); + + gf2m_to_mpz(G_ORDER, n); + gf2m_to_mpz(private_key, D); + gf2m_set_zero(r_out); + gf2m_set_zero(s_out); + + while(!mpz_cmp(r, zero) || !mpz_cmp(s, zero)) { + // Generate k in [1, n - 1] + generate_k(n, z, k); + element k_elem; + mpz_to_gf2m(k, k_elem); + + // Calculate P = kG + ec_point G, P; + gf2m_copy(G_X, G.x); + gf2m_copy(G_Y, G.y); + ec_point_mul(k_elem, &G, &P); + + // Calculate r = x_p mod n + gf2m_to_mpz(P.x, x_p); + mpz_mod(r, x_p, n); + + // Calculate s = k^-1(z + rD) mod n + if (mpz_invert(k_inv, k, n) == 0) { + fprintf(stderr, "An error occurred while calculating the inverse of k mod n.\n"); + fprintf(stderr, "The resulting signature will be invalid!\n"); + } + mpz_mul(med, r, D); + mpz_add(med, z, med); + mpz_mod(med, med, n); + + mpz_mul(s, k_inv, med); + mpz_mod(s, s, n); + } + mpz_to_gf2m(r, r_out); + mpz_to_gf2m(s, s_out); + + clear_mpz_list(9, r, s, n, D, zero, k, x_p, k_inv, med); +} + +int ecdsa_verify(const mpz_t z, const ec_point * public_key, const element r_input, const element s_input) { + ec_point Q, test; + ec_point_copy(public_key, &Q); + element zero = { 0 }; + + // If Q is the identity, Q is invalid + if (gf2m_is_equal(Q.x, zero) && gf2m_is_equal(Q.y, zero)) { + return 0; + } + // If Q is not a point on the curve, Q is invalid + if (!ec_point_on_curve(&Q)) { + return 0; + } + // If nQ is not the identity, Q is invalid (or n is messed up) + ec_point_mul(G_ORDER, &Q, &test); + if (!(gf2m_is_equal(test.x, zero) && gf2m_is_equal(test.y, zero))) { + return 0; + } + + // Public key is valid, now verify signature... + mpz_t r, s, n; + init_mpz_list(3, r, s, n); + gf2m_to_mpz(r_input, r); + gf2m_to_mpz(s_input, s); + gf2m_to_mpz(G_ORDER, n); + + // If r or s are not in [1, n - 1], sig is invalid + if ( (mpz_cmp_ui(r, 1) < 0 || mpz_cmp(r, n) > 0 || mpz_cmp(r, n) == 0) || + (mpz_cmp_ui(s, 1) < 0 || mpz_cmp(s, n) > 0 || mpz_cmp(s, n) == 0) ) { + clear_mpz_list(3, r, s, n); + return 0; + } + + // Calculate u_1 and u_2 + mpz_t s_inv, u_1, u_2; + init_mpz_list(3, s_inv, u_1, u_2); + + if (mpz_invert(s_inv, s, n) == 0) { + fprintf(stderr, "An error occurred while calculating the inverse of s mod n.\n"); + clear_mpz_list(6, r, s, n, s_inv, u_1, u_2); + return 0; + } + mpz_mul(u_1, z, s_inv); + mpz_mod(u_1, u_1, n); + mpz_mul(u_2, r, s_inv); + mpz_mod(u_2, u_2, n); + + // Calculate P3 = u_1G + u_2Q + element u_1_elem, u_2_elem; + mpz_to_gf2m(u_1, u_1_elem); + mpz_to_gf2m(u_2, u_2_elem); + ec_point G, P1, P2, P3; + gf2m_copy(G_X, G.x); + gf2m_copy(G_Y, G.y); + + ec_point_mul(u_1_elem, &G, &P1); + ec_point_mul(u_2_elem, &Q, &P2); + ec_point_add(&P1, &P2, &P3); + + // If P3 is the identity, sig is invalid + if (gf2m_is_equal(P3.x, zero) && gf2m_is_equal(P3.y, zero)) { + clear_mpz_list(6, r, s, n, s_inv, u_1, u_2); + return 0; + } + + // And finally, is r congruent to P3.x mod n? + mpz_t x_p; + mpz_init(x_p); + gf2m_to_mpz(P3.x, x_p); + + int is_congruent = mpz_congruent_p(r, x_p, n) != 0; + clear_mpz_list(7, r, s, n, s_inv, u_1, u_2, x_p); + return is_congruent; +} \ No newline at end of file diff --git a/src/ninty-233.cpp b/src/ninty-233.cpp deleted file mode 100644 index 9495e07..0000000 --- a/src/ninty-233.cpp +++ /dev/null @@ -1,204 +0,0 @@ -/* - ninty-233.cpp - - Copyright © 2018 Jbop (https://github.com/jbop1626); - - This file is a part of ninty-233. - - ninty-233 is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - ninty-233 is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . -*/ -#include "ninty-233.hpp" -#include -#include -#include - -/* - BigUnsigned <--> GF(2^m) element conversions -*/ -BigUnsigned gf2m_to_bigunsigned(const element src) { - BigUnsigned dst = src[0]; - for (int i = 1; i < 8; ++i) { - dst <<= 32; - dst += src[i]; - } - return dst; -} - -void bigunsigned_to_gf2m(const BigUnsigned & src, element dst) { - gf2m_set_zero(dst); - BigUnsigned temp = src; - for (int i = 7; i >= 0; --i) { - BigUnsigned low32 = temp & 0xFFFFFFFF; - dst[i] = low32.toUnsignedInt(); - temp >>= 32; - } -} - -/* - SHA-1 result as big (unsigned) integer -*/ -BigUnsigned sha1(uint8_t * input, int length) { - SHA1_HASH hash; - Sha1Context context; - - Sha1Initialise(&context); - Sha1Update(&context, input, length); -#if defined(IQUE_ECC) && (IQUE_ECC == 1) - uint8_t ique_magic[4] = { 0x06, 0x09, 0x19, 0x68 }; // iQue-specific magic - Sha1Update(&context, &ique_magic, 4); -#endif - Sha1Finalise(&context, &hash); - - BigUnsigned hash_bigint = hash.bytes[0]; - for (int i = 1; i < 20; ++i) { - hash_bigint <<= 8; - hash_bigint += hash.bytes[i]; - } - - return hash_bigint; -} - -/* - Generation of k, for signing -*/ -BigUnsigned generate_k(const BigUnsigned & n, const BigUnsigned & hash) { - BigUnsigned bu_temp = hash; - element elem_temp; - uint8_t os_temp[32]; - - std::random_device rd; - std::uniform_int_distribution dist(1, std::numeric_limits::max()); - unsigned int salt = dist(rd); - salt = salt * dist(rd); - bu_temp += salt; - - bigunsigned_to_gf2m(bu_temp, elem_temp); - elem_to_os(elem_temp, os_temp); - BigUnsigned k = sha1(os_temp, 32); - k *= k; - while(k >= n) { - k /= 7; - } - return k; -} - -/* - ECC algorithms -*/ -void ecdh(const uint8_t * private_key, const uint8_t * public_key, uint8_t * output) { - element private_copy; - ec_point public_copy; - ec_point shared_secret; - - os_to_elem(private_key, private_copy); - os_to_point(public_key, public_copy); - - ec_point_mul(private_copy, public_copy, shared_secret); - - point_to_os(shared_secret, output); -} - -void ecdsa_sign(const BigUnsigned z, const uint8_t * private_key, element r_out, element s_out) { - element private_copy; - os_to_elem(private_key, private_copy); - - BigUnsigned r = 0; - BigUnsigned s = 0; - BigUnsigned n = gf2m_to_bigunsigned(G_order); - BigUnsigned D = gf2m_to_bigunsigned(private_copy); - - while(r == 0 || s == 0) { - // Generate k in [1, n - 1] - BigUnsigned k = generate_k(n, z); - element k_elem; - bigunsigned_to_gf2m(k, k_elem); - - // Calculate P = kG - ec_point G; - gf2m_copy(G_x, G.x); - gf2m_copy(G_y, G.y); - ec_point P; - ec_point_mul(k_elem, G, P); - - // Calculate r = x_p mod n - BigUnsigned x_p = gf2m_to_bigunsigned(P.x); - r = x_p % n; - - // Calculate s = k^-1(z + rD) mod n - BigUnsigned k_inv = modinv(k, n); - BigUnsigned med = (z + (r * D)) % n; - s = (k_inv * med) % n; - } - bigunsigned_to_gf2m(r, r_out); - bigunsigned_to_gf2m(s, s_out); -} - -bool ecdsa_verify(const BigUnsigned z, const uint8_t * public_key, const element r_input, const element s_input) { - ec_point Q, test; - os_to_point(public_key, Q); - element zero = { 0 }; - - // If Q is the identity, Q is invalid - if (gf2m_is_equal(Q.x, zero) && gf2m_is_equal(Q.y, zero)) { - return false; - } - // If Q is not a point on the curve, Q is invalid - if (!ec_point_on_curve(Q)) { - return false; - } - // If nQ is not the identity, Q is invalid (or n is messed up) - ec_point_mul(G_order, Q, test); - if (!(gf2m_is_equal(test.x, zero) && gf2m_is_equal(test.y, zero))) { - return false; - } - // Public key is valid, now verify signature... - BigUnsigned r = gf2m_to_bigunsigned(r_input); - BigUnsigned s = gf2m_to_bigunsigned(s_input); - BigUnsigned n = gf2m_to_bigunsigned(G_order); - - // If r,s are not in [1, n - 1], sig is invalid - if (r < 1 || r >= n) { - return false; - } - if (s < 1 || s >= n) { - return false; - } - - // Calculate u_1 and u_2 - BigUnsigned s_inv = modinv(s, n); - BigUnsigned u_1 = (z * s_inv) % n; - BigUnsigned u_2 = (r * s_inv) % n; - - // Calculate P3 = u_1G + u_2Q - element u_1_elem, u_2_elem; - ec_point P1, P2, P3; - ec_point G; - gf2m_copy(G_x, G.x); - gf2m_copy(G_y, G.y); - bigunsigned_to_gf2m(u_1, u_1_elem); - bigunsigned_to_gf2m(u_2, u_2_elem); - - ec_point_mul(u_1_elem, G, P1); - ec_point_mul(u_2_elem, Q, P2); - ec_point_add(P1, P2, P3); - - // If P3 is the identity, sig is invalid - if (gf2m_is_equal(P3.x, zero) && gf2m_is_equal(P3.y, zero)) { - return false; - } - - // And finally, is r congruent to P3.x mod n? - BigUnsigned x_p = gf2m_to_bigunsigned(P3.x); - return r == (x_p % n); -} diff --git a/src/ninty-233.h b/src/ninty-233.h new file mode 100644 index 0000000..7733818 --- /dev/null +++ b/src/ninty-233.h @@ -0,0 +1,59 @@ +/* + ninty-233 + Library for ECC operations using keys defined with + sect233r1 / NIST B-233 -- the curve/domain parameters + used by Nintendo in the iQue Player and Wii. + + Copyright © 2018, 2019 Jbop (https://github.com/jbop1626) + + This file is a part of ninty-233. + + ninty-233 is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ninty-233 is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ +#ifndef NINTY_233_H +#define NINTY_233_H + +#include + +#include "ecc/ecc.h" +#include "mini-gmp/mini-gmp.h" + +#if defined (__cplusplus) +extern "C" { +#endif +/* + Multi-precision integer <--> GF(2^m) element conversions +*/ +void mpz_to_gf2m(const mpz_t src, element dst); +void gf2m_to_mpz(const element src, mpz_t dst); + +/* + SHA-1 result as multi-precision integer +*/ +#define NOT_IQUE_HASH 0 +#define IQUE_HASH 1 +void sha1(const uint8_t * input, uint32_t input_length, unsigned ique_flag, mpz_t hash_out); + +/* + ECC algorithms +*/ +void ecdh(const element private_key, const ec_point * public_key, ec_point * shared_secret_output); +void ecdsa_sign(const mpz_t hash, const element private_key, element r_out, element s_out); +int ecdsa_verify(const mpz_t hash, const ec_point * public_key, const element r_input, const element s_input); + +#if defined (__cplusplus) +} +#endif + +#endif diff --git a/src/ninty-233.hpp b/src/ninty-233.hpp deleted file mode 100644 index 329952b..0000000 --- a/src/ninty-233.hpp +++ /dev/null @@ -1,56 +0,0 @@ -/* - ninty-233 - Library for ECC operations using keys defined with - sect233r1 / NIST B-233 -- the curve/domain parameters - used by Nintendo in the iQue Player and Wii. - - Copyright © 2018 Jbop (https://github.com/jbop1626); - - This file is a part of ninty-233. - - ninty-233 is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - ninty-233 is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . -*/ -#ifndef NINTY_233 -#define NINTY_233 - -#include "bigint/include/BigIntegerLibrary.hpp" -#include "ecc/ecc.hpp" -#include "sha1/sha1.hpp" - -#define IQUE_ECC 1 - -/* - BigUnsigned <--> GF(2^m) element conversions -*/ -BigUnsigned gf2m_to_bigunsigned(const element src); -void bigunsigned_to_gf2m(const BigUnsigned & src, element dst); - -/* - SHA-1 result as big (unsigned) integer -*/ -BigUnsigned sha1(uint8_t * input, int length); - -/* - Generation of k, for signing -*/ -BigUnsigned generate_k(const BigUnsigned & n, const BigUnsigned & hash); - -/* - ECC algorithms -*/ -void ecdh(const uint8_t * private_key, const uint8_t * public_key, uint8_t * output); -void ecdsa_sign(const BigUnsigned hash, const uint8_t * private_key, element r_out, element s_out); -bool ecdsa_verify(const BigUnsigned hash, const uint8_t * public_key, const element r_input, const element s_input); - -#endif diff --git a/src/sha1/sha1.h b/src/sha1/sha1.h index 81b3a85..0985464 100644 --- a/src/sha1/sha1.h +++ b/src/sha1/sha1.h @@ -10,7 +10,8 @@ // This is free and unencumbered software released into the public domain - June 2013 waterjuice.org //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// -#pragma once +#ifndef NINTY_233_WJCRYPT_SHA1_H +#define NINTY_233_WJCRYPT_SHA1_H //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // IMPORTS @@ -22,6 +23,9 @@ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // TYPES //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// +#if defined (__cplusplus) +extern "C" { +#endif // Sha1Context - This must be initialised using Sha1Initialised. Do not modify the contents of this structure directly. typedef struct @@ -92,3 +96,9 @@ void uint32_t BufferSize, // [in] SHA1_HASH* Digest // [in] ); + +#if defined (__cplusplus) +} +#endif + +#endif diff --git a/src/sha1/sha1.hpp b/src/sha1/sha1.hpp deleted file mode 100644 index abeccee..0000000 --- a/src/sha1/sha1.hpp +++ /dev/null @@ -1,32 +0,0 @@ -/* - Copyright © 2018 Jbop (https://github.com/jbop1626); - - This file is a part of ninty-233. - - ninty-233 is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - ninty-233 is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see . -*/ - -#ifndef NINTY_233_SHA1_HPP -#define NINTY_233_SHA1_HPP - -#ifndef __cplusplus -#error Do not include this header in a C project -#endif - -extern "C" { -#include "sha1.h" -} - -#endif -