2015年9月

欧几里得算法

欧氏算法用来求出$ A, B $二数的最大公约数,其实就是辗转相除法
首先,下文标记的GCD = Greatest Common Divisor != CCP。
算法建立在一个显然的结论上:$$ GCD(A, B) = GCD(B, A \\, \mathbf{mod} \\, B) $$
即$$ GCD(A, B) | A \\, \mathbf{mod} \\, B $$
证明如下。
因为$ GCD(A, B) | A, B $,而$ A = \lfloor \frac{A}{B} \rfloor B + A \\, \mathbf{mod} \\, B $,等式左边整除$ GCD(A, B) $,加号左边也是,所以,加号右边整除$ GCD(A, B) $,即$ GCD(A, B) | A \\, \mathbf{mod} \\, B $。
以上咯。

template<class T>
T gcd(T a, T b)
{
    return b == 0 ? a : gcd(b, a % b);
}

拓欧

拓欧是来求解丢翻图方程$$ Ax+By=(A, B) $$
首先我们通过普通欧氏算法可以知道,算法必有一个终止态,即$ A_0 = (A, B), x_0=1, B_0=0, y_0=0 $。那我们将其视为特解,就得思考如何反推回原方程。
假设我们已经解出方程$ B_0x+(A_0 \\, \mathbf{mod} \\, B_0)y=(A, B) $的解,要反推回$ A_0x+B_0y=(A, B) $,则

$$ \begin{align} 记d &= (A, B) \\\\ 则d &= Bx_1 + (A \\, \mathbf{mod} \\, B)y_1 \\\\ &= Bx_1 + (A - \lfloor \frac{A}{B} \rfloor B)y_1 \\\\ &= Bx_1 + Ay_1 - \lfloor \frac{A}{B} \rfloor B y_1 \\\\ &= Ay_1 + B(x_1 - \lfloor \frac{a}{B} \rfloor y_1) \\\\ \end{align} $$

代码

template<class T>
void gcd(T a, T b, T& x, T& y, T& d)
{
    if(b == 0)
        d = a, x = 1, y = 0;
    else
        gcd(b, a % b, y, x, d), y -= x * (a / b);
}

快速幂

前几天有人想来这拿逆元的代码,然后没找到还鄙视我。被我回了一句你快速幂和拓欧都不会打么,然后他说不懂使用姿势……
其实我觉得我讲得够清楚了……一定是那家伙理解能力有问题……
那我就想为了方便这种懒人干脆就把这些基础全都写出来算了……
然后发现要我明确说出所以然还是挺困难的……还好想了一会儿想懂了……不过拓欧还没懂……


快速幂是能用来在$ O(\lg n) $的时间内求出$ a^n $的值。
首先我们知道,只用二的幂次方的数(以下简称二的幂)的和是能够表示所有正整数的,这即是二进制。
所以就能把这样的$ n $分解成二进制,或者说分成$ 2^0+2^1+2^2+ \cdots$(当然不是每一个幂都有)。
所以只需要让$ a $每次反复乘以自身($ a^n \to a^{2n} $),在需要当前幂时再把之乘入到最终答案里。
具体来说的话将$ n $二进制分解后从后往前看,当前位若为$ 1 $则乘,否则罢了。
复杂度显而易见是$ O(\lg n) $。

代码

template<class T>
T power(T a, T n, T P)
{
    T ans = 1;
    while(n > 0)
    {
        if(n & 1)
            ans = ans * a % P;
        a = a * a % P;
        n >>= 1;
    }
    return ans;
}

Baby Step Giant Step及其拓展算法

嗯没错,咱又来讲数学了QuQ。

BSGS

原始的BSGS算法是用来解决这样的一类问题:求关于$ x $的方程$$ A^x \equiv B \pmod P $$,其中$ 0 \le x \le P $,且$ P $为质数。

其实算法的思路很简单。

Giant Step

首先记$ m = \lceil \sqrt P \rceil $,则可将$ x $分解为$ x = i \times m + j $,即$ A^x=(A^m)^i \times A^j $,其中$ 0 \le i,j < m $。然后就可以枚举$ i $得到$ j $,$ O(\sqrt P) $级别的枚举。
令$ D=(A^m)^i $,即$ D A^j \equiv B \pmod P $,将$ A^j $视为一个整体解出来,即$ A^j \equiv BD^{-1} \pmod P $。$ D $的逆元可以很方便求出。线性筛逆元那一节我们提到过逆元可积,所以$ A^{-mi}=(A^{-m})^{i} $,所以只要一开始用$ O(\lg P) $计算$ A^{-m} $,接下来直接就有$ D_i^{-1}=D_{i-1}^{-1} \times A^{-m} $。

