HDU3589-雅可比符号

题意

给出a,n满足n是奇数，$1<a,nleqslant 1000000$，求雅可比符号$left(frac{a}{n} ight)$

分析

可以拿定义硬上，因子分解然后用勒让德符号

勒让德符号的欧拉判别法：设p是奇素数，a是不被p整除的正整数，则

$left(frac{a}{n} ight)=a^{frac{n-1}{2}} (mod\,p)$

雅可比符号是勒让德符号的推广

设n是正奇数，其素幂因子分解式为$n={p_1}^{t_1}{p_2}^{t_2}cdots{p_m}^{t_m}$

令a是与n互素的正整数，则雅可比符号定义是

$left(frac{a}{n} ight)=prodlimits_{i=0}^{m}left(frac{a}{p_i} ight)^{t_i}$

其中等式的右边是勒让德符号

这里给出另外一种方法，类比于欧几里得算法

设a和b是正整数，令$R_0=a,R_1=b$，利用带余除法，并提取出余数中2的最高次幂得

$R_0=R_1q_1+2^{s_1}R_2$

其中$s_1$是非负整数，$R_2$是小于$R_1$的正奇数

反复使用带余除法，并提取出余数中的2的最高次幂得

$R_1=R_2q_2+2^{s_2}R_3$

$R_2=R_3q_3+2^{s_3}R_4$

$cdots$

$R_{n-3}=R_{n-2}q_{n-2}+2^{s_{n-2}}R_{n-1}$

$R_{n-2}=R_{n-1}q_{n-1}+2^{s_{n-1}}cdot 1$

计算雅可比符号的定理：

$left(frac{a}{b} ight)=(-1)^{sumlimits_{i=1}^{n-1}(s_1frac{{R_i}^2-1}{8}+frac{(R_i-1)(R_{i+1}-1)}{4})}$

代码

#include <cstdio>
#include <cstring>
using namespace std;
int r[100],s[100];
int gcd(int a,int b){return b==0?a:gcd(b,a%b);}
int jcb(int a,int b){
    if(gcd(a,b)>1)return 0;
    r[0]=a;r[1]=b;
    memset(s,0,sizeof s);
    int n=1;
    while(r[n++]!=1){
        r[n]=r[n-2]%r[n-1];
        while(r[n]&1^1)r[n]>>=1,s[n-1]++;
    }
    int p=0;
    for(int i=1;i<n;i++)r[i]&=7,s[i]&=1;
    for(int i=1;i<n;i++){
        p^=(s[i]*((r[i]*r[i])-1)/8)&1;
        p^=((r[i]-1)*(r[i+1]-1)/4)&1;
    }
    return 1-2*p;
}
int main(){
    int a,n;
    while(~scanf("%d%d",&a,&n)){
        printf("%d
",jcb(a,n));
    }
    return 0;
}