From 709fe34e8e1ebcd97413aca46148250d3132ebb1 Mon Sep 17 00:00:00 2001 From: mobileflowllc Date: Thu, 29 Jun 2017 17:07:43 -0700 Subject: [PATCH] Add SHA512 compliant SRP client --- Homekit.xcodeproj/project.pbxproj | 28 + src/SRPClient.cpp | 191 +- src/SRPClient.h | 35 + src/md32_common.h | 428 +++ src/mini-gmp.c | 4381 +++++++++++++++++++++++++++++ src/mini-gmp.h | 298 ++ src/sha.h | 175 ++ src/sha1.c | 537 ++++ src/sha256.c | 399 +++ src/sha512.c | 408 +++ src/srp.c | 1047 +++++++ src/srp.h | 201 ++ test.cpp | 38 +- 13 files changed, 8161 insertions(+), 5 deletions(-) create mode 100644 src/md32_common.h create mode 100644 src/mini-gmp.c create mode 100644 src/mini-gmp.h create mode 100644 src/sha.h create mode 100644 src/sha1.c create mode 100644 src/sha256.c create mode 100644 src/sha512.c create mode 100644 src/srp.c create mode 100644 src/srp.h diff --git a/Homekit.xcodeproj/project.pbxproj b/Homekit.xcodeproj/project.pbxproj index d37992e..cc9175b 100644 --- a/Homekit.xcodeproj/project.pbxproj +++ b/Homekit.xcodeproj/project.pbxproj @@ -7,6 +7,11 @@ objects = { /* Begin PBXBuildFile section */ + 913ECAA21F05CD6000E910C5 /* srp.c in Sources */ = {isa = PBXBuildFile; fileRef = 913ECAA01F05CD6000E910C5 /* srp.c */; }; + 913ECAA71F05CD7000E910C5 /* sha1.c in Sources */ = {isa = PBXBuildFile; fileRef = 913ECAA41F05CD7000E910C5 /* sha1.c */; }; + 913ECAA81F05CD7000E910C5 /* sha256.c in Sources */ = {isa = PBXBuildFile; fileRef = 913ECAA51F05CD7000E910C5 /* sha256.c */; }; + 913ECAA91F05CD7000E910C5 /* sha512.c in Sources */ = {isa = PBXBuildFile; fileRef = 913ECAA61F05CD7000E910C5 /* sha512.c */; }; + 913ECAAC1F05CD7900E910C5 /* mini-gmp.c in Sources */ = {isa = PBXBuildFile; fileRef = 913ECAAA1F05CD7900E910C5 /* mini-gmp.c */; }; 91BBDB181F02C5A2009CBF0A /* test.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 91BBDB171F02C5A2009CBF0A /* test.cpp */; }; 91BBDB1C1F02C5B1009CBF0A /* homekit.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 91BBDB1A1F02C5B1009CBF0A /* homekit.cpp */; }; 91F870F21F0345BC00F77259 /* WebClient.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 91F870F01F0345BC00F77259 /* WebClient.cpp */; }; @@ -27,6 +32,15 @@ /* End PBXCopyFilesBuildPhase section */ /* Begin PBXFileReference section */ + 913ECAA01F05CD6000E910C5 /* srp.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = srp.c; sourceTree = ""; }; + 913ECAA11F05CD6000E910C5 /* srp.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = srp.h; sourceTree = ""; }; + 913ECAA31F05CD7000E910C5 /* sha.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = sha.h; sourceTree = ""; }; + 913ECAA41F05CD7000E910C5 /* sha1.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = sha1.c; sourceTree = ""; }; + 913ECAA51F05CD7000E910C5 /* sha256.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = sha256.c; sourceTree = ""; }; + 913ECAA61F05CD7000E910C5 /* sha512.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = sha512.c; sourceTree = ""; }; + 913ECAAA1F05CD7900E910C5 /* mini-gmp.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = "mini-gmp.c"; sourceTree = ""; }; + 913ECAAB1F05CD7900E910C5 /* mini-gmp.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = "mini-gmp.h"; sourceTree = ""; }; + 913ECAAD1F05CD8C00E910C5 /* md32_common.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = md32_common.h; sourceTree = ""; }; 9163EA621F02C5690098C4B8 /* Homekit */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = Homekit; sourceTree = BUILT_PRODUCTS_DIR; }; 91BBDB171F02C5A2009CBF0A /* test.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = test.cpp; sourceTree = SOURCE_ROOT; }; 91BBDB1A1F02C5B1009CBF0A /* homekit.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = homekit.cpp; sourceTree = ""; }; @@ -98,6 +112,15 @@ 91F870F51F04205E00F77259 /* TLV8.h */, 91F870F71F045B9700F77259 /* SRPClient.cpp */, 91F870F81F045B9700F77259 /* SRPClient.h */, + 913ECAA01F05CD6000E910C5 /* srp.c */, + 913ECAA11F05CD6000E910C5 /* srp.h */, + 913ECAA31F05CD7000E910C5 /* sha.h */, + 913ECAAD1F05CD8C00E910C5 /* md32_common.h */, + 913ECAA41F05CD7000E910C5 /* sha1.c */, + 913ECAA51F05CD7000E910C5 /* sha256.c */, + 913ECAA61F05CD7000E910C5 /* sha512.c */, + 913ECAAA1F05CD7900E910C5 /* mini-gmp.c */, + 913ECAAB1F05CD7900E910C5 /* mini-gmp.h */, ); path = src; sourceTree = ""; @@ -162,7 +185,12 @@ files = ( 91BBDB181F02C5A2009CBF0A /* test.cpp in Sources */, 91BBDB1C1F02C5B1009CBF0A /* homekit.cpp in Sources */, + 913ECAA81F05CD7000E910C5 /* sha256.c in Sources */, 91F870F91F045B9700F77259 /* SRPClient.cpp in Sources */, + 913ECAA91F05CD7000E910C5 /* sha512.c in Sources */, + 913ECAAC1F05CD7900E910C5 /* mini-gmp.c in Sources */, + 913ECAA21F05CD6000E910C5 /* srp.c in Sources */, + 913ECAA71F05CD7000E910C5 /* sha1.c in Sources */, 91F870F21F0345BC00F77259 /* WebClient.cpp in Sources */, 91F870F61F04205E00F77259 /* TLV8.cpp in Sources */, ); diff --git a/src/SRPClient.cpp b/src/SRPClient.cpp index 130d334..89851ca 100644 --- a/src/SRPClient.cpp +++ b/src/SRPClient.cpp @@ -17,17 +17,202 @@ SRPClient::~SRPClient() { } +int SRPClient::getChallenge(uint8_t **salt,uint16_t *salt_len, uint8_t **key, uint16_t *key_len) +{ + + + size_t len_s = 16; + size_t len_v = 0; + + const char *username = "alice"; + const char *password = "password123"; + + SRP_HashAlgorithm alg = SRP_SHA512; + SRP_NGType ng_type = SRP_NG_3072; // was SRP_NG_1024; // TEST_NG; + + // The test vectors from + // https://tools.ietf.org/html/rfc5054#appendix-B + + static const uint8_t srp_5054_salt[] = { + 0xBE, 0xB2, 0x53, 0x79, 0xD1, 0xA8, 0x58, 0x1E, 0xB5, 0xA7, 0x27, 0x67, 0x3A, 0x24, + 0x41, 0xEE + }; + + static const uint8_t srp_3072_salt[] = { + 0xBE, 0xB2, 0x53, 0x79, 0xD1, 0xA8, 0x58, 0x1E, 0xB5, 0xA7, 0x27, 0x67, 0x3A, 0x24, + 0x41, 0xEE + }; + + unsigned char *bytes_s = 0; + bytes_s = (unsigned char *)malloc(sizeof(srp_3072_salt)); + memcpy(bytes_s, srp_3072_salt, sizeof(srp_3072_salt)); + + unsigned char *bytes_v = 0; + + srp_create_salted_verification_key(alg, + ng_type, + username, + (const unsigned char *)password, + strlen(password), + &bytes_s, + &len_s, + &bytes_v, + &len_v, + NULL, + NULL); + + *salt = bytes_s; + *salt_len = sizeof(srp_3072_salt); + *key = bytes_v; + *key_len = len_v; + + return 0; +} + +int SRPClient::verifySession(uint8_t **serverKeyProof, uint16_t *proof_len, uint8_t *clientPublicKey, uint16_t client_key_len, uint8_t *clientKeyProof, uint16_t client_proof_len) +{ + struct SRPVerifier *ver; + struct SRPUser *usr; + + size_t len_s = 16; + size_t len_v = client_proof_len; + size_t len_A = 0; + size_t len_B = 0; + + const char *username = "alice"; + const char *password = "password123"; + + SRP_HashAlgorithm alg = SRP_SHA512; + SRP_NGType ng_type = SRP_NG_3072; // was SRP_NG_1024; // TEST_NG; + + // The test vectors from + // https://tools.ietf.org/html/rfc5054#appendix-B + + static const uint8_t srp_5054_salt[] = { + 0xBE, 0xB2, 0x53, 0x79, 0xD1, 0xA8, 0x58, 0x1E, 0xB5, 0xA7, 0x27, 0x67, 0x3A, 0x24, + 0x41, 0xEE + }; + + static const uint8_t srp_3072_salt[] = { + 0xBE, 0xB2, 0x53, 0x79, 0xD1, 0xA8, 0x58, 0x1E, 0xB5, 0xA7, 0x27, 0x67, 0x3A, 0x24, + 0x41, 0xEE + }; + + static const uint8_t srp_3072_a[] = { + 0x60, 0x97, 0x55, 0x27, 0x03, 0x5C, 0xF2, 0xAD, 0x19, 0x89, 0x80, 0x6F, 0x04, 0x07, + 0x21, 0x0B, 0xC8, 0x1E, 0xDC, 0x04, 0xE2, 0x76, 0x2A, 0x56, 0xAF, 0xD5, 0x29, 0xDD, + 0xDA, 0x2D, 0x43, 0x93 + }; + + static const uint8_t srp_5054_b[] = { + 0xE4, 0x87, 0xCB, 0x59, 0xD3, 0x1A, 0xC5, 0x50, 0x47, 0x1E, 0x81, 0xF0, 0x0F, 0x69, + 0x28, 0xE0, 0x1D, 0xDA, 0x08, 0xE9, 0x74, 0xA0, 0x04, 0xF4, 0x9E, 0x61, 0xF5, 0xD1, + 0x05, 0x28, 0x4D, 0x20, + }; + + + + unsigned char *bytes_s = 0; + bytes_s = (unsigned char *)malloc(sizeof(srp_3072_salt)); + memcpy(bytes_s, srp_3072_salt, sizeof(srp_3072_salt)); + + unsigned char *bytes_v = clientKeyProof; + + unsigned char *bytes_A = 0; + unsigned char *bytes_B = 0; + + usr = srp_user_new(alg, ng_type, username, username, (const unsigned char *)password, strlen(password), NULL, NULL); + + srp_user_start_authentication(usr, NULL, (unsigned char *)srp_3072_a, 32, &bytes_A, &len_A); + + + static const uint8_t srp_5054_A[] = { + 0x61, 0xD5, 0xE4, 0x90, 0xF6, 0xF1, 0xB7, 0x95, 0x47, 0xB0, 0x70, 0x4C, 0x43, 0x6F, + 0x52, 0x3D, 0xD0, 0xE5, 0x60, 0xF0, 0xC6, 0x41, 0x15, 0xBB, 0x72, 0x55, 0x7E, 0xC4, + 0x43, 0x52, 0xE8, 0x90, 0x32, 0x11, 0xC0, 0x46, 0x92, 0x27, 0x2D, 0x8B, 0x2D, 0x1A, + 0x53, 0x58, 0xA2, 0xCF, 0x1B, 0x6E, 0x0B, 0xFC, 0xF9, 0x9F, 0x92, 0x15, 0x30, 0xEC, + 0x8E, 0x39, 0x35, 0x61, 0x79, 0xEA, 0xE4, 0x5E, 0x42, 0xBA, 0x92, 0xAE, 0xAC, 0xED, + 0x82, 0x51, 0x71, 0xE1, 0xE8, 0xB9, 0xAF, 0x6D, 0x9C, 0x03, 0xE1, 0x32, 0x7F, 0x44, + 0xBE, 0x08, 0x7E, 0xF0, 0x65, 0x30, 0xE6, 0x9F, 0x66, 0x61, 0x52, 0x61, 0xEE, 0xF5, + 0x40, 0x73, 0xCA, 0x11, 0xCF, 0x58, 0x58, 0xF0, 0xED, 0xFD, 0xFE, 0x15, 0xEF, 0xEA, + 0xB3, 0x49, 0xEF, 0x5D, 0x76, 0x98, 0x8A, 0x36, 0x72, 0xFA, 0xC4, 0x7B, 0x07, 0x69, + 0x44, 0x7B, + }; + + static const uint8_t srp_3072_A[] = { + 0xFA, 0xB6, 0xF5, 0xD2, 0x61, 0x5D, 0x1E, 0x32, 0x35, 0x12, 0xE7, 0x99, 0x1C, 0xC3, + 0x74, 0x43, 0xF4, 0x87, 0xDA, 0x60, 0x4C, 0xA8, 0xC9, 0x23, 0x0F, 0xCB, 0x04, 0xE5, + 0x41, 0xDC, 0xE6, 0x28, 0x0B, 0x27, 0xCA, 0x46, 0x80, 0xB0, 0x37, 0x4F, 0x17, 0x9D, + 0xC3, 0xBD, 0xC7, 0x55, 0x3F, 0xE6, 0x24, 0x59, 0x79, 0x8C, 0x70, 0x1A, 0xD8, 0x64, + 0xA9, 0x13, 0x90, 0xA2, 0x8C, 0x93, 0xB6, 0x44, 0xAD, 0xBF, 0x9C, 0x00, 0x74, 0x5B, + 0x94, 0x2B, 0x79, 0xF9, 0x01, 0x2A, 0x21, 0xB9, 0xB7, 0x87, 0x82, 0x31, 0x9D, 0x83, + 0xA1, 0xF8, 0x36, 0x28, 0x66, 0xFB, 0xD6, 0xF4, 0x6B, 0xFC, 0x0D, 0xDB, 0x2E, 0x1A, + 0xB6, 0xE4, 0xB4, 0x5A, 0x99, 0x06, 0xB8, 0x2E, 0x37, 0xF0, 0x5D, 0x6F, 0x97, 0xF6, + 0xA3, 0xEB, 0x6E, 0x18, 0x20, 0x79, 0x75, 0x9C, 0x4F, 0x68, 0x47, 0x83, 0x7B, 0x62, + 0x32, 0x1A, 0xC1, 0xB4, 0xFA, 0x68, 0x64, 0x1F, 0xCB, 0x4B, 0xB9, 0x8D, 0xD6, 0x97, + 0xA0, 0xC7, 0x36, 0x41, 0x38, 0x5F, 0x4B, 0xAB, 0x25, 0xB7, 0x93, 0x58, 0x4C, 0xC3, + 0x9F, 0xC8, 0xD4, 0x8D, 0x4B, 0xD8, 0x67, 0xA9, 0xA3, 0xC1, 0x0F, 0x8E, 0xA1, 0x21, + 0x70, 0x26, 0x8E, 0x34, 0xFE, 0x3B, 0xBE, 0x6F, 0xF8, 0x99, 0x98, 0xD6, 0x0D, 0xA2, + 0xF3, 0xE4, 0x28, 0x3C, 0xBE, 0xC1, 0x39, 0x3D, 0x52, 0xAF, 0x72, 0x4A, 0x57, 0x23, + 0x0C, 0x60, 0x4E, 0x9F, 0xBC, 0xE5, 0x83, 0xD7, 0x61, 0x3E, 0x6B, 0xFF, 0xD6, 0x75, + 0x96, 0xAD, 0x12, 0x1A, 0x87, 0x07, 0xEE, 0xC4, 0x69, 0x44, 0x95, 0x70, 0x33, 0x68, + 0x6A, 0x15, 0x5F, 0x64, 0x4D, 0x5C, 0x58, 0x63, 0xB4, 0x8F, 0x61, 0xBD, 0xBF, 0x19, + 0xA5, 0x3E, 0xAB, 0x6D, 0xAD, 0x0A, 0x18, 0x6B, 0x8C, 0x15, 0x2E, 0x5F, 0x5D, 0x8C, + 0xAD, 0x4B, 0x0E, 0xF8, 0xAA, 0x4E, 0xA5, 0x00, 0x88, 0x34, 0xC3, 0xCD, 0x34, 0x2E, + 0x5E, 0x0F, 0x16, 0x7A, 0xD0, 0x45, 0x92, 0xCD, 0x8B, 0xD2, 0x79, 0x63, 0x93, 0x98, + 0xEF, 0x9E, 0x11, 0x4D, 0xFA, 0xAA, 0xB9, 0x19, 0xE1, 0x4E, 0x85, 0x09, 0x89, 0x22, + 0x4D, 0xDD, 0x98, 0x57, 0x6D, 0x79, 0x38, 0x5D, 0x22, 0x10, 0x90, 0x2E, 0x9F, 0x9B, + 0x1F, 0x2D, 0x86, 0xCF, 0xA4, 0x7E, 0xE2, 0x44, 0x63, 0x54, 0x65, 0xF7, 0x10, 0x58, + 0x42, 0x1A, 0x01, 0x84, 0xBE, 0x51, 0xDD, 0x10, 0xCC, 0x9D, 0x07, 0x9E, 0x6F, 0x16, + 0x04, 0xE7, 0xAA, 0x9B, 0x7C, 0xF7, 0x88, 0x3C, 0x7D, 0x4C, 0xE1, 0x2B, 0x06, 0xEB, + 0xE1, 0x60, 0x81, 0xE2, 0x3F, 0x27, 0xA2, 0x31, 0xD1, 0x84, 0x32, 0xD7, 0xD1, 0xBB, + 0x55, 0xC2, 0x8A, 0xE2, 0x1F, 0xFC, 0xF0, 0x05, 0xF5, 0x75, 0x28, 0xD1, 0x5A, 0x88, + 0x88, 0x1B, 0xB3, 0xBB, 0xB7, 0xFE + }; + + if (memcmp(&srp_3072_A, bytes_A, len_A) != 0) { + return -1; + } + + ver = srp_verifier_new(alg, + ng_type, + username, + (unsigned char *)srp_5054_salt, + len_s, + bytes_v, + len_v, + bytes_A, + len_A, + (unsigned char *)srp_5054_b, + 32, + &bytes_B, + &len_B, + NULL, + NULL); + + + *serverKeyProof = bytes_B; + *proof_len = len_B; + + return 0; + +} + void SRPClient::createSaltedVerificationKey(uint8_t * salt, uint8_t * verificationKey ) { // let salt = salt ?? Data(bytes: try! Random.generate(byteCount: 16)) - salt = (uint8_t *)malloc(16); + salt = (uint8_t *)malloc(16 + 20); uint8_t defaultSalt[16] = {0xBE, 0xB2, 0x53, 0x79, 0xD1, 0xA8, 0x58, 0x1E, 0xB5, 0xA7, 0x27, 0x67, 0x3A, 0x24, 0x41, 0xEE}; memcpy(salt, defaultSalt, 16); - uint8_t * key = (uint8_t *)malloc(64); + uint8_t * key = (uint8_t *)malloc(20); + static char namePass[18] = "alice:password123"; + crypto_hash_sha1((uint8_t *)key, (uint8_t *)namePass, 17); + + memcpy(salt+16, key, 16); + crypto_hash_sha1((uint8_t *)key, (uint8_t *)salt, 16 + 20); - crypto_hash_sha512((uint8_t *)key, _username, 5); // let x = calculate_x(algorithm: algorithm, salt: salt, username: username, password: password) diff --git a/src/SRPClient.h b/src/SRPClient.h index 4003873..36adbbc 100644 --- a/src/SRPClient.h +++ b/src/SRPClient.h @@ -13,20 +13,55 @@ # include "application.h" #endif +#include "srp.h" + #include #include #include #include +extern "C" SRP_Result srp_create_salted_verification_key(SRP_HashAlgorithm alg, + SRP_NGType ng_type, const char *username_for_verifier, + const unsigned char *password, size_t len_password, + unsigned char **bytes_s, size_t *len_s, + unsigned char **bytes_v, size_t *len_v, + const char *n_hex, const char *g_hex); + + +extern "C" struct SRPVerifier* srp_verifier_new(SRP_HashAlgorithm alg, SRP_NGType ng_type, + const char *username, + const unsigned char *bytes_s, size_t len_s, + const unsigned char *bytes_v, size_t len_v, + const unsigned char *bytes_A, size_t len_A, + const unsigned char *bytes_b, size_t len_b, + unsigned char** bytes_B, size_t *len_B, + const char* n_hex, const char* g_hex); + +extern "C" struct SRPUser *srp_user_new(SRP_HashAlgorithm alg, SRP_NGType ng_type, + const char *username, const char *username_for_verifier, + const unsigned char *bytes_password, size_t len_password, const char *n_hex, + const char *g_hex); + +extern "C" SRP_Result srp_user_start_authentication(struct SRPUser* usr, char **username, + const unsigned char *bytes_a, size_t len_a, + unsigned char **bytes_A, size_t* len_A); + + + class SRPClient { public: SRPClient(); ~SRPClient(); + int getChallenge(uint8_t **salt,uint16_t *salt_len, uint8_t **key, uint16_t *key_len); + + int verifySession(uint8_t **serverKeyProof, uint16_t *proof_len, uint8_t *clientPublicKey, uint16_t client_key_len, uint8_t *clientKeyProof, uint16_t client_proof_len); + void createSaltedVerificationKey(uint8_t * salt, uint8_t * verificationKey ); int crypto_hash_sha512(unsigned char *out, const unsigned char *in, unsigned long long inlen); int crypto_hash_sha1(unsigned char *out, const unsigned char *in, unsigned long long inlen); + private: uint8_t * _salt; diff --git a/src/md32_common.h b/src/md32_common.h new file mode 100644 index 0000000..96828d2 --- /dev/null +++ b/src/md32_common.h @@ -0,0 +1,428 @@ +/* crypto/md32_common.h */ +/* ==================================================================== + * Copyright (c) 1999-2007 The OpenSSL Project. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * + * 3. All advertising materials mentioning features or use of this + * software must display the following acknowledgment: + * "This product includes software developed by the OpenSSL Project + * for use in the OpenSSL Toolkit. (http://www.OpenSSL.org/)" + * + * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to + * endorse or promote products derived from this software without + * prior written permission. For written permission, please contact + * licensing@OpenSSL.org. + * + * 5. Products derived from this software may not be called "OpenSSL" + * nor may "OpenSSL" appear in their names without prior written + * permission of the OpenSSL Project. + * + * 6. Redistributions of any form whatsoever must retain the following + * acknowledgment: + * "This product includes software developed by the OpenSSL Project + * for use in the OpenSSL Toolkit (http://www.OpenSSL.org/)" + * + * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY + * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE OpenSSL PROJECT OR + * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, + * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED + * OF THE POSSIBILITY OF SUCH DAMAGE. + * ==================================================================== + * + */ + +/*- + * This is a generic 32 bit "collector" for message digest algorithms. + * Whenever needed it collects input character stream into chunks of + * 32 bit values and invokes a block function that performs actual hash + * calculations. + * + * Porting guide. + * + * Obligatory macros: + * + * DATA_ORDER_IS_BIG_ENDIAN or DATA_ORDER_IS_LITTLE_ENDIAN + * this macro defines byte order of input stream. + * HASH_CBLOCK + * size of a unit chunk HASH_BLOCK operates on. + * HASH_LONG + * has to be at lest 32 bit wide, if it's wider, then + * HASH_LONG_LOG2 *has to* be defined along + * HASH_CTX + * context structure that at least contains following + * members: + * typedef struct { + * ... + * HASH_LONG Nl,Nh; + * either { + * HASH_LONG data[HASH_LBLOCK]; + * unsigned char data[HASH_CBLOCK]; + * }; + * unsigned int num; + * ... + * } HASH_CTX; + * data[] vector is expected to be zeroed upon first call to + * HASH_UPDATE. + * HASH_UPDATE + * name of "Update" function, implemented here. + * HASH_TRANSFORM + * name of "Transform" function, implemented here. + * HASH_FINAL + * name of "Final" function, implemented here. + * HASH_BLOCK_DATA_ORDER + * name of "block" function capable of treating *unaligned* input + * message in original (data) byte order, implemented externally. + * HASH_MAKE_STRING + * macro convering context variables to an ASCII hash string. + * + * MD5 example: + * + * #define DATA_ORDER_IS_LITTLE_ENDIAN + * + * #define HASH_LONG MD5_LONG + * #define HASH_LONG_LOG2 MD5_LONG_LOG2 + * #define HASH_CTX MD5_CTX + * #define HASH_CBLOCK MD5_CBLOCK + * #define HASH_UPDATE MD5_Update + * #define HASH_TRANSFORM MD5_Transform + * #define HASH_FINAL MD5_Final + * #define HASH_BLOCK_DATA_ORDER md5_block_data_order + * + * + */ + +#if !defined(DATA_ORDER_IS_BIG_ENDIAN) && !defined(DATA_ORDER_IS_LITTLE_ENDIAN) +# error "DATA_ORDER must be defined!" +#endif + +#ifndef HASH_CBLOCK +# error "HASH_CBLOCK must be defined!" +#endif +#ifndef HASH_LONG +# error "HASH_LONG must be defined!" +#endif +#ifndef HASH_CTX +# error "HASH_CTX must be defined!" +#endif + +#ifndef HASH_UPDATE +# error "HASH_UPDATE must be defined!" +#endif +#ifndef HASH_TRANSFORM +# error "HASH_TRANSFORM must be defined!" +#endif +#ifndef HASH_FINAL +# error "HASH_FINAL must be defined!" +#endif + +#ifndef HASH_BLOCK_DATA_ORDER +# error "HASH_BLOCK_DATA_ORDER must be defined!" +#endif + +/* + * Engage compiler specific rotate intrinsic function if available. + */ +#undef ROTATE +#ifndef PEDANTIC +# if defined(_MSC_VER) +# define ROTATE(a,n) _lrotl(a,n) +# elif defined(__ICC) +# define ROTATE(a,n) _rotl(a,n) +# elif defined(__MWERKS__) +# if defined(__POWERPC__) +# define ROTATE(a,n) __rlwinm(a,n,0,31) +# elif defined(__MC68K__) + /* Motorola specific tweak. */ +# define ROTATE(a,n) ( n<24 ? __rol(a,n) : __ror(a,32-n) ) +# else +# define ROTATE(a,n) __rol(a,n) +# endif +# elif defined(__GNUC__) && __GNUC__>=2 && !defined(OPENSSL_NO_ASM) && !defined(OPENSSL_NO_INLINE_ASM) + /* + * Some GNU C inline assembler templates. Note that these are + * rotates by *constant* number of bits! But that's exactly + * what we need here... + * + */ +# if defined(__i386) || defined(__i386__) || defined(__x86_64) || defined(__x86_64__) +# define ROTATE(a,n) ({ register unsigned int ret; \ + asm ( \ + "roll %1,%0" \ + : "=r"(ret) \ + : "I"(n), "0"((unsigned int)(a)) \ + : "cc"); \ + ret; \ + }) +# elif defined(_ARCH_PPC) || defined(_ARCH_PPC64) || \ + defined(__powerpc) || defined(__ppc__) || defined(__powerpc64__) +# define ROTATE(a,n) ({ register unsigned int ret; \ + asm ( \ + "rlwinm %0,%1,%2,0,31" \ + : "=r"(ret) \ + : "r"(a), "I"(n)); \ + ret; \ + }) +# elif defined(__s390x__) +# define ROTATE(a,n) ({ register unsigned int ret; \ + asm ("rll %0,%1,%2" \ + : "=r"(ret) \ + : "r"(a), "I"(n)); \ + ret; \ + }) +# endif +# endif +#endif /* PEDANTIC */ + +#ifndef ROTATE +# define ROTATE(a,n) (((a)<<(n))|(((a)&0xffffffff)>>(32-(n)))) +#endif + +#if defined(DATA_ORDER_IS_BIG_ENDIAN) + +# ifndef PEDANTIC +# if defined(__GNUC__) && __GNUC__>=2 && !defined(OPENSSL_NO_ASM) && !defined(OPENSSL_NO_INLINE_ASM) +# if ((defined(__i386) || defined(__i386__)) && !defined(I386_ONLY)) || \ + (defined(__x86_64) || defined(__x86_64__)) +# if !defined(B_ENDIAN) + /* + * This gives ~30-40% performance improvement in SHA-256 compiled + * with gcc [on P4]. Well, first macro to be frank. We can pull + * this trick on x86* platforms only, because these CPUs can fetch + * unaligned data without raising an exception. + */ +# define HOST_c2l(c,l) ({ unsigned int r=*((const unsigned int *)(c)); \ + asm ("bswapl %0":"=r"(r):"0"(r)); \ + (c)+=4; (l)=r; }) +# define HOST_l2c(l,c) ({ unsigned int r=(l); \ + asm ("bswapl %0":"=r"(r):"0"(r)); \ + *((unsigned int *)(c))=r; (c)+=4; r; }) +# endif +# elif defined(__aarch64__) +# if defined(__BYTE_ORDER__) +# if defined(__ORDER_LITTLE_ENDIAN__) && __BYTE_ORDER__==__ORDER_LITTLE_ENDIAN__ +# define HOST_c2l(c,l) ({ unsigned int r; \ + asm ("rev %w0,%w1" \ + :"=r"(r) \ + :"r"(*((const unsigned int *)(c))));\ + (c)+=4; (l)=r; }) +# define HOST_l2c(l,c) ({ unsigned int r; \ + asm ("rev %w0,%w1" \ + :"=r"(r) \ + :"r"((unsigned int)(l)));\ + *((unsigned int *)(c))=r; (c)+=4; r; }) +# elif defined(__ORDER_BIG_ENDIAN__) && __BYTE_ORDER__==__ORDER_BIG_ENDIAN__ +# define HOST_c2l(c,l) ((l)=*((const unsigned int *)(c)), (c)+=4, (l)) +# define HOST_l2c(l,c) (*((unsigned int *)(c))=(l), (c)+=4, (l)) +# endif +# endif +# endif +# endif +# if defined(__s390__) || defined(__s390x__) +# define HOST_c2l(c,l) ((l)=*((const unsigned int *)(c)), (c)+=4, (l)) +# define HOST_l2c(l,c) (*((unsigned int *)(c))=(l), (c)+=4, (l)) +# endif +# endif + +# ifndef HOST_c2l +# define HOST_c2l(c,l) (l =(((unsigned long)(*((c)++)))<<24), \ + l|=(((unsigned long)(*((c)++)))<<16), \ + l|=(((unsigned long)(*((c)++)))<< 8), \ + l|=(((unsigned long)(*((c)++))) ) ) +# endif +# ifndef HOST_l2c +# define HOST_l2c(l,c) (*((c)++)=(unsigned char)(((l)>>24)&0xff), \ + *((c)++)=(unsigned char)(((l)>>16)&0xff), \ + *((c)++)=(unsigned char)(((l)>> 8)&0xff), \ + *((c)++)=(unsigned char)(((l) )&0xff), \ + l) +# endif + +#elif defined(DATA_ORDER_IS_LITTLE_ENDIAN) + +# ifndef PEDANTIC +# if defined(__GNUC__) && __GNUC__>=2 && !defined(OPENSSL_NO_ASM) && !defined(OPENSSL_NO_INLINE_ASM) +# if defined(__s390x__) +# define HOST_c2l(c,l) ({ asm ("lrv %0,%1" \ + :"=d"(l) :"m"(*(const unsigned int *)(c)));\ + (c)+=4; (l); }) +# define HOST_l2c(l,c) ({ asm ("strv %1,%0" \ + :"=m"(*(unsigned int *)(c)) :"d"(l));\ + (c)+=4; (l); }) +# endif +# endif +# if defined(__i386) || defined(__i386__) || defined(__x86_64) || defined(__x86_64__) +# ifndef B_ENDIAN + /* See comment in DATA_ORDER_IS_BIG_ENDIAN section. */ +# define HOST_c2l(c,l) ((l)=*((const unsigned int *)(c)), (c)+=4, l) +# define HOST_l2c(l,c) (*((unsigned int *)(c))=(l), (c)+=4, l) +# endif +# endif +# endif + +# ifndef HOST_c2l +# define HOST_c2l(c,l) (l =(((unsigned long)(*((c)++))) ), \ + l|=(((unsigned long)(*((c)++)))<< 8), \ + l|=(((unsigned long)(*((c)++)))<<16), \ + l|=(((unsigned long)(*((c)++)))<<24) ) +# endif +# ifndef HOST_l2c +# define HOST_l2c(l,c) (*((c)++)=(unsigned char)(((l) )&0xff), \ + *((c)++)=(unsigned char)(((l)>> 8)&0xff), \ + *((c)++)=(unsigned char)(((l)>>16)&0xff), \ + *((c)++)=(unsigned char)(((l)>>24)&0xff), \ + l) +# endif + +#endif + +/* + * Time for some action:-) + */ + +int HASH_UPDATE(HASH_CTX *c, const void *data_, size_t len) +{ + const unsigned char *data = data_; + unsigned char *p; + HASH_LONG l; + size_t n; + + if (len == 0) + return 1; + + l = (c->Nl + (((HASH_LONG) len) << 3)) & 0xffffffffUL; + /* + * 95-05-24 eay Fixed a bug with the overflow handling, thanks to Wei Dai + * for pointing it out. + */ + if (l < c->Nl) /* overflow */ + c->Nh++; + c->Nh += (HASH_LONG) (len >> 29); /* might cause compiler warning on + * 16-bit */ + c->Nl = l; + + n = c->num; + if (n != 0) { + p = (unsigned char *)c->data; + + if (len >= HASH_CBLOCK || len + n >= HASH_CBLOCK) { + memcpy(p + n, data, HASH_CBLOCK - n); + HASH_BLOCK_DATA_ORDER(c, p, 1); + n = HASH_CBLOCK - n; + data += n; + len -= n; + c->num = 0; + memset(p, 0, HASH_CBLOCK); /* keep it zeroed */ + } else { + memcpy(p + n, data, len); + c->num += (unsigned int)len; + return 1; + } + } + + n = len / HASH_CBLOCK; + if (n > 0) { + HASH_BLOCK_DATA_ORDER(c, data, n); + n *= HASH_CBLOCK; + data += n; + len -= n; + } + + if (len != 0) { + p = (unsigned char *)c->data; + c->num = (unsigned int)len; + memcpy(p, data, len); + } + return 1; +} + +void HASH_TRANSFORM(HASH_CTX *c, const unsigned char *data) +{ + HASH_BLOCK_DATA_ORDER(c, data, 1); +} + +int HASH_FINAL(unsigned char *md, HASH_CTX *c) +{ + unsigned char *p = (unsigned char *)c->data; + size_t n = c->num; + + p[n] = 0x80; /* there is always room for one */ + n++; + + if (n > (HASH_CBLOCK - 8)) { + memset(p + n, 0, HASH_CBLOCK - n); + n = 0; + HASH_BLOCK_DATA_ORDER(c, p, 1); + } + memset(p + n, 0, HASH_CBLOCK - 8 - n); + + p += HASH_CBLOCK - 8; +#if defined(DATA_ORDER_IS_BIG_ENDIAN) + (void)HOST_l2c(c->Nh, p); + (void)HOST_l2c(c->Nl, p); +#elif defined(DATA_ORDER_IS_LITTLE_ENDIAN) + (void)HOST_l2c(c->Nl, p); + (void)HOST_l2c(c->Nh, p); +#endif + p -= HASH_CBLOCK; + HASH_BLOCK_DATA_ORDER(c, p, 1); + c->num = 0; + memset(p, 0, HASH_CBLOCK); + +#ifndef HASH_MAKE_STRING +# error "HASH_MAKE_STRING must be defined!" +#else + HASH_MAKE_STRING(c, md); +#endif + + return 1; +} + +#ifndef MD32_REG_T +# if defined(__alpha) || defined(__sparcv9) || defined(__mips) +# define MD32_REG_T long +/* + * This comment was originaly written for MD5, which is why it + * discusses A-D. But it basically applies to all 32-bit digests, + * which is why it was moved to common header file. + * + * In case you wonder why A-D are declared as long and not + * as MD5_LONG. Doing so results in slight performance + * boost on LP64 architectures. The catch is we don't + * really care if 32 MSBs of a 64-bit register get polluted + * with eventual overflows as we *save* only 32 LSBs in + * *either* case. Now declaring 'em long excuses the compiler + * from keeping 32 MSBs zeroed resulting in 13% performance + * improvement under SPARC Solaris7/64 and 5% under AlphaLinux. + * Well, to be honest it should say that this *prevents* + * performance degradation. + * + */ +# else +/* + * Above is not absolute and there are LP64 compilers that + * generate better code if MD32_REG_T is defined int. The above + * pre-processor condition reflects the circumstances under which + * the conclusion was made and is subject to further extension. + * + */ +# define MD32_REG_T int +# endif +#endif diff --git a/src/mini-gmp.c b/src/mini-gmp.c new file mode 100644 index 0000000..9a62855 --- /dev/null +++ b/src/mini-gmp.c @@ -0,0 +1,4381 @@ +/* mini-gmp, a minimalistic implementation of a GNU GMP subset. + + Contributed to the GNU project by Niels Möller + +Copyright 1991-1997, 1999-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/. */ + +/* 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_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. */ +mp_limb_t +mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0) +{ + mp_limb_t r, p, m; + unsigned ul, uh; + unsigned ql, qh; + + /* First, do a 2/1 inverse. */ + /* The inverse m is defined as floor( (B^2 - 1 - u1)/u1 ), so that 0 < + * B^2 - (B + m) u1 <= u1 */ + assert (u1 >= GMP_LIMB_HIGHBIT); + + ul = u1 & GMP_LLIMB_MASK; + uh = u1 >> (GMP_LIMB_BITS / 2); + + qh = ~u1 / uh; + 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; + + /* Do a 3/2 division (with half limb size) */ + p = (r >> (GMP_LIMB_BITS / 2)) * qh + r; + ql = (p >> (GMP_LIMB_BITS / 2)) + 1; + + /* By the 3/2 method, 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; + } + + 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; +} + +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; + + k = 1 + (sn - 1) % info->exp; + + j = 0; + w = sp[j++]; + while (--k != 0) + w = w * b + sp[j++]; + + rp[0] = w; + + for (rn = (w > 0); 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) +{ + mp_size_t us = u->_mp_size; + + if (us > 0) + return (long) (u->_mp_d[0] & ~GMP_LIMB_HIGHBIT); + else if (us < 0) + return (long) (- u->_mp_d[0] | GMP_LIMB_HIGHBIT); + else + return 0; +} + +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) +{ + mp_size_t usize = u->_mp_size; + + return (usize > 0) - (usize < 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 */ + { + mp_limb_t ul = u->_mp_d[0]; + if ((mp_limb_t)GMP_NEG_CAST (unsigned long int, v) < ul) + return -1; + else + return (mp_limb_t)GMP_NEG_CAST (unsigned long int, v) > ul; + } +} + +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 + { + mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0; + return (ul > v) - (ul < 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) +{ + mp_size_t un = GMP_ABS (u->_mp_size); + mp_limb_t ul; + + if (un > 1) + return 1; + + ul = (un == 1) ? u->_mp_d[0] : 0; + + return (ul > v) - (ul < 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 sn; + 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') + { + sp++; + if (*sp == 'x' || *sp == 'X') + { + base = 16; + sp++; + } + else if (*sp == 'b' || *sp == 'B') + { + base = 2; + sp++; + } + else + base = 8; + } + else + base = 10; + } + + sn = strlen (sp); + dp = (unsigned char *) gmp_xalloc (sn + (sn == 0)); + + for (sn = 0; *sp; sp++) + { + unsigned digit; + + if (isspace ((unsigned char) *sp)) + continue; + 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[sn++] = digit; + } + + bits = mpn_base_power_of_two_p (base); + + if (bits > 0) + { + alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS; + rp = MPZ_REALLOC (r, alloc); + rn = mpn_set_str_bits (rp, dp, sn, bits); + } + else + { + struct mpn_base_info info; + mpn_get_base_info (&info, base); + alloc = (sn + info.exp - 1) / info.exp; + rp = MPZ_REALLOC (r, alloc); + rn = mpn_set_str_other (rp, dp, sn, base, &info); + } + 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.h b/src/mini-gmp.h new file mode 100644 index 0000000..bb5c637 --- /dev/null +++ b/src/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/sha.h b/src/sha.h new file mode 100644 index 0000000..0c1d748 --- /dev/null +++ b/src/sha.h @@ -0,0 +1,175 @@ +/* crypto/sha/sha.h */ +/* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com) + * All rights reserved. + * + * This package is an SSL implementation written + * by Eric Young (eay@cryptsoft.com). + * The implementation was written so as to conform with Netscapes SSL. + * + * This library is free for commercial and non-commercial use as long as + * the following conditions are aheared to. The following conditions + * apply to all code found in this distribution, be it the RC4, RSA, + * lhash, DES, etc., code; not just the SSL code. The SSL documentation + * included with this distribution is covered by the same copyright terms + * except that the holder is Tim Hudson (tjh@cryptsoft.com). + * + * Copyright remains Eric Young's, and as such any Copyright notices in + * the code are not to be removed. + * If this package is used in a product, Eric Young should be given attribution + * as the author of the parts of the library used. + * This can be in the form of a textual message at program startup or + * in documentation (online or textual) provided with the package. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * 3. All advertising materials mentioning features or use of this software + * must display the following acknowledgement: + * "This product includes cryptographic software written by + * Eric Young (eay@cryptsoft.com)" + * The word 'cryptographic' can be left out if the rouines from the library + * being used are not cryptographic related :-). + * 4. If you include any Windows specific code (or a derivative thereof) from + * the apps directory (application code) you must include an acknowledgement: + * "This product includes software written by Tim Hudson (tjh@cryptsoft.com)" + * + * THIS SOFTWARE IS PROVIDED BY ERIC YOUNG ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * + * The licence and distribution terms for any publically available version or + * derivative of this code cannot be changed. i.e. this code cannot simply be + * copied and put under another distribution licence + * [including the GNU Public Licence.] + */ + +#ifndef HEADER_SHA_H +# define HEADER_SHA_H + +# include + +#ifdef __cplusplus +extern "C" { +#endif + +# if defined(OPENSSL_NO_SHA) || (defined(OPENSSL_NO_SHA0) && defined(OPENSSL_NO_SHA1)) +# error SHA is disabled. +# endif + +# if defined(OPENSSL_FIPS) +# define FIPS_SHA_SIZE_T size_t +# endif + +/* + Compat stuff from OpenSSL land + */ + +/* crypto.h */ + +# define fips_md_init(alg) fips_md_init_ctx(alg, alg) + +# define fips_md_init_ctx(alg, cx) \ + int alg##_Init(cx##_CTX *c) +# define fips_cipher_abort(alg) while(0) + +/*- + * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + * ! SHA_LONG has to be at least 32 bits wide. If it's wider, then ! + * ! SHA_LONG_LOG2 has to be defined along. ! + * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + */ + +# if defined(__LP32__) +# define SHA_LONG unsigned long +# elif defined(__ILP64__) +# define SHA_LONG unsigned long +# define SHA_LONG_LOG2 3 +# else +# define SHA_LONG unsigned int +# endif + +# define SHA_LBLOCK 16 +# define SHA_CBLOCK (SHA_LBLOCK*4)/* SHA treats input data as a + * contiguous array of 32 bit wide + * big-endian values. */ +# define SHA_LAST_BLOCK (SHA_CBLOCK-8) +# define SHA_DIGEST_LENGTH 20 + +typedef struct SHAstate_st { + SHA_LONG h0, h1, h2, h3, h4; + SHA_LONG Nl, Nh; + SHA_LONG data[SHA_LBLOCK]; + unsigned int num; +} SHA_CTX; + +# define SHA256_CBLOCK (SHA_LBLOCK*4)/* SHA-256 treats input data as a + * contiguous array of 32 bit wide + * big-endian values. */ + +int SHA1_Init(SHA_CTX *c); +int SHA1_Update(SHA_CTX *c, const void *data, size_t len); +int SHA1_Final(unsigned char *md, SHA_CTX *c); +unsigned char *SHA1(const unsigned char *d, size_t n, unsigned char *md); +void SHA1_Transform(SHA_CTX *c, const unsigned char *data); + +# define SHA224_DIGEST_LENGTH 28 +# define SHA256_DIGEST_LENGTH 32 + +typedef struct SHA256state_st { + SHA_LONG h[8]; + SHA_LONG Nl, Nh; + SHA_LONG data[SHA_LBLOCK]; + unsigned int num, md_len; +} SHA256_CTX; + +# ifndef OPENSSL_NO_SHA256 +# ifdef OPENSSL_FIPS +int private_SHA224_Init(SHA256_CTX *c); +int private_SHA256_Init(SHA256_CTX *c); +# endif +int SHA224_Init(SHA256_CTX *c); +int SHA224_Update(SHA256_CTX *c, const void *data, size_t len); +int SHA224_Final(unsigned char *md, SHA256_CTX *c); +unsigned char *SHA224(const unsigned char *d, size_t n, unsigned char *md); +int SHA256_Init(SHA256_CTX *c); +int SHA256_Update(SHA256_CTX *c, const void *data, size_t len); +int SHA256_Final(unsigned char *md, SHA256_CTX *c); +unsigned char *SHA256(const unsigned char *d, size_t n, unsigned char *md); +void SHA256_Transform(SHA256_CTX *c, const unsigned char *data); +# endif + +# define SHA384_DIGEST_LENGTH 48 +# define SHA512_DIGEST_LENGTH 64 + +typedef struct SHA512state_st { + uint64_t state[8]; + uint64_t count[2]; + uint8_t buf[128]; + +} SHA512_CTX; + +int SHA512_Init(SHA512_CTX *c); +int SHA512_Update(SHA512_CTX *c, const void *data, size_t len); +int SHA512_Final(unsigned char *md, SHA512_CTX *c); +unsigned char *SHA512(const unsigned char *d, size_t n, unsigned char *md); +void SHA512_Transform(SHA512_CTX *c, const unsigned char *data); + + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/sha1.c b/src/sha1.c new file mode 100644 index 0000000..499719e --- /dev/null +++ b/src/sha1.c @@ -0,0 +1,537 @@ +/* crypto/sha/sha_locl.h */ +/* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com) + * All rights reserved. + * + * This package is an SSL implementation written + * by Eric Young (eay@cryptsoft.com). + * The implementation was written so as to conform with Netscapes SSL. + * + * This library is free for commercial and non-commercial use as long as + * the following conditions are aheared to. The following conditions + * apply to all code found in this distribution, be it the RC4, RSA, + * lhash, DES, etc., code; not just the SSL code. The SSL documentation + * included with this distribution is covered by the same copyright terms + * except that the holder is Tim Hudson (tjh@cryptsoft.com). + * + * Copyright remains Eric Young's, and as such any Copyright notices in + * the code are not to be removed. + * If this package is used in a product, Eric Young should be given attribution + * as the author of the parts of the library used. + * This can be in the form of a textual message at program startup or + * in documentation (online or textual) provided with the package. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * 3. All advertising materials mentioning features or use of this software + * must display the following acknowledgement: + * "This product includes cryptographic software written by + * Eric Young (eay@cryptsoft.com)" + * The word 'cryptographic' can be left out if the rouines from the library + * being used are not cryptographic related :-). + * 4. If you include any Windows specific code (or a derivative thereof) from + * the apps directory (application code) you must include an acknowledgement: + * "This product includes software written by Tim Hudson (tjh@cryptsoft.com)" + * + * THIS SOFTWARE IS PROVIDED BY ERIC YOUNG ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * + * The licence and distribution terms for any publically available version or + * derivative of this code cannot be changed. i.e. this code cannot simply be + * copied and put under another distribution licence + * [including the GNU Public Licence.] + */ + +#include +#include + +#include "sha.h" + +/* mem_clr.c */ +unsigned static char cleanse_ctr = 0; +static void OPENSSL_cleanse(void *ptr, size_t len) +{ + unsigned char *p = ptr; + size_t loop = len, ctr = cleanse_ctr; + while (loop--) { + *(p++) = (unsigned char)ctr; + ctr += (17 + ((size_t)p & 0xF)); + } + p = memchr(ptr, (unsigned char)ctr, len); + if (p) + ctr += (63 + (size_t)p); + cleanse_ctr = (unsigned char)ctr; +} + +#define DATA_ORDER_IS_BIG_ENDIAN + +#define SHA_1 + +#define HASH_LONG SHA_LONG +#define HASH_CTX SHA_CTX +#define HASH_CBLOCK SHA_CBLOCK +#define HASH_MAKE_STRING(c,s) do { \ + unsigned long ll; \ + ll=(c)->h0; (void)HOST_l2c(ll,(s)); \ + ll=(c)->h1; (void)HOST_l2c(ll,(s)); \ + ll=(c)->h2; (void)HOST_l2c(ll,(s)); \ + ll=(c)->h3; (void)HOST_l2c(ll,(s)); \ + ll=(c)->h4; (void)HOST_l2c(ll,(s)); \ + } while (0) + +#if defined(SHA_0) + +# define HASH_UPDATE SHA_Update +# define HASH_TRANSFORM SHA_Transform +# define HASH_FINAL SHA_Final +# define HASH_INIT SHA_Init +# define HASH_BLOCK_DATA_ORDER sha_block_data_order +# define Xupdate(a,ix,ia,ib,ic,id) (ix=(a)=(ia^ib^ic^id)) + +static void sha_block_data_order(SHA_CTX *c, const void *p, size_t num); + +#elif defined(SHA_1) + +# define HASH_UPDATE SHA1_Update +# define HASH_TRANSFORM SHA1_Transform +# define HASH_FINAL SHA1_Final +# define HASH_INIT SHA1_Init +# define HASH_BLOCK_DATA_ORDER sha1_block_data_order +# if defined(__MWERKS__) && defined(__MC68K__) + /* Metrowerks for Motorola fails otherwise:-( */ +# define Xupdate(a,ix,ia,ib,ic,id) do { (a)=(ia^ib^ic^id); \ + ix=(a)=ROTATE((a),1); \ + } while (0) +# else +# define Xupdate(a,ix,ia,ib,ic,id) ( (a)=(ia^ib^ic^id), \ + ix=(a)=ROTATE((a),1) \ + ) +# endif + +# ifndef SHA1_ASM +static +# endif +void sha1_block_data_order(SHA_CTX *c, const void *p, size_t num); + +#else +# error "Either SHA_0 or SHA_1 must be defined." +#endif + +#include "md32_common.h" + +#define INIT_DATA_h0 0x67452301UL +#define INIT_DATA_h1 0xefcdab89UL +#define INIT_DATA_h2 0x98badcfeUL +#define INIT_DATA_h3 0x10325476UL +#define INIT_DATA_h4 0xc3d2e1f0UL + +# define fips_md_init(alg) fips_md_init_ctx(alg, alg) +# define fips_md_init_ctx(alg, cx) \ + int alg##_Init(cx##_CTX *c) +# define fips_cipher_abort(alg) while(0) + +unsigned char *SHA1(const unsigned char *d, size_t n, unsigned char *md) +{ + SHA_CTX c; + static unsigned char m[SHA_DIGEST_LENGTH]; + + if (md == NULL) + md = m; + if (!SHA1_Init(&c)) + return NULL; + SHA1_Update(&c, d, n); + SHA1_Final(md, &c); + OPENSSL_cleanse(&c, sizeof(c)); + return (md); +} + +#ifdef SHA_0 +fips_md_init(SHA) +#else +fips_md_init_ctx(SHA1, SHA) +#endif +{ + memset(c, 0, sizeof(*c)); + c->h0 = INIT_DATA_h0; + c->h1 = INIT_DATA_h1; + c->h2 = INIT_DATA_h2; + c->h3 = INIT_DATA_h3; + c->h4 = INIT_DATA_h4; + return 1; +} + +#define K_00_19 0x5a827999UL +#define K_20_39 0x6ed9eba1UL +#define K_40_59 0x8f1bbcdcUL +#define K_60_79 0xca62c1d6UL + +/* + * As pointed out by Wei Dai , F() below can be simplified + * to the code in F_00_19. Wei attributes these optimisations to Peter + * Gutmann's SHS code, and he attributes it to Rich Schroeppel. #define + * F(x,y,z) (((x) & (y)) | ((~(x)) & (z))) I've just become aware of another + * tweak to be made, again from Wei Dai, in F_40_59, (x&a)|(y&a) -> (x|y)&a + */ +#define F_00_19(b,c,d) ((((c) ^ (d)) & (b)) ^ (d)) +#define F_20_39(b,c,d) ((b) ^ (c) ^ (d)) +#define F_40_59(b,c,d) (((b) & (c)) | (((b)|(c)) & (d))) +#define F_60_79(b,c,d) F_20_39(b,c,d) + +#ifndef OPENSSL_SMALL_FOOTPRINT + +# define BODY_00_15(i,a,b,c,d,e,f,xi) \ + (f)=xi+(e)+K_00_19+ROTATE((a),5)+F_00_19((b),(c),(d)); \ + (b)=ROTATE((b),30); + +# define BODY_16_19(i,a,b,c,d,e,f,xi,xa,xb,xc,xd) \ + Xupdate(f,xi,xa,xb,xc,xd); \ + (f)+=(e)+K_00_19+ROTATE((a),5)+F_00_19((b),(c),(d)); \ + (b)=ROTATE((b),30); + +# define BODY_20_31(i,a,b,c,d,e,f,xi,xa,xb,xc,xd) \ + Xupdate(f,xi,xa,xb,xc,xd); \ + (f)+=(e)+K_20_39+ROTATE((a),5)+F_20_39((b),(c),(d)); \ + (b)=ROTATE((b),30); + +# define BODY_32_39(i,a,b,c,d,e,f,xa,xb,xc,xd) \ + Xupdate(f,xa,xa,xb,xc,xd); \ + (f)+=(e)+K_20_39+ROTATE((a),5)+F_20_39((b),(c),(d)); \ + (b)=ROTATE((b),30); + +# define BODY_40_59(i,a,b,c,d,e,f,xa,xb,xc,xd) \ + Xupdate(f,xa,xa,xb,xc,xd); \ + (f)+=(e)+K_40_59+ROTATE((a),5)+F_40_59((b),(c),(d)); \ + (b)=ROTATE((b),30); + +# define BODY_60_79(i,a,b,c,d,e,f,xa,xb,xc,xd) \ + Xupdate(f,xa,xa,xb,xc,xd); \ + (f)=xa+(e)+K_60_79+ROTATE((a),5)+F_60_79((b),(c),(d)); \ + (b)=ROTATE((b),30); + +# ifdef X +# undef X +# endif +# ifndef MD32_XARRAY + /* + * Originally X was an array. As it's automatic it's natural + * to expect RISC compiler to accomodate at least part of it in + * the register bank, isn't it? Unfortunately not all compilers + * "find" this expectation reasonable:-( On order to make such + * compilers generate better code I replace X[] with a bunch of + * X0, X1, etc. See the function body below... + * + */ +# define X(i) XX##i +# else + /* + * However! Some compilers (most notably HP C) get overwhelmed by + * that many local variables so that we have to have the way to + * fall down to the original behavior. + */ +# define X(i) XX[i] +# endif + +# if !defined(SHA_1) || !defined(SHA1_ASM) +static void HASH_BLOCK_DATA_ORDER(SHA_CTX *c, const void *p, size_t num) +{ + const unsigned char *data = p; + register unsigned MD32_REG_T A, B, C, D, E, T, l; +# ifndef MD32_XARRAY + unsigned MD32_REG_T XX0, XX1, XX2, XX3, XX4, XX5, XX6, XX7, + XX8, XX9, XX10, XX11, XX12, XX13, XX14, XX15; +# else + SHA_LONG XX[16]; +# endif + + A = c->h0; + B = c->h1; + C = c->h2; + D = c->h3; + E = c->h4; + + for (;;) { + const union { + long one; + char little; + } is_endian = { + 1 + }; + + if (!is_endian.little && sizeof(SHA_LONG) == 4 + && ((size_t)p % 4) == 0) { + const SHA_LONG *W = (const SHA_LONG *)data; + + X(0) = W[0]; + X(1) = W[1]; + BODY_00_15(0, A, B, C, D, E, T, X(0)); + X(2) = W[2]; + BODY_00_15(1, T, A, B, C, D, E, X(1)); + X(3) = W[3]; + BODY_00_15(2, E, T, A, B, C, D, X(2)); + X(4) = W[4]; + BODY_00_15(3, D, E, T, A, B, C, X(3)); + X(5) = W[5]; + BODY_00_15(4, C, D, E, T, A, B, X(4)); + X(6) = W[6]; + BODY_00_15(5, B, C, D, E, T, A, X(5)); + X(7) = W[7]; + BODY_00_15(6, A, B, C, D, E, T, X(6)); + X(8) = W[8]; + BODY_00_15(7, T, A, B, C, D, E, X(7)); + X(9) = W[9]; + BODY_00_15(8, E, T, A, B, C, D, X(8)); + X(10) = W[10]; + BODY_00_15(9, D, E, T, A, B, C, X(9)); + X(11) = W[11]; + BODY_00_15(10, C, D, E, T, A, B, X(10)); + X(12) = W[12]; + BODY_00_15(11, B, C, D, E, T, A, X(11)); + X(13) = W[13]; + BODY_00_15(12, A, B, C, D, E, T, X(12)); + X(14) = W[14]; + BODY_00_15(13, T, A, B, C, D, E, X(13)); + X(15) = W[15]; + BODY_00_15(14, E, T, A, B, C, D, X(14)); + BODY_00_15(15, D, E, T, A, B, C, X(15)); + + data += SHA_CBLOCK; + } else { + (void)HOST_c2l(data, l); + X(0) = l; + (void)HOST_c2l(data, l); + X(1) = l; + BODY_00_15(0, A, B, C, D, E, T, X(0)); + (void)HOST_c2l(data, l); + X(2) = l; + BODY_00_15(1, T, A, B, C, D, E, X(1)); + (void)HOST_c2l(data, l); + X(3) = l; + BODY_00_15(2, E, T, A, B, C, D, X(2)); + (void)HOST_c2l(data, l); + X(4) = l; + BODY_00_15(3, D, E, T, A, B, C, X(3)); + (void)HOST_c2l(data, l); + X(5) = l; + BODY_00_15(4, C, D, E, T, A, B, X(4)); + (void)HOST_c2l(data, l); + X(6) = l; + BODY_00_15(5, B, C, D, E, T, A, X(5)); + (void)HOST_c2l(data, l); + X(7) = l; + BODY_00_15(6, A, B, C, D, E, T, X(6)); + (void)HOST_c2l(data, l); + X(8) = l; + BODY_00_15(7, T, A, B, C, D, E, X(7)); + (void)HOST_c2l(data, l); + X(9) = l; + BODY_00_15(8, E, T, A, B, C, D, X(8)); + (void)HOST_c2l(data, l); + X(10) = l; + BODY_00_15(9, D, E, T, A, B, C, X(9)); + (void)HOST_c2l(data, l); + X(11) = l; + BODY_00_15(10, C, D, E, T, A, B, X(10)); + (void)HOST_c2l(data, l); + X(12) = l; + BODY_00_15(11, B, C, D, E, T, A, X(11)); + (void)HOST_c2l(data, l); + X(13) = l; + BODY_00_15(12, A, B, C, D, E, T, X(12)); + (void)HOST_c2l(data, l); + X(14) = l; + BODY_00_15(13, T, A, B, C, D, E, X(13)); + (void)HOST_c2l(data, l); + X(15) = l; + BODY_00_15(14, E, T, A, B, C, D, X(14)); + BODY_00_15(15, D, E, T, A, B, C, X(15)); + } + + BODY_16_19(16, C, D, E, T, A, B, X(0), X(0), X(2), X(8), X(13)); + BODY_16_19(17, B, C, D, E, T, A, X(1), X(1), X(3), X(9), X(14)); + BODY_16_19(18, A, B, C, D, E, T, X(2), X(2), X(4), X(10), X(15)); + BODY_16_19(19, T, A, B, C, D, E, X(3), X(3), X(5), X(11), X(0)); + + BODY_20_31(20, E, T, A, B, C, D, X(4), X(4), X(6), X(12), X(1)); + BODY_20_31(21, D, E, T, A, B, C, X(5), X(5), X(7), X(13), X(2)); + BODY_20_31(22, C, D, E, T, A, B, X(6), X(6), X(8), X(14), X(3)); + BODY_20_31(23, B, C, D, E, T, A, X(7), X(7), X(9), X(15), X(4)); + BODY_20_31(24, A, B, C, D, E, T, X(8), X(8), X(10), X(0), X(5)); + BODY_20_31(25, T, A, B, C, D, E, X(9), X(9), X(11), X(1), X(6)); + BODY_20_31(26, E, T, A, B, C, D, X(10), X(10), X(12), X(2), X(7)); + BODY_20_31(27, D, E, T, A, B, C, X(11), X(11), X(13), X(3), X(8)); + BODY_20_31(28, C, D, E, T, A, B, X(12), X(12), X(14), X(4), X(9)); + BODY_20_31(29, B, C, D, E, T, A, X(13), X(13), X(15), X(5), X(10)); + BODY_20_31(30, A, B, C, D, E, T, X(14), X(14), X(0), X(6), X(11)); + BODY_20_31(31, T, A, B, C, D, E, X(15), X(15), X(1), X(7), X(12)); + + BODY_32_39(32, E, T, A, B, C, D, X(0), X(2), X(8), X(13)); + BODY_32_39(33, D, E, T, A, B, C, X(1), X(3), X(9), X(14)); + BODY_32_39(34, C, D, E, T, A, B, X(2), X(4), X(10), X(15)); + BODY_32_39(35, B, C, D, E, T, A, X(3), X(5), X(11), X(0)); + BODY_32_39(36, A, B, C, D, E, T, X(4), X(6), X(12), X(1)); + BODY_32_39(37, T, A, B, C, D, E, X(5), X(7), X(13), X(2)); + BODY_32_39(38, E, T, A, B, C, D, X(6), X(8), X(14), X(3)); + BODY_32_39(39, D, E, T, A, B, C, X(7), X(9), X(15), X(4)); + + BODY_40_59(40, C, D, E, T, A, B, X(8), X(10), X(0), X(5)); + BODY_40_59(41, B, C, D, E, T, A, X(9), X(11), X(1), X(6)); + BODY_40_59(42, A, B, C, D, E, T, X(10), X(12), X(2), X(7)); + BODY_40_59(43, T, A, B, C, D, E, X(11), X(13), X(3), X(8)); + BODY_40_59(44, E, T, A, B, C, D, X(12), X(14), X(4), X(9)); + BODY_40_59(45, D, E, T, A, B, C, X(13), X(15), X(5), X(10)); + BODY_40_59(46, C, D, E, T, A, B, X(14), X(0), X(6), X(11)); + BODY_40_59(47, B, C, D, E, T, A, X(15), X(1), X(7), X(12)); + BODY_40_59(48, A, B, C, D, E, T, X(0), X(2), X(8), X(13)); + BODY_40_59(49, T, A, B, C, D, E, X(1), X(3), X(9), X(14)); + BODY_40_59(50, E, T, A, B, C, D, X(2), X(4), X(10), X(15)); + BODY_40_59(51, D, E, T, A, B, C, X(3), X(5), X(11), X(0)); + BODY_40_59(52, C, D, E, T, A, B, X(4), X(6), X(12), X(1)); + BODY_40_59(53, B, C, D, E, T, A, X(5), X(7), X(13), X(2)); + BODY_40_59(54, A, B, C, D, E, T, X(6), X(8), X(14), X(3)); + BODY_40_59(55, T, A, B, C, D, E, X(7), X(9), X(15), X(4)); + BODY_40_59(56, E, T, A, B, C, D, X(8), X(10), X(0), X(5)); + BODY_40_59(57, D, E, T, A, B, C, X(9), X(11), X(1), X(6)); + BODY_40_59(58, C, D, E, T, A, B, X(10), X(12), X(2), X(7)); + BODY_40_59(59, B, C, D, E, T, A, X(11), X(13), X(3), X(8)); + + BODY_60_79(60, A, B, C, D, E, T, X(12), X(14), X(4), X(9)); + BODY_60_79(61, T, A, B, C, D, E, X(13), X(15), X(5), X(10)); + BODY_60_79(62, E, T, A, B, C, D, X(14), X(0), X(6), X(11)); + BODY_60_79(63, D, E, T, A, B, C, X(15), X(1), X(7), X(12)); + BODY_60_79(64, C, D, E, T, A, B, X(0), X(2), X(8), X(13)); + BODY_60_79(65, B, C, D, E, T, A, X(1), X(3), X(9), X(14)); + BODY_60_79(66, A, B, C, D, E, T, X(2), X(4), X(10), X(15)); + BODY_60_79(67, T, A, B, C, D, E, X(3), X(5), X(11), X(0)); + BODY_60_79(68, E, T, A, B, C, D, X(4), X(6), X(12), X(1)); + BODY_60_79(69, D, E, T, A, B, C, X(5), X(7), X(13), X(2)); + BODY_60_79(70, C, D, E, T, A, B, X(6), X(8), X(14), X(3)); + BODY_60_79(71, B, C, D, E, T, A, X(7), X(9), X(15), X(4)); + BODY_60_79(72, A, B, C, D, E, T, X(8), X(10), X(0), X(5)); + BODY_60_79(73, T, A, B, C, D, E, X(9), X(11), X(1), X(6)); + BODY_60_79(74, E, T, A, B, C, D, X(10), X(12), X(2), X(7)); + BODY_60_79(75, D, E, T, A, B, C, X(11), X(13), X(3), X(8)); + BODY_60_79(76, C, D, E, T, A, B, X(12), X(14), X(4), X(9)); + BODY_60_79(77, B, C, D, E, T, A, X(13), X(15), X(5), X(10)); + BODY_60_79(78, A, B, C, D, E, T, X(14), X(0), X(6), X(11)); + BODY_60_79(79, T, A, B, C, D, E, X(15), X(1), X(7), X(12)); + + c->h0 = (c->h0 + E) & 0xffffffffL; + c->h1 = (c->h1 + T) & 0xffffffffL; + c->h2 = (c->h2 + A) & 0xffffffffL; + c->h3 = (c->h3 + B) & 0xffffffffL; + c->h4 = (c->h4 + C) & 0xffffffffL; + + if (--num == 0) + break; + + A = c->h0; + B = c->h1; + C = c->h2; + D = c->h3; + E = c->h4; + + } +} +# endif + +#else /* OPENSSL_SMALL_FOOTPRINT */ + +# define BODY_00_15(xi) do { \ + T=E+K_00_19+F_00_19(B,C,D); \ + E=D, D=C, C=ROTATE(B,30), B=A; \ + A=ROTATE(A,5)+T+xi; } while(0) + +# define BODY_16_19(xa,xb,xc,xd) do { \ + Xupdate(T,xa,xa,xb,xc,xd); \ + T+=E+K_00_19+F_00_19(B,C,D); \ + E=D, D=C, C=ROTATE(B,30), B=A; \ + A=ROTATE(A,5)+T; } while(0) + +# define BODY_20_39(xa,xb,xc,xd) do { \ + Xupdate(T,xa,xa,xb,xc,xd); \ + T+=E+K_20_39+F_20_39(B,C,D); \ + E=D, D=C, C=ROTATE(B,30), B=A; \ + A=ROTATE(A,5)+T; } while(0) + +# define BODY_40_59(xa,xb,xc,xd) do { \ + Xupdate(T,xa,xa,xb,xc,xd); \ + T+=E+K_40_59+F_40_59(B,C,D); \ + E=D, D=C, C=ROTATE(B,30), B=A; \ + A=ROTATE(A,5)+T; } while(0) + +# define BODY_60_79(xa,xb,xc,xd) do { \ + Xupdate(T,xa,xa,xb,xc,xd); \ + T=E+K_60_79+F_60_79(B,C,D); \ + E=D, D=C, C=ROTATE(B,30), B=A; \ + A=ROTATE(A,5)+T+xa; } while(0) + +# if !defined(SHA_1) || !defined(SHA1_ASM) +static void HASH_BLOCK_DATA_ORDER(SHA_CTX *c, const void *p, size_t num) +{ + const unsigned char *data = p; + register unsigned MD32_REG_T A, B, C, D, E, T, l; + int i; + SHA_LONG X[16]; + + A = c->h0; + B = c->h1; + C = c->h2; + D = c->h3; + E = c->h4; + + for (;;) { + for (i = 0; i < 16; i++) { + HOST_c2l(data, l); + X[i] = l; + BODY_00_15(X[i]); + } + for (i = 0; i < 4; i++) { + BODY_16_19(X[i], X[i + 2], X[i + 8], X[(i + 13) & 15]); + } + for (; i < 24; i++) { + BODY_20_39(X[i & 15], X[(i + 2) & 15], X[(i + 8) & 15], + X[(i + 13) & 15]); + } + for (i = 0; i < 20; i++) { + BODY_40_59(X[(i + 8) & 15], X[(i + 10) & 15], X[i & 15], + X[(i + 5) & 15]); + } + for (i = 4; i < 24; i++) { + BODY_60_79(X[(i + 8) & 15], X[(i + 10) & 15], X[i & 15], + X[(i + 5) & 15]); + } + + c->h0 = (c->h0 + A) & 0xffffffffL; + c->h1 = (c->h1 + B) & 0xffffffffL; + c->h2 = (c->h2 + C) & 0xffffffffL; + c->h3 = (c->h3 + D) & 0xffffffffL; + c->h4 = (c->h4 + E) & 0xffffffffL; + + if (--num == 0) + break; + + A = c->h0; + B = c->h1; + C = c->h2; + D = c->h3; + E = c->h4; + + } +} +# endif + +#endif diff --git a/src/sha256.c b/src/sha256.c new file mode 100644 index 0000000..8c477e4 --- /dev/null +++ b/src/sha256.c @@ -0,0 +1,399 @@ +/* crypto/sha/sha256.c */ +/* ==================================================================== + * Copyright (c) 2004 The OpenSSL Project. All rights reserved + * according to the OpenSSL license [found in ../../LICENSE]. + * ==================================================================== + */ +# include +# include + +# include "sha.h" + +/* mem_clr.c */ +unsigned static char cleanse_ctr = 0; +static void OPENSSL_cleanse(void *ptr, size_t len) +{ + unsigned char *p = ptr; + size_t loop = len, ctr = cleanse_ctr; + while (loop--) { + *(p++) = (unsigned char)ctr; + ctr += (17 + ((size_t)p & 0xF)); + } + p = memchr(ptr, (unsigned char)ctr, len); + if (p) + ctr += (63 + (size_t)p); + cleanse_ctr = (unsigned char)ctr; +} + +# define fips_md_init(alg) fips_md_init_ctx(alg, alg) +# define fips_md_init_ctx(alg, cx) \ + int alg##_Init(cx##_CTX *c) +# define fips_cipher_abort(alg) while(0) + +fips_md_init_ctx(SHA224, SHA256) +{ + memset(c, 0, sizeof(*c)); + c->h[0] = 0xc1059ed8UL; + c->h[1] = 0x367cd507UL; + c->h[2] = 0x3070dd17UL; + c->h[3] = 0xf70e5939UL; + c->h[4] = 0xffc00b31UL; + c->h[5] = 0x68581511UL; + c->h[6] = 0x64f98fa7UL; + c->h[7] = 0xbefa4fa4UL; + c->md_len = SHA224_DIGEST_LENGTH; + return 1; +} + +fips_md_init(SHA256) +{ + memset(c, 0, sizeof(*c)); + c->h[0] = 0x6a09e667UL; + c->h[1] = 0xbb67ae85UL; + c->h[2] = 0x3c6ef372UL; + c->h[3] = 0xa54ff53aUL; + c->h[4] = 0x510e527fUL; + c->h[5] = 0x9b05688cUL; + c->h[6] = 0x1f83d9abUL; + c->h[7] = 0x5be0cd19UL; + c->md_len = SHA256_DIGEST_LENGTH; + return 1; +} + +unsigned char *SHA224(const unsigned char *d, size_t n, unsigned char *md) +{ + SHA256_CTX c; + static unsigned char m[SHA224_DIGEST_LENGTH]; + + if (md == NULL) + md = m; + SHA224_Init(&c); + SHA256_Update(&c, d, n); + SHA256_Final(md, &c); + OPENSSL_cleanse(&c, sizeof(c)); + return (md); +} + +unsigned char *SHA256(const unsigned char *d, size_t n, unsigned char *md) +{ + SHA256_CTX c; + static unsigned char m[SHA256_DIGEST_LENGTH]; + + if (md == NULL) + md = m; + SHA256_Init(&c); + SHA256_Update(&c, d, n); + SHA256_Final(md, &c); + OPENSSL_cleanse(&c, sizeof(c)); + return (md); +} + +int SHA224_Update(SHA256_CTX *c, const void *data, size_t len) +{ + return SHA256_Update(c, data, len); +} + +int SHA224_Final(unsigned char *md, SHA256_CTX *c) +{ + return SHA256_Final(md, c); +} + +# define DATA_ORDER_IS_BIG_ENDIAN + +# define HASH_LONG SHA_LONG +# define HASH_CTX SHA256_CTX +# define HASH_CBLOCK SHA_CBLOCK +/* + * Note that FIPS180-2 discusses "Truncation of the Hash Function Output." + * default: case below covers for it. It's not clear however if it's + * permitted to truncate to amount of bytes not divisible by 4. I bet not, + * but if it is, then default: case shall be extended. For reference. + * Idea behind separate cases for pre-defined lenghts is to let the + * compiler decide if it's appropriate to unroll small loops. + */ +# define HASH_MAKE_STRING(c,s) do { \ + unsigned long ll; \ + unsigned int nn; \ + switch ((c)->md_len) \ + { case SHA224_DIGEST_LENGTH: \ + for (nn=0;nnh[nn]; (void)HOST_l2c(ll,(s)); } \ + break; \ + case SHA256_DIGEST_LENGTH: \ + for (nn=0;nnh[nn]; (void)HOST_l2c(ll,(s)); } \ + break; \ + default: \ + if ((c)->md_len > SHA256_DIGEST_LENGTH) \ + return 0; \ + for (nn=0;nn<(c)->md_len/4;nn++) \ + { ll=(c)->h[nn]; (void)HOST_l2c(ll,(s)); } \ + break; \ + } \ + } while (0) + +# define HASH_UPDATE SHA256_Update +# define HASH_TRANSFORM SHA256_Transform +# define HASH_FINAL SHA256_Final +# define HASH_BLOCK_DATA_ORDER sha256_block_data_order +# ifndef SHA256_ASM +static +# endif +void sha256_block_data_order(SHA256_CTX *ctx, const void *in, size_t num); + +# include "md32_common.h" + +# ifndef SHA256_ASM +static const SHA_LONG K256[64] = { + 0x428a2f98UL, 0x71374491UL, 0xb5c0fbcfUL, 0xe9b5dba5UL, + 0x3956c25bUL, 0x59f111f1UL, 0x923f82a4UL, 0xab1c5ed5UL, + 0xd807aa98UL, 0x12835b01UL, 0x243185beUL, 0x550c7dc3UL, + 0x72be5d74UL, 0x80deb1feUL, 0x9bdc06a7UL, 0xc19bf174UL, + 0xe49b69c1UL, 0xefbe4786UL, 0x0fc19dc6UL, 0x240ca1ccUL, + 0x2de92c6fUL, 0x4a7484aaUL, 0x5cb0a9dcUL, 0x76f988daUL, + 0x983e5152UL, 0xa831c66dUL, 0xb00327c8UL, 0xbf597fc7UL, + 0xc6e00bf3UL, 0xd5a79147UL, 0x06ca6351UL, 0x14292967UL, + 0x27b70a85UL, 0x2e1b2138UL, 0x4d2c6dfcUL, 0x53380d13UL, + 0x650a7354UL, 0x766a0abbUL, 0x81c2c92eUL, 0x92722c85UL, + 0xa2bfe8a1UL, 0xa81a664bUL, 0xc24b8b70UL, 0xc76c51a3UL, + 0xd192e819UL, 0xd6990624UL, 0xf40e3585UL, 0x106aa070UL, + 0x19a4c116UL, 0x1e376c08UL, 0x2748774cUL, 0x34b0bcb5UL, + 0x391c0cb3UL, 0x4ed8aa4aUL, 0x5b9cca4fUL, 0x682e6ff3UL, + 0x748f82eeUL, 0x78a5636fUL, 0x84c87814UL, 0x8cc70208UL, + 0x90befffaUL, 0xa4506cebUL, 0xbef9a3f7UL, 0xc67178f2UL +}; + +/* + * FIPS specification refers to right rotations, while our ROTATE macro + * is left one. This is why you might notice that rotation coefficients + * differ from those observed in FIPS document by 32-N... + */ +# define Sigma0(x) (ROTATE((x),30) ^ ROTATE((x),19) ^ ROTATE((x),10)) +# define Sigma1(x) (ROTATE((x),26) ^ ROTATE((x),21) ^ ROTATE((x),7)) +# define sigma0(x) (ROTATE((x),25) ^ ROTATE((x),14) ^ ((x)>>3)) +# define sigma1(x) (ROTATE((x),15) ^ ROTATE((x),13) ^ ((x)>>10)) + +# define Ch(x,y,z) (((x) & (y)) ^ ((~(x)) & (z))) +# define Maj(x,y,z) (((x) & (y)) ^ ((x) & (z)) ^ ((y) & (z))) + +# ifdef OPENSSL_SMALL_FOOTPRINT + +static void sha256_block_data_order(SHA256_CTX *ctx, const void *in, + size_t num) +{ + unsigned MD32_REG_T a, b, c, d, e, f, g, h, s0, s1, T1, T2; + SHA_LONG X[16], l; + int i; + const unsigned char *data = in; + + while (num--) { + + a = ctx->h[0]; + b = ctx->h[1]; + c = ctx->h[2]; + d = ctx->h[3]; + e = ctx->h[4]; + f = ctx->h[5]; + g = ctx->h[6]; + h = ctx->h[7]; + + for (i = 0; i < 16; i++) { + HOST_c2l(data, l); + T1 = X[i] = l; + T1 += h + Sigma1(e) + Ch(e, f, g) + K256[i]; + T2 = Sigma0(a) + Maj(a, b, c); + h = g; + g = f; + f = e; + e = d + T1; + d = c; + c = b; + b = a; + a = T1 + T2; + } + + for (; i < 64; i++) { + s0 = X[(i + 1) & 0x0f]; + s0 = sigma0(s0); + s1 = X[(i + 14) & 0x0f]; + s1 = sigma1(s1); + + T1 = X[i & 0xf] += s0 + s1 + X[(i + 9) & 0xf]; + T1 += h + Sigma1(e) + Ch(e, f, g) + K256[i]; + T2 = Sigma0(a) + Maj(a, b, c); + h = g; + g = f; + f = e; + e = d + T1; + d = c; + c = b; + b = a; + a = T1 + T2; + } + + ctx->h[0] += a; + ctx->h[1] += b; + ctx->h[2] += c; + ctx->h[3] += d; + ctx->h[4] += e; + ctx->h[5] += f; + ctx->h[6] += g; + ctx->h[7] += h; + + } +} + +# else + +# define ROUND_00_15(i,a,b,c,d,e,f,g,h) do { \ + T1 += h + Sigma1(e) + Ch(e,f,g) + K256[i]; \ + h = Sigma0(a) + Maj(a,b,c); \ + d += T1; h += T1; } while (0) + +# define ROUND_16_63(i,a,b,c,d,e,f,g,h,X) do { \ + s0 = X[(i+1)&0x0f]; s0 = sigma0(s0); \ + s1 = X[(i+14)&0x0f]; s1 = sigma1(s1); \ + T1 = X[(i)&0x0f] += s0 + s1 + X[(i+9)&0x0f]; \ + ROUND_00_15(i,a,b,c,d,e,f,g,h); } while (0) + +static void sha256_block_data_order(SHA256_CTX *ctx, const void *in, + size_t num) +{ + unsigned MD32_REG_T a, b, c, d, e, f, g, h, s0, s1, T1; + SHA_LONG X[16]; + int i; + const unsigned char *data = in; + const union { + long one; + char little; + } is_endian = { + 1 + }; + + while (num--) { + + a = ctx->h[0]; + b = ctx->h[1]; + c = ctx->h[2]; + d = ctx->h[3]; + e = ctx->h[4]; + f = ctx->h[5]; + g = ctx->h[6]; + h = ctx->h[7]; + + if (!is_endian.little && sizeof(SHA_LONG) == 4 + && ((size_t)in % 4) == 0) { + const SHA_LONG *W = (const SHA_LONG *)data; + + T1 = X[0] = W[0]; + ROUND_00_15(0, a, b, c, d, e, f, g, h); + T1 = X[1] = W[1]; + ROUND_00_15(1, h, a, b, c, d, e, f, g); + T1 = X[2] = W[2]; + ROUND_00_15(2, g, h, a, b, c, d, e, f); + T1 = X[3] = W[3]; + ROUND_00_15(3, f, g, h, a, b, c, d, e); + T1 = X[4] = W[4]; + ROUND_00_15(4, e, f, g, h, a, b, c, d); + T1 = X[5] = W[5]; + ROUND_00_15(5, d, e, f, g, h, a, b, c); + T1 = X[6] = W[6]; + ROUND_00_15(6, c, d, e, f, g, h, a, b); + T1 = X[7] = W[7]; + ROUND_00_15(7, b, c, d, e, f, g, h, a); + T1 = X[8] = W[8]; + ROUND_00_15(8, a, b, c, d, e, f, g, h); + T1 = X[9] = W[9]; + ROUND_00_15(9, h, a, b, c, d, e, f, g); + T1 = X[10] = W[10]; + ROUND_00_15(10, g, h, a, b, c, d, e, f); + T1 = X[11] = W[11]; + ROUND_00_15(11, f, g, h, a, b, c, d, e); + T1 = X[12] = W[12]; + ROUND_00_15(12, e, f, g, h, a, b, c, d); + T1 = X[13] = W[13]; + ROUND_00_15(13, d, e, f, g, h, a, b, c); + T1 = X[14] = W[14]; + ROUND_00_15(14, c, d, e, f, g, h, a, b); + T1 = X[15] = W[15]; + ROUND_00_15(15, b, c, d, e, f, g, h, a); + + data += SHA256_CBLOCK; + } else { + SHA_LONG l; + + HOST_c2l(data, l); + T1 = X[0] = l; + ROUND_00_15(0, a, b, c, d, e, f, g, h); + HOST_c2l(data, l); + T1 = X[1] = l; + ROUND_00_15(1, h, a, b, c, d, e, f, g); + HOST_c2l(data, l); + T1 = X[2] = l; + ROUND_00_15(2, g, h, a, b, c, d, e, f); + HOST_c2l(data, l); + T1 = X[3] = l; + ROUND_00_15(3, f, g, h, a, b, c, d, e); + HOST_c2l(data, l); + T1 = X[4] = l; + ROUND_00_15(4, e, f, g, h, a, b, c, d); + HOST_c2l(data, l); + T1 = X[5] = l; + ROUND_00_15(5, d, e, f, g, h, a, b, c); + HOST_c2l(data, l); + T1 = X[6] = l; + ROUND_00_15(6, c, d, e, f, g, h, a, b); + HOST_c2l(data, l); + T1 = X[7] = l; + ROUND_00_15(7, b, c, d, e, f, g, h, a); + HOST_c2l(data, l); + T1 = X[8] = l; + ROUND_00_15(8, a, b, c, d, e, f, g, h); + HOST_c2l(data, l); + T1 = X[9] = l; + ROUND_00_15(9, h, a, b, c, d, e, f, g); + HOST_c2l(data, l); + T1 = X[10] = l; + ROUND_00_15(10, g, h, a, b, c, d, e, f); + HOST_c2l(data, l); + T1 = X[11] = l; + ROUND_00_15(11, f, g, h, a, b, c, d, e); + HOST_c2l(data, l); + T1 = X[12] = l; + ROUND_00_15(12, e, f, g, h, a, b, c, d); + HOST_c2l(data, l); + T1 = X[13] = l; + ROUND_00_15(13, d, e, f, g, h, a, b, c); + HOST_c2l(data, l); + T1 = X[14] = l; + ROUND_00_15(14, c, d, e, f, g, h, a, b); + HOST_c2l(data, l); + T1 = X[15] = l; + ROUND_00_15(15, b, c, d, e, f, g, h, a); + } + + for (i = 16; i < 64; i += 8) { + ROUND_16_63(i + 0, a, b, c, d, e, f, g, h, X); + ROUND_16_63(i + 1, h, a, b, c, d, e, f, g, X); + ROUND_16_63(i + 2, g, h, a, b, c, d, e, f, X); + ROUND_16_63(i + 3, f, g, h, a, b, c, d, e, X); + ROUND_16_63(i + 4, e, f, g, h, a, b, c, d, X); + ROUND_16_63(i + 5, d, e, f, g, h, a, b, c, X); + ROUND_16_63(i + 6, c, d, e, f, g, h, a, b, X); + ROUND_16_63(i + 7, b, c, d, e, f, g, h, a, X); + } + + ctx->h[0] += a; + ctx->h[1] += b; + ctx->h[2] += c; + ctx->h[3] += d; + ctx->h[4] += e; + ctx->h[5] += f; + ctx->h[6] += g; + ctx->h[7] += h; + + } +} + +# endif +# endif /* SHA256_ASM */ diff --git a/src/sha512.c b/src/sha512.c new file mode 100644 index 0000000..f3f4ed0 --- /dev/null +++ b/src/sha512.c @@ -0,0 +1,408 @@ +// +// sha512.c +// HomeKitSRPGMP +// +// Created by d. nye on 6/29/17. +// Copyright © 2017 Mobile Flow LLC. All rights reserved. +// + +# include +# include + +# include "sha.h" + +// https://github.com/jedisct1/libsodium/src/libsodium/crypto_hash/sha512/cp/hash_sha512_cp.c + +/*- + * Copyright 2005,2007,2009 Colin Percival + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * + */ + +#define crypto_hash_sha512_BYTES 64 +//#define NATIVE_BIG_ENDIAN 1 +typedef struct crypto_hash_sha512_state { + uint64_t state[8]; + uint64_t count[2]; + uint8_t buf[128]; +} crypto_hash_sha512_state_t; + +#define ROTR64(X, B) rotr64((X), (B)) +static inline uint64_t +rotr64(const uint64_t x, const int b) +{ + return (x >> b) | (x << (64 - b)); +} + +#define STORE64_BE(DST, W) store64_be((DST), (W)) +static inline void +store64_be(uint8_t dst[8], uint64_t w) +{ +#ifdef NATIVE_BIG_ENDIAN + memcpy(dst, &w, sizeof w); +#else + dst[7] = (uint8_t) w; w >>= 8; + dst[6] = (uint8_t) w; w >>= 8; + dst[5] = (uint8_t) w; w >>= 8; + dst[4] = (uint8_t) w; w >>= 8; + dst[3] = (uint8_t) w; w >>= 8; + dst[2] = (uint8_t) w; w >>= 8; + dst[1] = (uint8_t) w; w >>= 8; + dst[0] = (uint8_t) w; +#endif +} + +#define LOAD64_BE(SRC) load64_be(SRC) +static inline uint64_t +load64_be(const uint8_t src[8]) +{ +#ifdef NATIVE_BIG_ENDIAN + uint64_t w; + memcpy(&w, src, sizeof w); + return w; +#else + uint64_t w = (uint64_t) src[7]; + w |= (uint64_t) src[6] << 8; + w |= (uint64_t) src[5] << 16; + w |= (uint64_t) src[4] << 24; + w |= (uint64_t) src[3] << 32; + w |= (uint64_t) src[2] << 40; + w |= (uint64_t) src[1] << 48; + w |= (uint64_t) src[0] << 56; + return w; +#endif +} + +void +sodium_memzero(void *const pnt, const size_t len) +{ + volatile unsigned char *volatile pnt_ = + (volatile unsigned char *volatile) pnt; + size_t i = (size_t) 0U; + + while (i < len) { + pnt_[i++] = 0U; + } +} + +static void +be64enc_vect(unsigned char *dst, const uint64_t *src, size_t len) +{ + size_t i; + + for (i = 0; i < len / 8; i++) { + STORE64_BE(dst + i * 8, src[i]); + } +} + +static void +be64dec_vect(uint64_t *dst, const unsigned char *src, size_t len) +{ + size_t i; + uint64_t w; + + for (i = 0; i < len / 8; i++) { + dst[i] = LOAD64_BE(src + i * 8); + } +} + +static const uint64_t Krnd[80] = { + 0x428a2f98d728ae22ULL, 0x7137449123ef65cdULL, 0xb5c0fbcfec4d3b2fULL, + 0xe9b5dba58189dbbcULL, 0x3956c25bf348b538ULL, 0x59f111f1b605d019ULL, + 0x923f82a4af194f9bULL, 0xab1c5ed5da6d8118ULL, 0xd807aa98a3030242ULL, + 0x12835b0145706fbeULL, 0x243185be4ee4b28cULL, 0x550c7dc3d5ffb4e2ULL, + 0x72be5d74f27b896fULL, 0x80deb1fe3b1696b1ULL, 0x9bdc06a725c71235ULL, + 0xc19bf174cf692694ULL, 0xe49b69c19ef14ad2ULL, 0xefbe4786384f25e3ULL, + 0x0fc19dc68b8cd5b5ULL, 0x240ca1cc77ac9c65ULL, 0x2de92c6f592b0275ULL, + 0x4a7484aa6ea6e483ULL, 0x5cb0a9dcbd41fbd4ULL, 0x76f988da831153b5ULL, + 0x983e5152ee66dfabULL, 0xa831c66d2db43210ULL, 0xb00327c898fb213fULL, + 0xbf597fc7beef0ee4ULL, 0xc6e00bf33da88fc2ULL, 0xd5a79147930aa725ULL, + 0x06ca6351e003826fULL, 0x142929670a0e6e70ULL, 0x27b70a8546d22ffcULL, + 0x2e1b21385c26c926ULL, 0x4d2c6dfc5ac42aedULL, 0x53380d139d95b3dfULL, + 0x650a73548baf63deULL, 0x766a0abb3c77b2a8ULL, 0x81c2c92e47edaee6ULL, + 0x92722c851482353bULL, 0xa2bfe8a14cf10364ULL, 0xa81a664bbc423001ULL, + 0xc24b8b70d0f89791ULL, 0xc76c51a30654be30ULL, 0xd192e819d6ef5218ULL, + 0xd69906245565a910ULL, 0xf40e35855771202aULL, 0x106aa07032bbd1b8ULL, + 0x19a4c116b8d2d0c8ULL, 0x1e376c085141ab53ULL, 0x2748774cdf8eeb99ULL, + 0x34b0bcb5e19b48a8ULL, 0x391c0cb3c5c95a63ULL, 0x4ed8aa4ae3418acbULL, + 0x5b9cca4f7763e373ULL, 0x682e6ff3d6b2b8a3ULL, 0x748f82ee5defb2fcULL, + 0x78a5636f43172f60ULL, 0x84c87814a1f0ab72ULL, 0x8cc702081a6439ecULL, + 0x90befffa23631e28ULL, 0xa4506cebde82bde9ULL, 0xbef9a3f7b2c67915ULL, + 0xc67178f2e372532bULL, 0xca273eceea26619cULL, 0xd186b8c721c0c207ULL, + 0xeada7dd6cde0eb1eULL, 0xf57d4f7fee6ed178ULL, 0x06f067aa72176fbaULL, + 0x0a637dc5a2c898a6ULL, 0x113f9804bef90daeULL, 0x1b710b35131c471bULL, + 0x28db77f523047d84ULL, 0x32caab7b40c72493ULL, 0x3c9ebe0a15c9bebcULL, + 0x431d67c49c100d4cULL, 0x4cc5d4becb3e42b6ULL, 0x597f299cfc657e2aULL, + 0x5fcb6fab3ad6faecULL, 0x6c44198c4a475817ULL +}; + +#define Ch(x, y, z) ((x & (y ^ z)) ^ z) +#define Maj(x, y, z) ((x & (y | z)) | (y & z)) +#define SHR(x, n) (x >> n) +#define ROTR(x, n) ROTR64(x, n) +#define S0(x) (ROTR(x, 28) ^ ROTR(x, 34) ^ ROTR(x, 39)) +#define S1(x) (ROTR(x, 14) ^ ROTR(x, 18) ^ ROTR(x, 41)) +#define s0(x) (ROTR(x, 1) ^ ROTR(x, 8) ^ SHR(x, 7)) +#define s1(x) (ROTR(x, 19) ^ ROTR(x, 61) ^ SHR(x, 6)) + +#define RND(a, b, c, d, e, f, g, h, k) \ +h += S1(e) + Ch(e, f, g) + k; \ +d += h; \ +h += S0(a) + Maj(a, b, c); + +#define RNDr(S, W, i, ii) \ +RND(S[(80 - i) % 8], S[(81 - i) % 8], S[(82 - i) % 8], S[(83 - i) % 8], \ +S[(84 - i) % 8], S[(85 - i) % 8], S[(86 - i) % 8], S[(87 - i) % 8], \ +W[i + ii] + Krnd[i + ii]) + +#define MSCH(W, ii, i) \ +W[i + ii + 16] = \ +s1(W[i + ii + 14]) + W[i + ii + 9] + s0(W[i + ii + 1]) + W[i + ii] + + +static void +SHA512_Transform_Internal(uint64_t *state, const uint8_t block[128], uint64_t W[80], + uint64_t S[8]) +{ + int i; + + be64dec_vect(W, block, 128); + memcpy(S, state, 64); + for (i = 0; i < 80; i += 16) { + RNDr(S, W, 0, i); + RNDr(S, W, 1, i); + RNDr(S, W, 2, i); + RNDr(S, W, 3, i); + RNDr(S, W, 4, i); + RNDr(S, W, 5, i); + RNDr(S, W, 6, i); + RNDr(S, W, 7, i); + RNDr(S, W, 8, i); + RNDr(S, W, 9, i); + RNDr(S, W, 10, i); + RNDr(S, W, 11, i); + RNDr(S, W, 12, i); + RNDr(S, W, 13, i); + RNDr(S, W, 14, i); + RNDr(S, W, 15, i); + if (i == 64) { + break; + } + MSCH(W, 0, i); + MSCH(W, 1, i); + MSCH(W, 2, i); + MSCH(W, 3, i); + MSCH(W, 4, i); + MSCH(W, 5, i); + MSCH(W, 6, i); + MSCH(W, 7, i); + MSCH(W, 8, i); + MSCH(W, 9, i); + MSCH(W, 10, i); + MSCH(W, 11, i); + MSCH(W, 12, i); + MSCH(W, 13, i); + MSCH(W, 14, i); + MSCH(W, 15, i); + } + for (i = 0; i < 8; i++) { + state[i] += S[i]; + } +} + +static const uint8_t PAD[128] = { + 0x80, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 +}; + +static void +SHA512_Pad(crypto_hash_sha512_state_t *state, uint64_t tmp64[80 + 8]) +{ + uint64_t r; + uint64_t i; + + r = (state->count[1] >> 3) & 0x7f; + if (r < 112) { + for (i = 0; i < 112 - r; i++) { + state->buf[r + i] = PAD[i]; + } + } else { + for (i = 0; i < 128 - r; i++) { + state->buf[r + i] = PAD[i]; + } + SHA512_Transform_Internal(state->state, state->buf, &tmp64[0], &tmp64[80]); + memset(&state->buf[0], 0, 112); + } + be64enc_vect(&state->buf[112], state->count, 16); + SHA512_Transform_Internal(state->state, state->buf, &tmp64[0], &tmp64[80]); +} + +int +crypto_hash_sha512_init(crypto_hash_sha512_state_t *state) +{ + static const uint64_t sha512_initial_state[8] = { + 0x6a09e667f3bcc908ULL, 0xbb67ae8584caa73bULL, 0x3c6ef372fe94f82bULL, + 0xa54ff53a5f1d36f1ULL, 0x510e527fade682d1ULL, 0x9b05688c2b3e6c1fULL, + 0x1f83d9abfb41bd6bULL, 0x5be0cd19137e2179ULL + }; + + state->count[0] = state->count[1] = (uint64_t) 0U; + memcpy(state->state, sha512_initial_state, sizeof sha512_initial_state); + + return 0; +} + + +int +crypto_hash_sha512_update(crypto_hash_sha512_state_t *state, + const unsigned char *in, unsigned long long inlen) +{ + uint64_t tmp64[80 + 8]; + uint64_t bitlen[2]; + unsigned long long i; + unsigned long long r; + + if (inlen <= 0U) { + return 0; + } + r = (unsigned long long) ((state->count[1] >> 3) & 0x7f); + + bitlen[1] = ((uint64_t) inlen) << 3; + bitlen[0] = ((uint64_t) inlen) >> 61; + if ((state->count[1] += bitlen[1]) < bitlen[1]) { + state->count[0]++; + } + state->count[0] += bitlen[0]; + if (inlen < 128 - r) { + for (i = 0; i < inlen; i++) { + state->buf[r + i] = in[i]; + } + return 0; + } + for (i = 0; i < 128 - r; i++) { + state->buf[r + i] = in[i]; + } + SHA512_Transform_Internal(state->state, state->buf, &tmp64[0], &tmp64[80]); + in += 128 - r; + inlen -= 128 - r; + + while (inlen >= 128) { + SHA512_Transform_Internal(state->state, in, &tmp64[0], &tmp64[80]); + in += 128; + inlen -= 128; + } + inlen &= 127; + for (i = 0; i < inlen; i++) { + state->buf[i] = in[i]; + } + sodium_memzero((void *) tmp64, sizeof tmp64); + + return 0; +} + +int +crypto_hash_sha512_final(crypto_hash_sha512_state_t *state, unsigned char *out) +{ + uint64_t tmp64[80 + 8]; + + SHA512_Pad(state, tmp64); + be64enc_vect(out, state->state, 64); + sodium_memzero((void *) tmp64, sizeof tmp64); + sodium_memzero((void *) state, sizeof *state); + + return 0; +} + +// Interface + +int +SHA512_Init(SHA512_CTX *c) +{ + + crypto_hash_sha512_state_t *state = (crypto_hash_sha512_state_t *)c; + + crypto_hash_sha512_init(state); + + return 0; +} + +int +SHA512_Update(SHA512_CTX *c, const void *data, size_t len) +{ + crypto_hash_sha512_state_t *state = (crypto_hash_sha512_state_t *)c; + + return crypto_hash_sha512_update(state, (const unsigned char *)data, (unsigned long long)len); + +} + +int +SHA512_Final(unsigned char *md, SHA512_CTX *c) +{ + crypto_hash_sha512_state_t *state = (crypto_hash_sha512_state_t *)c; + + return crypto_hash_sha512_final(state, md); + +} + +unsigned char*SHA512(const unsigned char *d, size_t n, unsigned char *md) +{ + SHA512_CTX c; + static unsigned char m[SHA512_DIGEST_LENGTH]; + + if (md == NULL) + md = m; + SHA512_Init(&c); + SHA512_Update(&c, d, n); + SHA512_Final(md, &c); +// OPENSSL_cleanse(&c, sizeof(c)); + return (md); + +} + +void +SHA512_Transform(SHA512_CTX *c, const unsigned char *data) +{ + + // Change block order... + +} + +/* +int SRPClient::crypto_hash_sha512(unsigned char *out, const unsigned char *in, unsigned long long inlen) +{ + + crypto_hash_sha512_state state; + + crypto_hash_sha512_init(&state); + crypto_hash_sha512_update(&state, in, inlen); + crypto_hash_sha512_final(&state, out); + + return 0; + +} +*/ diff --git a/src/srp.c b/src/srp.c new file mode 100644 index 0000000..2bf818c --- /dev/null +++ b/src/srp.c @@ -0,0 +1,1047 @@ +/* + * Secure Remote Password 6a implementation + * https://github.com/est31/csrp-gmp + * + * The MIT License (MIT) + * + * Copyright (c) 2010, 2013 Tom Cocagne, 2015 est31 + * + * Permission is hereby granted, free of charge, to any person obtaining a copy of + * this software and associated documentation files (the "Software"), to deal in + * the Software without restriction, including without limitation the rights to + * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies + * of the Software, and to permit persons to whom the Software is furnished to do + * so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ + +// clang-format off +#ifdef WIN32 + #include + #include +#else + #include +#endif +// clang-format on + +#include +#include +#include + +#include "mini-gmp.h" +#include "sha.h" + +#include "srp.h" +#define CSRP_USE_SHA1 +#define CSRP_USE_SHA256 + +#define srp_dbg_data(data, datalen, prevtext) ; +/*void srp_dbg_data(unsigned char * data, size_t datalen, char * prevtext) +{ + printf(prevtext); + size_t i; + for (i = 0; i < datalen; i++) + { + printf("%02X", data[i]); + } + printf("\n"); +}*/ + +static int g_initialized = 0; + +#define RAND_BUFF_MAX 128 +static unsigned int g_rand_idx; +static unsigned char g_rand_buff[RAND_BUFF_MAX]; + +void *(*srp_alloc)(size_t) = &malloc; +void *(*srp_realloc)(void *, size_t) = &realloc; +void (*srp_free)(void *) = &free; + +// clang-format off +void srp_set_memory_functions( + void *(*new_srp_alloc)(size_t), + void *(*new_srp_realloc)(void *, size_t), + void (*new_srp_free)(void *)) +{ + srp_alloc = new_srp_alloc; + srp_realloc = new_srp_realloc; + srp_free = new_srp_free; +} +// clang-format on + +typedef struct { + mpz_t N; + mpz_t g; +} NGConstant; + +struct NGHex { + const char *n_hex; + const char *g_hex; +}; + +/* All constants here were pulled from Appendix A of RFC 5054 */ +static struct NGHex global_Ng_constants[] = { + {/* 1024 */ + "EEAF0AB9ADB38DD69C33F80AFA8FC5E86072618775FF3C0B9EA2314C" + "9C256576D674DF7496EA81D3383B4813D692C6E0E0D5D8E250B98BE4" + "8E495C1D6089DAD15DC7D7B46154D6B6CE8EF4AD69B15D4982559B29" + "7BCF1885C529F566660E57EC68EDBC3C05726CC02FD4CBF4976EAA9A" + "FD5138FE8376435B9FC61D2FC0EB06E3", + "2"}, + {/* 2048 */ + "AC6BDB41324A9A9BF166DE5E1389582FAF72B6651987EE07FC319294" + "3DB56050A37329CBB4A099ED8193E0757767A13DD52312AB4B03310D" + "CD7F48A9DA04FD50E8083969EDB767B0CF6095179A163AB3661A05FB" + "D5FAAAE82918A9962F0B93B855F97993EC975EEAA80D740ADBF4FF74" + "7359D041D5C33EA71D281E446B14773BCA97B43A23FB801676BD207A" + "436C6481F1D2B9078717461A5B9D32E688F87748544523B524B0D57D" + "5EA77A2775D2ECFA032CFBDBF52FB3786160279004E57AE6AF874E73" + "03CE53299CCC041C7BC308D82A5698F3A8D0C38271AE35F8E9DBFBB6" + "94B5C803D89F7AE435DE236D525F54759B65E372FCD68EF20FA7111F" + "9E4AFF73", + "2"}, + {/* 3072 */ + "FFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74" + "020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F1437" + "4FE1356D6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7ED" + "EE386BFB5A899FA5AE9F24117C4B1FE649286651ECE45B3DC2007CB8A163BF05" + "98DA48361C55D39A69163FA8FD24CF5F83655D23DCA3AD961C62F356208552BB" + "9ED529077096966D670C354E4ABC9804F1746C08CA18217C32905E462E36CE3B" + "E39E772C180E86039B2783A2EC07A28FB5C55DF06F4C52C9DE2BCBF695581718" + "3995497CEA956AE515D2261898FA051015728E5A8AAAC42DAD33170D04507A33" + "A85521ABDF1CBA64ECFB850458DBEF0A8AEA71575D060C7DB3970F85A6E1E4C7" + "ABF5AE8CDB0933D71E8C94E04A25619DCEE3D2261AD2EE6BF12FFA06D98A0864" + "D87602733EC86A64521F2B18177B200CBBE117577A615D6C770988C0BAD946E2" + "08E24FA074E5AB3143DB5BFCE0FD108E4B82D120A93AD2CAFFFFFFFFFFFFFFFF", + "5"}, + {/* 4096 */ + "FFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD129024E08" + "8A67CC74020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B" + "302B0A6DF25F14374FE1356D6D51C245E485B576625E7EC6F44C42E9" + "A637ED6B0BFF5CB6F406B7EDEE386BFB5A899FA5AE9F24117C4B1FE6" + "49286651ECE45B3DC2007CB8A163BF0598DA48361C55D39A69163FA8" + "FD24CF5F83655D23DCA3AD961C62F356208552BB9ED529077096966D" + "670C354E4ABC9804F1746C08CA18217C32905E462E36CE3BE39E772C" + "180E86039B2783A2EC07A28FB5C55DF06F4C52C9DE2BCBF695581718" + "3995497CEA956AE515D2261898FA051015728E5A8AAAC42DAD33170D" + "04507A33A85521ABDF1CBA64ECFB850458DBEF0A8AEA71575D060C7D" + "B3970F85A6E1E4C7ABF5AE8CDB0933D71E8C94E04A25619DCEE3D226" + "1AD2EE6BF12FFA06D98A0864D87602733EC86A64521F2B18177B200C" + "BBE117577A615D6C770988C0BAD946E208E24FA074E5AB3143DB5BFC" + "E0FD108E4B82D120A92108011A723C12A787E6D788719A10BDBA5B26" + "99C327186AF4E23C1A946834B6150BDA2583E9CA2AD44CE8DBBBC2DB" + "04DE8EF92E8EFC141FBECAA6287C59474E6BC05D99B2964FA090C3A2" + "233BA186515BE7ED1F612970CEE2D7AFB81BDD762170481CD0069127" + "D5B05AA993B4EA988D8FDDC186FFB7DC90A6C08F4DF435C934063199" + "FFFFFFFFFFFFFFFF", + "5"}, + {/* 8192 */ + "FFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD129024E08" + "8A67CC74020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B" + "302B0A6DF25F14374FE1356D6D51C245E485B576625E7EC6F44C42E9" + "A637ED6B0BFF5CB6F406B7EDEE386BFB5A899FA5AE9F24117C4B1FE6" + "49286651ECE45B3DC2007CB8A163BF0598DA48361C55D39A69163FA8" + "FD24CF5F83655D23DCA3AD961C62F356208552BB9ED529077096966D" + "670C354E4ABC9804F1746C08CA18217C32905E462E36CE3BE39E772C" + "180E86039B2783A2EC07A28FB5C55DF06F4C52C9DE2BCBF695581718" + "3995497CEA956AE515D2261898FA051015728E5A8AAAC42DAD33170D" + "04507A33A85521ABDF1CBA64ECFB850458DBEF0A8AEA71575D060C7D" + "B3970F85A6E1E4C7ABF5AE8CDB0933D71E8C94E04A25619DCEE3D226" + "1AD2EE6BF12FFA06D98A0864D87602733EC86A64521F2B18177B200C" + "BBE117577A615D6C770988C0BAD946E208E24FA074E5AB3143DB5BFC" + "E0FD108E4B82D120A92108011A723C12A787E6D788719A10BDBA5B26" + "99C327186AF4E23C1A946834B6150BDA2583E9CA2AD44CE8DBBBC2DB" + "04DE8EF92E8EFC141FBECAA6287C59474E6BC05D99B2964FA090C3A2" + "233BA186515BE7ED1F612970CEE2D7AFB81BDD762170481CD0069127" + "D5B05AA993B4EA988D8FDDC186FFB7DC90A6C08F4DF435C934028492" + "36C3FAB4D27C7026C1D4DCB2602646DEC9751E763DBA37BDF8FF9406" + "AD9E530EE5DB382F413001AEB06A53ED9027D831179727B0865A8918" + "DA3EDBEBCF9B14ED44CE6CBACED4BB1BDB7F1447E6CC254B33205151" + "2BD7AF426FB8F401378CD2BF5983CA01C64B92ECF032EA15D1721D03" + "F482D7CE6E74FEF6D55E702F46980C82B5A84031900B1C9E59E7C97F" + "BEC7E8F323A97A7E36CC88BE0F1D45B7FF585AC54BD407B22B4154AA" + "CC8F6D7EBF48E1D814CC5ED20F8037E0A79715EEF29BE32806A1D58B" + "B7C5DA76F550AA3D8A1FBFF0EB19CCB1A313D55CDA56C9EC2EF29632" + "387FE8D76E3C0468043E8F663F4860EE12BF2D5B0B7474D6E694F91E" + "6DBE115974A3926F12FEE5E438777CB6A932DF8CD8BEC4D073B931BA" + "3BC832B68D9DD300741FA7BF8AFC47ED2576F6936BA424663AAB639C" + "5AE4F5683423B4742BF1C978238F16CBE39D652DE3FDB8BEFC848AD9" + "22222E04A4037C0713EB57A81A23F0C73473FC646CEA306B4BCBC886" + "2F8385DDFA9D4B7FA2C087E879683303ED5BDD3A062B3CF5B3A278A6" + "6D2A13F83F44F82DDF310EE074AB6A364597E899A0255DC164F31CC5" + "0846851DF9AB48195DED7EA1B1D510BD7EE74D73FAF36BC31ECFA268" + "359046F4EB879F924009438B481C6CD7889A002ED5EE382BC9190DA6" + "FC026E479558E4475677E9AA9E3050E2765694DFC81F56E880B96E71" + "60C980DD98EDD3DFFFFFFFFFFFFFFFFF", + "13"}, + {0, 0} /* null sentinel */ +}; + +static void delete_ng(NGConstant *ng) +{ + if (ng) { + mpz_clear(ng->N); + mpz_clear(ng->g); + srp_free(ng); + } +} + +static NGConstant *new_ng(SRP_NGType ng_type, const char *n_hex, const char *g_hex) +{ + NGConstant *ng = (NGConstant *)srp_alloc(sizeof(NGConstant)); + + if (!ng) return 0; + + mpz_init(ng->N); + mpz_init(ng->g); + + if (ng_type != SRP_NG_CUSTOM) { + n_hex = global_Ng_constants[ng_type].n_hex; + g_hex = global_Ng_constants[ng_type].g_hex; + } + + int rv = 0; + rv = mpz_set_str(ng->N, n_hex, 16); + rv = rv | mpz_set_str(ng->g, g_hex, 16); + + if (rv) { + delete_ng(ng); + return 0; + } + + return ng; +} + +typedef union { + SHA_CTX sha; + SHA256_CTX sha256; + SHA512_CTX sha512; +} HashCTX; + +struct SRPVerifier { + SRP_HashAlgorithm hash_alg; + NGConstant *ng; + + char *username; + unsigned char *bytes_B; + int authenticated; + + unsigned char M[SHA512_DIGEST_LENGTH]; + unsigned char H_AMK[SHA512_DIGEST_LENGTH]; + unsigned char session_key[SHA512_DIGEST_LENGTH]; +}; + +struct SRPUser { + SRP_HashAlgorithm hash_alg; + NGConstant *ng; + + mpz_t a; + mpz_t A; + mpz_t S; + + unsigned char *bytes_A; + int authenticated; + + char *username; + char *username_verifier; + unsigned char *password; + size_t password_len; + + unsigned char M[SHA512_DIGEST_LENGTH]; + unsigned char H_AMK[SHA512_DIGEST_LENGTH]; + unsigned char session_key[SHA512_DIGEST_LENGTH]; +}; + +// clang-format off +static int hash_init(SRP_HashAlgorithm alg, HashCTX *c) +{ + switch (alg) { +#ifdef CSRP_USE_SHA1 + case SRP_SHA1: return SHA1_Init(&c->sha); +#endif + /* + case SRP_SHA224: return SHA224_Init(&c->sha256); + */ +#ifdef CSRP_USE_SHA256 + case SRP_SHA256: return SHA256_Init(&c->sha256); +#endif + /* + case SRP_SHA384: return SHA384_Init(&c->sha512); + */ + case SRP_SHA512: return SHA512_Init(&c->sha512); + + default: return -1; + }; +} +static int hash_update( SRP_HashAlgorithm alg, HashCTX *c, const void *data, size_t len ) +{ + switch (alg) { +#ifdef CSRP_USE_SHA1 + case SRP_SHA1: return SHA1_Update(&c->sha, data, len); +#endif + /* + case SRP_SHA224: return SHA224_Update(&c->sha256, data, len); + */ +#ifdef CSRP_USE_SHA256 + case SRP_SHA256: return SHA256_Update(&c->sha256, data, len); +#endif + /* + case SRP_SHA384: return SHA384_Update(&c->sha512, data, len); + */ + case SRP_SHA512: return SHA512_Update(&c->sha512, data, len); + + default: return -1; + }; +} +static int hash_final( SRP_HashAlgorithm alg, HashCTX *c, unsigned char *md ) +{ + switch (alg) { +#ifdef CSRP_USE_SHA1 + case SRP_SHA1: return SHA1_Final(md, &c->sha); +#endif + /* + case SRP_SHA224: return SHA224_Final(md, &c->sha256); + */ +#ifdef CSRP_USE_SHA256 + case SRP_SHA256: return SHA256_Final(md, &c->sha256); +#endif + /* + case SRP_SHA384: return SHA384_Final(md, &c->sha512); + */ + case SRP_SHA512: return SHA512_Final(md, &c->sha512); + + default: return -1; + }; +} +static unsigned char *hash(SRP_HashAlgorithm alg, const unsigned char *d, size_t n, unsigned char *md) +{ + switch (alg) { +#ifdef CSRP_USE_SHA1 + case SRP_SHA1: return SHA1(d, n, md); +#endif + /* + case SRP_SHA224: return SHA224( d, n, md ); + */ +#ifdef CSRP_USE_SHA256 + case SRP_SHA256: return SHA256(d, n, md); +#endif + /* + case SRP_SHA384: return SHA384( d, n, md ); + */ + case SRP_SHA512: return SHA512( d, n, md ); + + default: return 0; + }; +} +static size_t hash_length(SRP_HashAlgorithm alg) +{ + switch (alg) { +#ifdef CSRP_USE_SHA1 + case SRP_SHA1: return SHA_DIGEST_LENGTH; +#endif + /* + case SRP_SHA224: return SHA224_DIGEST_LENGTH; + */ +#ifdef CSRP_USE_SHA256 + case SRP_SHA256: return SHA256_DIGEST_LENGTH; +#endif + /* + case SRP_SHA384: return SHA384_DIGEST_LENGTH; + */ + case SRP_SHA512: return SHA512_DIGEST_LENGTH; + + default: return -1; + }; +} +// clang-format on + +inline static int mpz_num_bytes(const mpz_t op) +{ + return (mpz_sizeinbase(op, 2) + 7) / 8; +} + +inline static void mpz_to_bin(const mpz_t op, unsigned char *to) +{ + mpz_export(to, NULL, 1, 1, 1, 0, op); +} + +inline static void mpz_from_bin(const unsigned char *s, size_t len, mpz_t ret) +{ + mpz_import(ret, len, 1, 1, 1, 0, s); +} + +// set op to (op1 * op2) mod d, using tmp for the calculation +inline static void mpz_mulm( + mpz_t op, const mpz_t op1, const mpz_t op2, const mpz_t d, mpz_t tmp) +{ + mpz_mul(tmp, op1, op2); + mpz_mod(op, tmp, d); +} + +// set op to (op1 + op2) mod d, using tmp for the calculation +inline static void mpz_addm( + mpz_t op, const mpz_t op1, const mpz_t op2, const mpz_t d, mpz_t tmp) +{ + mpz_add(tmp, op1, op2); + mpz_mod(op, tmp, d); +} + +// set op to (op1 - op2) mod d, using tmp for the calculation +inline static void mpz_subm( + mpz_t op, const mpz_t op1, const mpz_t op2, const mpz_t d, mpz_t tmp) +{ + mpz_sub(tmp, op1, op2); + mpz_mod(op, tmp, d); +} + +static SRP_Result H_nn( + mpz_t result, SRP_HashAlgorithm alg, const mpz_t N, const mpz_t n1, const mpz_t n2) +{ + unsigned char buff[SHA512_DIGEST_LENGTH]; + size_t len_N = mpz_num_bytes(N); + size_t len_n1 = mpz_num_bytes(n1); + size_t len_n2 = mpz_num_bytes(n2); + size_t nbytes = len_N + len_N; + unsigned char *bin = (unsigned char *)srp_alloc(nbytes); + if (!bin) return SRP_ERR; + if (len_n1 > len_N || len_n2 > len_N) { + srp_free(bin); + return SRP_ERR; + } + memset(bin, 0, nbytes); + mpz_to_bin(n1, bin + (len_N - len_n1)); + mpz_to_bin(n2, bin + (len_N + len_N - len_n2)); + hash(alg, bin, nbytes, buff); + srp_free(bin); + mpz_from_bin(buff, hash_length(alg), result); + return SRP_OK; +} + +static SRP_Result H_ns(mpz_t result, SRP_HashAlgorithm alg, const unsigned char *n, + size_t len_n, const unsigned char *bytes, size_t len_bytes) +{ + unsigned char buff[SHA512_DIGEST_LENGTH]; + size_t nbytes = len_n + len_bytes; + unsigned char *bin = (unsigned char *)srp_alloc(nbytes); + if (!bin) return SRP_ERR; + memcpy(bin, n, len_n); + memcpy(bin + len_n, bytes, len_bytes); + hash(alg, bin, nbytes, buff); + srp_free(bin); + mpz_from_bin(buff, hash_length(alg), result); + return SRP_OK; +} + +static int calculate_x(mpz_t result, SRP_HashAlgorithm alg, const unsigned char *salt, + size_t salt_len, const char *username, const unsigned char *password, + size_t password_len) +{ + unsigned char ucp_hash[SHA512_DIGEST_LENGTH]; + HashCTX ctx; + hash_init(alg, &ctx); + + srp_dbg_data((char *)username, strlen(username), "Username for x: "); + srp_dbg_data((char *)password, password_len, "Password for x: "); + hash_update(alg, &ctx, username, strlen(username)); + hash_update(alg, &ctx, ":", 1); + hash_update(alg, &ctx, password, password_len); + + hash_final(alg, &ctx, ucp_hash); + + return H_ns(result, alg, salt, salt_len, ucp_hash, hash_length(alg)); +} + +static SRP_Result update_hash_n(SRP_HashAlgorithm alg, HashCTX *ctx, const mpz_t n) +{ + size_t len = mpz_num_bytes(n); + unsigned char *n_bytes = (unsigned char *)srp_alloc(len); + if (!n_bytes) return SRP_ERR; + mpz_to_bin(n, n_bytes); + hash_update(alg, ctx, n_bytes, len); + srp_free(n_bytes); + return SRP_OK; +} + +static SRP_Result hash_num(SRP_HashAlgorithm alg, const mpz_t n, unsigned char *dest) +{ + int nbytes = mpz_num_bytes(n); + unsigned char *bin = (unsigned char *)srp_alloc(nbytes); + if (!bin) return SRP_ERR; + mpz_to_bin(n, bin); + hash(alg, bin, nbytes, dest); + srp_free(bin); + return SRP_OK; +} + +static SRP_Result calculate_M(SRP_HashAlgorithm alg, NGConstant *ng, unsigned char *dest, + const char *I, const unsigned char *s_bytes, size_t s_len, const mpz_t A, + const mpz_t B, const unsigned char *K) +{ + unsigned char H_N[SHA512_DIGEST_LENGTH]; + unsigned char H_g[SHA512_DIGEST_LENGTH]; + unsigned char H_I[SHA512_DIGEST_LENGTH]; + unsigned char H_xor[SHA512_DIGEST_LENGTH]; + HashCTX ctx; + size_t i = 0; + size_t hash_len = hash_length(alg); + + if (!hash_num(alg, ng->N, H_N)) return SRP_ERR; + if (!hash_num(alg, ng->g, H_g)) return SRP_ERR; + + hash(alg, (const unsigned char *)I, strlen(I), H_I); + + for (i = 0; i < hash_len; i++) + H_xor[i] = H_N[i] ^ H_g[i]; + + hash_init(alg, &ctx); + + hash_update(alg, &ctx, H_xor, hash_len); + hash_update(alg, &ctx, H_I, hash_len); + hash_update(alg, &ctx, s_bytes, s_len); + if (!update_hash_n(alg, &ctx, A)) return SRP_ERR; + if (!update_hash_n(alg, &ctx, B)) return SRP_ERR; + hash_update(alg, &ctx, K, hash_len); + + hash_final(alg, &ctx, dest); + return SRP_OK; +} + +static SRP_Result calculate_H_AMK(SRP_HashAlgorithm alg, unsigned char *dest, + const mpz_t A, const unsigned char *M, const unsigned char *K) +{ + HashCTX ctx; + + hash_init(alg, &ctx); + + if (!update_hash_n(alg, &ctx, A)) return SRP_ERR; + hash_update(alg, &ctx, M, hash_length(alg)); + hash_update(alg, &ctx, K, hash_length(alg)); + + hash_final(alg, &ctx, dest); + return SRP_OK; +} + +static SRP_Result fill_buff() +{ + g_rand_idx = 0; + +#ifdef WIN32 + HCRYPTPROV wctx; +#else + FILE *fp = 0; +#endif + +#ifdef WIN32 + + if (!CryptAcquireContext(&wctx, NULL, NULL, PROV_RSA_FULL, CRYPT_VERIFYCONTEXT)) + return SRP_ERR; + if (!CryptGenRandom(wctx, sizeof(g_rand_buff), (BYTE *)g_rand_buff)) return SRP_ERR; + if (!CryptReleaseContext(wctx, 0)) return SRP_ERR; + +#else + fp = fopen("/dev/urandom", "r"); + + if (!fp) return SRP_ERR; + + if (fread(g_rand_buff, sizeof(g_rand_buff), 1, fp) != 1) return SRP_ERR; + if (fclose(fp)) return SRP_ERR; +#endif + return SRP_OK; +} + +static SRP_Result mpz_fill_random(mpz_t num) +{ + // was call: BN_rand(num, 256, -1, 0); + if (RAND_BUFF_MAX - g_rand_idx < 32) + if (fill_buff() != SRP_OK) return SRP_ERR; + mpz_from_bin((const unsigned char *)(&g_rand_buff[g_rand_idx]), 32, num); + g_rand_idx += 32; + return SRP_OK; +} + +static SRP_Result init_random() +{ + if (g_initialized) return SRP_OK; + SRP_Result ret = fill_buff(); + g_initialized = (ret == SRP_OK); + return ret; +} + +#define srp_dbg_num(num, text) ; +/*void srp_dbg_num(mpz_t num, char * prevtext) +{ + int len_num = mpz_num_bytes(num); + char *bytes_num = (char*) srp_alloc(len_num); + mpz_to_bin(num, (unsigned char *) bytes_num); + srp_dbg_data(bytes_num, len_num, prevtext); + srp_free(bytes_num); + +}*/ + +/*********************************************************************************************************** + * + * Exported Functions + * + ***********************************************************************************************************/ + +// clang-format off +SRP_Result srp_create_salted_verification_key( SRP_HashAlgorithm alg, + SRP_NGType ng_type, const char *username_for_verifier, + const unsigned char *password, size_t len_password, + unsigned char **bytes_s, size_t *len_s, + unsigned char **bytes_v, size_t *len_v, + const char *n_hex, const char *g_hex ) +{ + SRP_Result ret = SRP_OK; + + mpz_t v; mpz_init(v); + mpz_t x; mpz_init(x); + // clang-format on + + NGConstant *ng = new_ng(ng_type, n_hex, g_hex); + + if (!ng) goto error_and_exit; + + if (init_random() != SRP_OK) /* Only happens once */ + goto error_and_exit; + + if (*bytes_s == NULL) { + size_t size_to_fill = 16; + *len_s = size_to_fill; + if (RAND_BUFF_MAX - g_rand_idx < size_to_fill) + if (fill_buff() != SRP_OK) goto error_and_exit; + *bytes_s = (unsigned char *)srp_alloc(size_to_fill); + if (!*bytes_s) goto error_and_exit; + memcpy(*bytes_s, &g_rand_buff + g_rand_idx, size_to_fill); + g_rand_idx += size_to_fill; + } + + if (!calculate_x( + x, alg, *bytes_s, *len_s, username_for_verifier, password, len_password)) + goto error_and_exit; + + srp_dbg_num(x, "Server calculated x: "); + + mpz_powm(v, ng->g, x, ng->N); + + *len_v = mpz_num_bytes(v); + + *bytes_v = (unsigned char *)srp_alloc(*len_v); + + if (!*bytes_v) goto error_and_exit; + + mpz_to_bin(v, *bytes_v); + +cleanup_and_exit: + delete_ng(ng); + mpz_clear(v); + mpz_clear(x); + return ret; +error_and_exit: + ret = SRP_ERR; + goto cleanup_and_exit; +} + +// clang-format off + +/* Out: bytes_B, len_B. + * + * On failure, bytes_B will be set to NULL and len_B will be set to 0 + */ +struct SRPVerifier *srp_verifier_new(SRP_HashAlgorithm alg, + SRP_NGType ng_type, const char *username, + const unsigned char *bytes_s, size_t len_s, + const unsigned char *bytes_v, size_t len_v, + const unsigned char *bytes_A, size_t len_A, + const unsigned char *bytes_b, size_t len_b, + unsigned char **bytes_B, size_t *len_B, + const char *n_hex, const char *g_hex ) +{ + mpz_t v; mpz_init(v); mpz_from_bin(bytes_v, len_v, v); + mpz_t A; mpz_init(A); mpz_from_bin(bytes_A, len_A, A); + mpz_t u; mpz_init(u); + mpz_t B; mpz_init(B); + mpz_t S; mpz_init(S); + mpz_t b; mpz_init(b); + mpz_t k; mpz_init(k); + mpz_t tmp1; mpz_init(tmp1); + mpz_t tmp2; mpz_init(tmp2); + mpz_t tmp3; mpz_init(tmp3); + // clang-format on + size_t ulen = strlen(username) + 1; + NGConstant *ng = new_ng(ng_type, n_hex, g_hex); + struct SRPVerifier *ver = 0; + + *len_B = 0; + *bytes_B = 0; + + if (!ng) goto cleanup_and_exit; + + ver = (struct SRPVerifier *)srp_alloc(sizeof(struct SRPVerifier)); + + if (!ver) goto cleanup_and_exit; + + if (init_random() != SRP_OK) { /* Only happens once */ + srp_free(ver); + ver = 0; + goto cleanup_and_exit; + } + + ver->username = (char *)srp_alloc(ulen); + ver->hash_alg = alg; + ver->ng = ng; + + if (!ver->username) { + srp_free(ver); + ver = 0; + goto cleanup_and_exit; + } + + memcpy((char *)ver->username, username, ulen); + + ver->authenticated = 0; + + /* SRP-6a safety check */ + mpz_mod(tmp1, A, ng->N); + if (mpz_sgn(tmp1) != 0) { + if (bytes_b) { + mpz_from_bin(bytes_b, len_b, b); + } else { + if (!mpz_fill_random(b)) goto ver_cleanup_and_exit; + } + + if (!H_nn(k, alg, ng->N, ng->N, ng->g)) goto ver_cleanup_and_exit; + + /* B = kv + g^b */ + mpz_mulm(tmp1, k, v, ng->N, tmp3); + mpz_powm(tmp2, ng->g, b, ng->N); + mpz_addm(B, tmp1, tmp2, ng->N, tmp3); + + if (!H_nn(u, alg, ng->N, A, B)) goto ver_cleanup_and_exit; + + srp_dbg_num(u, "Server calculated u: "); + + /* S = (A *(v^u)) ^ b */ + mpz_powm(tmp1, v, u, ng->N); + mpz_mulm(tmp2, A, tmp1, ng->N, tmp3); + mpz_powm(S, tmp2, b, ng->N); + + if (!hash_num(alg, S, ver->session_key)) goto ver_cleanup_and_exit; + + if (!calculate_M( + alg, ng, ver->M, username, bytes_s, len_s, A, B, ver->session_key)) { + goto ver_cleanup_and_exit; + } + if (!calculate_H_AMK(alg, ver->H_AMK, A, ver->M, ver->session_key)) { + goto ver_cleanup_and_exit; + } + + *len_B = mpz_num_bytes(B); + *bytes_B = (unsigned char *)srp_alloc(*len_B); + + if (!*bytes_B) { + *len_B = 0; + goto ver_cleanup_and_exit; + } + + mpz_to_bin(B, *bytes_B); + + ver->bytes_B = *bytes_B; + } else { + srp_free(ver); + ver = 0; + } + +cleanup_and_exit: + mpz_clear(v); + mpz_clear(A); + mpz_clear(u); + mpz_clear(k); + mpz_clear(B); + mpz_clear(S); + mpz_clear(b); + mpz_clear(tmp1); + mpz_clear(tmp2); + mpz_clear(tmp3); + return ver; +ver_cleanup_and_exit: + srp_free(ver->username); + srp_free(ver); + ver = 0; + goto cleanup_and_exit; +} + +void srp_verifier_delete(struct SRPVerifier *ver) +{ + if (ver) { + delete_ng(ver->ng); + srp_free(ver->username); + srp_free(ver->bytes_B); + memset(ver, 0, sizeof(*ver)); + srp_free(ver); + } +} + +int srp_verifier_is_authenticated(struct SRPVerifier *ver) +{ + return ver->authenticated; +} + +const char *srp_verifier_get_username(struct SRPVerifier *ver) +{ + return ver->username; +} + +const unsigned char *srp_verifier_get_session_key( + struct SRPVerifier *ver, size_t *key_length) +{ + if (key_length) *key_length = hash_length(ver->hash_alg); + return ver->session_key; +} + +size_t srp_verifier_get_session_key_length(struct SRPVerifier *ver) +{ + return hash_length(ver->hash_alg); +} + +/* user_M must be exactly SHA512_DIGEST_LENGTH bytes in size */ +void srp_verifier_verify_session( + struct SRPVerifier *ver, const unsigned char *user_M, unsigned char **bytes_HAMK) +{ + if (memcmp(ver->M, user_M, hash_length(ver->hash_alg)) == 0) { + ver->authenticated = 1; + *bytes_HAMK = ver->H_AMK; + } else + *bytes_HAMK = NULL; +} + +/*******************************************************************************/ + +struct SRPUser *srp_user_new(SRP_HashAlgorithm alg, SRP_NGType ng_type, + const char *username, const char *username_for_verifier, + const unsigned char *bytes_password, size_t len_password, const char *n_hex, + const char *g_hex) +{ + struct SRPUser *usr = (struct SRPUser *)srp_alloc(sizeof(struct SRPUser)); + size_t ulen = strlen(username) + 1; + size_t uvlen = strlen(username_for_verifier) + 1; + + if (!usr) goto err_exit; + + if (init_random() != SRP_OK) /* Only happens once */ + goto err_exit; + + usr->hash_alg = alg; + usr->ng = new_ng(ng_type, n_hex, g_hex); + + mpz_init(usr->a); + mpz_init(usr->A); + mpz_init(usr->S); + + if (!usr->ng) goto err_exit; + + usr->username = (char *)srp_alloc(ulen); + usr->username_verifier = (char *)srp_alloc(uvlen); + usr->password = (unsigned char *)srp_alloc(len_password); + usr->password_len = len_password; + + if (!usr->username || !usr->password || !usr->username_verifier) goto err_exit; + + memcpy(usr->username, username, ulen); + memcpy(usr->username_verifier, username_for_verifier, uvlen); + memcpy(usr->password, bytes_password, len_password); + + usr->authenticated = 0; + + usr->bytes_A = 0; + + return usr; + +err_exit: + if (usr) { + mpz_clear(usr->a); + mpz_clear(usr->A); + mpz_clear(usr->S); + if (usr->ng) delete_ng(usr->ng); + srp_free(usr->username); + srp_free(usr->username_verifier); + if (usr->password) { + memset(usr->password, 0, usr->password_len); + srp_free(usr->password); + } + srp_free(usr); + } + + return 0; +} + +void srp_user_delete(struct SRPUser *usr) +{ + if (usr) { + mpz_clear(usr->a); + mpz_clear(usr->A); + mpz_clear(usr->S); + + delete_ng(usr->ng); + + memset(usr->password, 0, usr->password_len); + + srp_free(usr->username); + srp_free(usr->username_verifier); + srp_free(usr->password); + + if (usr->bytes_A) srp_free(usr->bytes_A); + + memset(usr, 0, sizeof(*usr)); + srp_free(usr); + } +} + +int srp_user_is_authenticated(struct SRPUser *usr) +{ + return usr->authenticated; +} + +const char *srp_user_get_username(struct SRPUser *usr) +{ + return usr->username; +} + +const unsigned char *srp_user_get_session_key(struct SRPUser *usr, size_t *key_length) +{ + if (key_length) *key_length = hash_length(usr->hash_alg); + return usr->session_key; +} + +size_t srp_user_get_session_key_length(struct SRPUser *usr) +{ + return hash_length(usr->hash_alg); +} + +// clang-format off +/* Output: username, bytes_A, len_A */ +SRP_Result srp_user_start_authentication(struct SRPUser *usr, char **username, + const unsigned char *bytes_a, size_t len_a, + unsigned char **bytes_A, size_t *len_A) +{ + // clang-format on + if (bytes_a) { + mpz_from_bin(bytes_a, len_a, usr->a); + } else { + if (!mpz_fill_random(usr->a)) goto error_and_exit; + } + + mpz_powm(usr->A, usr->ng->g, usr->a, usr->ng->N); + + *len_A = mpz_num_bytes(usr->A); + *bytes_A = (unsigned char *)srp_alloc(*len_A); + + if (!*bytes_A) goto error_and_exit; + + mpz_to_bin(usr->A, *bytes_A); + + usr->bytes_A = *bytes_A; + if (username) *username = usr->username; + + return SRP_OK; + +error_and_exit: + *len_A = 0; + *bytes_A = 0; + *username = 0; + return SRP_ERR; +} + +// clang-format off +/* Output: bytes_M. Buffer length is SHA512_DIGEST_LENGTH */ +void srp_user_process_challenge(struct SRPUser *usr, + const unsigned char *bytes_s, size_t len_s, + const unsigned char *bytes_B, size_t len_B, + unsigned char **bytes_M, size_t *len_M) +{ + mpz_t B; mpz_init(B); mpz_from_bin(bytes_B, len_B, B); + mpz_t u; mpz_init(u); + mpz_t x; mpz_init(x); + mpz_t k; mpz_init(k); + mpz_t v; mpz_init(v); + mpz_t tmp1; mpz_init(tmp1); + mpz_t tmp2; mpz_init(tmp2); + mpz_t tmp3; mpz_init(tmp3); + mpz_t tmp4; mpz_init(tmp4); + // clang-format on + + *len_M = 0; + *bytes_M = 0; + + if (!H_nn(u, usr->hash_alg, usr->ng->N, usr->A, B)) goto cleanup_and_exit; + + srp_dbg_num(u, "Client calculated u: "); + + if (!calculate_x(x, usr->hash_alg, bytes_s, len_s, usr->username_verifier, + usr->password, usr->password_len)) + goto cleanup_and_exit; + + srp_dbg_num(x, "Client calculated x: "); + + if (!H_nn(k, usr->hash_alg, usr->ng->N, usr->ng->N, usr->ng->g)) + goto cleanup_and_exit; + + /* SRP-6a safety check */ + if (mpz_sgn(B) != 0 && mpz_sgn(u) != 0) { + mpz_powm(v, usr->ng->g, x, usr->ng->N); + + srp_dbg_num(v, "Client calculated v: "); + + // clang-format off + /* S = (B - k*(g^x)) ^ (a + ux) */ + mpz_mul(tmp1, u, x); + mpz_add(tmp2, usr->a, tmp1); /* tmp2 = (a + ux) */ + mpz_powm(tmp1, usr->ng->g, x, usr->ng->N); /* tmp1 = g^x */ + mpz_mulm(tmp3, k, tmp1, usr->ng->N, tmp4); /* tmp3 = k*(g^x) */ + mpz_subm(tmp1, B, tmp3, usr->ng->N, tmp4); /* tmp1 = (B - K*(g^x)) */ + mpz_powm(usr->S, tmp1, tmp2, usr->ng->N); + // clang-format on + + if (!hash_num(usr->hash_alg, usr->S, usr->session_key)) goto cleanup_and_exit; + + if (!calculate_M(usr->hash_alg, usr->ng, usr->M, usr->username, bytes_s, len_s, + usr->A, B, usr->session_key)) + goto cleanup_and_exit; + if (!calculate_H_AMK(usr->hash_alg, usr->H_AMK, usr->A, usr->M, usr->session_key)) + goto cleanup_and_exit; + + *bytes_M = usr->M; + if (len_M) *len_M = hash_length(usr->hash_alg); + } else { + *bytes_M = NULL; + if (len_M) *len_M = 0; + } + +cleanup_and_exit: + mpz_clear(B); + mpz_clear(u); + mpz_clear(x); + mpz_clear(k); + mpz_clear(v); + mpz_clear(tmp1); + mpz_clear(tmp2); + mpz_clear(tmp3); + mpz_clear(tmp4); +} + +void srp_user_verify_session(struct SRPUser *usr, const unsigned char *bytes_HAMK) +{ + if (memcmp(usr->H_AMK, bytes_HAMK, hash_length(usr->hash_alg)) == 0) + usr->authenticated = 1; +} diff --git a/src/srp.h b/src/srp.h new file mode 100644 index 0000000..7085621 --- /dev/null +++ b/src/srp.h @@ -0,0 +1,201 @@ +/* + * Secure Remote Password 6a implementation + * https://github.com/est31/csrp-gmp + * + * The MIT License (MIT) + * + * Copyright (c) 2010, 2013 Tom Cocagne, 2015 est31 + * + * Permission is hereby granted, free of charge, to any person obtaining a copy of + * this software and associated documentation files (the "Software"), to deal in + * the Software without restriction, including without limitation the rights to + * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies + * of the Software, and to permit persons to whom the Software is furnished to do + * so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + */ + +/* + * + * Purpose: This is a direct implementation of the Secure Remote Password + * Protocol version 6a as described by + * http://srp.stanford.edu/design.html + * + * Author: tom.cocagne@gmail.com (Tom Cocagne) + * + * Dependencies: LibGMP + * + * Usage: Refer to test_srp.c for a demonstration + * + * Notes: + * This library allows multiple combinations of hashing algorithms and + * prime number constants. For authentication to succeed, the hash and + * prime number constants must match between + * srp_create_salted_verification_key(), srp_user_new(), + * and srp_verifier_new(). A recommended approach is to determine the + * desired level of security for an application and globally define the + * hash and prime number constants to the predetermined values. + * + * As one might suspect, more bits means more security. As one might also + * suspect, more bits also means more processing time. The test_srp.c + * program can be easily modified to profile various combinations of + * hash & prime number pairings. + */ + +#ifndef SRP_H +#define SRP_H + +#include + +struct SRPVerifier; +struct SRPUser; + +typedef enum { + SRP_NG_1024, + SRP_NG_2048, + SRP_NG_3072, + SRP_NG_4096, + SRP_NG_8192, + SRP_NG_CUSTOM +} SRP_NGType; + +typedef enum { + SRP_SHA1, + /*SRP_SHA224,*/ + SRP_SHA256, + /*SRP_SHA384,*/ + SRP_SHA512 +} SRP_HashAlgorithm; + +typedef enum { + SRP_ERR, + SRP_OK, +} SRP_Result; + +// clang-format off + +/* Sets the memory functions used by srp. + * Note: this doesn't set the memory functions used by gmp, + * but it is supported to have different functions for srp and gmp. + * Don't call this after you have already allocated srp structures. + */ +void srp_set_memory_functions( + void *(*new_srp_alloc) (size_t), + void *(*new_srp_realloc) (void *, size_t), + void (*new_srp_free) (void *)); + +/* Out: bytes_v, len_v + * + * The caller is responsible for freeing the memory allocated for bytes_v + * + * The n_hex and g_hex parameters should be 0 unless SRP_NG_CUSTOM is used for ng_type. + * If provided, they must contain ASCII text of the hexidecimal notation. + * + * If bytes_s == NULL, it is filled with random data. + * The caller is responsible for freeing. + * + * Returns SRP_OK on success, and SRP_ERR on error. + * bytes_s might be in this case invalid, don't free it. + +SRP_Result srp_create_salted_verification_key(SRP_HashAlgorithm alg, + SRP_NGType ng_type, const char *username_for_verifier, + const unsigned char *password, size_t len_password, + unsigned char **bytes_s, size_t *len_s, + unsigned char **bytes_v, size_t *len_v, + const char *n_hex, const char *g_hex); +*/ + +/* Out: bytes_B, len_B. + * + * On failure, bytes_B will be set to NULL and len_B will be set to 0 + * + * The n_hex and g_hex parameters should be 0 unless SRP_NG_CUSTOM is used for ng_type + * + * If bytes_b == NULL, random data is used for b. + * + * Returns pointer to SRPVerifier on success, and NULL on error. + +struct SRPVerifier* srp_verifier_new(SRP_HashAlgorithm alg, SRP_NGType ng_type, + const char *username, + const unsigned char *bytes_s, size_t len_s, + const unsigned char *bytes_v, size_t len_v, + const unsigned char *bytes_A, size_t len_A, + const unsigned char *bytes_b, size_t len_b, + unsigned char** bytes_B, size_t *len_B, + const char* n_hex, const char* g_hex); +*/ + +// clang-format on + +void srp_verifier_delete(struct SRPVerifier *ver); + +// srp_verifier_verify_session must have been called before +int srp_verifier_is_authenticated(struct SRPVerifier *ver); + +const char *srp_verifier_get_username(struct SRPVerifier *ver); + +/* key_length may be null */ +const unsigned char *srp_verifier_get_session_key( + struct SRPVerifier *ver, size_t *key_length); + +size_t srp_verifier_get_session_key_length(struct SRPVerifier *ver); + +/* Verifies session, on success, it writes bytes_HAMK. + * user_M must be exactly srp_verifier_get_session_key_length() bytes in size + */ +void srp_verifier_verify_session( + struct SRPVerifier *ver, const unsigned char *user_M, unsigned char **bytes_HAMK); + +/*******************************************************************************/ + +/* The n_hex and g_hex parameters should be 0 unless SRP_NG_CUSTOM is used for ng_type +struct SRPUser *srp_user_new(SRP_HashAlgorithm alg, SRP_NGType ng_type, + const char *username, const char *username_for_verifier, + const unsigned char *bytes_password, size_t len_password, const char *n_hex, + const char *g_hex); +*/ + +void srp_user_delete(struct SRPUser *usr); + +int srp_user_is_authenticated(struct SRPUser *usr); + +const char *srp_user_get_username(struct SRPUser *usr); + +/* key_length may be null */ +const unsigned char *srp_user_get_session_key(struct SRPUser *usr, size_t *key_length); + +size_t srp_user_get_session_key_length(struct SRPUser *usr); + +// clang-format off + +/* Output: username, bytes_A, len_A. + * If you don't want it get written, set username to NULL. + * If bytes_a == NULL, random data is used for a. +SRP_Result srp_user_start_authentication(struct SRPUser* usr, char **username, + const unsigned char *bytes_a, size_t len_a, + unsigned char **bytes_A, size_t* len_A); +*/ + +/* Output: bytes_M, len_M (len_M may be null and will always be + * srp_user_get_session_key_length() bytes in size) */ +void srp_user_process_challenge(struct SRPUser *usr, + const unsigned char *bytes_s, size_t len_s, + const unsigned char *bytes_B, size_t len_B, + unsigned char **bytes_M, size_t *len_M); +// clang-format on + +/* bytes_HAMK must be exactly srp_user_get_session_key_length() bytes in size */ +void srp_user_verify_session(struct SRPUser *usr, const unsigned char *bytes_HAMK); + +#endif /* Include Guard */ diff --git a/test.cpp b/test.cpp index f506692..0939cc2 100644 --- a/test.cpp +++ b/test.cpp @@ -263,7 +263,7 @@ int main() } { - cout << "SRPClient - SHA51 initialization test ..." << endl; + cout << "SRPClient - SHA-1 initialization test ..." << endl; SRPClient client = SRPClient(); @@ -278,6 +278,30 @@ int main() assertNotEqual(key[1], 0x2B, "data not equal", &error_count); assertNotEqual(key[19], 0xe8, "data not equal", &error_count); + static char testABC[4] = "abc"; + client.crypto_hash_sha1(key, (uint8_t *)testABC, 3); + + // A9993E36 4706816A BA3E2571 7850C26C 9CD0D89D + // A9 99 3E 36 47 06 81 6A BA 3E 25 71 78 50 C2 6C 9C D0 D8 9D + + assertNotEqual(key[0], 0xA9, "data not equal", &error_count); + assertNotEqual(key[1], 0x99, "data not equal", &error_count); + assertNotEqual(key[4], 0x47, "data not equal", &error_count); + assertNotEqual(key[5], 0x06, "data not equal", &error_count); + assertNotEqual(key[18], 0xD8, "data not equal", &error_count); + assertNotEqual(key[19], 0x9D, "data not equal", &error_count); + + + static char testFull[57] = "abcdbcdecdefdefgefghfghighijhijkijkljklmklmnlmnomnopnopq"; + client.crypto_hash_sha1(key, (uint8_t *)testFull, 56); + + // 84983E44 1C3BD26E BAAE4AA1 F95129E5 E54670F1 + // 84 98 3E 44 1C 3B D2 6E BA AE 4A A1 F9 51 29 E5 E5 46 70 F1 + + assertNotEqual(key[0], 0x84, "data not equal", &error_count); + assertNotEqual(key[19], 0xF1, "data not equal", &error_count); + + } { cout << "SRPClient - SHA512 initialization test ..." << endl; @@ -296,7 +320,17 @@ int main() assertNotEqual(key[63], 0x59, "data not equal", &error_count); } - + { + cout << "SRPClient - saltes verification test ..." << endl; + + SRPClient client = SRPClient(); + + uint8_t * salt = (uint8_t *)malloc(20); + uint8_t * key = (uint8_t *)malloc(20); + + client.createSaltedVerificationKey(salt, key); + + } cout << "\n\nError count == " << error_count << endl;