Tonelli–Shanks Algorithm 二次剩余系解法 (Ural 1132. Square Root)

wikipedia上的解释和证明：http://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm

The Tonelli–Shanks algorithm (referred to by Shanks as the RESSOL algorithm) is used within modular arithmetic to solve a congruence of the form

x^2 \equiv n \pmod p

where n is a quadratic residue (mod p), and p is an odd prime.

Tonelli–Shanks cannot be used for composite moduli; finding square roots modulo composite numbers is a computational problem equivalent to integer factorization.[1]

An equivalent, but slightly more redundant version of this algorithm was developed by Alberto Tonelli in 1891. The version discussed here was developed independently by Daniel Shanksin 1973, who explained:

"My tardiness in learning of these historical references was because I had lent Volume 1 of Dickson's History to a friend and it was never returned."[2]


(Note: All \equiv are taken to mean \pmod p, unless indicated otherwise).[edit]The algorithm

Inputs: p, an odd prime. n, an integer which is a quadratic residue (mod p), meaning that the Legendre symbol \bigl(\tfrac{n}{p}\bigr)=1.

Outputs: R, an integer satisfying R^2 \equiv n.

  1. Factor out powers of 2 from p − 1, defining Q and S as: p-1 = Q2^S with Q odd. Note that if S = 1, i.e. p \equiv 3 \pmod 4, then solutions are given directly by R \equiv \pm n^{\frac{p+1}{4}}.
  2. Select a z such that the Legendre symbol \bigl(\tfrac{z}{p}\bigr)=-1 (that is, z should be a quadratic non-residue modulo p), and set c \equiv z^Q.
  3. Let R \equiv n^{\frac{Q+1}{2}}, t\equiv n^Q, M = S.
  4. Loop:
    1. If t \equiv 1, return R.
    2. Otherwise, find the lowest i, 0 < i < M, such that t^{2^i} \equiv 1; e.g. via repeated squaring.
    3. Let b \equiv c^{2^{M-i-1}}, and set R \equiv Rb, \; t \equiv tb^2, c \equiv b^2 and M =\; i.

Once you have solved the congruence with R the second solution is p − R.

Example

