Calcul de la racine carrée entière

Il était nécessaire de vérifier si l'entier est un carré et, dans l'affirmative, de calculer la racine. Et je veux le faire en arithmétique entière. Il est clair que la méthode Newton peut être implémentée en nombres entiers, mais elle nécessite une division à chaque étape. Mais peut-il en être autrement? Trouvez la racine carrée modulo la puissance de deux, et vérifiez si ce sera une racine carrée ordinaire.

Vous pouvez vous limiter aux nombres impairs: pour un nombre pair, si le nombre de chiffres zéro les moins significatifs est impair, alors il n'y a pas de racine, mais s'il est pair, alors vous pouvez déplacer le nombre vers la droite, calculer la racine de l'impaire et revenir à la moitié gauche du nombre d'origine de zéro bits.

Pour N impairs et 2 ^k , k> 3, si N ≡ 1 mod 8, alors il y a 4 racines différentes modulo 2 ^k , sinon il n'y a pas de racines. Nous avons besoin de la plus petite des quatre racines de x. Les trois autres racines sont 2 ^k - x, 2 ^k-1 + x et 2 ^k - 2 ^k-1 - x

Je voudrais quelque chose de similaire au calcul du modulo inverse 2 ^k - doublant le nombre de bits valides par itération.

Supposons que nous ayons déjà une racine x ₀ de N modulo 2 ^k : N - x ₀ ² = 2 ^k a
Et nous voulons trouver x ₁ = x ₀ + 2 ^k-1 y , de telle sorte que dans N - x ₁ ^2, il y ait plus de bits zéro inférieurs.
N - (x ₀ + 2 ^k-1 y) ² = 2 ^k a - 2 ^k x ₀ * y - 2 ^2k-2 y ²
Diviser par 2 ^k : a - x ₀ * y - 2 ^k-2 y ²
Et équivaut à 0 modulo 2 ^k-2 : y = a * x ₀ ^-1 mod 2 ^k-2
Reçu x ₁ = x ₀ + 2 ^k-1 a * (x ₀ ^-1 mod 2 ^k-2 )
Et enfin x ₁ = x ₀ + (N - x ₀ ² ) / 2 * (x ₀ ^-1 mod 2 ^k-2 )

Des k bits de l'itération suivante, vous obtenez 2 (k-1) bits. En même temps, nous considérons l'opposé de la racine à chaque itération.

Code de test:

uint8_t sqr16(uint16_t n) { if (n % 8 != 1) return 0; uint16_t sqr = (n + 1) / 2; //4 bit uint16_t inv = 2 - sqr; sqr += inv * (n-sqr*sqr)/2; //6 bit inv *= 2 - sqr * inv; sqr += inv * (n-sqr*sqr)/2; //10 bit //inv *= 2 - sqr * inv; if (sqr & 256) sqr = 0u - sqr; sqr = (uint8_t)sqr; // lowest root if (n == sqr*sqr) return sqr; return 0; } 

En ajoutant quelques itérations, nous obtenons la racine de uint_64