@@ -115,6 +115,36 @@ macro_rules! unsigned_fn {
115
115
if n <= <$HalfBitsT>:: MAX as $UnsignedT {
116
116
$HalfBitsT( n as $HalfBitsT) as $UnsignedT
117
117
} else {
118
+ // The normalization shift satisfies the Karatsuba square root
119
+ // algorithm precondition "a₃ ≥ b/4" where a₃ is the most
120
+ // significant quarter of the bits and b is the number of
121
+ // values in `$UnsignedT`.
122
+ //
123
+ // In binary, b would thus be 1 followed by a count of
124
+ // `<$UnsignedT>::BITS` 0s. b/4 would then be 1 followed by
125
+ // `<$UnsignedT>::BITS - 2` 0s. a₃ must then be at least b/4,
126
+ // which means that the most significant bit or its neighbor
127
+ // must be a 1.
128
+ //
129
+ // The reason to shift by an even number of bits is because an
130
+ // even number of bits produces the square root shifted to the
131
+ // left by half of the normalization shift:
132
+ //
133
+ // sqrt(n << (2 * p))
134
+ // sqrt(2.pow(2 * p) * n)
135
+ // sqrt(2.pow(2 * p)) * sqrt(n)
136
+ // 2.pow(p) * sqrt(n)
137
+ // sqrt(n) << p
138
+ //
139
+ // Shifting by an odd number of bits leaves an ugly sqrt(2)
140
+ // multiplied in:
141
+ //
142
+ // sqrt(n << (2 * p + 1))
143
+ // sqrt(2.pow(2 * p + 1) * n)
144
+ // sqrt(2 * 2.pow(2 * p) * n)
145
+ // sqrt(2) * sqrt(2.pow(2 * p)) * sqrt(n)
146
+ // sqrt(2) * 2.pow(p) * sqrt(n)
147
+ // sqrt(2) * (sqrt(n) << p)
118
148
const EVEN_MAKING_BITMASK : u32 = !1 ;
119
149
let normalization_shift = n. leading_zeros( ) & EVEN_MAKING_BITMASK ;
120
150
n <<= normalization_shift;
0 commit comments