Когда то давно, во время ковыряния «в» и изучения «как» очень хорошего и полезного пакета OpenSSL и как всегда неожиданно возникла одна простая идея и как все такие очень неожиданные идеи канула в лету.

Но сухой остаток остался — была найдена ошибка в OpenSSL, в умножении большого числа на BN_ULONG и небольшая программа извлечения квадратного корня побитно. Сообщение об ошибке ушло на багтрекинг и было поправлено (пользуясь случаем извиняюсь за свою излишнюю эмоциональность тогда, не каждый день в OpenSSL ошибки находишь), а вот ту самую небольшую программу нахождения квадратного корня побитово по модулю 2^n, где n это количество бит\разрядность и предлагаю вашему вниманию.

Если рассмотреть алгоритм умножения двух чисел побитно, в столбик, то по значению i-х битов сомножителей, зная перенос, можно быстро определить бит результата — так и зная результат и предполагая распределение битов в сомножителях можно быстро эти биты вычислить. Как всегда ничего путного не получилось, но при равенстве сомножителей, т.е. при извлечении квадратного корня, соответствующий бит сомножителя (он теперь один) можно быстро получить,

Очевидно, что корень это или нет станет ясно когда будут получены n/2 битов, или если по модулю степени двойки, то полученные два числа и будут корнями.

Небольшие допущения — считаем, что извлекаем корень из нечетного числа. Если последние биты нули и их четное количество, то можно их отбросить.

Теперь конечно главная часть — код.

Для проверки и испытания взят тот самый OpenSSL для получения больших простых чисел. После число умножается на себя с помощью BN_mul и в функции square_root дважды вычисляется корень. Вычисления проводятся дважды, так как последние биты сомножителя или 11 или 01 неразличимы для этого алгоритма. Для хранения чисел используются или BIGNUM пакета OpenSSL или bitset одинаковой длины.

Итак код

#include <bitset>
#include <stdio.h>
#include <openssl/rsa.h>
#include <openssl/bn.h>

using namespace std;
#define DIM 512

int ret = 0;
RSA *r = NULL;
BIGNUM *bne = NULL;
BIO *bp_public = NULL, *bp_private = NULL;

int bits = DIM;
unsigned long e = RSA_F4;

bool generate_key() {
	r = RSA_new();
	ret = RSA_generate_key_ex(r, bits, bne, NULL);
	return (ret );
}

bitset<DIM> square_root(bitset<DIM> key, int prim_1) {
	bitset<DIM> prim;
	int carry = prim_1;
	int i, j, ie;
	prim[0] = 1;
	prim[1] = prim_1;

	ie = DIM / 2;
	for (i = 2; i < ie; i++) {
		for (j = 1; j < i; j++)
			carry = carry + (int) (prim[j] * prim[i - j]);

		bool i1 = i & 1;
		int q2 = (carry / 2) & 1;
		int key1 = (int) key[i + 1];
		if (!i1 && q2 != key1)
			prim[i] = 1;
		if (!i1 && q2 == key1)
			prim[i] = 0;
		if (i1 && q2 == key1)
			prim[i] = prim[(i + 1) / 2];
		if (i1 && q2 != key1)
			prim[i] = 1 - (int) prim[(i + 1) / 2];

		carry += 2 * (int) prim[i];
		carry /= 2;
	}
	return prim;
}

int main() {
	bitset<DIM> bit0_sqrt(0);
	bitset<DIM> bit1_sqrt(0);
	bitset<DIM> bit_p2(1);
	bne = BN_new();
	ret = BN_set_word(bne, e);
	char *pc, *qc, *rc;

	if (generate_key() == 1) {
		BIGNUM *rez = NULL;
		BIGNUM *p2 = NULL;
		BIGNUM *tmp = NULL;
		const BIGNUM *n = NULL;
		const BIGNUM *e = NULL;
		const BIGNUM *d = NULL;
		const BIGNUM *p = NULL;
		const BIGNUM *q = NULL;

		RSA_get0_key(r, &n, &e, &d);
		RSA_get0_factors(r, &p, &q);

		pc = BN_bn2hex(p);
		printf("  p = 0x%s\n", pc);

		p2 = BN_new();
		tmp = BN_new();
		BN_CTX *ctx;
		ctx = BN_CTX_new();
		BN_CTX_start(ctx);
		BN_mul(p2, p, p, ctx);
		pc = BN_bn2hex(p2);
		printf("p^2 = 0x%s\n", pc);

		bit_p2 = 0;
		for (int i = 0; i < BN_num_bits(p2); i++)
			if (BN_is_bit_set(p2, i))
				bit_p2[i] = 1;

		rez = BN_new();
		bit0_sqrt = square_root(bit_p2, 0);
		for (int i = 0; i < DIM; i++)
			if (bit0_sqrt[i])
				BN_set_bit(rez, i);
			else
				BN_clear_bit(rez, i);
		rc = BN_bn2hex(rez);
		printf(" 0r = 0x%s\n", rc);

		BN_sqr(tmp, rez, ctx);
		rc = BN_bn2hex(tmp);
		printf(" p2 = 0x%s\n", rc);

		bit1_sqrt = square_root(bit_p2, 1);
		for (int i = 0; i < DIM; i++)
			if (bit1_sqrt[i])
				BN_set_bit(rez, i);
			else
				BN_clear_bit(rez, i);
		rc = BN_bn2hex(rez);
		printf(" 1r = 0x%s\n", rc);

		BN_sqr(tmp, rez, ctx);
		rc = BN_bn2hex(tmp);
		printf(" p2 = 0x%s\n", rc);
	}

	BIO_free_all(bp_public);
	BIO_free_all(bp_private);
	RSA_free(r);
	BN_free(bne);
}

Код простой и извиняюсь, что одним куском — но там большая часть main это служебные вызовы OpenSSL, в программе вычисления корня оценка и вычисление битов всего в 10 строк, Ясно, что DIM — это размерность.

   p = 0xE5CBB3DF3D2E3F56A3DEEAE03204A37995970BFD98FE7242EB37BFFC4935BFD3

p^2 = 0xCE4611E425E1B6E3C6594862F0C53A61E3D2460CD86A1709B992806B3B920C89
F2F33861A38225ABFB4A95E65852BEB5930F7968120D65F039A83417531A87E9

    0r = 0x1A344C20C2D1C0A95C21151FCDFB5C866A68F40267018DBD14C84003B6CA402D
    p2 = 0x02AEAA25AB8538367E9B72A28CBBF36EB8A42E11A66D3283E3230072A9268CE3
F2F33861A38225ABFB4A95E65852BEB5930F7968120D65F039A83417531A87E9

    1r = 0xE5CBB3DF3D2E3F56A3DEEAE03204A37995970BFD98FE7242EB37BFFC4935BFD3
    p2 = 0xCE4611E425E1B6E3C6594862F0C53A61E3D2460CD86A1709B992806B3B920C89
F2F33861A38225ABFB4A95E65852BEB5930F7968120D65F039A83417531A87E9

вывод наглядно демонстрирует, что оба числа и 0r и 1r квадратные корни числа p2 по модулю 2 в степени DIM. Т.е. 0r * 0r == p2 mod 2^DIM и 1r * 1r == p2 mod 2^DIM. И к тому же 1r * 1r == p2.

Алгоритм Герона быстр и нет нужды наверно в возне с битами, но как этюд программа хороша.

Этюд в битовых тонах