はじめに

rustの勉強を兼ねて、フェルマーの因数分解法のプログラムを作って、それを少しだけ性能向上をしてみました
タイトルの100倍以上というのは、ナイーブな実装に比べた数値です
フェルマーの因数分解法を知っている人ならわかると思いますが、RSAのような巨大な素数が対象の場合には天文学的なパターンがあるため、100倍はもとより10億倍ぐらい高速化しても誤差みたいなものです
なので実用的な内容ではなく、頭の体操的なものだと思ってください

環境

rustc/cargo のバージョンは 1.85.1 を使用しました。

また関連する依存関係は（すべてのソースで必要になるわけではありませんが）以下の４つがあれば動くと思います。

Cargo.toml

[dependencies]
num-bigint = "0.4.6"
num-integer = "0.1.46"
num-traits = "0.2.19"
once_cell = "1.21.3"

フェルマーの因数分解法とは

ご存じの方も多いと思いますが簡単に説明すると、フェルマーの因数分解法（Fermat's factorization method; 単にフェルマー法と呼ばれることも）とは、同じぐらいの大きさの因数に分解するのが非常に効率的な因数分解の方法です。

$N = pq$ という関係（ここで $p$ および $q$ はともに奇数、かつ $p > q$ とします）のとき

$a$ を「 $p$ と $q$ の平均」： $a = (p + b) / 2$
$b$ を「 $a$ と $p$ （あるいは $q$ ）の差」： $b = (p - q)/2$

とします（RSAのような奇素数×奇素数を想定しているので $p$, $q$ がともに奇数であることを前提にしていますが、実際にはともに偶数でも可。奇数と偶数の組み合わせだと平均が整数にならないのでダメ）。
言い換えると

$p = a + b$
$q = a - b$

の関係になっています。
そうすると、

$N = pq = (a + b) (a - b) = a^2 - b^2$

と表せます。

フェルマーの因数分解法では、

$a_1 = \lceil \sqrt{N} \rceil$
$a_n = a_{n-1} + 1$

として $a_n^2 - N$ が平方数になれば $a = a_n$ 、 $b = \sqrt{ a_n^2 - N }$ となり、因数分解ができるというものです。

RSA の鍵生成の場合には、通常、素数 $p$ と $q$ がある程度以上近い場合には作り直す等の対応をしているので（例えば FIPS.186-5では上位100bitに違いがなければ作り直すことになっている）、実装上の誤りがなければフェルマー法で素因数分解することはできないはずです。何年か前に実装上の問題で解けるRSA鍵が見つかったというニュースがありましたが、そのような実装の誤りはめったにないはずです。

ナイーブな実装

アルゴリズムをそのまま実装してみます。いわゆる「ナイーブな実装」（アルゴリズムや定義をそのまま、最適化や工夫なしで実装したコード）というやつです。

ナイーブな実装

use num_bigint::BigUint;
use num_integer::Integer;

pub fn factorize(n: &BigUint, round: usize) -> Option<(BigUint, BigUint)> {
    if n.is_even() {
        return Some((BigUint::from(2u8), n >> 1));
    }

    let mut a = n.sqrt();
    if &a * &a == *n {
        return Some((a.clone(), a));
    }

    for _ in 0..round {
        a += 1u8;
        let b2 = &a * &a - n;
        let b = b2.sqrt();
        if &b * &b == b2 {
            return Some((&a - &b, &a + &b));
        }
    }
    None
}

引数 n として偶数が渡された場合には None を返してもよいのかもしれませんが、とりあえず「 2 」と「 2 で割った値」を返すことにします。

また、念のため n が平方数の場合もチェックをしておきます。

ちなみに、round はどこまで頑張って探索するかの値です。

性能を測るため、同じぐらいの大きさの素数を掛け合わせた数字を渡して処理時間を測ってみます。

毎回ランダムで生成すると性能の評価がしにくいので、適当なプログラムを作って1024bitの乱数値と、それの下から指定のbitを反転させた値を作って、それらの値の次の素数を探して、２つの素数のペアを生成しました。
これらの事前に用意したペアを掛け合わせたものを引数 $n$ に渡して、処理時間を測ってみます。

性能測定用のソースはちょっと長いので折りたたんでおきます。

性能測定用ソース

性能測定用ソース

use num_bigint::BigUint;
use std::time::Instant;
use yoshi389111_fermat_factor::factorize;