Baby Step

那得到$ A^j $怎么得到$ j $ 呢?开始预处理的时候只须$ O(\sqrt P) $枚举$ A^j $加入到一个hash表里,然后$ O(1) $查询即可。
总复杂度为$ O(\sqrt P + \lg P) $。

代码(POJ2417,9124K,79ms)

#include <iostream>
#include <cmath>
#include <cstring>
using namespace std;
template<int hashmod = 1 << 20 >
class HashTable
{
private:
    int hash[hashmod + 10], stack[1 << 16], stacktop;
    int value[hashmod + 10];
    int locate(int x)
    {
        int h = x % hashmod;
        while(hash[h] != -1 && hash[h] != x)
            h++, h %= hashmod;
        return h;
    }
public:
    HashTable() : stacktop(0){memset(hash, -1, sizeof hash);}
    void insert(int x, int v)
    {
        int h = locate(x);
        if(hash[h] == -1)
            hash[h] = x, value[h] = v, stack[stacktop++] = h;
    }
    int get(int x)
    {
        int h = locate(x);
        if(hash[h] != -1)
            return value[h];
        return -1;
    }
    void clear(){while(stacktop)hash[stack[--stacktop]] = -1;}
};
HashTable<> hash;

long long power(long long a, long long b, long long c)
{
    long long ans = 1;
    a %= c;
    while(b > 0)
    {
        if(b & 1)
            ans *= a, ans %= c;
        a *= a, a %= c;
        b >>= 1;
    }
    return ans;
}

long long BSGS(long long A, long long B, long long C)
{
    long long sqrtC = ceil(sqrt(C));
    long long base = 1;
    hash.clear();
    for(int i = 0; i < sqrtC; i++, base *= A, base %= C) //Baby Step
        hash.insert(base, i);
    long long inv = power(base, C - 2, C);
    for(int i = 0, invD = 1, j; i < sqrtC; i++, invD = invD * inv % C) //Giant Step
        if((j = hash.get(B * invD % C)) != -1)
            return i * sqrtC + j;
    return -1;
}

int main()
{
    long long A, B, C;
    while(cin>>C>>A>>B)
    {
        long long ans = BSGS(A, B, C);
        if(ans == -1)
            cout<<"no solution"<<endl;
        else
            cout<<ans<<endl;
    }
    return 0;
}

拓展BSGS

拓展后的BSGS,可解决模数不为质数的情况,下文标记为$ C $,即$$ A^x \equiv B \pmod C $$
首先有这样一个性质:$$ ad \equiv bd \pmod {cd} \Leftrightarrow a \equiv b \pmod c $$易证,略。
所以我们可以先执行消因子:

while GCD(A, C) != 1
    if B mod GCD(A, C) != 0 then
        No Solution
    B /= GCD(A, C)
    C /= GCD(A, C)
    D *= A / GCD(A, C)
    cnt ++

注意到其实$ A $是有$ x $个,所以才要执行若干次消元。
若$ A, C $不互质,则不存在$ A $不存在关于$ C $的逆元(证明见),显然以下无解。
然后问题就变成了$ D \cdot A^{j-cnt} \equiv B \pmod C $。即$ x = i \times m + j + cnt $忽略掉$ cnt $后和上述BSGS就一模一样了。
但还是有一个问题,因为我们默认了$ x \ge cnt $。所以直接$ O(\lg n) $枚举小于cnt的情况,即直接验证是否$ A^i \equiv B \pmod C $,特判即可。
总复杂度仍为$ O(\sqrt C \lg C) $。【好像还是哪里有问题】

代码(POJ3243)

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

template<class T = int, int hashmod = 3894229>
class HashTable
{
private:
    int hash[hashmod + 10], stack[1 << 16], stacktop;
    T value[hashmod + 10];

    int locate(int x)
    {
        int h = x % hashmod;
        while(hash[h] != -1 && hash[h] != x)
            h++;
        return h;
    }
public:
    HashTable() : stacktop(0){memset(hash, -1, sizeof hash);}
    void insert(int x, int v)
    {
        int h = locate(x);
        if(hash[h] == -1)
            hash[h] = x, value[h] = v, stack[stacktop++] = h;
    }
    int get(int x)
    {
        int h = locate(x);
        return hash[h] == x ? value[h] : -1;
    }
    void clear(){while(stacktop) hash[stack[--stacktop]] = -1;}
};
HashTable<> hash;

