关闭

哪些素数可以表示成两平方数之和? hdoj 3542 费马降阶

标签: vbnull
3300人阅读 评论(0) 收藏 举报
分类:

摘自数论概论的内容:

素数的两平方数之和定理:设p是素数,则p是两平方数之和的充要条件是p= 1(mod 4) (或 p = 2).

两平方数之和定理实际上由两个陈述组成:

陈述1:如果p是两平方数之和,则p = 1(mod 4).

证明:设p = a^2 + b^2,p是奇数,所以a,b为一奇一偶,设a为奇数, b为偶数.比如 a= 2*n+1 b = 2*m.p = a^2 + b^2 = 4n^2+4n+1+4m^2 = 1 (mod 4).

陈述2:如果= 1(mod 4),则p是两平方数之和. 这个的证明很麻烦,主要依据费马降阶法,可以参考数论概论第26章。

简单的说,如果= 1(mod 4),不直接获得p是两平方数之和,而是将p的某个倍数表示成两个平方数之和。由二次互反律知x^2=-1(mod p)有一解,令x = a,b = 1,

a*a + b*b = Mp.利用费马降阶不断减小p的倍数使其可以表示两平方数之和,最终使p变成两平方数之和。如何利用已知的a, b, M来产生新的a, b, M.有恒等式:

(v^2+v^2)(a^2+b^2) = (ua+vb)^2 + (va-ub)^2.降阶程序有5个断言,只列出内容:1)a^2 + b^2 = Mp; 应用恒等式,我们选取的u,v满足u=a(mid M), v= b(mod M)

-M/2<= u, v, <= M/2. 于是有,u^2 + v^2 = a^2 + b^2= 0 (mod M),u^2 + v^2能被M整除,设u^2 + v^2 = Mr.其余四个断言陈述:2)r>=1; 3)r < M; 4)ua + vb能被M整除,

5)va-ub能被M整除。

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <cmath>
using namespace std;

typedef __int64 lint;

lint pow_mod(lint r, lint x, lint p) {
    lint pm = 1;
    while (x) {
        if (x&1)
            pm = (pm*r)%p;
        r = r*r%p;
        x >>= 1;
    }
    return pm;
}

int main()
{
    lint p, a, b, r, x, s, M, u, v, k, y, z, pm;
    while (scanf("%I64d", &p) != EOF) {
        if ((p-1)%4)
            printf("Illegal\n");
        else {
            b = 1;
            srand(NULL);
            r = rand()%(p-2)+2;
            x = (p-1)>>2;
            pm = pow_mod(r, x, p);
            while ((pm*pm)%p != p - 1) {
                r = rand()%(p-1)+1;
                pm = pow_mod(r, x, p);
            }
            a = pm;
            s = a*a + b*b;
            while (s != p) {
                M = s/p;
                k = M>>1;
                u = (a%M + M)%M;
                v = (b%M + M)%M;
                if (u > k)
                    u = M - u;
                if (v > k)
                    v = M - v;
                if ((u*a + v*b)%M)
                    swap(a, b);
                y = (u*a + v*b)/M;
                z = (v*a - u*b)/M;
                s = y*y + z*z;
                a = y;
                b = z;
            }
            if (a < 0)
                a = -a;
            if (b < 0)
                b = -b;
            if (a > b)
                swap(a, b);
            printf("Legal %I64d %I64d\n", a, b);
        }
    }
    return 0;
}


1
0

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:103982次
    • 积分:1889
    • 等级:
    • 排名:千里之外
    • 原创:83篇
    • 转载:6篇
    • 译文:0篇
    • 评论:13条