Perhitungan akar kuadrat integer

Ada kebutuhan untuk memeriksa apakah bilangan bulat adalah persegi, dan jika demikian, hitung root. Dan saya ingin melakukan ini dalam hitung bilangan bulat. Jelas bahwa metode Newton dapat diimplementasikan dalam bilangan bulat, tetapi membutuhkan pembagian di setiap langkah. Tapi bisakah sebaliknya? Temukan modulo root kuadrat kekuatan dua, dan periksa apakah itu akan menjadi root kuadrat biasa.

Anda dapat membatasi diri Anda ke bilangan ganjil: untuk bilangan genap, jika jumlah nol nol paling signifikan adalah ganjil, maka tidak ada root, tetapi jika bilangan genap, maka Anda dapat menggeser bilangan ke kanan, menghitung akar ganjil, dan menggeser kembali ke setengah kiri jumlah awal nol bit.

Untuk N aneh dan 2 ^k , k> 3, jika N ≡ 1 mod 8, maka ada 4 akar modulo 2 ^{k yang berbeda} , jika tidak maka tidak ada akar. Kita membutuhkan yang terkecil dari empat akar x. Tiga akar lainnya adalah 2 ^k - x, 2 ^k-1 + x dan 2 ^k - 2 ^k-1 - x

Saya ingin sesuatu yang mirip dengan menghitung inversi modulo 2 ^k - menggandakan jumlah bit yang valid per iterasi.

Misalkan kita sudah memiliki root x ₀ dari N modulo 2 ^k : N - x ₀ ² = 2 ^k a
Dan kami ingin menemukan x ₁ = x ₀ + 2 ^k-1 y sedemikian rupa sehingga dalam N - x ₁ ² ada lebih banyak bit nol yang lebih rendah.
N - (x ₀ + 2 ^k-1 y) ² = 2 ^k a - 2 ^k x ₀ * y - 2 ^2k-2 y ²
Bagi dengan 2 ^k : a - x ₀ * y - 2 ^k-2 y ²
Dan sama dengan 0 modulo 2 ^k-2 : y = a * x ₀ ^-1 mod 2 ^k-2
Diterima x ₁ = x ₀ + 2 ^k-1 a * (x ₀ ^-1 mod 2 ^k-2 )
Dan akhirnya x ₁ = x ₀ + (N - x ₀ ² ) / 2 * (x ₀ ^-1 mod 2 ^k-2 )

Dari bit k dalam iterasi berikutnya, Anda mendapatkan 2 (k-1) bit. Pada saat yang sama, kami mempertimbangkan kebalikan dari akar pada setiap iterasi.

Kode uji:

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

Menambahkan beberapa iterasi, kami mendapatkan root dari uint_64