jacobi  平方剰余記号 (a/p) の値を計算する

int jacobi(long a, long p);

a  平方剰余記号の分子
p  平方剰余記号の分母

平方剰余記号の値 +1 または -1。引数により計算不能の場合は 0。

int jacobi(long a, long p)
{
    int  x = 1;
    long tmp;
    
    if (p <= 0) return 0;
    if (p > 2 && (p & 1) == 0) return 0;
    if ((a = a % p) <= 0) return 0;
   
    while (a > 1) {
        if (a & 1) {
            if ((a & 3) == 3 && (p & 3) == 3) x = -x;
            tmp = p;
            p = a;
            a = tmp;
            a = a % p;
        } else {
            a >>= 1;
            if ((p & 7) == 3 || (p & 7) == 5) x = -x;
        }
    }
    if (a == 0) return 0;
    return x;
}

　　　　x²≡a　(mod p)

なる x は存在するときと、存在しないときがある。たとえば、p = 7 のとき
　　　　1²≡1、 2²≡4、3²≡2、 4²≡2、 5²≡4、6²≡1
であるから、a≡1 または 2 または 4 のときのみ解がある。解があるとき、 a を mod p での平方剰余という。a がある数を平方して p で割った余りと なるからである。解がないときは、平方非剰余という。このとき、

　　　　(a/p)= 1 　　解があり
　　　　(a/p)=-1 　　解がなし

と定義される。また、この (a/p) を平方剰余記号という。

平方剰余記号の値を計算するのに、以下の性質が利用される。

　　　　(1/p) = 1
　　　　(ab/p) = (a/p)(b/p)

　　　　(2/p) = 1、p≡1 または 7 (mod 8) 時
　　　　(2/p) = -1、p≡3 または 5 (mod 8) 時

　　　　奇数 a と互いに素な p に対して、ガウス定理より
　　　　(a/p) = (p/a)、p≡1 または a≡1 (mod 4) 時
　　　　(a/p) = -(p/a)、p≡3 かつ a≡3 (mod 4) 時

上の性質から、a の値が奇数か偶数かで異なる処理をする。偶数の場合は、 (2/p)の値をかけておき、a 自身の値を半分にする。奇数の場合は、 分母と分子を入れ替え、ユークリッド互除法を適用する。この繰り返しを a が 1 になるまで行う。

したがって、本アルゴリズムでは平方剰余記号の値を log(p) に比例時間で 計算される。