void gcd(long long a, long long b, long long& x, long long& y, long long& d)
{
    if(b == 0)
        d = a, x = 1, y = 0;
    else
        gcd(b, a % b, y, x, d), y -= x * (a / b);
}

long long gcd(long long a, long long b)
{
    return b == 0 ? a : gcd(b, a % b);
}

long long BSGS(long long A, long long B, long long C)
{
    for(long long tmp = 1, i = 0; i < 32; i++, tmp = tmp * A % C)
        if(tmp == B)
            return i;
    long long cnt = 0, D = 1;
    for(int i; (i = gcd(A, C)) != 1; cnt++)
    {
        if(B % i != 0)
            return -1;
        B /= i, C /= i;
        D = D * A / i % C;
    }

    int sqrtC = static_cast<int>(ceil(sqrt(C)));
    long long base = 1;
    hash.clear();
    for(int i = 0; i < sqrtC; i++, base = base * A % C)
        hash.insert(base, i);
    for(long long i = 0; i < sqrtC; i++)
    {
        long long x, y, d;
        gcd(D, C, x, y, d);
           long long j = hash.get((x + C) % C * B % C);
        if(j != -1)
            return i * sqrtC + j + cnt;
        D = D * base % C;
    }
    return -1;
}

int main()
{
    long long A, B, C;
    while(cin>>A>>C>>B && (A || B || C))
    {
        B %= C;
        long long ans = BSGS(A, B, C);
        if(ans == -1)
            cout<<"No Solution"<<endl;
        else
            cout<<ans<<endl;
    }
    return 0;
}

POJ3243另一份代码

#include <iostream>
#include <cmath>
#include <cstring>
using namespace std;
template<int hashmod = 1 << 21>
class HashMap
{
private:
    int hash[hashmod + 10], stack[1 << 22], stacktop;
    int value[hashmod + 10];
    int locate(int x)
    {
        int h = x % hashmod;
        while(hash[h] != x && hash[h] != -1)
            h++, h %= hashmod;
        return h;
    }
public:
    HashMap():stacktop(0){memset(hash, -1, sizeof hash);}
    void insert(int x, int v)
    {
        int h = locate(x);
        if(hash[h] == -1)
            hash[h] = x, value[h] = v, stack[stacktop++] = h;
    }
    int get(int x)
    {
        int h = locate(x);
        if(hash[h] == x)
            return value[h];
        return -1;
    }
    void clear(){while(stacktop){hash[stack[--stacktop]] = -1;}}
};
HashMap<> hash;

template<class T>
T gcd(T a, T b, T& x, T& y)
{
    if(b == 0)
    {
        x = 1, y = 0; return a;
    }
    T ans = gcd(b, a % b, y, x);
    y -= x * (a / b);
    return ans;
}

template<class T>
T BSGS(T A, T B, T C)
{
    T n1, n2;
    for(T tmp = 1, i = 0; i < 32; i++, tmp = tmp * A % C)
        if(tmp == B)
            return i;
    T cnt = 0, D = 1;
    for(T i; (i = gcd(A, C, n1, n2)) != 1; cnt++)
    {
        if(B % i != 0)
            return -1;
        B /= i, C /= i;
        D = D * A / i % C;
    }

    T sqrtC = ceil(sqrt(C));
    T base = 1;
    hash.clear();
    for(T i = 0; i < sqrtC; i++, base = base * A % C)
        hash.insert(base, i);
    T inv, invD;
    gcd(base, C, inv, n1); inv = (inv + C) % C;
    gcd(D, C, invD, n1); invD = (invD + C) % C;
    for(T i = 0, j; i < sqrtC; i++, invD = invD * inv % C)
        if((j = hash.get(invD * B % C)) != -1)
            return i * sqrtC + j + cnt;
    return -1;
}

int main()
{
    long long A, B, C;
    while(cin>>A>>C>>B && (A || B || C))
    {
        B %= C;
        long long ans = BSGS(A, B, C);
        if(ans != -1)
            cout<<ans<<endl;
        else
            cout<<"No Solution"<<endl;
    }
    return 0;
}

其它运用

很容易就可以猜到的。其实就是计算对数。求$ \log_a b $的值$ x $。变形即得$ a^x=b $。