From 40573dd79cb51db6d7110acccde8446e5c8088d0 Mon Sep 17 00:00:00 2001 From: SChernykh Date: Fri, 3 Aug 2018 11:41:41 +0200 Subject: [PATCH] Cryptonight variant 2 Contains two modifications to improve ASIC resistance: shuffle and integer math. Shuffle makes use of the whole 64-byte cache line instead of 16 bytes only, making Cryptonight 4 times more demanding for memory bandwidth. Integer math adds 64:32 bit integer division followed by 64 bit integer square root, adding large and unavoidable computational latency to the main loop. More details and performance numbers: https://github.com/SChernykh/xmr-stak-cpu/blob/master/README.md --- src/crypto/slow-hash.c | 250 ++++++++++++++++++++++++++------- src/crypto/variant2_int_sqrt.h | 163 +++++++++++++++++++++ tests/hash/CMakeLists.txt | 6 +- tests/hash/main.cpp | 202 +++++++++++++++++++++++++- tests/hash/tests-slow-2.txt | 10 ++ 5 files changed, 576 insertions(+), 55 deletions(-) create mode 100644 src/crypto/variant2_int_sqrt.h create mode 100644 tests/hash/tests-slow-2.txt diff --git a/src/crypto/slow-hash.c b/src/crypto/slow-hash.c index d3c4da428..6ae0b11c9 100644 --- a/src/crypto/slow-hash.c +++ b/src/crypto/slow-hash.c @@ -38,6 +38,7 @@ #include "common/int-util.h" #include "hash-ops.h" #include "oaes_lib.h" +#include "variant2_int_sqrt.h" #define MEMORY (1 << 21) // 2MB scratchpad #define ITER (1 << 20) @@ -50,7 +51,7 @@ extern int aesb_single_round(const uint8_t *in, uint8_t*out, const uint8_t *expa extern int aesb_pseudo_round(const uint8_t *in, uint8_t *out, const uint8_t *expandedKey); #define VARIANT1_1(p) \ - do if (variant > 0) \ + do if (variant == 1) \ { \ const uint8_t tmp = ((const uint8_t*)(p))[11]; \ static const uint32_t table = 0x75310; \ @@ -59,7 +60,7 @@ extern int aesb_pseudo_round(const uint8_t *in, uint8_t *out, const uint8_t *exp } while(0) #define VARIANT1_2(p) \ - do if (variant > 0) \ + do if (variant == 1) \ { \ xor64(p, tweak1_2); \ } while(0) @@ -67,7 +68,7 @@ extern int aesb_pseudo_round(const uint8_t *in, uint8_t *out, const uint8_t *exp #define VARIANT1_CHECK() \ do if (length < 43) \ { \ - fprintf(stderr, "Cryptonight variants need at least 43 bytes of data"); \ + fprintf(stderr, "Cryptonight variant 1 needs at least 43 bytes of data"); \ _exit(1); \ } while(0) @@ -75,7 +76,7 @@ extern int aesb_pseudo_round(const uint8_t *in, uint8_t *out, const uint8_t *exp #define VARIANT1_PORTABLE_INIT() \ uint8_t tweak1_2[8]; \ - do if (variant > 0) \ + do if (variant == 1) \ { \ VARIANT1_CHECK(); \ memcpy(&tweak1_2, &state.hs.b[192], sizeof(tweak1_2)); \ @@ -83,11 +84,119 @@ extern int aesb_pseudo_round(const uint8_t *in, uint8_t *out, const uint8_t *exp } while(0) #define VARIANT1_INIT64() \ - if (variant > 0) \ + if (variant == 1) \ { \ VARIANT1_CHECK(); \ } \ - const uint64_t tweak1_2 = variant > 0 ? (state.hs.w[24] ^ (*((const uint64_t*)NONCE_POINTER))) : 0 + const uint64_t tweak1_2 = (variant == 1) ? (state.hs.w[24] ^ (*((const uint64_t*)NONCE_POINTER))) : 0 + +#define VARIANT2_INIT64() \ + uint64_t division_result = 0; \ + uint64_t sqrt_result = 0; \ + do if (variant >= 2) \ + { \ + U64(b)[2] = state.hs.w[8] ^ state.hs.w[10]; \ + U64(b)[3] = state.hs.w[9] ^ state.hs.w[11]; \ + division_result = state.hs.w[12]; \ + sqrt_result = state.hs.w[13]; \ + } while (0) + +#define VARIANT2_PORTABLE_INIT() \ + uint64_t division_result = 0; \ + uint64_t sqrt_result = 0; \ + do if (variant >= 2) \ + { \ + memcpy(b + AES_BLOCK_SIZE, state.hs.b + 64, AES_BLOCK_SIZE); \ + xor64(b + AES_BLOCK_SIZE, state.hs.b + 80); \ + xor64(b + AES_BLOCK_SIZE + 8, state.hs.b + 88); \ + division_result = state.hs.w[12]; \ + sqrt_result = state.hs.w[13]; \ + } while (0) + +#define VARIANT2_SHUFFLE_ADD_SSE2(base_ptr, offset) \ + do if (variant >= 2) \ + { \ + const __m128i chunk1 = _mm_load_si128((__m128i *)((base_ptr) + ((offset) ^ 0x10))); \ + const __m128i chunk2 = _mm_load_si128((__m128i *)((base_ptr) + ((offset) ^ 0x20))); \ + const __m128i chunk3 = _mm_load_si128((__m128i *)((base_ptr) + ((offset) ^ 0x30))); \ + _mm_store_si128((__m128i *)((base_ptr) + ((offset) ^ 0x10)), _mm_add_epi64(chunk3, _b1)); \ + _mm_store_si128((__m128i *)((base_ptr) + ((offset) ^ 0x20)), _mm_add_epi64(chunk1, _b)); \ + _mm_store_si128((__m128i *)((base_ptr) + ((offset) ^ 0x30)), _mm_add_epi64(chunk2, _a)); \ + } while (0) + +#define VARIANT2_SHUFFLE_ADD_NEON(base_ptr, offset) \ + do if (variant >= 2) \ + { \ + const uint64x2_t chunk1 = vld1q_u64(U64((base_ptr) + ((offset) ^ 0x10))); \ + const uint64x2_t chunk2 = vld1q_u64(U64((base_ptr) + ((offset) ^ 0x20))); \ + const uint64x2_t chunk3 = vld1q_u64(U64((base_ptr) + ((offset) ^ 0x30))); \ + vst1q_u64(U64((base_ptr) + ((offset) ^ 0x10)), vaddq_u64(chunk3, vreinterpretq_u64_u8(_b1))); \ + vst1q_u64(U64((base_ptr) + ((offset) ^ 0x20)), vaddq_u64(chunk1, vreinterpretq_u64_u8(_b))); \ + vst1q_u64(U64((base_ptr) + ((offset) ^ 0x30)), vaddq_u64(chunk2, vreinterpretq_u64_u8(_a))); \ + } while (0) + +#define VARIANT2_PORTABLE_SHUFFLE_ADD(base_ptr, offset) \ + do if (variant >= 2) \ + { \ + uint64_t* chunk1 = U64((base_ptr) + ((offset) ^ 0x10)); \ + uint64_t* chunk2 = U64((base_ptr) + ((offset) ^ 0x20)); \ + uint64_t* chunk3 = U64((base_ptr) + ((offset) ^ 0x30)); \ + \ + const uint64_t chunk1_old[2] = { chunk1[0], chunk1[1] }; \ + \ + uint64_t b1[2]; \ + memcpy(b1, b + 16, 16); \ + chunk1[0] = chunk3[0] + b1[0]; \ + chunk1[1] = chunk3[1] + b1[1]; \ + \ + uint64_t a0[2]; \ + memcpy(a0, a, 16); \ + chunk3[0] = chunk2[0] + a0[0]; \ + chunk3[1] = chunk2[1] + a0[1]; \ + \ + uint64_t b0[2]; \ + memcpy(b0, b, 16); \ + chunk2[0] = chunk1_old[0] + b0[0]; \ + chunk2[1] = chunk1_old[1] + b0[1]; \ + } while (0) + +#define VARIANT2_INTEGER_MATH_DIVISION_STEP(b, ptr) \ + ((uint64_t*)(b))[0] ^= division_result ^ (sqrt_result << 32); \ + { \ + const uint64_t dividend = ((uint64_t*)(ptr))[1]; \ + const uint32_t divisor = (((uint64_t*)(ptr))[0] + (uint32_t)(sqrt_result << 1)) | 0x80000001UL; \ + division_result = ((uint32_t)(dividend / divisor)) + \ + (((uint64_t)(dividend % divisor)) << 32); \ + } \ + const uint64_t sqrt_input = ((uint64_t*)(ptr))[0] + division_result + +#define VARIANT2_INTEGER_MATH_SSE2(b, ptr) \ + do if (variant >= 2) \ + { \ + VARIANT2_INTEGER_MATH_DIVISION_STEP(b, ptr); \ + VARIANT2_INTEGER_MATH_SQRT_STEP_SSE2(); \ + VARIANT2_INTEGER_MATH_SQRT_FIXUP(sqrt_result); \ + } while(0) + +#if defined DBL_MANT_DIG && (DBL_MANT_DIG >= 50) + // double precision floating point type has enough bits of precision on current platform + #define VARIANT2_PORTABLE_INTEGER_MATH(b, ptr) \ + do if (variant >= 2) \ + { \ + VARIANT2_INTEGER_MATH_DIVISION_STEP(b, ptr); \ + VARIANT2_INTEGER_MATH_SQRT_STEP_FP64(); \ + VARIANT2_INTEGER_MATH_SQRT_FIXUP(sqrt_result); \ + } while (0) +#else + // double precision floating point type is not good enough on current platform + // fall back to the reference code (integer only) + #define VARIANT2_PORTABLE_INTEGER_MATH(b, ptr) \ + do if (variant >= 2) \ + { \ + VARIANT2_INTEGER_MATH_DIVISION_STEP(b, ptr); \ + VARIANT2_INTEGER_MATH_SQRT_STEP_REF(); \ + } while (0) +#endif #if !defined NO_AES && (defined(__x86_64__) || (defined(_MSC_VER) && defined(_WIN64))) // Optimised code below, uses x86-specific intrinsics, SSE2, AES-NI @@ -164,19 +273,22 @@ extern int aesb_pseudo_round(const uint8_t *in, uint8_t *out, const uint8_t *exp * This code is based upon an optimized implementation by dga. */ #define post_aes() \ + VARIANT2_SHUFFLE_ADD_SSE2(hp_state, j); \ _mm_store_si128(R128(c), _c); \ - _b = _mm_xor_si128(_b, _c); \ - _mm_store_si128(R128(&hp_state[j]), _b); \ + _mm_store_si128(R128(&hp_state[j]), _mm_xor_si128(_b, _c)); \ VARIANT1_1(&hp_state[j]); \ j = state_index(c); \ p = U64(&hp_state[j]); \ b[0] = p[0]; b[1] = p[1]; \ + VARIANT2_INTEGER_MATH_SSE2(b, c); \ __mul(); \ + VARIANT2_SHUFFLE_ADD_SSE2(hp_state, j); \ a[0] += hi; a[1] += lo; \ p = U64(&hp_state[j]); \ p[0] = a[0]; p[1] = a[1]; \ a[0] ^= b[0]; a[1] ^= b[1]; \ VARIANT1_2(p + 1); \ + _b1 = _b; \ _b = _c; \ #if defined(_MSC_VER) @@ -570,10 +682,10 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int uint8_t text[INIT_SIZE_BYTE]; RDATA_ALIGN16 uint64_t a[2]; - RDATA_ALIGN16 uint64_t b[2]; + RDATA_ALIGN16 uint64_t b[4]; RDATA_ALIGN16 uint64_t c[2]; union cn_slow_hash_state state; - __m128i _a, _b, _c; + __m128i _a, _b, _b1, _c; uint64_t hi, lo; size_t i, j; @@ -599,6 +711,7 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int memcpy(text, state.init, INIT_SIZE_BYTE); VARIANT1_INIT64(); + VARIANT2_INIT64(); /* CryptoNight Step 2: Iteratively encrypt the results from Keccak to fill * the 2MB large random access buffer. @@ -637,6 +750,7 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int */ _b = _mm_load_si128(R128(b)); + _b1 = _mm_load_si128(R128(b) + 1); // Two independent versions, one with AES, one without, to ensure that // the useAes test is only performed once, not every iteration. if(useAes) @@ -761,19 +875,22 @@ union cn_slow_hash_state _a = vld1q_u8((const uint8_t *)a); \ #define post_aes() \ + VARIANT2_SHUFFLE_ADD_NEON(hp_state, j); \ vst1q_u8((uint8_t *)c, _c); \ - _b = veorq_u8(_b, _c); \ - vst1q_u8(&hp_state[j], _b); \ + vst1q_u8(&hp_state[j], veorq_u8(_b, _c)); \ VARIANT1_1(&hp_state[j]); \ j = state_index(c); \ p = U64(&hp_state[j]); \ b[0] = p[0]; b[1] = p[1]; \ + VARIANT2_PORTABLE_INTEGER_MATH(b, c); \ __mul(); \ + VARIANT2_SHUFFLE_ADD_NEON(hp_state, j); \ a[0] += hi; a[1] += lo; \ p = U64(&hp_state[j]); \ p[0] = a[0]; p[1] = a[1]; \ a[0] ^= b[0]; a[1] ^= b[1]; \ VARIANT1_2(p + 1); \ + _b1 = _b; \ _b = _c; \ @@ -937,10 +1054,10 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int uint8_t text[INIT_SIZE_BYTE]; RDATA_ALIGN16 uint64_t a[2]; - RDATA_ALIGN16 uint64_t b[2]; + RDATA_ALIGN16 uint64_t b[4]; RDATA_ALIGN16 uint64_t c[2]; union cn_slow_hash_state state; - uint8x16_t _a, _b, _c, zero = {0}; + uint8x16_t _a, _b, _b1, _c, zero = {0}; uint64_t hi, lo; size_t i, j; @@ -961,6 +1078,7 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int memcpy(text, state.init, INIT_SIZE_BYTE); VARIANT1_INIT64(); + VARIANT2_INIT64(); /* CryptoNight Step 2: Iteratively encrypt the results from Keccak to fill * the 2MB large random access buffer. @@ -984,7 +1102,7 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int */ _b = vld1q_u8((const uint8_t *)b); - + _b1 = vld1q_u8(((const uint8_t *)b) + AES_BLOCK_SIZE); for(i = 0; i < ITER / 2; i++) { @@ -1104,6 +1222,11 @@ __asm__ __volatile__( #endif /* !aarch64 */ #endif // NO_OPTIMIZED_MULTIPLY_ON_ARM +STATIC INLINE void copy_block(uint8_t* dst, const uint8_t* src) +{ + memcpy(dst, src, AES_BLOCK_SIZE); +} + STATIC INLINE void sum_half_blocks(uint8_t* a, const uint8_t* b) { uint64_t a0, a1, b0, b1; @@ -1138,7 +1261,9 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int { uint8_t text[INIT_SIZE_BYTE]; uint8_t a[AES_BLOCK_SIZE]; - uint8_t b[AES_BLOCK_SIZE]; + uint8_t b[AES_BLOCK_SIZE * 2]; + uint8_t c[AES_BLOCK_SIZE]; + uint8_t c1[AES_BLOCK_SIZE]; uint8_t d[AES_BLOCK_SIZE]; uint8_t aes_key[AES_KEY_SIZE]; RDATA_ALIGN16 uint8_t expandedKey[256]; @@ -1166,11 +1291,12 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int } memcpy(text, state.init, INIT_SIZE_BYTE); - VARIANT1_INIT64(); - aes_ctx = (oaes_ctx *) oaes_alloc(); oaes_key_import_data(aes_ctx, state.hs.b, AES_KEY_SIZE); + VARIANT1_INIT64(); + VARIANT2_INIT64(); + // use aligned data memcpy(expandedKey, aes_ctx->key->exp_data, aes_ctx->key->exp_data_len); for(i = 0; i < MEMORY / INIT_SIZE_BYTE; i++) @@ -1191,23 +1317,33 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int #define state_index(x) ((*(uint32_t *) x) & MASK) // Iteration 1 - p = &long_state[state_index(a)]; + j = state_index(a); + p = &long_state[j]; aesb_single_round(p, p, a); + copy_block(c1, p); - xor_blocks(b, p); - swap_blocks(b, p); - swap_blocks(a, b); + VARIANT2_PORTABLE_SHUFFLE_ADD(long_state, j); + xor_blocks(p, b); VARIANT1_1(p); // Iteration 2 - p = &long_state[state_index(a)]; + j = state_index(c1); + p = &long_state[j]; + copy_block(c, p); - mul(a, p, d); - sum_half_blocks(b, d); - swap_blocks(b, p); - xor_blocks(b, p); - swap_blocks(a, b); - VARIANT1_2(U64(p) + 1); + VARIANT2_PORTABLE_INTEGER_MATH(c, c1); + mul(c1, c, d); + VARIANT2_PORTABLE_SHUFFLE_ADD(long_state, j); + sum_half_blocks(a, d); + swap_blocks(a, c); + xor_blocks(a, c); + VARIANT1_2(U64(c) + 1); + copy_block(p, c); + + if (variant >= 2) { + copy_block(b + AES_BLOCK_SIZE, b); + } + copy_block(b, c1); } memcpy(text, state.init, INIT_SIZE_BYTE); @@ -1331,8 +1467,9 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int union cn_slow_hash_state state; uint8_t text[INIT_SIZE_BYTE]; uint8_t a[AES_BLOCK_SIZE]; - uint8_t b[AES_BLOCK_SIZE]; - uint8_t c[AES_BLOCK_SIZE]; + uint8_t b[AES_BLOCK_SIZE * 2]; + uint8_t c1[AES_BLOCK_SIZE]; + uint8_t c2[AES_BLOCK_SIZE]; uint8_t d[AES_BLOCK_SIZE]; size_t i, j; uint8_t aes_key[AES_KEY_SIZE]; @@ -1348,6 +1485,7 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int aes_ctx = (oaes_ctx *) oaes_alloc(); VARIANT1_PORTABLE_INIT(); + VARIANT2_PORTABLE_INIT(); oaes_key_import_data(aes_ctx, aes_key, AES_KEY_SIZE); for (i = 0; i < MEMORY / INIT_SIZE_BYTE; i++) { @@ -1357,9 +1495,9 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int memcpy(&long_state[i * INIT_SIZE_BYTE], text, INIT_SIZE_BYTE); } - for (i = 0; i < 16; i++) { - a[i] = state.k[ i] ^ state.k[32 + i]; - b[i] = state.k[16 + i] ^ state.k[48 + i]; + for (i = 0; i < AES_BLOCK_SIZE; i++) { + a[i] = state.k[ i] ^ state.k[AES_BLOCK_SIZE * 2 + i]; + b[i] = state.k[AES_BLOCK_SIZE + i] ^ state.k[AES_BLOCK_SIZE * 3 + i]; } for (i = 0; i < ITER / 2; i++) { @@ -1368,26 +1506,32 @@ void cn_slow_hash(const void *data, size_t length, char *hash, int variant, int * next address <-+ */ /* Iteration 1 */ - j = e2i(a, MEMORY / AES_BLOCK_SIZE); - copy_block(c, &long_state[j * AES_BLOCK_SIZE]); - aesb_single_round(c, c, a); - xor_blocks(b, c); - swap_blocks(b, c); - copy_block(&long_state[j * AES_BLOCK_SIZE], c); - assert(j == e2i(a, MEMORY / AES_BLOCK_SIZE)); - swap_blocks(a, b); - VARIANT1_1(&long_state[j * AES_BLOCK_SIZE]); + j = e2i(a, MEMORY / AES_BLOCK_SIZE) * AES_BLOCK_SIZE; + copy_block(c1, &long_state[j]); + aesb_single_round(c1, c1, a); + VARIANT2_PORTABLE_SHUFFLE_ADD(long_state, j); + copy_block(&long_state[j], c1); + xor_blocks(&long_state[j], b); + assert(j == e2i(a, MEMORY / AES_BLOCK_SIZE) * AES_BLOCK_SIZE); + VARIANT1_1(&long_state[j]); /* Iteration 2 */ - j = e2i(a, MEMORY / AES_BLOCK_SIZE); - copy_block(c, &long_state[j * AES_BLOCK_SIZE]); - mul(a, c, d); - sum_half_blocks(b, d); - swap_blocks(b, c); - xor_blocks(b, c); - VARIANT1_2(c + 8); - copy_block(&long_state[j * AES_BLOCK_SIZE], c); - assert(j == e2i(a, MEMORY / AES_BLOCK_SIZE)); - swap_blocks(a, b); + j = e2i(c1, MEMORY / AES_BLOCK_SIZE) * AES_BLOCK_SIZE; + copy_block(c2, &long_state[j]); + VARIANT2_PORTABLE_INTEGER_MATH(c2, c1); + mul(c1, c2, d); + VARIANT2_PORTABLE_SHUFFLE_ADD(long_state, j); + swap_blocks(a, c1); + sum_half_blocks(c1, d); + swap_blocks(c1, c2); + xor_blocks(c1, c2); + VARIANT1_2(c2 + 8); + copy_block(&long_state[j], c2); + assert(j == e2i(a, MEMORY / AES_BLOCK_SIZE) * AES_BLOCK_SIZE); + if (variant >= 2) { + copy_block(b + AES_BLOCK_SIZE, b); + } + copy_block(b, a); + copy_block(a, c1); } memcpy(text, state.init, INIT_SIZE_BYTE); diff --git a/src/crypto/variant2_int_sqrt.h b/src/crypto/variant2_int_sqrt.h new file mode 100644 index 000000000..b405bb798 --- /dev/null +++ b/src/crypto/variant2_int_sqrt.h @@ -0,0 +1,163 @@ +#ifndef VARIANT2_INT_SQRT_H +#define VARIANT2_INT_SQRT_H + +#include +#include + +#define VARIANT2_INTEGER_MATH_SQRT_STEP_SSE2() \ + do { \ + const __m128i exp_double_bias = _mm_set_epi64x(0, 1023ULL << 52); \ + __m128d x = _mm_castsi128_pd(_mm_add_epi64(_mm_cvtsi64_si128(sqrt_input >> 12), exp_double_bias)); \ + x = _mm_sqrt_sd(_mm_setzero_pd(), x); \ + sqrt_result = (uint64_t)(_mm_cvtsi128_si64(_mm_sub_epi64(_mm_castpd_si128(x), exp_double_bias))) >> 19; \ + } while(0) + +#define VARIANT2_INTEGER_MATH_SQRT_STEP_FP64() \ + do { \ + sqrt_result = sqrt(sqrt_input + 18446744073709551616.0) * 2.0 - 8589934592.0; \ + } while(0) + +#define VARIANT2_INTEGER_MATH_SQRT_STEP_REF() \ + sqrt_result = integer_square_root_v2(sqrt_input) + +// Reference implementation of the integer square root for Cryptonight variant 2 +// Computes integer part of "sqrt(2^64 + n) * 2 - 2^33" +// +// In other words, given 64-bit unsigned integer n: +// 1) Write it as x = 1.NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN000... in binary (1 <= x < 2, all 64 bits of n are used) +// 2) Calculate sqrt(x) = 1.0RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR... (1 <= sqrt(x) < sqrt(2), so it will always start with "1.0" in binary) +// 3) Take 32 bits that come after "1.0" and return them as a 32-bit unsigned integer, discard all remaining bits +// +// Some sample inputs and outputs: +// +// Input | Output | Exact value of "sqrt(2^64 + n) * 2 - 2^33" +// -----------------|------------|------------------------------------------- +// 0 | 0 | 0 +// 2^32 | 0 | 0.99999999994179233909330885695244... +// 2^32 + 1 | 1 | 1.0000000001746229827200734316305... +// 2^50 | 262140 | 262140.00012206565608606978175873... +// 2^55 + 20963331 | 8384515 | 8384515.9999999997673963974959744... +// 2^55 + 20963332 | 8384516 | 8384516 +// 2^62 + 26599786 | 1013904242 | 1013904242.9999999999479374853545... +// 2^62 + 26599787 | 1013904243 | 1013904243.0000000001561875439364... +// 2^64 - 1 | 3558067407 | 3558067407.9041987696409179931096... + +// The reference implementation as it is now uses only unsigned int64 arithmetic, so it can't have undefined behavior +// It was tested once for all edge cases and confirmed correct +static inline uint32_t integer_square_root_v2(uint64_t n) +{ + uint64_t r = 1ULL << 63; + + for (uint64_t bit = 1ULL << 60; bit; bit >>= 2) + { + const bool b = (n < r + bit); + const uint64_t n_next = n - (r + bit); + const uint64_t r_next = r + bit * 2; + n = b ? n : n_next; + r = b ? r : r_next; + r >>= 1; + } + + return r * 2 + ((n > r) ? 1 : 0); +} + +/* +VARIANT2_INTEGER_MATH_SQRT_FIXUP checks that "r" is an integer part of "sqrt(2^64 + sqrt_input) * 2 - 2^33" and adds or subtracts 1 if needed +It's hard to understand how it works, so here is a full calculation of formulas used in VARIANT2_INTEGER_MATH_SQRT_FIXUP + +The following inequalities must hold for r if it's an integer part of "sqrt(2^64 + sqrt_input) * 2 - 2^33": +1) r <= sqrt(2^64 + sqrt_input) * 2 - 2^33 +2) r + 1 > sqrt(2^64 + sqrt_input) * 2 - 2^33 + +We need to check them using only unsigned integer arithmetic to avoid rounding errors and undefined behavior + +First inequality: r <= sqrt(2^64 + sqrt_input) * 2 - 2^33 +----------------------------------------------------------------------------------- +r <= sqrt(2^64 + sqrt_input) * 2 - 2^33 +r + 2^33 <= sqrt(2^64 + sqrt_input) * 2 +r/2 + 2^32 <= sqrt(2^64 + sqrt_input) +(r/2 + 2^32)^2 <= 2^64 + sqrt_input + +Rewrite r as r = s * 2 + b (s = trunc(r/2), b is 0 or 1) + +((s*2+b)/2 + 2^32)^2 <= 2^64 + sqrt_input +(s*2+b)^2/4 + 2*2^32*(s*2+b)/2 + 2^64 <= 2^64 + sqrt_input +(s*2+b)^2/4 + 2*2^32*(s*2+b)/2 <= sqrt_input +(s*2+b)^2/4 + 2^32*r <= sqrt_input +(s^2*4+2*s*2*b+b^2)/4 + 2^32*r <= sqrt_input +s^2+s*b+b^2/4 + 2^32*r <= sqrt_input +s*(s+b) + b^2/4 + 2^32*r <= sqrt_input + +Let r2 = s*(s+b) + r*2^32 +r2 + b^2/4 <= sqrt_input + +If this inequality doesn't hold, then we must decrement r: IF "r2 + b^2/4 > sqrt_input" THEN r = r - 1 + +b can be 0 or 1 +If b is 0 then we need to compare "r2 > sqrt_input" +If b is 1 then b^2/4 = 0.25, so we need to compare "r2 + 0.25 > sqrt_input" +Since both r2 and sqrt_input are integers, we can safely replace it with "r2 + 1 > sqrt_input" +----------------------------------------------------------------------------------- +Both cases can be merged to a single expression "r2 + b > sqrt_input" +----------------------------------------------------------------------------------- +There will be no overflow when calculating "r2 + b", so it's safe to compare with sqrt_input: +r2 + b = s*(s+b) + r*2^32 + b +The largest value s, b and r can have is s = 1779033703, b = 1, r = 3558067407 when sqrt_input = 2^64 - 1 +r2 + b <= 1779033703*1779033704 + 3558067407*2^32 + 1 = 18446744068217447385 < 2^64 + +Second inequality: r + 1 > sqrt(2^64 + sqrt_input) * 2 - 2^33 +----------------------------------------------------------------------------------- +r + 1 > sqrt(2^64 + sqrt_input) * 2 - 2^33 +r + 1 + 2^33 > sqrt(2^64 + sqrt_input) * 2 +((r+1)/2 + 2^32)^2 > 2^64 + sqrt_input + +Rewrite r as r = s * 2 + b (s = trunc(r/2), b is 0 or 1) + +((s*2+b+1)/2 + 2^32)^2 > 2^64 + sqrt_input +(s*2+b+1)^2/4 + 2*(s*2+b+1)/2*2^32 + 2^64 > 2^64 + sqrt_input +(s*2+b+1)^2/4 + (s*2+b+1)*2^32 > sqrt_input +(s*2+b+1)^2/4 + (r+1)*2^32 > sqrt_input +(s*2+(b+1))^2/4 + r*2^32 + 2^32 > sqrt_input +(s^2*4+2*s*2*(b+1)+(b+1)^2)/4 + r*2^32 + 2^32 > sqrt_input +s^2+s*(b+1)+(b+1)^2/4 + r*2^32 + 2^32 > sqrt_input +s*(s+b) + s + (b+1)^2/4 + r*2^32 + 2^32 > sqrt_input + +Let r2 = s*(s+b) + r*2^32 + +r2 + s + (b+1)^2/4 + 2^32 > sqrt_input +r2 + 2^32 + (b+1)^2/4 > sqrt_input - s + +If this inequality doesn't hold, then we must decrement r: IF "r2 + 2^32 + (b+1)^2/4 <= sqrt_input - s" THEN r = r - 1 +b can be 0 or 1 +If b is 0 then we need to compare "r2 + 2^32 + 1/4 <= sqrt_input - s" which is equal to "r2 + 2^32 < sqrt_input - s" because all numbers here are integers +If b is 1 then (b+1)^2/4 = 1, so we need to compare "r2 + 2^32 + 1 <= sqrt_input - s" which is also equal to "r2 + 2^32 < sqrt_input - s" +----------------------------------------------------------------------------------- +Both cases can be merged to a single expression "r2 + 2^32 < sqrt_input - s" +----------------------------------------------------------------------------------- +There will be no overflow when calculating "r2 + 2^32": +r2 + 2^32 = s*(s+b) + r*2^32 + 2^32 = s*(s+b) + (r+1)*2^32 +The largest value s, b and r can have is s = 1779033703, b = 1, r = 3558067407 when sqrt_input = 2^64 - 1 +r2 + b <= 1779033703*1779033704 + 3558067408*2^32 = 18446744072512414680 < 2^64 + +There will be no integer overflow when calculating "sqrt_input - s", i.e. "sqrt_input >= s" at all times: +s = trunc(r/2) = trunc(sqrt(2^64 + sqrt_input) - 2^32) < sqrt(2^64 + sqrt_input) - 2^32 + 1 +sqrt_input > sqrt(2^64 + sqrt_input) - 2^32 + 1 +sqrt_input + 2^32 - 1 > sqrt(2^64 + sqrt_input) +(sqrt_input + 2^32 - 1)^2 > sqrt_input + 2^64 +sqrt_input^2 + 2*sqrt_input*(2^32 - 1) + (2^32-1)^2 > sqrt_input + 2^64 +sqrt_input^2 + sqrt_input*(2^33 - 2) + (2^32-1)^2 > sqrt_input + 2^64 +sqrt_input^2 + sqrt_input*(2^33 - 3) + (2^32-1)^2 > 2^64 +sqrt_input^2 + sqrt_input*(2^33 - 3) + 2^64-2^33+1 > 2^64 +sqrt_input^2 + sqrt_input*(2^33 - 3) - 2^33 + 1 > 0 +This inequality is true if sqrt_input > 1 and it's easy to check that s = 0 if sqrt_input is 0 or 1, so there will be no integer overflow +*/ + +#define VARIANT2_INTEGER_MATH_SQRT_FIXUP(r) \ + do { \ + const uint64_t s = r >> 1; \ + const uint64_t b = r & 1; \ + const uint64_t r2 = (uint64_t)(s) * (s + b) + (r << 32); \ + r += ((r2 + b > sqrt_input) ? -1 : 0) + ((r2 + (1ULL << 32) < sqrt_input - s) ? 1 : 0); \ + } while(0) + +#endif diff --git a/tests/hash/CMakeLists.txt b/tests/hash/CMakeLists.txt index 92abeca20..433cf94e9 100644 --- a/tests/hash/CMakeLists.txt +++ b/tests/hash/CMakeLists.txt @@ -43,8 +43,12 @@ set_property(TARGET hash-tests PROPERTY FOLDER "tests") -foreach (hash IN ITEMS fast slow slow-1 tree extra-blake extra-groestl extra-jh extra-skein) +foreach (hash IN ITEMS fast slow slow-1 slow-2 tree extra-blake extra-groestl extra-jh extra-skein) add_test( NAME "hash-${hash}" COMMAND hash-tests "${hash}" "${CMAKE_CURRENT_SOURCE_DIR}/tests-${hash}.txt") endforeach () + +add_test( + NAME "hash-variant2-int-sqrt" + COMMAND hash-tests "variant2_int_sqrt") diff --git a/tests/hash/main.cpp b/tests/hash/main.cpp index cc5b9ba66..7767d0d3b 100644 --- a/tests/hash/main.cpp +++ b/tests/hash/main.cpp @@ -33,9 +33,11 @@ #include #include #include +#include #include "warnings.h" #include "crypto/hash.h" +#include "crypto/variant2_int_sqrt.h" #include "../io.h" using namespace std; @@ -57,6 +59,9 @@ extern "C" { static void cn_slow_hash_1(const void *data, size_t length, char *hash) { return cn_slow_hash(data, length, hash, 1/*variant*/, 0/*prehashed*/); } + static void cn_slow_hash_2(const void *data, size_t length, char *hash) { + return cn_slow_hash(data, length, hash, 2/*variant*/, 0/*prehashed*/); + } } POP_WARNINGS @@ -67,7 +72,10 @@ struct hash_func { } hashes[] = {{"fast", cn_fast_hash}, {"slow", cn_slow_hash_0}, {"tree", hash_tree}, {"extra-blake", hash_extra_blake}, {"extra-groestl", hash_extra_groestl}, {"extra-jh", hash_extra_jh}, {"extra-skein", hash_extra_skein}, - {"slow-1", cn_slow_hash_1}}; + {"slow-1", cn_slow_hash_1}, {"slow-2", cn_slow_hash_2}}; + +int test_variant2_int_sqrt(); +int test_variant2_int_sqrt_ref(); int main(int argc, char *argv[]) { hash_f *f; @@ -78,6 +86,36 @@ int main(int argc, char *argv[]) { size_t test = 0; bool error = false; if (argc != 3) { + if ((argc == 2) && (strcmp(argv[1], "variant2_int_sqrt") == 0)) { + if (test_variant2_int_sqrt_ref() != 0) { + return 1; + } + const int round_modes[3] = { FE_DOWNWARD, FE_TONEAREST, FE_UPWARD }; + for (int i = 0; i < 3; ++i) { + std::fesetround(round_modes[i]); + const int result = test_variant2_int_sqrt(); + if (result != 0) { + cerr << "FPU round mode was set to "; + switch (round_modes[i]) { + case FE_DOWNWARD: + cerr << "FE_DOWNWARD"; + break; + case FE_TONEAREST: + cerr << "FE_TONEAREST"; + break; + case FE_UPWARD: + cerr << "FE_UPWARD"; + break; + default: + cerr << "unknown"; + break; + } + cerr << endl; + return result; + } + } + return 0; + } cerr << "Wrong number of arguments" << endl; return 1; } @@ -127,3 +165,165 @@ int main(int argc, char *argv[]) { } return error ? 1 : 0; } + +#if defined(__x86_64__) || (defined(_MSC_VER) && defined(_WIN64)) + +#include + +#if defined(_MSC_VER) || defined(__MINGW32__) + #include +#else + #include +#endif + +#endif + +static inline bool test_variant2_int_sqrt_sse(const uint64_t sqrt_input, const uint64_t correct_result) +{ +#if defined(__x86_64__) || (defined(_MSC_VER) && defined(_WIN64)) + uint64_t sqrt_result; + VARIANT2_INTEGER_MATH_SQRT_STEP_SSE2(); + VARIANT2_INTEGER_MATH_SQRT_FIXUP(sqrt_result); + if (sqrt_result != correct_result) { + cerr << "Integer sqrt (SSE2 version) returned incorrect result for N = " << sqrt_input << endl; + cerr << "Expected result: " << correct_result << endl; + cerr << "Returned result: " << sqrt_result << endl; + return false; + } +#endif + + return true; +} + +static inline bool test_variant2_int_sqrt_fp64(const uint64_t sqrt_input, const uint64_t correct_result) +{ +#if defined DBL_MANT_DIG && (DBL_MANT_DIG >= 50) + uint64_t sqrt_result; + VARIANT2_INTEGER_MATH_SQRT_STEP_FP64(); + VARIANT2_INTEGER_MATH_SQRT_FIXUP(sqrt_result); + if (sqrt_result != correct_result) { + cerr << "Integer sqrt (FP64 version) returned incorrect result for N = " << sqrt_input << endl; + cerr << "Expected result: " << correct_result << endl; + cerr << "Returned result: " << sqrt_result << endl; + return false; + } +#endif + + return true; +} + +static inline bool test_variant2_int_sqrt_ref(const uint64_t sqrt_input, const uint64_t correct_result) +{ + uint64_t sqrt_result; + VARIANT2_INTEGER_MATH_SQRT_STEP_REF(); + if (sqrt_result != correct_result) { + cerr << "Integer sqrt (reference version) returned incorrect result for N = " << sqrt_input << endl; + cerr << "Expected result: " << correct_result << endl; + cerr << "Returned result: " << sqrt_result << endl; + return false; + } + + return true; +} + +static inline bool test_variant2_int_sqrt(const uint64_t sqrt_input, const uint64_t correct_result) +{ + if (!test_variant2_int_sqrt_sse(sqrt_input, correct_result)) { + return false; + } + if (!test_variant2_int_sqrt_fp64(sqrt_input, correct_result)) { + return false; + } + + return true; +} + +int test_variant2_int_sqrt() +{ + if (!test_variant2_int_sqrt(0, 0)) { + return 1; + } + if (!test_variant2_int_sqrt(1ULL << 63, 1930543745UL)) { + return 1; + } + if (!test_variant2_int_sqrt(uint64_t(-1), 3558067407UL)) { + return 1; + } + + for (uint64_t i = 1; i <= 3558067407UL; ++i) { + // "i" is integer part of "sqrt(2^64 + n) * 2 - 2^33" + // n = (i/2 + 2^32)^2 - 2^64 + + const uint64_t i0 = i >> 1; + uint64_t n1; + if ((i & 1) == 0) { + // n = (i/2 + 2^32)^2 - 2^64 + // n = i^2/4 + 2*2^32*i/2 + 2^64 - 2^64 + // n = i^2/4 + 2^32*i + // i is even, so i^2 is divisible by 4: + // n = (i^2 >> 2) + (i << 32) + + // int_sqrt_v2(i^2/4 + 2^32*i - 1) must be equal to i - 1 + // int_sqrt_v2(i^2/4 + 2^32*i) must be equal to i + n1 = i0 * i0 + (i << 32) - 1; + } + else { + // n = (i/2 + 2^32)^2 - 2^64 + // n = i^2/4 + 2*2^32*i/2 + 2^64 - 2^64 + // n = i^2/4 + 2^32*i + // i is odd, so i = i0*2+1 (i0 = i >> 1) + // n = (i0*2+1)^2/4 + 2^32*i + // n = (i0^2*4+i0*4+1)/4 + 2^32*i + // n = i0^2+i0+1/4 + 2^32*i + // i0^2+i0 + 2^32*i < n < i0^2+i0+1 + 2^32*i + + // int_sqrt_v2(i0^2+i0 + 2^32*i) must be equal to i - 1 + // int_sqrt_v2(i0^2+i0+1 + 2^32*i) must be equal to i + n1 = i0 * i0 + i0 + (i << 32); + } + + if (!test_variant2_int_sqrt(n1, i - 1)) { + return 1; + } + if (!test_variant2_int_sqrt(n1 + 1, i)) { + return 1; + } + } + + return 0; +} + +int test_variant2_int_sqrt_ref() +{ + if (!test_variant2_int_sqrt_ref(0, 0)) { + return 1; + } + if (!test_variant2_int_sqrt_ref(1ULL << 63, 1930543745UL)) { + return 1; + } + if (!test_variant2_int_sqrt_ref(uint64_t(-1), 3558067407UL)) { + return 1; + } + + // Reference version is slow, so we test only every 83th edge case + // "i += 83" because 1 + 83 * 42868282 = 3558067407 + for (uint64_t i = 1; i <= 3558067407UL; i += 83) { + const uint64_t i0 = i >> 1; + uint64_t n1; + if ((i & 1) == 0) { + n1 = i0 * i0 + (i << 32) - 1; + } + else { + n1 = i0 * i0 + i0 + (i << 32); + } + + if (!test_variant2_int_sqrt_ref(n1, i - 1)) { + return 1; + } + if (!test_variant2_int_sqrt_ref(n1 + 1, i)) { + return 1; + } + } + + return 0; +} diff --git a/tests/hash/tests-slow-2.txt b/tests/hash/tests-slow-2.txt new file mode 100644 index 000000000..8f90d05c9 --- /dev/null +++ b/tests/hash/tests-slow-2.txt @@ -0,0 +1,10 @@ +4cf1ff9ca46eb433b36cd9f70e02b14cc06bfd18ca77fa9ccaafd1fd96c674b0 5468697320697320612074657374205468697320697320612074657374205468697320697320612074657374 +7d292e43f4751714ec07dbcb0e4bbffe2a7afb6066420960684ff57d7474c871 4c6f72656d20697073756d20646f6c6f722073697420616d65742c20636f6e73656374657475722061646970697363696e67 +335563425256edebf1d92dc342369c2f4770ebb4112ba975659bd8a0f210abd0 656c69742c2073656420646f20656975736d6f642074656d706f7220696e6369646964756e74207574206c61626f7265 +47758e86d2f57210366cec36fff26f9464d89efd116fe6ef28b718b5da120801 657420646f6c6f7265206d61676e6120616c697175612e20557420656e696d206164206d696e696d2076656e69616d2c +48787b48d5c68f0c1dd825c32580af741cc0ee314f08133135c1e86d87a24a95 71756973206e6f737472756420657865726369746174696f6e20756c6c616d636f206c61626f726973206e697369 +93bdf47495854f7cfaaca1af8c0f39ef4a3024c10eb0dea23726b0e06ef29e84 757420616c697175697020657820656120636f6d6d6f646f20636f6e7365717561742e20447569732061757465 +a375a71d0541057ccc96719150dfe10b6e6f486b19cf4a0835e19605413a8417 697275726520646f6c6f7220696e20726570726568656e646572697420696e20766f6c7570746174652076656c6974 +163478a76f8f1432533fbdd1284d65c89f37479e54f20841c6ce4eba56c73854 657373652063696c6c756d20646f6c6f726520657520667567696174206e756c6c612070617269617475722e +356b0470c6eea75cad7a108179e232905b23bdaf03c2824c6e619d503ee93677 4578636570746575722073696e74206f6363616563617420637570696461746174206e6f6e2070726f6964656e742c +a47e2b007dc25bb279e197a1b91f67ecebe2ddd8791cd32dd2cb76dd21ed943f 73756e7420696e2063756c706120717569206f666669636961206465736572756e74206d6f6c6c697420616e696d20696420657374206c61626f72756d2e