fn main() {
    let prime_pairs = vec![
        (
            512,
            "98e4402433270625811dbfcc1aabe2c04a1a62c8a198e6c70f2d214f9b7c07f5965694745a2da3305f0d7bfa5b50631873bd26121e4d907011c4d8b206dc0b6b58d2a98773e00e895748caa94e1268c4b5a5e80c886f725925b92fa5facbc9943f31d4e632cb478fb930136a4cf9fa209be20e0d51afc7f14b6bac0861f9a1db",
            "98e4402433270625811dbfcc1aabe2c04a1a62c8a198e6c70f2d214f9b7c07f5965694745a2da3305f0d7bfa5b50631873bd26121e4d907011c4d8b206dc0b6bd8d2a98773e00e895748caa94e1268c4b5a5e80c886f725925b92fa5facbc9943f31d4e632cb478fb930136a4cf9fa209be20e0d51afc7f14b6bac0861f99475",
        ),
        (
            513,
            "8e759c64f4252646977c0148cb11d08600569187cc5abb81da46390e6214550287e0d268e63456304cf1a594c2caf274af6f3b37e64f610dd9ecba1a9055e9eedbfabe1e63ec29e60074e924c0246442adbb9605994092297c47a4d21a128fe4c8530d11c3a3c826b58a55a8d102222eb8427646dd3b11e49c4e9845031e3d4b",
            "8e759c64f4252646977c0148cb11d08600569187cc5abb81da46390e6214550287e0d268e63456304cf1a594c2caf274af6f3b37e64f610dd9ecba1a9055e9efdbfabe1e63ec29e60074e924c0246442adbb9605994092297c47a4d21a128fe4c8530d11c3a3c826b58a55a8d102222eb8427646dd3b11e49c4e9845031e401f",
        ),
        (
            514,
            "854c507107603d393712713d9f0ff8db7a6033da925430f3e8976e333ece2fbb916a15ce5df41a82b66abf060500c2467e93ec4414430269c87a4c7133c795dda631c0e07458a0bbb06b071c35466345ea9ade9f46826e4b9c05c7e81c4b4c2f3980da7d7ff8c584802da03df0609fa5be2557e18f7106cf7ac31c616e92b9e3",
            "854c507107603d393712713d9f0ff8db7a6033da925430f3e8976e333ece2fbb916a15ce5df41a82b66abf060500c2467e93ec4414430269c87a4c7133c795dfa631c0e07458a0bbb06b071c35466345ea9ade9f46826e4b9c05c7e81c4b4c2f3980da7d7ff8c584802da03df0609fa5be2557e18f7106cf7ac31c616e92b9f5",
        ),
        (
            515,
            "85eadb05968a80ac944af5c2062f4bc571516719a7a58a99583b06f4d615b062544c29fdad50b42c7463483e3ade1963d2e57adc98a1a2aa54cb1b97026b3dc369640a6506c0d2126e0147323328cc3948e8db46d98e37e30b01b4057f4e231424f1557ec5fb3bd7280762bf295c45fcb9382c523caf737f94a0419a3d39e86d",
            "85eadb05968a80ac944af5c2062f4bc571516719a7a58a99583b06f4d615b062544c29fdad50b42c7463483e3ade1963d2e57adc98a1a2aa54cb1b97026b3dc769640a6506c0d2126e0147323328cc3948e8db46d98e37e30b01b4057f4e231424f1557ec5fb3bd7280762bf295c45fcb9382c523caf737f94a0419a3d39df93",
        ),
        (
            516,
            "d18087cafc53f26732d5152664cc19690e87bd467e8c42164be37d2df5204e945287a1f37ec2fc0d676d844db2ad4b1737d7e77d4dba8701400b927f46a286323c06b565396a09318e1058b65e146f2477baaa3ab2de1a29b84d4636e1df501e6d0cbdbf49f5b22a89cffd665caca33536867ca786687df9a8b921e1e7a9f053",
            "d18087cafc53f26732d5152664cc19690e87bd467e8c42164be37d2df5204e945287a1f37ec2fc0d676d844db2ad4b1737d7e77d4dba8701400b927f46a2863a3c06b565396a09318e1058b65e146f2477baaa3ab2de1a29b84d4636e1df501e6d0cbdbf49f5b22a89cffd665caca33536867ca786687df9a8b921e1e7a9f3b7",
        ),
        (
            517,
            "a7bf1b79eef33164908faeee76c1870defde7e82d3fad366fcd522abac288501a83e5d76379862c0e34e134cd21c50dfaabb4f4bb871e6241662adcfcff767039ef4d128de0fb3ad93c3d7a6eb6a1a921635e9fc3a3965f8a629946bceb06581055dc26709304316209ebe536d34ca6cb617625f9e5ab3fc52fbf00b424ca0ff",
            "a7bf1b79eef33164908faeee76c1870defde7e82d3fad366fcd522abac288501a83e5d76379862c0e34e134cd21c50dfaabb4f4bb871e6241662adcfcff767139ef4d128de0fb3ad93c3d7a6eb6a1a921635e9fc3a3965f8a629946bceb06581055dc26709304316209ebe536d34ca6cb617625f9e5ab3fc52fbf00b424c9da7",
        ),
        (
            518,
            "f1917aa83aed0b8acf4c952c5735f10c2e61a7d20589627fec428d797916f9c84f350ed37ffdb932135ee39cd6bcc711e8cf7cc9a1627dc56c312b02c3b9194e43fc205daefebab72bcbc1571d08724e6f02e51b63f9287c35e034ff2e243d276b745bdf3cb9fba633cf70c274799b4b3182e667873a89452f8e3a565ac422ef",
            "f1917aa83aed0b8acf4c952c5735f10c2e61a7d20589627fec428d797916f9c84f350ed37ffdb932135ee39cd6bcc711e8cf7cc9a1627dc56c312b02c3b9196e43fc205daefebab72bcbc1571d08724e6f02e51b63f9287c35e034ff2e243d276b745bdf3cb9fba633cf70c274799b4b3182e667873a89452f8e3a565ac41bdb",
        ),
        (
            519,
            "b01a3d84d107319843fddece002c8df85cdeff1dcd64e55abce6e046ec63e8bd4536fd6f027dc04d78b50c2733cd6b6a642e77d9f5908799b5c0e7e215a0d8bd8e5de2b4ba27e2b3c5a54febd380cb2abc617ca0bc75abd5d576eec3fb0e684006f538af9c56a5e24d99f5cd5d98fea5bbfab6f8a0be7e9272c64e68ebc55113",
            "b01a3d84d107319843fddece002c8df85cdeff1dcd64e55abce6e046ec63e8bd4536fd6f027dc04d78b50c2733cd6b6a642e77d9f5908799b5c0e7e215a0d8fd8e5de2b4ba27e2b3c5a54febd380cb2abc617ca0bc75abd5d576eec3fb0e684006f538af9c56a5e24d99f5cd5d98fea5bbfab6f8a0be7e9272c64e68ebc55107",
        ),
        (
            520,
            "83735706cb0ba451bc0d72545658ea17f55c664f4f1c478ae78023f57a797d723c2f93eb59783739f2eb5770c52f48a24993021cc6378096824f62b1835c9701b97fe4ecb181226c6774a1811da21a733319af5bb587e35745e7e851657f45762e9771e862b9c39542a5583e199b949b5597ae2a9dd19afe7f69cbcc7e9ab193",
            "83735706cb0ba451bc0d72545658ea17f55c664f4f1c478ae78023f57a797d723c2f93eb59783739f2eb5770c52f48a24993021cc6378096824f62b1835c9781b97fe4ecb181226c6774a1811da21a733319af5bb587e35745e7e851657f45762e9771e862b9c39542a5583e199b949b5597ae2a9dd19afe7f69cbcc7e9ab195",
        ),
        (
            521,
            "fa87d133ac13b1812ab02fc2103a19c8fa808d631fb2277705b6795963265051a34eb7ac2581dda91f6afab779d1782831d81e881dd96299b1c0522caa28ced3eb1002887398e557933701bf79dde5dd9814ebd31e233a53a306f6fa823f5bf3a583ea32272e1ba06c5ba233eecd255ac486b9d025c332952ca152e0855c7b15",
            "fa87d133ac13b1812ab02fc2103a19c8fa808d631fb2277705b6795963265051a34eb7ac2581dda91f6afab779d1782831d81e881dd96299b1c0522caa28cfd3eb1002887398e557933701bf79dde5dd9814ebd31e233a53a306f6fa823f5bf3a583ea32272e1ba06c5ba233eecd255ac486b9d025c332952ca152e0855c810d",
        ),
        (
            522,
            "bf30899cd145cb8a6dd788b47b5daa8976bab3d29975bbbd9e3a713a3bde3b58d8c40302a3ee586863efc5970292102bf80d340680068fc574c989e24eb394f3f5c3bc7150bf6c744ab647f7d3321f7ddc251b217f17c703b848713d9b9fef7a43d10dad42a4cf17a0e2300ba0322681c8107659de43f5f376a7122206115f7f",
            "bf30899cd145cb8a6dd788b47b5daa8976bab3d29975bbbd9e3a713a3bde3b58d8c40302a3ee586863efc5970292102bf80d340680068fc574c989e24eb396f3f5c3bc7150bf6c744ab647f7d3321f7ddc251b217f17c703b848713d9b9fef7a43d10dad42a4cf17a0e2300ba0322681c8107659de43f5f376a71222061162c7",
        ),
        (
            523,
            "b0a05c0cd5b72e6b2033cc5ae6804031ee92e3daa95935b5aa1e4d70b603160f2a6d8200306fc4a91d7834b15cc66c2395ee53375df16ea85e0ce4c37e8451839d7bb49595adfbf319b4590534b151e10e7e35cce1fce6814070b795c89f94ae963d49b60d26faf37535cf544ca765b5263e0c606f87d729cfb8742ad375b98b",
            "b0a05c0cd5b72e6b2033cc5ae6804031ee92e3daa95935b5aa1e4d70b603160f2a6d8200306fc4a91d7834b15cc66c2395ee53375df16ea85e0ce4c37e8455839d7bb49595adfbf319b4590534b151e10e7e35cce1fce6814070b795c89f94ae963d49b60d26faf37535cf544ca765b5263e0c606f87d729cfb8742ad375b9b3",
        ),
        (
            524,
            "baa1e3b27e57e6bd91ca46725c0cfc060e00a20fbc45649ce6cc8357ef93c3ed472c7c36452e324c2c812ad005f7cee5b7847520df18a1fccdc94a06fd9334f725515761d738becaf4361333a019d93b8be2a6d1efbbfeaf6cfd4b34a19b974602df42f58ba0fba6c9a5d494807a9d8acbf47e740c7242c0cc8dc5447c0c4c85",
            "baa1e3b27e57e6bd91ca46725c0cfc060e00a20fbc45649ce6cc8357ef93c3ed472c7c36452e324c2c812ad005f7cee5b7847520df18a1fccdc94a06fd933cf725515761d738becaf4361333a019d93b8be2a6d1efbbfeaf6cfd4b34a19b974602df42f58ba0fba6c9a5d494807a9d8acbf47e740c7242c0cc8dc5447c0c512f",
        ),
        (
            525,
            "adf92b8952fd423ec67c96a6c5004e1953bd2157b26536d372434dfe299d751043cfc63375133f8e33daecfcaf1fa12939330fe38f6705edbe296d9ad5fe49b3c82eccb6df41cd57afbb61bf6de391483c3aaae1350d1bd10159327095f88671dd81d3ef6c18826b82d6978dede0fed1316e58a39536de87a7126e96f79eeb7d",
            "adf92b8952fd423ec67c96a6c5004e1953bd2157b26536d372434dfe299d751043cfc63375133f8e33daecfcaf1fa12939330fe38f6705edbe296d9ad5fe59b3c82eccb6df41cd57afbb61bf6de391483c3aaae1350d1bd10159327095f88671dd81d3ef6c18826b82d6978dede0fed1316e58a39536de87a7126e96f79ee9b5",
        ),
        (
            526,
            "d1fd64f105f09d0e715477aebdbdf7d8a547e11dc04ee63ae8c8a1bae6162bd3ba8e84005289dc1b21d46fccda1e581d56aea1aaa9e73f6d67de66e5883f4a4ebc996478d460da71ef039ec223d9747d82a41e157f5d13c0638c0ae53de4cb5285e089bc15c242c74e7721aec4889abbd469cbc497dc4895746f4764738930ed",
            "d1fd64f105f09d0e715477aebdbdf7d8a547e11dc04ee63ae8c8a1bae6162bd3ba8e84005289dc1b21d46fccda1e581d56aea1aaa9e73f6d67de66e5883f6a4ebc996478d460da71ef039ec223d9747d82a41e157f5d13c0638c0ae53de4cb5285e089bc15c242c74e7721aec4889abbd469cbc497dc4895746f476473893165",
        ),
        (
            527,
            "8dd7e43e89e9f5482d249737c0ef7770a2e8a2be2637a68f7b7fae366336b77c2f9c18666eeca835a96803de68510d731c149fd390c4b622b363db33e5aabcd0357e94eeb467cd8a27a388739ed1c04e3b8cbc5184a647eb5b3be70bc08c7942ced86ddc2059443c8e82682f6a0e40e4ec0ed25c99732d15cb2b47bba226acdb",
            "8dd7e43e89e9f5482d249737c0ef7770a2e8a2be2637a68f7b7fae366336b77c2f9c18666eeca835a96803de68510d731c149fd390c4b622b363db33e5aafcd0357e94eeb467cd8a27a388739ed1c04e3b8cbc5184a647eb5b3be70bc08c7942ced86ddc2059443c8e82682f6a0e40e4ec0ed25c99732d15cb2b47bba226a8b3",
        ),
        (
            528,
            "91f11e6ae3250cf19a75706073950acec6b5e2b1c16b49901dfab39251672fc74350cb9cbcf41645043d9071e2e8bd343004e8f491fc1ecf6dc984ace8a66281d339bb9b1f95c7b4fa2b1d664524176494cd2014d659246c1cc70a392764565cffaba0b77e6829b84d37e5d8c709097aa3ab9345ab44c4dc9c00219ec14909f3",
            "91f11e6ae3250cf19a75706073950acec6b5e2b1c16b49901dfab39251672fc74350cb9cbcf41645043d9071e2e8bd343004e8f491fc1ecf6dc984ace8a6e281d339bb9b1f95c7b4fa2b1d664524176494cd2014d659246c1cc70a392764565cffaba0b77e6829b84d37e5d8c709097aa3ab9345ab44c4dc9c00219ec149100d",
        ),
        (
            529,
            "cb5ba17b8fa5447279d871936674b6e6da7655752ae7c2f4057cae437bbfd577c903bb485f66895d430242862fe50ab9fae1ff370b1a78431734387ab27a4122d13c6bbf7bd997707f7d4720249026977ad8875b0b5cb260cf464037523d1705607addf814876d038d0e0bb9a7a0c6e3999f19c2b318bdf6e19df55f3e1d00af",
            "cb5ba17b8fa5447279d871936674b6e6da7655752ae7c2f4057cae437bbfd577c903bb485f66895d430242862fe50ab9fae1ff370b1a78431734387ab27b4122d13c6bbf7bd997707f7d4720249026977ad8875b0b5cb260cf464037523d1705607addf814876d038d0e0bb9a7a0c6e3999f19c2b318bdf6e19df55f3e1d035d",
        ),
    ];
    for (x, p, q) in prime_pairs.iter() {
        let p = BigUint::parse_bytes(p.as_bytes(), 16).unwrap();
        let q = BigUint::parse_bytes(q.as_bytes(), 16).unwrap();

        let n = &p * &q;
        let start = Instant::now();
        let result = factorize(&n, 200_000_000);
        let duration = start.elapsed();
        println!("Fermat factorization took {:?} for x = {}", duration, x);
        assert_eq!(result, Some((p, q)));
    }
}