Solving the congruence  x^2 \equiv 10 \pmod {13} . It is clear that 13 is odd, and since 10^{\frac{13-1}{2}} = 10^6 \equiv 1 \pmod {13}, 10 is a quadratic residue (by Euler's criterion).

  • Step 1: Observe p-1 = 12 = 3 \cdot 2^2  so Q=3, S=2.
  • Step 2: Take z=2 as the quadratic nonresidue (2 is a quadratic nonresidue since 2^{\frac{13-1}{2}} = -1 \pmod {13} (again, Euler's criterion)). Set  c = 2^3 \equiv 8 \pmod {13}.
  • Step 3: R=10^2 \equiv -4, \; t\equiv 10^3 \equiv -1 \pmod {13}, M = 2.
  • Step 4: Now we start the loop:  t \not\equiv 1 \pmod {13} so 0 < i <\; 2; i.e. i = \;1.
    • Let  b \equiv 8^{2^{2-1-1}} \equiv 8 \pmod {13}, so b^2 \equiv 8^2 \equiv -1 \pmod {13}.
    • Set R=-4\cdot8 \equiv 7 \pmod {13} . Set t \equiv -1 \cdot -1 \equiv 1 \pmod {13}, and M =\;1.
    • We restart the loop, and since t \equiv 1 \pmod{13} we are done, returning R\equiv7 \pmod {13}.

Indeed, observe that 7^2 = 49 \equiv 10 \pmod {13}  and naturally also (-7)^2 \equiv 6^2 \equiv 10 \pmod {13} . So the algorithm yields two solutions to our congruence.

Proof

First write p-1=Q2^S. Now write r \equiv n^{\frac{Q+1}{2}}\pmod p and t \equiv n^Q \pmod p, observing that r^2 \equiv nt \pmod p. This latter congruence will be true after every iteration of the algorithm's main loop. If at any point, t \equiv 1 \pmod p  then r^2 \equiv n \pmod p  and the algorithm terminates with R \equiv \pm r \pmod p.

If t \not\equiv 1 \pmod p , then consider z, a quadratic non-residue of p. Let c \equiv z^Q \pmod p. Then c^{2^S} \equiv (z^Q)^{2^S} \equiv z^{2^SQ}\equiv z^{p-1} \equiv 1 \pmod p and  c^{2^{S-1}} \equiv z^\frac{p-1}{2}\equiv -1 \pmod p, which shows that the order of c is 2^S.

Similarly we have t^{2^S} \equiv 1 \pmod p, so the order of t divides 2^S. Suppose the order of t is 2^{S'}. Since n is a square modulo p, t \equiv n^Q \pmod p is also a square, and hence S'\leq S-1 .

Now we set  b \equiv c^{2^{S-S'-1}} \pmod p and with this r' \equiv br \pmod p, c' \equiv b^2 \pmod p and  t' \equiv c't \pmod p. As before, r'^2 \equiv nt' \pmod pholds; however with this construction both t and  c' have order 2^{S'}. This implies that t' has order 2^{S''} with  S'' < S' .

If S'' \equiv 0 \pmod p then t' \equiv 1 \pmod p, and the algorithm stops, returning R \equiv \pm r' \pmod p. Else, we restart the loop with analogous definitions of b', r'', c''and t'' until we arrive at an S^{(j)'} that equals 0. Since the sequence of S is strictly decreasing the algorithm terminates.

 
 
算法的大体过程已经说的很清楚，然后模拟那个过程就可以了，不过对于Ural上的这个题，要对n = 2特判一下。详见代码：
 
  
//#pragma comment(linker,"/STACK:327680000,327680000")
#include <iostream>
#include <cstdio>
#include <cmath>
#include <vector>
#include <cstring>
#include <algorithm>
#include <string>
#include <set>
#include <functional>
#include <numeric>
#include <sstream>
#include <stack>
#include <map>
#include <queue>

#define CL(arr, val)    memset(arr, val, sizeof(arr))
#define REP(i, n)       for((i) = 0; (i) < (n); ++(i))
#define FOR(i, l, h)    for((i) = (l); (i) <= (h); ++(i))
#define FORD(i, h, l)   for((i) = (h); (i) >= (l); --(i))
#define L(x)    (x) << 1
#define R(x)    (x) << 1 | 1
#define MID(l, r)   (l + r) >> 1
#define Min(x, y)   (x) < (y) ? (x) : (y)
#define Max(x, y)   (x) < (y) ? (y) : (x)
#define E(x)        (1 << (x))
#define iabs(x)     (x) < 0 ? -(x) : (x)
#define OUT(x)  printf("%I64d\n", x)
#define Read()  freopen("data.in", "r", stdin)
#define Write() freopen("data.out", "w", stdout);

typedef long long LL;
const double eps = 1e-8;
const double pi = acos(-1.0);
const double inf = ~0u>>2;


using namespace std;

LL MOD;

LL mod_exp(LL a, LL b) {
    LL res = 1;
    while(b > 0) {
        if(b&1)    res = (res*a)%MOD;
        a = (a*a)%MOD;
        b >>= 1;
    }
    return res;
}

LL solve(LL n, LL p) {
    int Q = p - 1, S = 0;
    while(Q%2 == 0) {
        Q >>= 1;
        S++;
    }
    if(S == 1) {return mod_exp(n, (p + 1)/4);}
    int z;
    while(1) {
        z = 1 + rand()%(p - 1);
        if(mod_exp(z, (p - 1)/2) != 1)   break;
    }
    LL c = mod_exp(z, Q);
    LL R = mod_exp(n, (Q + 1)/2);
    LL t = mod_exp(n, Q);
    LL M = S, b, i;
    while(1) {
        if(t%p == 1)  break;
        for(i = 1; i < M; ++i) {
            if(mod_exp(t, 1<<i) == 1)    break;
        }
        b = mod_exp(c, 1<<(M-i-1));
        R = (R*b)%p;
        t = (t*b*b)%p;
        c = (b*b)%p;
        M = i;
    }
    return (R%p + p)%p;
}

int main() {
    //Read();

    int T, a, n, i;
    scanf("%d", &T);
    while(T--) {
        scanf("%d%d", &a, &n);
        if(n == 2) {    //***
            if(a%n == 1)    printf("1\n");
            else    puts("No root");
            continue;
        }
        MOD = n;
        if(mod_exp(a, (n-1)/2) != 1)    {puts("No root"); continue; }

        i = solve(a, n);
        if(i == n - i)  printf("%d\n", i);
        else    printf("%d %d\n", min(i, n - i), max(i, n - i));
    }
    return 0;
}

 

 
 POJ 1808 勒让德符号(Legendre symbol)判定
 
View Code 
//#pragma comment(linker,"/STACK:327680000,327680000")
#include <iostream>
#include <cstdio>
#include <cmath>
#include <vector>
#include <cstring>
#include <algorithm>
#include <string>
#include <set>
#include <functional>
#include <numeric>
#include <sstream>
#include <stack>
#include <map>
#include <queue>

#define CL(arr, val)    memset(arr, val, sizeof(arr))
#define REP(i, n)       for((i) = 0; (i) < (n); ++(i))
#define FOR(i, l, h)    for((i) = (l); (i) <= (h); ++(i))
#define FORD(i, h, l)   for((i) = (h); (i) >= (l); --(i))
#define L(x)    (x) << 1
#define R(x)    (x) << 1 | 1
#define MID(l, r)   (l + r) >> 1
#define Min(x, y)   (x) < (y) ? (x) : (y)
#define Max(x, y)   (x) < (y) ? (y) : (x)
#define E(x)        (1 << (x))
#define iabs(x)     (x) < 0 ? -(x) : (x)
#define OUT(x)  printf("%I64d\n", x)
#define Read()  freopen("data.in", "r", stdin)
#define Write() freopen("data.out", "w", stdout);

typedef long long LL;
const double eps = 1e-8;
const double pi = acos(-1.0);
const double inf = ~0u>>2;


using namespace std;

LL MOD;

LL mod_exp(LL a, LL b) {
    LL res = 1;
    while(b > 0) {
        if(b&1)    res = (res*a)%MOD;
        a = (a*a)%MOD;
        b >>= 1;
    }
    return res;
}

int main() {
    //Read();

    LL a, p;
    int T, cas = 0;
    scanf("%d", &T);
    while(T--) {
        scanf("%lld%lld", &a, &p);
        MOD = p;
        printf("Scenario #%d:\n", ++cas);
        if(mod_exp((a%p+p)%p, (p - 1)/2) == 1)  puts("1");    //注意a有可能是负数
        else    puts("-1");
        cout << endl;
    }
}

 

 
 
 
 
 
posted @ 2012-11-01 16:30  AC_Von  阅读(2929)  评论(0)  编辑  收藏  举报