[POJ1845&POJ1061]扩展欧几里得应用两例

扩展欧几里得是用于求解不定方程、线性同余方程和乘法逆元的常用算法。

下面是代码：

function Euclid(a,b:int64;var x,y:int64):int64;
var t:int64;
begin
    if b=0 then
    begin
        x:=1;y:=0;exit(a);
    end else
    begin
        Euclid:=Euclid(b,a mod b,x,y);
        t:=x;x:=y;y:=t-(a div b)*y;
    end;
end;

 

下面出现了其中后两个应用。（虽然个人认为不定方程和同余方程可以互相转换）

---

POJ1061 

　　线性同余简单应用。

program poj1061;
var x0,y0,m,n,L,a,b,d,x,y,tem,c:int64;

function f(a:int64):int64;
begin
    if a<0 then a:=a+(-a) div L*L;
    while a<=0 do a:=a+L;
    exit(a);
end;

function Euclid(a,b:int64;var x,y:int64):int64;
var t:int64;
begin
    if b=0 then
    begin
        x:=1;y:=0;exit(a);
    end else
    begin
        Euclid:=Euclid(b,a mod b,x,y);
        t:=x;x:=y;y:=t-(a div b)*y;
    end;
end;

begin
    readln(x0,y0,m,n,L);
    a:=f(m-n);b:=L;    
    d:=Euclid(a,b,x,y);
    c:=f(y0-x0);
    // 由于扩展欧几里得算法里要求最大公约数，所以要处理成正数
    // 根据同余方程的性质，a,b,c加上L对答案不会有影响。
    if c mod d<>0 then 
    begin
        writeln('Impossible');
        // 同余方程有整数解的条件。
        halt;
    end;
    x:=x*c div d;
    if x<0 then x:=x+(-x) div (b div d)*(b div d);
    while x<0 do x:=x+(b div d);
    //这里是保证跳的步数为自然数
    if x-(b div d)>=0 then x:=x-x div (b div d)*(b div d);
    while x-(b div d)>=0 do x:=x-(b div d);
    //这里是保证跳的步数是最小步数
    writeln(x);
end.

---

POJ1845

　　由于姿势不大对...交了很多遍处理了很多细节才过。

　　因此觉得这道题挺好的...因为本来觉得只是一道求逆元的模板题，但实际上想要A掉它并不容易。

　　首先先讲一下算法部分。

　　求a^b的约数和。

　　= (p1^0+p1^1+...p1^k1)*(p2^0+p2^1...+p2^k2)...*(pn^0+pn^1+...pn*kn) (p为a的质因子，ki为第i个质因子的个数)

　　我们可以用素数拆分，然后对于每一部分，用等比数列求和公式。

　　分子部分可以用取模的快速幂解决。

　　分母部分乘上其逆元。

　　大体的思路已经构建好了。

　　但是细节超级超级多.....>_< 已经被搞疯啦...直接在代码里贴出来好了。

---

program poj1845;
const tt=9901;
var a,b,ans,k:int64;
      i:longint;

function Euclid(a,b:int64;var x,y:int64):int64;
var t:int64;
begin
    if b=0 then
    begin
        x:=1;y:=0;exit(a);
    end else
    begin
        Euclid:=Euclid(b,a mod b,x,y);
        t:=x;x:=y;y:=t-(a div b)*y;
    end;
end;

function mul(a,b:int64):int64;
var ans,w:int64;
begin
    ans:=1;w:=a mod tt;
    while b>0 do 
    begin
        if b and 1=1 then ans:=(ans*w)mod tt;
        w:=(w*w) mod tt;
        b:=b >> 1;
    end;
    exit(ans);
end;

function inverse(p:int64):int64;
var x,y:int64;
begin
    if Euclid(p,tt,x,y)<>1 then exit(-1) else
    begin
        while x<tt do inc(x,tt);            
        x:=x mod tt;
        // 我们需要找到一个正的逆元，根据x= x0+b/d*t可以得到，d=1,b=tt。所以只要不停加上tt即可。
        // 虽然刚开始的直觉就是不停加上tt直到找到一个正的为止...但是实际上是有理论依据的。
        exit(x);
    end;
end;
    
begin
    while not eof do 
    // 由于被多组数据坑过...所以不管怎样还是加上为好
    begin
        readln(a,b);
        if a=0 then 
        begin
            writeln(0);
            // 这是一个细节需要特判
            continue;
        end;        
        ans:=1;
        for i:=2 to trunc(sqrt(a)) do if a mod i=0 then 
        begin
            k:=0;
            while a mod i=0 do 
            begin
                inc(k);a:=a div i;
            end; 
            k:=k*b;
            ans:=((ans*((mul(i,k+1)+tt-1) mod tt)) mod tt*inverse(i-1))mod tt;
            // 本来这里的inverse(i-1）也需要处理i mod tt=1的情况，但是由于数据范围i的上限正好是7000多不用考虑。
    end;
    // 刚开始下面的一块都忘了QAQ
        if a<>1 then 
        begin
            if a mod tt=1 then 
            begin
                ans:=ans*(b+1) mod tt;            
            // 这是最容易忽略的一个细节。当a mod tt=1的时候，inverse(i)是算不出来的，这个时候思考一下便发现，a mod tt=1
            // a^k mod tt=1，所以最后mod tt余数就是b+1.
            // 刚开始由于直接输出了（b+1） WA了一次...后来发现剩下的a不一定是a本身。
        end else
            begin
                i:=a;k:=b;        
                ans:=((ans*((mul(i,k+1)+tt-1) mod tt)) mod tt*inverse(i-1)) mod tt;
            end;
        end;    
        writeln(ans);
    end;
end.