今回の例では $p$ と $q$ が 1024bit なので、$N$ は 2048bit 程度の RSA 鍵レベルを想定しています。

また下から 512ビット目が違うものから、529ビット目が違うものまでを用意して実測してみました。

ナイーブな実装での実測結果は以下の通りです（時間がかかるので、525ビット目までで打ち切りました）

Fermat factorization took 187.839µs for x = 512
Fermat factorization took 135.823µs for x = 513
Fermat factorization took 165.65µs for x = 514
Fermat factorization took 447.306µs for x = 515
Fermat factorization took 1.07017ms for x = 516
Fermat factorization took 4.258076ms for x = 517
Fermat factorization took 17.636901ms for x = 518
Fermat factorization took 66.855215ms for x = 519
Fermat factorization took 283.505197ms for x = 520
Fermat factorization took 580.439379ms for x = 521
Fermat factorization took 2.895906006s for x = 522
Fermat factorization took 12.281638067s for x = 523
Fermat factorization took 47.286195122s for x = 524
Fermat factorization took 200.997451656s for x = 525
^C

下から521ビットが違う程度まではほぼ一瞬、下から525ビットが違う素因数だと3分以上かかりました。

通常、性能実測するには複数回計測するものだと思いますが、そこまでしっかりしたものではないので、この１ショットを基準値としてしまいます（実際には複数回実測していますが、あまり大きな違いはありませんでした）。

