Electroneum
variant2_int_sqrt.h
Go to the documentation of this file.
1 #ifndef VARIANT2_INT_SQRT_H
2 #define VARIANT2_INT_SQRT_H
3 
4 #include <math.h>
5 #include <float.h>
6 
7 #define VARIANT2_INTEGER_MATH_SQRT_STEP_SSE2() \
8  do { \
9  const __m128i exp_double_bias = _mm_set_epi64x(0, 1023ULL << 52); \
10  __m128d x = _mm_castsi128_pd(_mm_add_epi64(_mm_cvtsi64_si128(sqrt_input >> 12), exp_double_bias)); \
11  x = _mm_sqrt_sd(_mm_setzero_pd(), x); \
12  sqrt_result = (uint64_t)(_mm_cvtsi128_si64(_mm_sub_epi64(_mm_castpd_si128(x), exp_double_bias))) >> 19; \
13  } while(0)
14 
15 #define VARIANT2_INTEGER_MATH_SQRT_STEP_FP64() \
16  do { \
17  sqrt_result = sqrt(sqrt_input + 18446744073709551616.0) * 2.0 - 8589934592.0; \
18  } while(0)
19 
20 #define VARIANT2_INTEGER_MATH_SQRT_STEP_REF() \
21  sqrt_result = integer_square_root_v2(sqrt_input)
22 
23 // Reference implementation of the integer square root for Cryptonight variant 2
24 // Computes integer part of "sqrt(2^64 + n) * 2 - 2^33"
25 //
26 // In other words, given 64-bit unsigned integer n:
27 // 1) Write it as x = 1.NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN000... in binary (1 <= x < 2, all 64 bits of n are used)
28 // 2) Calculate sqrt(x) = 1.0RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR... (1 <= sqrt(x) < sqrt(2), so it will always start with "1.0" in binary)
29 // 3) Take 32 bits that come after "1.0" and return them as a 32-bit unsigned integer, discard all remaining bits
30 //
31 // Some sample inputs and outputs:
32 //
33 // Input | Output | Exact value of "sqrt(2^64 + n) * 2 - 2^33"
34 // -----------------|------------|-------------------------------------------
35 // 0 | 0 | 0
36 // 2^32 | 0 | 0.99999999994179233909330885695244...
37 // 2^32 + 1 | 1 | 1.0000000001746229827200734316305...
38 // 2^50 | 262140 | 262140.00012206565608606978175873...
39 // 2^55 + 20963331 | 8384515 | 8384515.9999999997673963974959744...
40 // 2^55 + 20963332 | 8384516 | 8384516
41 // 2^62 + 26599786 | 1013904242 | 1013904242.9999999999479374853545...
42 // 2^62 + 26599787 | 1013904243 | 1013904243.0000000001561875439364...
43 // 2^64 - 1 | 3558067407 | 3558067407.9041987696409179931096...
44 
45 // The reference implementation as it is now uses only unsigned int64 arithmetic, so it can't have undefined behavior
46 // It was tested once for all edge cases and confirmed correct
47 static inline uint32_t integer_square_root_v2(uint64_t n)
48 {
49  uint64_t r = 1ULL << 63;
50 
51  for (uint64_t bit = 1ULL << 60; bit; bit >>= 2)
52  {
53  const bool b = (n < r + bit);
54  const uint64_t n_next = n - (r + bit);
55  const uint64_t r_next = r + bit * 2;
56  n = b ? n : n_next;
57  r = b ? r : r_next;
58  r >>= 1;
59  }
60 
61  return r * 2 + ((n > r) ? 1 : 0);
62 }
63 
64 /*
65 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
66 It's hard to understand how it works, so here is a full calculation of formulas used in VARIANT2_INTEGER_MATH_SQRT_FIXUP
67 
68 The following inequalities must hold for r if it's an integer part of "sqrt(2^64 + sqrt_input) * 2 - 2^33":
69 1) r <= sqrt(2^64 + sqrt_input) * 2 - 2^33
70 2) r + 1 > sqrt(2^64 + sqrt_input) * 2 - 2^33
71 
72 We need to check them using only unsigned integer arithmetic to avoid rounding errors and undefined behavior
73 
74 First inequality: r <= sqrt(2^64 + sqrt_input) * 2 - 2^33
75 -----------------------------------------------------------------------------------
76 r <= sqrt(2^64 + sqrt_input) * 2 - 2^33
77 r + 2^33 <= sqrt(2^64 + sqrt_input) * 2
78 r/2 + 2^32 <= sqrt(2^64 + sqrt_input)
79 (r/2 + 2^32)^2 <= 2^64 + sqrt_input
80 
81 Rewrite r as r = s * 2 + b (s = trunc(r/2), b is 0 or 1)
82 
83 ((s*2+b)/2 + 2^32)^2 <= 2^64 + sqrt_input
84 (s*2+b)^2/4 + 2*2^32*(s*2+b)/2 + 2^64 <= 2^64 + sqrt_input
85 (s*2+b)^2/4 + 2*2^32*(s*2+b)/2 <= sqrt_input
86 (s*2+b)^2/4 + 2^32*r <= sqrt_input
87 (s^2*4+2*s*2*b+b^2)/4 + 2^32*r <= sqrt_input
88 s^2+s*b+b^2/4 + 2^32*r <= sqrt_input
89 s*(s+b) + b^2/4 + 2^32*r <= sqrt_input
90 
91 Let r2 = s*(s+b) + r*2^32
92 r2 + b^2/4 <= sqrt_input
93 
94 If this inequality doesn't hold, then we must decrement r: IF "r2 + b^2/4 > sqrt_input" THEN r = r - 1
95 
96 b can be 0 or 1
97 If b is 0 then we need to compare "r2 > sqrt_input"
98 If b is 1 then b^2/4 = 0.25, so we need to compare "r2 + 0.25 > sqrt_input"
99 Since both r2 and sqrt_input are integers, we can safely replace it with "r2 + 1 > sqrt_input"
100 -----------------------------------------------------------------------------------
101 Both cases can be merged to a single expression "r2 + b > sqrt_input"
102 -----------------------------------------------------------------------------------
103 There will be no overflow when calculating "r2 + b", so it's safe to compare with sqrt_input:
104 r2 + b = s*(s+b) + r*2^32 + b
105 The largest value s, b and r can have is s = 1779033703, b = 1, r = 3558067407 when sqrt_input = 2^64 - 1
106 r2 + b <= 1779033703*1779033704 + 3558067407*2^32 + 1 = 18446744068217447385 < 2^64
107 
108 Second inequality: r + 1 > sqrt(2^64 + sqrt_input) * 2 - 2^33
109 -----------------------------------------------------------------------------------
110 r + 1 > sqrt(2^64 + sqrt_input) * 2 - 2^33
111 r + 1 + 2^33 > sqrt(2^64 + sqrt_input) * 2
112 ((r+1)/2 + 2^32)^2 > 2^64 + sqrt_input
113 
114 Rewrite r as r = s * 2 + b (s = trunc(r/2), b is 0 or 1)
115 
116 ((s*2+b+1)/2 + 2^32)^2 > 2^64 + sqrt_input
117 (s*2+b+1)^2/4 + 2*(s*2+b+1)/2*2^32 + 2^64 > 2^64 + sqrt_input
118 (s*2+b+1)^2/4 + (s*2+b+1)*2^32 > sqrt_input
119 (s*2+b+1)^2/4 + (r+1)*2^32 > sqrt_input
120 (s*2+(b+1))^2/4 + r*2^32 + 2^32 > sqrt_input
121 (s^2*4+2*s*2*(b+1)+(b+1)^2)/4 + r*2^32 + 2^32 > sqrt_input
122 s^2+s*(b+1)+(b+1)^2/4 + r*2^32 + 2^32 > sqrt_input
123 s*(s+b) + s + (b+1)^2/4 + r*2^32 + 2^32 > sqrt_input
124 
125 Let r2 = s*(s+b) + r*2^32
126 
127 r2 + s + (b+1)^2/4 + 2^32 > sqrt_input
128 r2 + 2^32 + (b+1)^2/4 > sqrt_input - s
129 
130 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
131 b can be 0 or 1
132 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
133 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"
134 -----------------------------------------------------------------------------------
135 Both cases can be merged to a single expression "r2 + 2^32 < sqrt_input - s"
136 -----------------------------------------------------------------------------------
137 There will be no overflow when calculating "r2 + 2^32":
138 r2 + 2^32 = s*(s+b) + r*2^32 + 2^32 = s*(s+b) + (r+1)*2^32
139 The largest value s, b and r can have is s = 1779033703, b = 1, r = 3558067407 when sqrt_input = 2^64 - 1
140 r2 + b <= 1779033703*1779033704 + 3558067408*2^32 = 18446744072512414680 < 2^64
141 
142 There will be no integer overflow when calculating "sqrt_input - s", i.e. "sqrt_input >= s" at all times:
143 s = trunc(r/2) = trunc(sqrt(2^64 + sqrt_input) - 2^32) < sqrt(2^64 + sqrt_input) - 2^32 + 1
144 sqrt_input > sqrt(2^64 + sqrt_input) - 2^32 + 1
145 sqrt_input + 2^32 - 1 > sqrt(2^64 + sqrt_input)
146 (sqrt_input + 2^32 - 1)^2 > sqrt_input + 2^64
147 sqrt_input^2 + 2*sqrt_input*(2^32 - 1) + (2^32-1)^2 > sqrt_input + 2^64
148 sqrt_input^2 + sqrt_input*(2^33 - 2) + (2^32-1)^2 > sqrt_input + 2^64
149 sqrt_input^2 + sqrt_input*(2^33 - 3) + (2^32-1)^2 > 2^64
150 sqrt_input^2 + sqrt_input*(2^33 - 3) + 2^64-2^33+1 > 2^64
151 sqrt_input^2 + sqrt_input*(2^33 - 3) - 2^33 + 1 > 0
152 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
153 */
154 
155 #define VARIANT2_INTEGER_MATH_SQRT_FIXUP(r) \
156  do { \
157  const uint64_t s = r >> 1; \
158  const uint64_t b = r & 1; \
159  const uint64_t r2 = (uint64_t)(s) * (s + b) + (r << 32); \
160  r += ((r2 + b > sqrt_input) ? -1 : 0) + ((r2 + (1ULL << 32) < sqrt_input - s) ? 1 : 0); \
161  } while(0)
162 
163 #endif
unsigned int uint32_t
Definition: stdint.h:126
unsigned __int64 uint64_t
Definition: stdint.h:136