素数のまとめとか

SPOJのPRICにて２位取りました。
　
SPOJ - Prime Checker - Ranking
https://www.spoj.pl/ranks/PRIC/
　
が……１位が取れませんので、日本勢頑張ってねということで、一旦バトンタッチします。
（ネタバレが嫌いな人は、読まないでください）
　
また４月ぐらいになったら時間取れるので、それまでに他の方々にアイデアを出しておいて欲しかったりとか、あるいは２位の僕だけでなく、１位も抜いちゃって欲しいですとの願いをこめて、日記にまとめておきます。
　
■主だったアイデア
　
僕の２位実装は、独自のアイデアは何一つありません。
ネットに出回ってるアイデアを繋ぎ合わせただけです。
　
一番大切な記事は、以下の２つ。
　
PRoxy Diary - 2007-12-27 [ICPC] Prime Checker
http://mono.kmc.gr.jp/~oxy/d/?date=20071227
　
に出てくる26418法と、
　
菊やんの雑記帳 - 2010-02-28 ヒープソートによる素数判定？
http://d.hatena.ne.jp/kikx/20100228
　
に出てくる等差数列上でもエラストテネスのふるいが使えるという事実。
（ヒープソートは重要じゃないので割愛）
　
それから、そこまで大切じゃないけど参考にした記事として、sheさんの以下のMontgomery Reductionの記事。
　
兼雑記 - 2007-12-29 [Program] PRIC 21:20
http://d.hatena.ne.jp/shinichiro_h/20071229
　
■26418法の亜種
　
26418法は亜種が沢山存在しており、oxyさんの132090法より良いものも、探せば見つかります。
　
2008年9月21日に8.35秒を出した時は、26418法（亜種）＋フェルマーテスト（＋Montgomery Reduction）の組み合わせを最適にチューニングしたものでした。（別途、試し割も使ってますが、これは基本技なので割愛）
このときの26418法亜種には、6515670法（＝2*3*5*7*19*23*71）を用いています。
　
菊やんさんの等差数列上のエラストテネスのふるいを使うにあたっては、310270法（＝2*5*19*23*71）を使っています。
　
ちなみに、26418法の超亜種として2法が存在していて、こいつはpricで問題となる数列の逆関数化として使います。
　
■26418法の最適亜種の求め方
　
iを1から65536ぐらいまで変化させつつ、(i*1234567890) mod Rの値をしらみつぶしに見ていきます。（ただし、R=2^{31}）
単純に試し割り目的で使う場合は、小さな素数を沢山約数として持っている方が有利（3の倍数であることはほぼ必須）で、結果として6515670法（＝2*3*5*7*19*23*71）が最も有利という結論でした。
しかし、等差数列上のエラストテネスのふるいでは、等差数列がいかに長続きするかが重要（33,333,333個の制限と、2^31の制限の２つがある）なので、必ずしも沢山の素数を含んでいなくても良いわけです。
でも幸運なことに、そうやって見つけた一番良い条件のものが、5と19と23と71の倍数だったので、これらについては26418法の、まとめて試し割り的な恩恵を受けております。
　
■等差数列上のエラストテネスのふるい
　
考え方は菊やんさんのサイトに書かれているので割愛します。
公差をDとしたときの、各素数を法とする公差の逆数（D^{-1} mod p)は、拡張ユークリッドの互除法を使って求めます。
このD^{-1}をかけることで、現在の位置が求まり、そこから次の倍数地点が分かります。
　
ここまでの実装で、だいたい４秒前後が目安です。
　
■D^{-1}とRとR^{-1}
　
Montgomery Reductionの処理自体でR^{-1} mod pがかかってしまうので、これにR mod pをかける処理を加えなきゃいけません。
なので、(D^{-1}・R) mod pを、前処理時点で求めておけばokです。
　
これで3.2秒前後になると思います。
　
■その他の工夫
　
・３３Ｍの出力を高速化する試み
・作業メモリを効率的に使う試み（ページ切り替えのコスト削減）
　
やってない・できなかった試み
　
・mmap　（stdoutに対して適用できなかった…適用する方法を知ってる人は情報下さい）
・前処理のプリプロセッサ化　（0.1秒切ってるので、後回しにしてます）
・Barrett Reduction
・SSEとMMXによる高速化　（spojサイトはSSE2に対応していません。SSE1とMMXになら対応しているので、この辺を使って高速化できないかなぁと……思ったのですが、思いつかないので、まだ試してません。ちなみにSSE2にさえ対応してもらえたら、フェルマーテストが２倍早くなる実装を既に作ってあったりします……もうフェルマーテストはしてないので、何の意味もないですが）
　
■spojのcpuid
　
pricでは、１回のプログラム実行につき0〜33,333,333の値を返し、その結果を見ることができます。
（全ての答えを求めておいて、最後に出力する文字数を変える）
私は利便上、0〜65535の値を返す様にして、１回の実行につき16bitの情報をサーバーから抜き出す様にして、内部情報を得る様にしています。
　
と言っても、spojから抜くべき情報はcpuidぐらいですね。
　
2011年元旦時点のspojのcpuidの結果を、以下にまとめます。
　
eax=0でcpuidを呼び出した時の結果：

	上位	下位
eax	0	2
ebx	30062	25927
ecx	27749	29806
edx	18789	28265

　
ebx->edx->ecxの順に下位からASCIIコードとして読んでいくと、GenuineIntelと読めます。
　
eax=1でcpuidを呼び出した時の結果：

	上位	下位
eax	0	1697
ebx	0	3
ecx	0	0
edx	899	64511

　
eax=1の時のedxの値から、SSE1とMMXは使えます。SSE2は使えません。
おそらくPentium3です。
　
（なお、2008年9月にやってみた時も、ほぼ同じ結果だったと記憶しています）