ちなみに実測しているのは、しょぼいノートパソコンです。またデバッグビルドで動かしています。

毎回二乗するのをやめる

ナイーブな実装をみて、時間のかかっていそうな個所は二乗と、平方根を計算するところではないかと思います。

まず最初は a の二乗をループ内で行うのをやめてみました。

$a_{n+1}^2 = (a_n + 1)^2 = a_n^2 + 2 a_n + 1$

なので、事前に計算した値 diff から 2 * a + 1 を引くようにすれば、ループ内での $a^2$ を削除できます。

aの二乗をやめる

use num_bigint::BigUint;
use num_integer::Integer;

pub fn factorize(n: &BigUint, round: usize) -> Option<(BigUint, BigUint)> {
    if n.is_even() {
        return Some((BigUint::from(2u8), n >> 1));
    }

    let mut a = n.sqrt();
    if &a * &a == *n {
        return Some((a.clone(), a));
    }

    a += 1u8;
    let mut diff = &a * &a - n;

    for _ in 0..round {
        let b = diff.sqrt();
        if &b * &b == diff {
            return Some((&a - &b, &a + &b));
        }
        diff += &a * 2u8 + 1u8; // a^2 を高速化
        a += 1u8;
    }
    None
}

性能を実測してみました。

Fermat factorization took 230.305µs for x = 512
Fermat factorization took 187.24µs for x = 513
Fermat factorization took 191.4µs for x = 514
Fermat factorization took 550.693µs for x = 515
Fermat factorization took 1.030859ms for x = 516
Fermat factorization took 3.94046ms for x = 517
Fermat factorization took 9.560979ms for x = 518
Fermat factorization took 48.838874ms for x = 519
Fermat factorization took 222.162604ms for x = 520
Fermat factorization took 418.151074ms for x = 521
Fermat factorization took 2.193511162s for x = 522
Fermat factorization took 9.709374758s for x = 523
Fermat factorization took 35.328254102s for x = 524
Fermat factorization took 154.255474859s for x = 525
^C

時間のかかっていた524ビットや525ビットなどでは少し早くなった気がします。
逆に515ビットまでは少し遅くなっているかもしれませんが、一瞬で終わっているのであまり気にしないこととします。

二乗と平方根をやめる（ダメ）

では、二乗を平方根をすべてやめたらよいかというとそうでもないようです。

試しに、二乗や平方根を使わないフェルマー法の実装をしてみました。

簡単に説明すると

s は 2 * a + 1 を保持していて、a に加算する差分の値
t は 2 * b + 1 を保持していて、b に加算する差分の値

となっています。

二乗と平方根をやめる

use num_bigint::BigUint;
use num_integer::Integer;

pub fn factorize(n: &BigUint, round: usize) -> Option<(BigUint, BigUint)> {
    if n.is_even() {
        return Some((BigUint::from(2u8), n >> 1));
    }

    let a = n.sqrt();
    if &a * &a == *n {
        return Some((a.clone(), a));
    }
    let a = a + 1u8;
    let mut s = &a * 2u8 + 1u8;
    let mut t = BigUint::from(3u8); // b = 1; t = 2b + 1;
    let mut m = &a * &a - 1u8; // m = a^2 - b^2

    for _ in 0..round {
        while &m > n {
            m -= &t;
            t += 2u8;
        }
        if &m == n {
            let p = (&s + &t - 2u8) / 2u8;
            let q = (&s - &t) / 2u8;
            return Some((p, q));
        }
        m += &s;
        s += 2u8;
    }
    None
}

しかしこれの性能実測すると、１つも結果が返ってきません。

少し考えてみればわかりますが、下から512bit目が違うということは、1ラウンド目ですら b の計算をするのに $2^{512}$ ぐらいの値を足し算で求める必要があるということになります。

RSAのような巨大な素数を扱っていると $2^{512}$ が小さく感じられるかもしれませんが、$2^{500} \approx 10^{150}$ ぐらいなので、かなり巨大な数です（1無量大数が $10^{68}$ から $10^{88}$ ぐらい。1グーゴルが $10^{100}$）。

つまり二乗や平方根は必要で、最小限必要な回数だけ実行するようにする必要があると考えます。

（ちなみに、1ラウンド目で $a_1$ が $\lceil \sqrt {N} \rceil$ とすると、$a_1^2 - N$ は $1$ から $2a$ の範囲になります。そうすると対応する $b_1$ はおおよそ $\sqrt{a}$ 、$p$ と $q$ の差は $2b$ なので、$p$ 、 $q$ のビット長の約半分（$N$ のビット長の約1/4）＋1～2ビットぐらいまでは1ラウンドで解ける可能性があることになります。つまり $N$ が 2048ビットであれば、下から513～514ビットが違うぐらいだと1ラウンドで解けるかもしれません）

4を法とする平方剰余をチェックする

b を計算する（ diff が平方数であることを確認する）のに、平方根や二乗をするのが必須だとして、性能を向上するのにその回数をできるだけ少なくしたいです。

そこで思いついたのが、平方数は $\bmod$ をとると、あるパターンが出てくる、というものです。いわゆる平方剰余（quadratic residue）というやつですね。

具体的には平方数を $4$ で割った余りは、$0$ か $1$ しかあえりません。なので diff を 4 で割った余りを計算して、0 か 1 の時だけ、平方数であるか確認すればよいはずです。

4を法とする平方剰余をチェックする

use num_bigint::BigUint;
use num_integer::Integer;
use num_traits::ToPrimitive;

pub fn factorize(n: &BigUint, round: usize) -> Option<(BigUint, BigUint)> {
    if n.is_even() {
        return Some((BigUint::from(2u8), n >> 1));
    }

    let mut a = n.sqrt();
    if &a * &a == *n {
        return Some((a.clone(), a));
    }

    a += 1u8;
    let mut diff = &a * &a - n;
    for _ in 0..round {
        // 法を4とした平方剰余をチェック
        let m = (&diff % 4u8).to_u64().unwrap();
        if m == 0 || m == 1 {
            let b = diff.sqrt();
            if &b * &b == diff {
                return Some((&a - &b, &a + &b));
            }
        }
        diff += &a * 2u8 + 1u8;
        a += 1u8;
    }
    None
}

実測してみます。

Fermat factorization took 366.959µs for x = 512
Fermat factorization took 190.483µs for x = 513
Fermat factorization took 317.786µs for x = 514
Fermat factorization took 380.121µs for x = 515
Fermat factorization took 867.015µs for x = 516
Fermat factorization took 1.756878ms for x = 517
Fermat factorization took 4.221878ms for x = 518
Fermat factorization took 25.084245ms for x = 519
Fermat factorization took 126.906119ms for x = 520
Fermat factorization took 235.793942ms for x = 521
Fermat factorization took 1.21905644s for x = 522
Fermat factorization took 5.746704873s for x = 523
Fermat factorization took 22.353793219s for x = 524
Fermat factorization took 89.659819985s for x = 525
^C

