PRIC

http://www.spoj.pl/ranks/PRIC/

oxyさんとこ で見たので解いてたんだけどお話になりませんでした。無念。

細かい話は oxyさんの方に出てるので略。

私は Wikipedia 調べたら出てきたミラーラビンとかいうヤツを使うことにしてやってました。フェルマーのヤツよりちょっぴり除算少なめになる感じかなぁと思ってたのですがあんまり変わらんかったみたい。

でまぁとにかく一個素数計算するごとに除算の数が平均 90 回程度とかとても大変なので、それを減らすのをアレコレ考えたり調べたりしてたのは似たような感じなんですが、 oxy さんがうまくいかなかったという Montgomery のなんちゃらというのは効果がありました。

めんどくさいからコードはって終わり。まぁつまりなんか 32bit 整数 k,n,r に対して (k^n)%r という計算が除算 1 回でできるらしいです。なんかすごいですね。

uint ndash(uint r1) {
    uint xn,t;
    uint d=r1;
    xn = r1;
loop:
    t=d*xn;
    if(t==1)return -xn;
    xn=xn*(2-t);
    goto loop;
}

uint mpow(uint e,uint n){
    //uint b=(1LL<<32)%n;
    uint b;
    __asm__("div %1;"
            : "=d"(b) : "r"(n), "d"(1), "a"(0));
    uint a = b<<1;
    uint nd=ndash(n);
    for(int i=31;i;){
        --i;
        __asm__("mul %1;"            // t=b*b
                "mov %%edx, %%ebx;"
                "mul %3;"            // m=t*nd
                "mul %2;"            // t+=m*n
                "adc %%ebx, %%edx;"
                : "=d"(b) : "a"(b), "S"(n), "D"(nd));
        if (b>=n)b-=n;
        if((e>>i)&1){
            __asm__("mov %1, %%eax;"
                    "mul %2;"            // t=a*b
                    "mov %%edx, %%ebx;"
                    "mul %4;"            // m=t*nd
                    "mul %3;"            // t+=m*n
                    "adc %%ebx, %%edx;"
                    : "=d"(b) : "g"(a), "b"(b), "S"(n), "D"(nd)
                    : "%eax");
        }
    }
    __asm__("mov %%eax, %%ebx;"  // t=b
            "mul %3;"            // m=t*nd
            "mul %2;"            // t+=m*n
            "add %%ebx, %%eax;"
            "adc $0, %%edx;"
            : "=d"(b) : "a"(b), "S"(n), "D"(nd));
    return b;
}

n*n'=-1 % r となる数字がいるらしく普通は拡張エラトステネスのふるいかなんかで求めるらしいのですが、それだと除算がそれなりに必要でつまらないなぁと思っていたら、なんかニュートン法でできるとハッカーのたのしみに書いてあったのでそれでやったら速くなりました、と。

これで除算はほとんど影響してないはずなので頑張って SSE2 でやったらもっと速くなるかもしれません。一応 SPOJ のサーバで SSE2 命令でコケなかった気がします。

でもまぁ oxy さんとこに書いてあった 26418 法とかいうのは思いつかなかったのでそこで最後まで差がつきまくってたのですが、 oxy さんのネタバレを見て適当にこれを使ってやると 18 秒付近までは速くなりました。

まぁ今となっては絶望的な差がついてて、こんだけ差が出るってことは sqrt(2**31) までの素数の乗法元求めておいて、それつかって計算とかかなーと実装してみたのですが遅くて話になりませんでした。 26418 的な周期を動的に計算したりとかできるのかなぁ。