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
This commit is contained in:
SChernykh 2018-08-03 11:41:41 +02:00 committed by wowario
parent 1f8b0ddc83
commit 40573dd79c
No known key found for this signature in database
GPG Key ID: 24DCBE762DE9C111
5 changed files with 576 additions and 55 deletions

View File

@ -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);

View File

@ -0,0 +1,163 @@
#ifndef VARIANT2_INT_SQRT_H
#define VARIANT2_INT_SQRT_H
#include <math.h>
#include <float.h>
#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

View File

@ -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")

View File

@ -33,9 +33,11 @@
#include <iomanip>
#include <ios>
#include <string>
#include <cfenv>
#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 <emmintrin.h>
#if defined(_MSC_VER) || defined(__MINGW32__)
#include <intrin.h>
#else
#include <wmmintrin.h>
#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;
}

View File

@ -0,0 +1,10 @@
4cf1ff9ca46eb433b36cd9f70e02b14cc06bfd18ca77fa9ccaafd1fd96c674b0 5468697320697320612074657374205468697320697320612074657374205468697320697320612074657374
7d292e43f4751714ec07dbcb0e4bbffe2a7afb6066420960684ff57d7474c871 4c6f72656d20697073756d20646f6c6f722073697420616d65742c20636f6e73656374657475722061646970697363696e67
335563425256edebf1d92dc342369c2f4770ebb4112ba975659bd8a0f210abd0 656c69742c2073656420646f20656975736d6f642074656d706f7220696e6369646964756e74207574206c61626f7265
47758e86d2f57210366cec36fff26f9464d89efd116fe6ef28b718b5da120801 657420646f6c6f7265206d61676e6120616c697175612e20557420656e696d206164206d696e696d2076656e69616d2c
48787b48d5c68f0c1dd825c32580af741cc0ee314f08133135c1e86d87a24a95 71756973206e6f737472756420657865726369746174696f6e20756c6c616d636f206c61626f726973206e697369
93bdf47495854f7cfaaca1af8c0f39ef4a3024c10eb0dea23726b0e06ef29e84 757420616c697175697020657820656120636f6d6d6f646f20636f6e7365717561742e20447569732061757465
a375a71d0541057ccc96719150dfe10b6e6f486b19cf4a0835e19605413a8417 697275726520646f6c6f7220696e20726570726568656e646572697420696e20766f6c7570746174652076656c6974
163478a76f8f1432533fbdd1284d65c89f37479e54f20841c6ce4eba56c73854 657373652063696c6c756d20646f6c6f726520657520667567696174206e756c6c612070617269617475722e
356b0470c6eea75cad7a108179e232905b23bdaf03c2824c6e619d503ee93677 4578636570746575722073696e74206f6363616563617420637570696461746174206e6f6e2070726f6964656e742c
a47e2b007dc25bb279e197a1b91f67ecebe2ddd8791cd32dd2cb76dd21ed943f 73756e7420696e2063756c706120717569206f666669636961206465736572756e74206d6f6c6c697420616e696d20696420657374206c61626f72756d2e