525ビットで約1.5分まで短くなりました（ナイーブな実装では3分以上かかっていた）。

256を法とする平方剰余をチェックする

先ほどは4を法とした平方剰余をチェックしましたが、別に4である必要はないです。
もっと大きい値を法としてチェックしたほうが、不要な演算が削除できるかもしれません。

そこで、256を法としてチェックしてみることにしました。256の平方剰余の集合をハードコーディングするのは面倒なので、平方剰余かどうかをチェックする機能を作ってみます。

quadratic_residue.rs

use std::sync::Arc;

/// Represents a set of quadratic residues modulo a given modulus.
pub(crate) struct QuadraticResidueSet {
    /// The modulus for the quadratic residues.
    modulus: usize,
    /// A boolean vector indicating which numbers are quadratic residues.
    residues: Vec<bool>,
}

impl QuadraticResidueSet {
    /// Creates a new `QuadraticResidueSet` for a given modulus.
    pub(crate) fn new(modulus: usize) -> Arc<Self> {
        let mut residues = vec![false; modulus];
        for i in 0..modulus {
            residues[(i * i) % modulus] = true;
        }
        Arc::new(Self { modulus, residues })
    }

    /// Returns the modulus of the quadratic residue set.
    #[allow(dead_code)]
    pub(crate) fn modulus(self: &Arc<Self>) -> usize {
        self.modulus
    }

    /// Checks if a given number is a quadratic residue modulo the set's modulus.
    pub(crate) fn is_quadratic_residue(self: &Arc<Self>, a: usize) -> bool {
        self.residues[a % self.modulus]
    }
}

まずは平方剰余の集合を計算するロジックです。
平方剰余の求め方は、実際に二乗して $\bmod$ をとればよいので簡単です。
前掲のwikipediaでは、平方剰余は全パターンを計算しなくても、半分まで見ればよいと書いてありましたが、あまり気にせず全部をみています。

あわせて、平方剰余かどうかをチェックするロジックも用意しました。これは毎回巨大な数字 diff の剰余を求めるのではなく、事前に指定の値を法とした a および diff を求めておけば、平方剰余かどうかの判断もその小さい値で計算できると考えたためです。

quadratic_residue_checker.rs

use crate::quadratic_residue::QuadraticResidueSet;
use num_bigint::BigUint;
use num_traits::ToPrimitive;
use std::sync::Arc;

/// This is a checker that determines whether the `diff` in Fermat's method could be a perfect square, and it also has a function to advance the state.
pub(crate) struct QuadraticResidueChecker {
    /// The set of quadratic residues.
    residues: Arc<QuadraticResidueSet>,
    /// Current value of 'a' modulo the modulus.
    a: usize,
    /// Current value of 'diff' modulo the modulus.
    diff: usize,
}

impl QuadraticResidueChecker {
    /// Creates a new `QuadraticResidueChecker` with given initial values.
    pub(crate) fn new(residues: &Arc<QuadraticResidueSet>, a: &BigUint, diff: &BigUint) -> Self {
        let residues = Arc::clone(residues);
        let modulus = residues.modulus();
        let a = (a % modulus).to_usize().unwrap();
        let diff = (diff % modulus).to_usize().unwrap();
        Self { residues, a, diff }
    }

    /// Advances the state of the checker by a given number of steps.
    pub(crate) fn advance(&mut self, value: usize) {
        let modulus = self.residues.modulus();
        let value = value % modulus;
        self.diff = (self.diff + self.a * 2 * value + value * value) % modulus;
        self.a = (self.a + value) % modulus;
    }

    /// Checks if the current status of `diff` is a possible square number.
    pub(crate) fn is_possible(&self) -> bool {
        self.residues.is_quadratic_residue(self.diff)
    }

    /// Returns the current value of 'a' modulo the modulus.
    #[allow(dead_code)]
    pub(crate) fn get_a(&self) -> usize {
        self.a
    }

    /// Returns the current value of 'diff' modulo the modulus.
    #[allow(dead_code)]
    pub(crate) fn get_diff(&self) -> usize {
        self.diff
    }
}

pub(crate) trait CheckerCreater {
    fn checker(&self, initial_a: &BigUint, initial_diff: &BigUint) -> QuadraticResidueChecker;
}

impl CheckerCreater for Arc<QuadraticResidueSet> {
    fn checker(&self, initial_a: &BigUint, initial_diff: &BigUint) -> QuadraticResidueChecker {
        QuadraticResidueChecker::new(self, initial_a, initial_diff)
    }
}

では、実際に256を法としてチェックしてみます。

256を法として平方剰余をチェックする

mod quadratic_residue;
mod quadratic_residue_checker;

use crate::quadratic_residue::QuadraticResidueSet;
use crate::quadratic_residue_checker::CheckerCreater;
use num_bigint::BigUint;
use num_integer::Integer;
use once_cell::sync::Lazy;
use std::sync::Arc;

static QUADRATIC_RESIDUE_SET: Lazy<Arc<QuadraticResidueSet>> =
    Lazy::new(|| QuadraticResidueSet::new(256));

pub fn factorize(n: &BigUint, round: usize) -> Option<(BigUint, BigUint)> {
    if n.is_even() {
        return Some((BigUint::from(2u8), n >> 1));
    }

    let mut a = n.sqrt();
    if &a * &a == *n {
        return Some((a.clone(), a));
    }

    a += 1u8;
    let mut diff = &a * &a - n;
    // 法256を用いた平方剰余をチェック
    let mut checker = QUADRATIC_RESIDUE_SET.checker(&a, &diff);
    for _ in 0..round {
        if checker.is_possible() {
            let b = diff.sqrt();
            if &b * &b == diff {
                return Some((&a - &b, &a + &b));
            }
        }
        checker.advance(1);
        diff += &a * 2u8 + 1u8;
        a += 1u8;
    }
    None
}

早速実測してみます。

Fermat factorization took 146.191µs for x = 512
Fermat factorization took 102.165µs for x = 513
Fermat factorization took 117.752µs for x = 514
Fermat factorization took 129.549µs for x = 515
Fermat factorization took 136.817µs for x = 516
Fermat factorization took 414.136µs for x = 517
Fermat factorization took 1.844548ms for x = 518
Fermat factorization took 8.15412ms for x = 519
Fermat factorization took 69.080568ms for x = 520
Fermat factorization took 49.682107ms for x = 521
Fermat factorization took 265.722728ms for x = 522
Fermat factorization took 947.872525ms for x = 523
Fermat factorization took 10.107649321s for x = 524
Fermat factorization took 14.865193126s for x = 525
Fermat factorization took 56.773441601s for x = 526
^C

平方剰余の集合は初回のみ計算しているので、おそらく512ビットの時にその分のオーバーヘッドが乗ってしまっていますが、さほど時間がかからないので無視してしまいます。

525ビットで15秒ぐらいなので、明らかに性能向上が出来ていると思います（ナイーブな実装に比べて10倍ちょっと速くなった）。

256を法として、次の平方剰余までのステップを計算

さらなる性能向上を考えてみます。

前の処理では、a を毎回 +1 していますが、よく考えてみれば a と diff の剰余がわかっていれば、次の平方剰余まで a を何回足せばよいかは事前に計算できるように思います。計算してみます。

