Связь и интернет Архив Программирование
   
Сделать стартовойСделать закладку            
   ПОИСК  
   
Главная / Алгоритмы / Математика / Теория чисел /
8  Perl
8  PHP
8  JavaScript
8  HTML
8  DHTML
8  XML
8  CSS
8  C / C++
8  Pascal и Delphi
8  Турбо Ассемблер
8  MySQL
8  CASE-технологии
8  Алгоритмы
8  Python
8  Обратная связь
8  Гостевая книга
Новости о мире


Квадратный корень по простому модулю - Программирование от RIN.RU
Квадратный корень по простому модулю

Алгоритм решения уравнения x2=a(mod p), где s - const, p - простое, x - из Zp.


  1. Вычислить символ Лежандра (a,p): (a,p)=-1 - решений нет: a - квадр невычет

  2. Найти случайное b, такое что (b,p)=-1 (b - квадр невычет)

  3. Представить p-1 в виде p-1=2s*t, где t - нечетное

  4. Вычислить (a-1)mod p - обратный элемент в кольце Zp

  5. Положить c=(bt)(mod p), r=(a(t+1)/2)(mod p)

  6. Цикл с i=1 до s-1:
    6.1 Вычислить
    d=[(r2*a-1)^(2s-i-1)][mod p]
    6.2 Если d=-1(mod p), то r=(r*c)mod p
    6.3 Положить c=(c*c)mod p

  7. Ответ: r, -r. Можно выбрать один из них, например, положительный.


В соответствующих статьях даны алгоритм вычисления обратного элемента в кольце p (т.е обратного элемента по модулю p) и алгоритм возведения в целую степень по модулю.


Общая сложность алгоритма - log4(n).


Исходник на Си



/* Author: Pate Williams (c) 1998 */
#include <stdio.h>
#include <stdlib.h>
int JACOBI(long a, long n)
{
int s;
long a1, b = a, e = 0, m, n1;
if (a == 0) return 0;
if (a == 1) return 1;
while ((b & 1) == 0)
b >>= 1, e++;
a1 = b;
m = n % 8;
if (!(e & 1)) s = 1;
else if (m == 1 || m == 7) s = + 1;
else if (m == 3 || m == 5) s = - 1;
if (n % 4 == 3 && a1 % 4 == 3) s = - s;
if (a1 != 1) n1 = n % a1; else n1 = 1;
return s * JACOBI(n1, a1);
}


void extended_euclid(long a, long b, long *x,
long *y, long *d)
/* calculates a * *x + b * *y = gcd(a, b) = *d */
{
long q, r, x1, x2, y1, y2;


if (b == 0) {
*d = a, *x = 1, *y = 0;
return;
}
x2 = 1, x1 = 0, y2 = 0, y1 = 1;
#ifdef DEBUG
printf("------------------------------");
printf("-------------------\n");
printf("q r x y a b ");
printf("x2 x1 y2 y1\n");
printf("------------------------------");
printf("-------------------\n");
#endif
while (b > 0) {
q = a / b, r = a - q * b;
*x = x2 - q * x1, *y = y2 - q * y1;
a = b, b = r;
x2 = x1, x1 = *x, y2 = y1, y1 = *y;
#ifdef DEBUG
printf("%4ld %4ld %4ld %4ld ", q, r, *x, *y);
printf("%4ld %4ld %4ld %4ld ", a, b, x2, x1);
printf("%4ld %4ld\n", y2, y1);
#endif
}
*d = a, *x = x2, *y = y2;
#ifdef DEBUG
printf("------------------------------");
printf("-------------------\n");
#endif
}


long inverse(long a, long b)
/* returns the inverse of a modulo b if it exists 0 otherwise */
{
long d, x, y;
extended_euclid(a, b, &x, &y, &d);
if (d == 1) return x;
return 0;
}
long exp_mod(long x, long b, long n)
/* returns x ^ b mod n */
{
long a = 1, s = x;
while (b != 0) {
if (b & 1) a = (a * s) % n;
b >>= 1;
if (b != 0) s = (s * s) % n;
}
return a;
}


long square_root_mod(long a, long p)
/* returns the square root of a modulo an odd prime p
if it exists 0 otherwise */
{
long ai, b, c, d, e, i, r, s = 0, t = p - 1;
/* is a quadratic nonresidue */
if (JACOBI(a, p) == - 1) return 0;
/* find quadratic nonresidue */
do
do b = rand() % p; while (b == 0);
while (JACOBI(b, p) != - 1);
/* write p - 1 = 2 ^ s * t for odd t */
while (!(t & 1)) s++, t >>= 1;
ai = inverse(a, p);
c = exp_mod(b, t, p);
r = exp_mod(a, (t + 1) / 2, p);
for (i = 1; i < s; i++) {
e = exp_mod(2, s - i - 1, p);
d = exp_mod((r * r % p) * ai % p, e, p);
if (d == p - 1) r = r * c % p;
c = c * c % p;
}
return r;
}


int main(void)
{
long a, p;
printf("x ^ 2 = a mod p\n");
for (;;) {
printf("a or 0 to quit = ");
scanf("%ld", &a);
if (a == 0) break;
printf("p = ");
scanf("%ld", &p);
printf("square_root_mod(%ld, %ld) = %ld\n", a, p, square_root_mod(a, p));
}
return 0;
}




 8  Комментарии к статье  8 8  Обсудить в чате

 
  
  
    Copyright ©  RIN 2003 - 2004      * Обратная связь