Cálculo da raiz quadrada inteira

Era necessário verificar se o número inteiro é um quadrado e, se sim, calcular a raiz. E eu quero fazer isso em aritmética inteira. É claro que o método de Newton pode ser implementado em números inteiros, mas requer divisão a cada passo. Mas pode ser de outra forma? Encontre o módulo da raiz quadrada com a potência de dois e verifique se será uma raiz quadrada comum.

Você pode restringir-se a números ímpares: para um número par, se o número de dígitos zero menos significativos for ímpar, não haverá raiz, mas se for par, poderá mudar o número para a direita, calcular a raiz do ímpar e voltar à metade esquerda do número original de zero bits.

Para N e 2 ^k ímpares, k> 3, se N ≡ 1 mod 8, existem 4 raízes diferentes no módulo 2 ^k , caso contrário não há raízes. Precisamos da menor das quatro raízes de x. As outras três raízes são 2 ^k - x, 2 ^k-1 + x e 2 ^k - 2 ^k-1 - x

Eu gostaria de algo semelhante ao cálculo do módulo inverso 2 ^k - dobrando o número de bits válidos por iteração.

Suponha que já tenhamos uma raiz x ₀ de N módulo 2 ^k : N - x ₀ ² = 2 ^k a
E queremos encontrar x ₁ = x ₀ + 2 ^k-1 y de modo que em N - x ₁ ² haja mais bits zero mais baixos.
N - (x ₀ + 2 ^k-1 y) ² = 2 ^k a - 2 ^k x ₀ * y - 2 ^2k-2 y ²
Divida por 2 ^k : a - x ₀ * y - 2 ^k-2 y ²
E equacione 0 módulo 2 ^k-2 : y = a * x ₀ ^-1 mod 2 ^k-2
Recebido x ₁ = x ₀ + 2 ^k-1 a * (x ₀ ^-1 mod 2 ^k-2 )
E finalmente x ₁ = x ₀ + (N - x ₀ ² ) / 2 * (x ₀ ^-1 mod 2 ^k-2 )

Dos k bits na próxima iteração, você obtém 2 (k-1) bits. Ao mesmo tempo, consideramos o oposto da raiz em cada iteração.

Código do teste:

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; } 

Adicionando algumas iterações, obtemos a raiz de uint_64