quadratic_residue_stepper.rs

use crate::quadratic_residue::QuadraticResidueSet;
use crate::quadratic_residue_checker::{CheckerCreater, QuadraticResidueChecker};
use num_bigint::BigUint;
use std::sync::Arc;

/// Precomputed step table for quadratic residues.
pub(crate) struct QuadraticResidueStepTable {
    /// The modulus for the quadratic residues.
    modulus: usize,
    /// A table mapping (a, diff) pairs to their step values.
    table: Vec<usize>,
}

impl QuadraticResidueStepTable {
    /// Creates a new `QuadraticResidueStepTable` for a given set of quadratic residues.
    pub(crate) fn new(residues: &Arc<QuadraticResidueSet>) -> Self {
        let modulus = residues.modulus();
        let mut table = vec![0usize; modulus * modulus];
        for a in 0..modulus {
            for diff in 0..modulus {
                table[a * modulus + diff] = Self::calc_step(residues, a, diff);
            }
        }
        Self { modulus, table }
    }

    /// Calculates the step value for a given (a, diff) pair.
    fn calc_step(residues: &Arc<QuadraticResidueSet>, a: usize, diff: usize) -> usize {
        let modulus = residues.modulus();
        if modulus % 2 == 0 && a % 2 == diff % 2 {
            return 0;
        }
        let a = BigUint::from(a);
        let diff = BigUint::from(diff);
        let mut checker = residues.checker(&a, &diff);
        for x in 1..(modulus * modulus) {
            checker.advance(1);
            if checker.is_possible() {
                return x;
            }
        }
        panic!("Not found");
    }

    /// Returns the step value for a given (a, diff) pair.
    pub(crate) fn get_step(&self, a: usize, diff: usize) -> usize {
        self.table[(a % self.modulus) * self.modulus + (diff % self.modulus)]
    }
}

ここでちょっと気を付けないといけないのは、単にすべてのパターンでテーブルを作ろうと思うと、無限ループになってしまうということです。
具体的には calc_step に渡される a と diff がどちらも偶数の場合に、いくつ足しても平方剰余にならない、ということがあり得ます。これは、最初の説明であった $p$ と $q$ がともに奇数であればあり得ないパターンです。そこで、a と diff がともに偶数の場合にはステップを 0 としています。
（もし $p$ と $q$ がともに偶数の場合なども解けるようにしたいなら、もう少し条件を考えたほうが良いかもしれませんが、ここでは $p$ と $q$ はともに奇数の前提で進めます）

法256をもとに、次の平方剰余までのステップを計算

mod quadratic_residue;
mod quadratic_residue_checker;
mod quadratic_residue_stepper;

use crate::quadratic_residue::QuadraticResidueSet;
use crate::quadratic_residue_checker::CheckerCreater;
use crate::quadratic_residue_stepper::QuadraticResidueStepTable;
use num_bigint::BigUint;
use num_integer::Integer;
use once_cell::sync::Lazy;
use std::sync::Arc;

static QUADRATIC_RESIDUE_SET: Lazy<Arc<QuadraticResidueSet>> =
    Lazy::new(|| QuadraticResidueSet::new(256));

static QUADRATIC_RESIDUE_STEP_TABLE: Lazy<QuadraticResidueStepTable> =
    Lazy::new(|| QuadraticResidueStepTable::new(&QUADRATIC_RESIDUE_SET));

pub fn factorize(n: &BigUint, round: usize) -> Option<(BigUint, BigUint)> {
    if n.is_even() {
        return Some((BigUint::from(2u8), n >> 1));
    }

    let mut a = n.sqrt();
    if &a * &a == *n {
        return Some((a.clone(), a));
    }

    a += 1u8;
    let mut diff = &a * &a - n;
    let mut checker = QUADRATIC_RESIDUE_SET.checker(&a, &diff);
    for _ in 0..round {
        let b = diff.sqrt();
        if &b * &b == diff {
            return Some((&a - &b, &a + &b));
        }
        // 法256をもとに、次の平方剰余までのステップを計算
        let step = QUADRATIC_RESIDUE_STEP_TABLE.get_step(checker.get_a(), checker.get_diff());
        checker.advance(step);
        diff += &a * 2u8 * step + step * step;
        a += step;
    }
    None
}

性能実測してみます。

Fermat factorization took 165.872µs for x = 512
Fermat factorization took 101.361µs for x = 513
Fermat factorization took 116.353µs for x = 514
Fermat factorization took 71.919885ms for x = 515
Fermat factorization took 166.084µs for x = 516
Fermat factorization took 383.587µs for x = 517
Fermat factorization took 1.015508ms for x = 518
Fermat factorization took 9.508589ms for x = 519
Fermat factorization took 71.426702ms for x = 520
Fermat factorization took 31.871831ms for x = 521
Fermat factorization took 154.296398ms for x = 522
Fermat factorization took 637.631111ms for x = 523
Fermat factorization took 9.696094606s for x = 524
Fermat factorization took 9.927856895s for x = 525
Fermat factorization took 34.134617547s for x = 526
Fermat factorization took 200.221531595s for x = 527
^C

525ビットで10秒ぐらいなので、少し速くなった気がします（ナイーブな実装に比べて20倍ぐらい）。

複数の法で平方剰余をチェックする

現在は 256 を法としてチェックしていますが、これって16進数でいえば、一番下の1バイトだけで平方数かどうかをチェックしているようなものだと思います。
もっと全体的にチェックしたほうが、効率が（つまり性能が）上がるんじゃないかと考えました。

なので、256以外の値でも剰余を求めてみることにします。

次の平方剰余までのステップは今まで通り 256 を法として求めておき、それとは別に256とは互いに素な値でも平方剰余チェックを行うことにしてみます。

256は $2^8$ なので、2以外の小さな素数の積で同じぐらいの数字を適当に選んでみました。

255 = 3 * 5 * 17
253 = 11 * 23
247 = 13 * 19
203 = 7 * 29

これで、素数の小さい方から10個はチェックしていることになると思います。

複数の法で平方剰余をチェックする

mod quadratic_residue;
mod quadratic_residue_checker;
mod quadratic_residue_stepper;

use crate::quadratic_residue::QuadraticResidueSet;
use crate::quadratic_residue_checker::CheckerCreater;
use crate::quadratic_residue_stepper::QuadraticResidueStepTable;
use num_bigint::BigUint;
use num_integer::Integer;
use once_cell::sync::Lazy;
use std::sync::Arc;

static QUADRATIC_RESIDUE_SET: Lazy<Arc<QuadraticResidueSet>> =
    Lazy::new(|| QuadraticResidueSet::new(256));

static QUADRATIC_RESIDUE_STEP_TABLE: Lazy<QuadraticResidueStepTable> =
    Lazy::new(|| QuadraticResidueStepTable::new(&QUADRATIC_RESIDUE_SET));

static AUX_QUADRATIC_RESIDUE_SETS: Lazy<Vec<Arc<QuadraticResidueSet>>> = Lazy::new(|| {
    vec![
        QuadraticResidueSet::new(255),
        QuadraticResidueSet::new(253),
        QuadraticResidueSet::new(247),
        QuadraticResidueSet::new(203),
    ]
});

