@@ -115,6 +115,36 @@ macro_rules! unsigned_fn {
115115 if n <= <$HalfBitsT>:: MAX as $UnsignedT {
116116 $HalfBitsT( n as $HalfBitsT) as $UnsignedT
117117 } 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)
118148 const EVEN_MAKING_BITMASK : u32 = !1 ;
119149 let normalization_shift = n. leading_zeros( ) & EVEN_MAKING_BITMASK ;
120150 n <<= normalization_shift;
0 commit comments