pub fn factorize(n: &BigUint, round: usize) -> Option<(BigUint, BigUint)> {
    if n.is_even() {
        return Some((BigUint::from(2u8), n >> 1));
    }

    let mut a = n.sqrt();
    if &a * &a == *n {
        return Some((a.clone(), a));
    }

    a += 1u8;
    let mut diff = &a * &a - n;
    let mut checker = QUADRATIC_RESIDUE_SET.checker(&a, &diff);
    let mut aux_checkers: Vec<_> = AUX_QUADRATIC_RESIDUE_SETS
        .iter()
        .map(|set| set.checker(&a, &diff))
        .collect();
    for _ in 0..round {
        if aux_checkers.iter().all(|c| c.is_possible()) {
            let b = diff.sqrt();
            if &b * &b == diff {
                return Some((&a - &b, &a + &b));
            }
        }
        let step = QUADRATIC_RESIDUE_STEP_TABLE.get_step(checker.get_a(), checker.get_diff());
        checker.advance(step);
        aux_checkers.iter_mut().for_each(|c| c.advance(step));
        diff += &a * 2u8 * step + step * step;
        a += step;
    }
    None
}

Fermat factorization took 362.551µs for x = 512
Fermat factorization took 217.108µs for x = 513
Fermat factorization took 235.157µs for x = 514
Fermat factorization took 90.979835ms for x = 515
Fermat factorization took 197.477µs for x = 516
Fermat factorization took 184.494µs for x = 517
Fermat factorization took 541.98µs for x = 518
Fermat factorization took 486.146µs for x = 519
Fermat factorization took 2.919195ms for x = 520
Fermat factorization took 1.603563ms for x = 521
Fermat factorization took 8.861655ms for x = 522
Fermat factorization took 34.22904ms for x = 523
Fermat factorization took 495.501513ms for x = 524
Fermat factorization took 463.836189ms for x = 525
Fermat factorization took 1.481605613s for x = 526
Fermat factorization took 7.966876889s for x = 527
Fermat factorization took 139.260361398s for x = 528
^C

525ビットが0.5秒ぐらい。思ったよりも速くなりました。ナイーブな実装に比べて400倍ぐらい？
（524ビットの方が、処理時間がかかっているのは実行時間の揺れですかね。処理時間も短くなってきたので、揺れ幅が相対的に目立つようになったのかな。今までもあまり違わないぐらいの時間だったので、テスト用データが偏っているのかもしれません）

チェックする法を厳選する

いままでチェックしている法は、いわば適当に選んだだけの値ですが、もっと数学的により良い法( $\bmod$ )があるかもしれません。確認してみました。

3～256 の範囲で一番絞り込める、つまり平方剰余の集合の数が少ないものを調べてみました。

この範囲だと 240 が一番良いようです（1024までの範囲で調べると 1008 の方がより絞り込めそうですが）。

というわけで、今まで次の平方剰余までのステップを計算するのに 256 を法として使っていましたが、240に切り替えてみます。

他の併用する法も、互いに素になる値でより絞り込めるものを選定しました。

チェックする法を厳選する

mod quadratic_residue;
mod quadratic_residue_checker;
mod quadratic_residue_stepper;

use crate::quadratic_residue::QuadraticResidueSet;
use crate::quadratic_residue_checker::CheckerCreater;
use crate::quadratic_residue_stepper::QuadraticResidueStepTable;
use num_bigint::BigUint;
use num_integer::Integer;
use once_cell::sync::Lazy;
use std::sync::Arc;

static QUADRATIC_RESIDUE_SET: Lazy<Arc<QuadraticResidueSet>> =
    Lazy::new(|| QuadraticResidueSet::new(240));

static QUADRATIC_RESIDUE_STEP_TABLE: Lazy<QuadraticResidueStepTable> =
    Lazy::new(|| QuadraticResidueStepTable::new(&QUADRATIC_RESIDUE_SET));

static AUX_QUADRATIC_RESIDUE_SETS: Lazy<Vec<Arc<QuadraticResidueSet>>> = Lazy::new(|| {
    vec![
        QuadraticResidueSet::new(247),
        QuadraticResidueSet::new(253),
        QuadraticResidueSet::new(217),
        QuadraticResidueSet::new(251),
        QuadraticResidueSet::new(241),
    ]
});

pub fn factorize(n: &BigUint, round: usize) -> Option<(BigUint, BigUint)> {
    if n.is_even() {
        return Some((BigUint::from(2u8), n >> 1));
    }

    let mut a = n.sqrt();
    if &a * &a == *n {
        return Some((a.clone(), a));
    }

    a += 1u8;
    let mut diff = &a * &a - n;
    let mut checker = QUADRATIC_RESIDUE_SET.checker(&a, &diff);
    let mut aux_checkers: Vec<_> = AUX_QUADRATIC_RESIDUE_SETS
        .iter()
        .map(|set| set.checker(&a, &diff))
        .collect();
    for _ in 0..round {
        if aux_checkers.iter().all(|c| c.is_possible()) {
            let b = diff.sqrt();
            if &b * &b == diff {
                return Some((&a - &b, &a + &b));
            }
        }
        let step = QUADRATIC_RESIDUE_STEP_TABLE.get_step(checker.get_a(), checker.get_diff());
        checker.advance(step);
        aux_checkers.iter_mut().for_each(|c| c.advance(step));
        diff += &a * 2u8 * step + step * step;
        a += step;
    }
    None
}

実測値です。

Fermat factorization took 322.73µs for x = 512
Fermat factorization took 130.888µs for x = 513
Fermat factorization took 147.79µs for x = 514
Fermat factorization took 99.445238ms for x = 515
Fermat factorization took 321.145µs for x = 516
Fermat factorization took 341.322µs for x = 517
Fermat factorization took 440.94µs for x = 518
Fermat factorization took 385.902µs for x = 519
Fermat factorization took 661.727µs for x = 520
Fermat factorization took 1.086505ms for x = 521
Fermat factorization took 6.927326ms for x = 522
Fermat factorization took 30.080282ms for x = 523
Fermat factorization took 109.303327ms for x = 524
Fermat factorization took 412.114258ms for x = 525
Fermat factorization took 874.033006ms for x = 526
Fermat factorization took 7.225859417s for x = 527
Fermat factorization took 27.797995773s for x = 528
Fermat factorization took 106.865781055s for x = 529

事前に用意した下から529ビット目が違う値まで処理できました。
525ビットで400ミリ秒ちょっと。ナイーブな実装に比べて500倍弱ぐらいですかね。

BigUintの加算を最小限にする

平方剰余のチェックをする法をさらに厳選すればもう少し速くなる気がしますが、いったんそちらは終了して別の方法を考えます。

前掲のロジックでは、平方剰余チェッカーの加算( checker や aux_checkers の advance() )と、diff の加算を同期して実行していますが、 BigUint の加算を毎回するのは遅いような気がしてきました。

そこで、平方剰余チェッカーがすべて平方剰余である場合のみ、diff をまとめて計算するようにしてみます。

まずは、前掲の quadratic_residue_stepper.rs に機能を追加します。

quadratic_residue_stepper.rsの差分

// 前掲のロジック(QuadraticResidueStepTable)は省略

/// Stepper to efficiently advance through quadratic residues.
#[allow(dead_code)]
pub(crate) struct QuadraticResidueStepper {
    residues: Arc<QuadraticResidueSet>,
    aux_residues: Vec<Arc<QuadraticResidueSet>>,
    step_table: QuadraticResidueStepTable,
}

#[allow(dead_code)]
impl QuadraticResidueStepper {
    /// Creates a new `QuadraticResidueStepper` for a given set of quadratic residues.
    pub(crate) fn new(modulus: usize, aux_moduluses: &[usize]) -> Arc<Self> {
        let residues = QuadraticResidueSet::new(modulus);
        let aux_residues = aux_moduluses
            .iter()
            .map(|&m| QuadraticResidueSet::new(m))
            .collect();
        let step_table = QuadraticResidueStepTable::new(&residues);
        Arc::new(Self {
            residues,
            aux_residues,
            step_table,
        })
    }

    /// Creates an iterator to step through quadratic residues starting from given `a` and `diff`.
    pub(crate) fn iter(
        self: &Arc<Self>,
        a: &BigUint,
        diff: &BigUint,
    ) -> impl Iterator<Item = usize> + use<> {
        let checker = self.residues.checker(a, diff);
        let aux_checkers = self
            .aux_residues
            .iter()
            .map(|r| r.checker(a, diff))
            .collect();
        QuadraticResidueStepIter {
            stepper: Arc::clone(self),
            checker,
            aux_checkers,
        }
    }
}

/// Iterator to step through quadratic residues.
pub(crate) struct QuadraticResidueStepIter {
    /// The stepper used for stepping.
    stepper: Arc<QuadraticResidueStepper>,
    /// The main checker for quadratic residues.
    checker: QuadraticResidueChecker,
    /// Auxiliary checkers for additional residue classes.
    aux_checkers: Vec<QuadraticResidueChecker>,
}

impl Iterator for QuadraticResidueStepIter {
    type Item = usize;

    fn next(&mut self) -> Option<Self::Item> {
        let mut total_step = 0;

        loop {
            let step = self
                .stepper
                .step_table
                .get_step(self.checker.get_a(), self.checker.get_diff());
            if step == 0 {
                panic!("Step is zero");
            }

            self.checker.advance(step);
            self.aux_checkers
                .iter_mut()
                .for_each(|aux_checker| aux_checker.advance(step));
            total_step += step;

            if self
                .aux_checkers
                .iter()
                .all(QuadraticResidueChecker::is_possible)
            {
                break;
            }
        }
        Some(total_step)
    }
}

mod quadratic_residue;
mod quadratic_residue_checker;
mod quadratic_residue_stepper;

use crate::quadratic_residue_stepper::QuadraticResidueStepper;
use num_bigint::BigUint;
use num_integer::Integer;
use once_cell::sync::Lazy;
use std::sync::Arc;

static STEPPER: Lazy<Arc<QuadraticResidueStepper>> =
    Lazy::new(|| QuadraticResidueStepper::new(240, &[247, 253, 217, 251, 241]));

pub fn factorize(n: &BigUint, round: usize) -> Option<(BigUint, BigUint)> {
    if n.is_even() {
        return Some((BigUint::from(2u8), n >> 1));
    }

    let mut a = n.sqrt();
    if &a * &a == *n {
        return Some((a.clone(), a));
    }

    a += 1u8;
    let mut diff = &a * &a - n;

    let mut step_iter = STEPPER.iter(&a, &diff);
    for _ in 0..round {
        let b = diff.sqrt();
        if &b * &b == diff {
            return Some((&a - &b, &a + &b));
        }
        let step = step_iter.next()? as u64;
        diff += &a * 2u8 * step + step * step;
        a += step;
    }
    None
}

Fermat factorization took 85.130328ms for x = 512
Fermat factorization took 123.331µs for x = 513
Fermat factorization took 142.354µs for x = 514
Fermat factorization took 191.073µs for x = 515
Fermat factorization took 216.225µs for x = 516
Fermat factorization took 460.879µs for x = 517
Fermat factorization took 465.528µs for x = 518
Fermat factorization took 297.494µs for x = 519
Fermat factorization took 540.14µs for x = 520
Fermat factorization took 549.44µs for x = 521
Fermat factorization took 1.923304ms for x = 522
Fermat factorization took 7.05437ms for x = 523
Fermat factorization took 35.178108ms for x = 524
Fermat factorization took 129.491239ms for x = 525
Fermat factorization took 198.968605ms for x = 526
Fermat factorization took 1.83574098s for x = 527
Fermat factorization took 6.710756374s for x = 528
Fermat factorization took 28.894852058s for x = 529

525ビットで130ミリ秒ぐらい。ナイーブな実装に比べて約1500倍ぐらいでしょうか。

相変わらず、ワンショットしかデータを載せていないので1500倍になったと堂々とは言いにくいですが、タイトルの100倍以上になったとは言っても大丈夫ではないかと思います。

リリースビルドで実行する

前述のようにここまでの実行はデバッグビルドで動かしていました。
念のため、ナイーブな実装と最終的な実装をリリースビルドでも実行してみました。

まずは、ナイーブな実装のリリースビルド版です。

Fermat factorization took 17.346µs for x = 512
Fermat factorization took 10.154µs for x = 513
Fermat factorization took 14.596µs for x = 514
Fermat factorization took 26.23µs for x = 515
Fermat factorization took 109.15µs for x = 516
Fermat factorization took 262.402µs for x = 517
Fermat factorization took 764.89µs for x = 518
Fermat factorization took 5.596754ms for x = 519
Fermat factorization took 29.992451ms for x = 520
Fermat factorization took 56.150718ms for x = 521
Fermat factorization took 263.951167ms for x = 522
Fermat factorization took 938.070446ms for x = 523
Fermat factorization took 3.446330091s for x = 524
Fermat factorization took 14.45814434s for x = 525
Fermat factorization took 45.68594813s for x = 526
Fermat factorization took 304.282836035s for x = 527
^C

次に最終的な実装のリリースビルド版です。

Fermat factorization took 10.398677ms for x = 512
Fermat factorization took 20.292µs for x = 513
Fermat factorization took 32.023µs for x = 514
Fermat factorization took 26.951µs for x = 515
Fermat factorization took 37.625µs for x = 516
Fermat factorization took 37.731µs for x = 517
Fermat factorization took 34.031µs for x = 518
Fermat factorization took 45.234µs for x = 519
Fermat factorization took 86.981µs for x = 520
Fermat factorization took 104.631µs for x = 521
Fermat factorization took 327.948µs for x = 522
Fermat factorization took 833.978µs for x = 523
Fermat factorization took 4.13047ms for x = 524
Fermat factorization took 22.880555ms for x = 525
Fermat factorization took 45.999753ms for x = 526
Fermat factorization took 307.211774ms for x = 527
Fermat factorization took 1.048897491s for x = 528
Fermat factorization took 3.985943878s for x = 529

525ビットベースで、約600倍ぐらいですかね。デバッグビルドよりも差が縮まったけど、それでも100倍以上になったといってもよさそう？
初回(512ビット)が遅いのは、いろいろなオーバーヘッドが乗ってしまっているのでしょうね。