🐶

RustでMD5を実装（ついでにsin64を計算）

2021/09/09に公開

Rust

math

tech

概要

MD5ハッシュアルゴリズムをRFC1321に記載された仕様に従ってRustで実装します。付録として、その中で利用される $4294967296|\sin i|$ の整数部分という値の計算方法も示します。

背景

しばらく前から主に公式ドキュメントでRustを学び始めました。Zennの以下の記事にもお世話になっています。

学ぶうちにちょっとしたサイズのアルゴリズムを実装してみたくなり、今春放送されたゴジラSP^[1]で重要な役割を果たしたMD5ハッシュの実装を試してみることにしました。

MD5とは

MD5は任意長のビット列を128ビットという固定長のハッシュ値に変換するアルゴリズムです。したがって理想としては任意長のビット列を引数に取るべきですが、本記事では文字列に対するMD5ハッシュ関数をRustで実装します。

MD5 - Wikipedia

考案者はRon Rivest^[2]。当時既に脆弱性の発見されていたMD4の後継として1991年に設計され、翌1992年にアルゴリズムの仕様がRFC1321として公開されました。この文書の付録にはC言語による実装が例示されています。

RFC1321（正誤表付き）

主にこの仕様書に沿って実装を進めます。手本として以下のクレートや実装例も挙げます。実際に利用する場合にはよく整備されたこれらのクレートのいずれかをインポートするべきでしょう。

MD5ハッシュ化はmd5sumコマンドでも手軽に行えます。

本記事で目標とするところの文字列に対するMD5ハッシュはたとえば、

$ printf '%s' "Rust" | md5sum
f5e265d607cb720058fc166e00083fe8 *-

のようにコマンドラインから試すことができます^[3]。

この結果は文字列Rustに対するMD5ハッシュがf5e265d607cb720058fc166e00083fe8であることを示します。ハッシュ値は16進で32桁、すなわち128ビットであることが確かめられます。

バイト順の問題

MD5のアルゴリズム自体は単純で、ビット演算さえ分かれば特段難しいところはありません。唯一引っ掛かるとすればバイト順（エンディアン）の問題です。

RFC1321の文中、"low-order bytes first"かそれと同等な表現が何度か登場します。これはハッシュ化の対象になる元のビット列に対して、リトルエンディアンの順で語を読み取ることを意味します。

リトルエンディアンが何であるか（というかビッグとリトルのどっちがどっちか）は実例を見るのが手っ取り早いでしょう。Rustではbyteorderクレートをdependenciesに加えることで利用できます。

byteorder - crates.io: Rust Package Registry

例として、バイト列を32ビットの語として読み取ってバッファへ書き込む関数read_u32_intoを見ます。

Cargo.tomlに

[dependencies]
byteorder = "1.4.3"

を書き加えて以下を実行：

byteorder_demo/src/main.rs

use byteorder::{BigEndian, ByteOrder, LittleEndian};

fn main() {
    let bytes: [u8; 12] = [
        0x01, 0x23, 0x45, 0x67, 0x89, 0xab, 0xcd, 0xef, 0x11, 0x22, 0x33, 0x44,
    ];
    let mut words: [u32; 3] = [0; 3];
    LittleEndian::read_u32_into(&bytes, &mut words);
    println!("{:x?}", words); //[67452301, efcdab89, 44332211]
    assert_eq!(words[0], 0x_67_45_23_01);
    assert_eq!(words[1], 0x_ef_cd_ab_89);
    assert_eq!(words[2], 0x_44_33_22_11);

    BigEndian::read_u32_into(&bytes, &mut words);
    println!("{:x?}", words); //[1234567, 89abcdef, 11223344]
    assert_eq!(words[0], 0x_01_23_45_67);
    assert_eq!(words[1], 0x_89_ab_cd_ef);
    assert_eq!(words[2], 0x_11_22_33_44);
}

diy_md5/main.rs at main · roiban1344/diy_md5

ソースとなる長さ12のバイト列bytesに対して、wordsへとLittleEndian, BigEndianそれぞれの順で32ビットの語3つが読み取られています。この場合、LittleEndianのほうが"low-order bytes first"であるというのは、「4バイト=32ビットを単位とする語の下位に来るビットがバイトの列上では前方に来る」という意味になります。「語」の長さは場合によりますが、RFC1321ではここに挙げた例と同様、32ビットを単位とする取り決めになっています。

Rustで実装

RFC1321の3. MD5 Algorithm Description に沿って実装を進めていきます。アルゴリズムはStep 1.からStep 5.に分けて説明されています。

準備

Rustにおいて、文字列：Stringやstr型はUTF-8による符号化で有効なビット列であって、charやu8の配列とは明確に区別されています。そのため、atや[]で特定位置のバイトを自由に取り出すことはできません^[4]。

文字列をバイトのベクターVec<u8>として扱うため、まず以下の変換を行います。

let mut message = message.as_bytes().to_vec();

後述する手順でパディング（延長）を行うため、mutで可変変数として定義しました。

また、所有権を奪わないため、パラメータの型は&strにとります。つまり以下のような関数の中身を実装することになります。

fn md5(message: &str) -> u128{
    todo!();
}

なお、UTF-8で文字列を符号化してバイトの列で表現する方法は一意であるため、as_bytesにバイト順の問題は生じません^[5]。「語」（この場合単一の文字の符号位置のこと）が固定長である影響でバイト順が問題になるUTF-16やUTF-32とは対照的です。

Step 1. Append Padding Bits

mod 512で448に合同なビット長まで、メッセージに対してパディングを行います。

具体的には、まず1を追加して、その後ろに0を続けます^[6]。パディングはメッセージのビット長が初めからmod 512で448に等しくても必ず行い、したがって付加されるビット列の長さは最小1（メッセージのビット長 $\equiv 447\mod 512$ のとき）、最大512（メッセージのビット長 $\equiv 448\mod 512$ のとき）になります。

例：

\begin{align*} \underbrace{1110010}_{元のビット長\,8}\rightarrow\underbrace{\underbrace{1110010}_{元のビット長\,8}\underbrace{1}_{1\,個}\underbrace{00\cdots 0}_{439\,個}}_{パディング後のビット長\,448} \end{align*}

いま扱うのはバイトの列限定です。したがって、メッセージの長さは8の倍数ビットです。メッセージのバイト長 $m$ の64に対する剰余を $m\equiv r\,(0\leq 64\leq 63)\,\mod64$ としましょう。 $8r$ が $mod 512$ でのビット長の剰余です。

パディングのバイト長 $p$ はより具体的には、

\begin{align*} p &= \left\{\begin{array}{cc} 56 - r & r < 56 \\ 120 - r & r \geq 56 \end{array}\right. \end{align*}

で与えられます。以下に示すRustコードではこの式中の変数はそれぞれ

$m$ : message_len_in_bytes
$r$ : mod_64
$p$ : padding_len_in_bytes

と対応付けます。

    let message_len_in_bytes = message.len();
    let mod_64 = message_len_in_bytes & 0b111111; //下位6ビットが64に対する剰余
    let padding_len_in_bytes = if mod_64 < 56 {
        56 - mod_64
    } else {
        120 - mod_64
    };

このステップでありうる最大長さのパディングを定数配列 PADDING として持っておいて、この配列のスライスによってメッセージを延長します。

const PADDING: [u8; 64] = [
    0b10000000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0,
];

message.extend(&PADDING[..padding_len_in_bytes]);

Step 2. Append Length

元のメッセージのビット長の64ビット二進表現（つまりu64）を32ビットの語2つとみなし、下位の語が先に来るように、Step1.でパディングされたメッセージの末尾に追加します。

？？？

These bits are appended as two 32-bit words and appended low-order word first in accordance with the previous conventions.

......というのはRFC1321に書かれている通りの説明ですが、要は64ビットの語として要素の各バイトをリトルエンディアンで並べます。推測するに、2. Terminology amd Notationの項で"word"を32ビットと規定しているため、こういうやや回りくどい説明になっているのでしょう。

例：元のメッセージが $776379490800$ ビット $=97047436350$ バイトの場合

\begin{align*} (776379490800)_{10}= & ({\rm b4c3d2e1f0})_{16}\\ & \phantom{00000} \hspace{-7pt}\searrow\hspace{-12pt}\swarrow\\ \underbrace{\underbrace{{\rm 1a2b3c4d5e\cdots ff}}_{97047436350バイト} \underbrace{8000\cdots 0}_{58バイト}}_{\rm Step 1でパディング後のメッセージ} & \underbrace{{\rm f0e1d2c3b4}}_{8バイト} \end{align*}

ビット長の64ビット表現 8 * message_len_in_bytes as u64 を、u8 の長さ8の配列 buf へリトルエンディアンで書き込み、それをStep 1でパディングされたメッセージ末尾に更に添加します。

    let mut buf: [u8; 8] = [0; 8];
    LittleEndian::write_u64(&mut buf, 8 * message_len_in_bytes as u64);
    message.extend(&buf);

Step 3. Initialize MD buffer

4つの32ビット（=4バイト）の語（バッファー）を以下の値で初期化します。

word A: 01 23 45 67
word B: 89 ab cd ef
word C: fe dc ba 98
word D: 76 54 32 10

Numeric separatorを挟んでRustのリテラルで表現するとこうなります：

let mut A: u32 = 0x_67_45_23_01;
let mut B: u32 = 0x_ef_cd_ab_89;
let mut C: u32 = 0x_98_ba_dc_fe;
let mut D: u32 = 0x_10_32_54_76;

ただし、Rustのデフォルトの変数の命名に関する規定では、変数は小文字のスネークケースであることが求められます。#[allow(non_snake_case)] 属性を付けることで警告を抑制することは可能^[7]ですが、ここでは命名指針に従うことにして小文字で定義します。

let mut a = 0x67452301u32;
let mut b = 0xefcdab89u32;
let mut c = 0x98badcfeu32;
let mut d = 0x10325476u32;

ちなみに a=1732584193 だけ素数です。

Step 4.Process Message in 16-word Blocks

パディングされた512の倍数のビット長のメッセージを、先頭から32ビットの語を単位として、Step 3で定義された A,B,C,D へ作用させていきます。RFC1321では"Step 4."で一括りにされていますが、ここでは4.1から4.3に分けて見ていきます。

4.1 補助関数

32ビットの語3つをとって32ビットの語1つを返す、4つの補助的な関数を定義します。

F(X,Y,Z) = XY v not(X) Z
G(X,Y,Z) = XZ v Y not(Z)
H(X,Y,Z) = X xor Y xor Z
I(X,Y,Z) = Y xor (X v not(Z))

ビット演算に関するRFC1321中の上の表現とRustでの表現は下表で対応付けられます。

RFC1321	Rust
`XY`	`X & Y`
`X v Y`	`X
`not(X)`	`!X`
`X xor Y`	`X^Y`

演算子の評価順序はRFC132とRustで一致しています。

この関数も小文字で定義します。i をループのダミー変数以外の用途で用いるのが微妙ですが、まあ登場しないので問題ないでしょう。

fn f(x: u32, y: u32, z: u32) -> u32 {
    x & y | !x & z
}

fn g(x: u32, y: u32, z: u32) -> u32 {
    x & z | y & !z
}

fn h(x: u32, y: u32, z: u32) -> u32 {
    x ^ y ^ z
}

fn i(x: u32, y: u32, z: u32) -> u32 {
    y ^ (x | !z)
}

ところで、このRFC1321のこの箇所には以下のような記述があります。

if the bits of X, Y, and Z are independent and unbiased, the each bit of F(X,Y,Z) will be independent and unbiased.

if the corresponding bits of X, Y, and Z are independent and unbiased, then each bit of G(X,Y,Z), H(X,Y,Z), and I(X,Y,Z) will be independent and unbiased.

どういうことでしょうか。

F,G,H,I は全てビット毎（bit-wise）な演算から成る関数であるため、引数と返り値の関係は (x,y,z)=(0,0,0),(0,0,1),...,(1,1,1) の8通りで本質的には尽きています。実際に計算してみると、以下のような関係になります。

    for xyz in 0..8 {
        let (x, y, z) = ((xyz >> 2) & 1, (xyz >> 1) & 1, xyz & 1);
        let f = f(x, y, z) & 1;
        let g = g(x, y, z) & 1;
        let h = h(x, y, z) & 1;
        let i = i(x, y, z) & 1;
        println!("{:03b} {} {} {} {}", xyz, f, g, h, i);
    }

出力：

xyz F G H I
000 0 0 0 1
001 1 0 1 0
010 0 1 1 0
011 1 0 0 1
100 0 0 1 1
101 0 1 0 1
110 1 1 0 0
111 1 1 1 0

各関数で返り値の0と1がちょうど4ずつ現れていることが観察できます。おそらくこのことを言っているのでしょう^[8]。

4.2 512ビットのブロックへの切り出し

パディングされて512の倍数のビット長になったメッセージの先頭から512ビット（=64バイト）の「ブロック」を切り出し、さらにそのブロックを16個の32ビット（=4バイト）の語に分割します。この16個の語を使ってA, B, C, Dをラウンド1から4までの操作（それぞれF, G, H, Iを利用）1セット（4.3で後述）で変換します。ブロックが全てなくなるまで同じセットを繰り返します。

「ブロックの切り出し」は std::io::Cursor の read によって実現します。read メソッドは Read トレイトに含まれるため、これもスコープに入れます。

use std::io::{Cursor, Read};

Cursor in std::io - Rust

Cursor は read で読み取った分次回の読み取り位置を移動するため、ループ変数による読み取り位置のコントロールは不要です。また、read の返り値は Result 型であるため、unwrap で値を取り出します。全て読み終わった後で尚読み取ろうとした場合 Err 型が返るためこれによってループを停止することができますが、読み取り回数は (メッセージのバイト長)/64 で決まっているため、この回数だけループすることにします。

    let mut cursor = Cursor::new(&message);
    let n = message.len() >> 6;

    for _ in 0..n {
        let mut block = [0; 64];
        cursor.read(&mut block).unwrap();
        todo!();　//4.3で実装
    }

512ビットのブロックを32ビットの語16個に分割します。既に登場しましたが、LittleEndian の read_u32_into のシグネチャは

fn read_u32_into(src: &[u8], dst: &mut [u32])

となっており、バイト列 src を32ビットの語の配列 dst に書き込むことができます。

let mut x: [0u32; 16] = [0; 16];
LittleEndian::read_u32_into(&block, &mut x);

こうして x に32ビットの語16個が書き込まれました。

4.3 バッファへ作用

A, B, C, D に対して語を作用させていきます。

処理はラウンド1から4まであり、ラウンド1は以下の疑似コードで表されます。

     /* Round 1. */
     /* Let [abcd k s i] denote the operation
          a = b + ((a + F(b,c,d) + X[k] + T[i]) <<< s). */
     /* Do the following 16 operations. */
     [ABCD  0  7  1]  [DABC  1 12  2]  [CDAB  2 17  3]  [BCDA  3 22  4]
     [ABCD  4  7  5]  [DABC  5 12  6]  [CDAB  6 17  7]  [BCDA  7 22  8]
     [ABCD  8  7  9]  [DABC  9 12 10]  [CDAB 10 17 11]  [BCDA 11 22 12]
     [ABCD 12  7 13]  [DABC 13 12 14]  [CDAB 14 17 15]  [BCDA 15 22 16]

ラウンド2, 3, 4では同様の操作を、F を G, H, I に置き換えて行います。

+はmod 32による和です。一見わざわざ説明するまでもないことのように思われますが、Rustはデバッグモードによるコンパイル時には、+演算子によるオーバーフローを検知してpanic!を起こすため注意が必要です。

最も荒っぽくこれを解消する方法はリリースモード(Cargoなら--releaseオプション)でコンパイルすることです。#![allow(arithmetic_overflow)]属性で明示的にチェックを抑制することもできます。

真っ当な解決方法は、数値型に定義されている wrapping_add を使うことです^[9]。少々読みにくくはなりますが、u32なら $\mod 2^{32}$ による計算を保証された挙動として実行できます。

T[i] は $2^{32} |\sin(i)|$ （引数はラジアン単位）の整数部分です（付録参照）。値は定数のテーブルとして用意しておきます。

const T: [u32; 65] = [
    0, 0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, 0xa8304613,
    0xfd469501, 0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e,
    0x49b40821, 0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x2441453, 0xd8a1e681,
    0xe7d3fbc8, 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8, 0x676f02d9,
    0x8d2a4c8a, 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60,
    0xbebfbc70, 0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x4881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8,
    0xc4ac5665, 0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, 0xffeff47d,
    0x85845dd1, 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1, 0xf7537e82, 0xbd3af235, 0x2ad7d2bb,
    0xeb86d391,
];

（T[0] は不使用だが記法を一致させるために保持）

<<< は左方向への巡回シフトです。rotate_left メソッドがu32に備わっています^[10]。

Rustによるコードでも略記法 [abcd k s i] を使えると嬉しそうです。関数やクロージャとして定義することもできますが、x のスコープを正しく考慮したり、F,G,H,I の部分を入れ替えることを考えるとなかなか面倒です。

「マクロを定義するマクロ」で実現してしまうことにします。

macro_rules! define_operation {
    ($macro_name:ident, $f:expr) => {
        macro_rules! $macro_name {
            ($a:expr, $b:expr, $c:expr, $d:expr, $k:expr, $s:expr, $i:expr) => {
                $a = $b.wrapping_add(
                    $a.wrapping_add(
                        $f($b, $c, $d).wrapping_add(x[$k]).wrapping_add(T[$i]),
                    )
                    .rotate_left($s),
                );
            };
        }
    };
}

このマクロを使うと、ラウンド1が以下のように書けます。角括弧を使ってRFC1321と見た目を揃えられる点も些細ながら嬉しいです。

define_operation!(ff, f);
ff![a, b, c, d, 0, 7, 1];
ff![d, a, b, c, 1, 12, 2];
ff![c, d, a, b, 2, 17, 3];
ff![b, c, d, a, 3, 22, 4];
ff![a, b, c, d, 4, 7, 5];
ff![d, a, b, c, 5, 12, 6];
ff![c, d, a, b, 6, 17, 7];
ff![b, c, d, a, 7, 22, 8];
ff![a, b, c, d, 8, 7, 9];
ff![d, a, b, c, 9, 12, 10];
ff![c, d, a, b, 10, 17, 11];
ff![b, c, d, a, 11, 22, 12];
ff![a, b, c, d, 12, 7, 13];
ff![d, a, b, c, 13, 12, 14];
ff![c, d, a, b, 14, 17, 15];
ff![b, c, d, a, 15, 22, 16];

同様に、マクロ定義+操作をラウンド2から4まで書きます。

ラウンド1から4終了後に、ラウンド1開始前の値をA, B, C, Dに加えます。これが512ビットのブロックに対する操作1セットです。

let (aa, bb, cc, dd) = (a, b, c, d);
/* Round 1 */
/* Round 2 */
/* Round 3 */
/* Round 4 */
a = a.wrapping_add(aa);
b = b.wrapping_add(bb);
c = c.wrapping_add(cc);
d = d.wrapping_add(dd);

この操作1セットを、ブロックがパディングされたメッセージの末尾に達するまで繰り返します。

Step 5. Output

Step 4の操作が完了した後、A,B,C,D をリトルエンディアンでビット列として結合します。例として、

A = 0xaeb0434b;
B = 0xcd2456e3;
C = 0x1810b995;
D = 0x31c23d9b;

なら、

4b43b0aee35624cd95b910189b3dc231

が最終的なハッシュ値になります。実装としては、長さ16のバイトの配列にリトルエンディアンで書き込んだ後、ビッグエンディアンで読み取ればよいです。

    let mut buf: [u8;16] = [0; 16];
    LittleEndian::write_u32_into(&[a, b, c, d], &mut buf);
    BigEndian::read_u128(&buf)

以下にコードの全体をまとめます。RFC1321 A.5に付記されているテストスイートを含めました。

コード全体

diy_md5/lib.rs

use byteorder::{ByteOrder, BE, LE};
use std::io::{Cursor, Read};

const PADDING: [u8; 64] = [
    0b10000000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0,
];

const T: [u32; 65] = [
    0, 0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, 0xa8304613,
    0xfd469501, 0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e,
    0x49b40821, 0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x2441453, 0xd8a1e681,
    0xe7d3fbc8, 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8, 0x676f02d9,
    0x8d2a4c8a, 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60,
    0xbebfbc70, 0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x4881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8,
    0xc4ac5665, 0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, 0xffeff47d,
    0x85845dd1, 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1, 0xf7537e82, 0xbd3af235, 0x2ad7d2bb,
    0xeb86d391,
];

fn f(x: u32, y: u32, z: u32) -> u32 {
    x & y | !x & z
}

fn g(x: u32, y: u32, z: u32) -> u32 {
    x & z | y & !z
}

fn h(x: u32, y: u32, z: u32) -> u32 {
    x ^ y ^ z
}

fn i(x: u32, y: u32, z: u32) -> u32 {
    y ^ (x | !z)
}

pub fn md5(message: &str) -> u128 {
    let mut message = message.as_bytes().to_vec();

    /* Step 1. Append Padding Bits */
    let message_len_in_bytes = message.len();
    let mod_64 = message_len_in_bytes & 0b111111;
    let padding_len_in_bytes = if mod_64 < 56 {
        56 - mod_64
    } else {
        120 - mod_64
    };
    message.extend(&PADDING[..padding_len_in_bytes]);

    /* Step 2. Append Length */
    let mut buf = [0; 8];
    LE::write_u64(&mut buf, 8 * message_len_in_bytes as u64);
    message.extend(&buf);

    /* Step 3. Initialize MD Buffer */
    let mut a = 0x67452301u32;
    let mut b = 0xefcdab89u32;
    let mut c = 0x98badcfeu32;
    let mut d = 0x10325476u32;

    /* Step 4. Process Message in 16-Word Blocks */
    let mut cursor = Cursor::new(&message);
    let n = message.len() >> 6;
    for _ in 0..n {
        let mut block = [0; 64];
        cursor.read(&mut block).unwrap();
        let mut x = [0; 16];
        LE::read_u32_into(&block, &mut x);
        macro_rules! define_operation {
            ($macro_name: ident, $func: ident) => {
                macro_rules! $macro_name {
                    ($a:expr, $b:expr, $c:expr, $d: expr, $k:expr, $s:expr, $i:expr) => {
                        $a = $b.wrapping_add(
                            $a.wrapping_add($func($b, $c, $d))
                                .wrapping_add(x[$k])
                                .wrapping_add(T[$i])
                                .rotate_left($s),
                        )
                    };
                }
            };
        }
        let (aa, bb, cc, dd) = (a, b, c, d);

        /* Round 1. */
        define_operation!(ff, f);
        ff![a, b, c, d, 0, 7, 1];
        ff![d, a, b, c, 1, 12, 2];
        ff![c, d, a, b, 2, 17, 3];
        ff![b, c, d, a, 3, 22, 4];
        ff![a, b, c, d, 4, 7, 5];
        ff![d, a, b, c, 5, 12, 6];
        ff![c, d, a, b, 6, 17, 7];
        ff![b, c, d, a, 7, 22, 8];
        ff![a, b, c, d, 8, 7, 9];
        ff![d, a, b, c, 9, 12, 10];
        ff![c, d, a, b, 10, 17, 11];
        ff![b, c, d, a, 11, 22, 12];
        ff![a, b, c, d, 12, 7, 13];
        ff![d, a, b, c, 13, 12, 14];
        ff![c, d, a, b, 14, 17, 15];
        ff![b, c, d, a, 15, 22, 16];

        /* Round 2. */
        define_operation!(gg, g);
        gg![a, b, c, d, 1, 5, 17];
        gg![d, a, b, c, 6, 9, 18];
        gg![c, d, a, b, 11, 14, 19];
        gg![b, c, d, a, 0, 20, 20];
        gg![a, b, c, d, 5, 5, 21];
        gg![d, a, b, c, 10, 9, 22];
        gg![c, d, a, b, 15, 14, 23];
        gg![b, c, d, a, 4, 20, 24];
        gg![a, b, c, d, 9, 5, 25];
        gg![d, a, b, c, 14, 9, 26];
        gg![c, d, a, b, 3, 14, 27];
        gg![b, c, d, a, 8, 20, 28];
        gg![a, b, c, d, 13, 5, 29];
        gg![d, a, b, c, 2, 9, 30];
        gg![c, d, a, b, 7, 14, 31];
        gg![b, c, d, a, 12, 20, 32];

        /* Round 3. */
        define_operation!(hh, h);
        hh![a, b, c, d, 5, 4, 33];
        hh![d, a, b, c, 8, 11, 34];
        hh![c, d, a, b, 11, 16, 35];
        hh![b, c, d, a, 14, 23, 36];
        hh![a, b, c, d, 1, 4, 37];
        hh![d, a, b, c, 4, 11, 38];
        hh![c, d, a, b, 7, 16, 39];
        hh![b, c, d, a, 10, 23, 40];
        hh![a, b, c, d, 13, 4, 41];
        hh![d, a, b, c, 0, 11, 42];
        hh![c, d, a, b, 3, 16, 43];
        hh![b, c, d, a, 6, 23, 44];
        hh![a, b, c, d, 9, 4, 45];
        hh![d, a, b, c, 12, 11, 46];
        hh![c, d, a, b, 15, 16, 47];
        hh![b, c, d, a, 2, 23, 48];

        /* Round 4. */
        define_operation!(ii, i);
        ii![a, b, c, d, 0, 6, 49];
        ii![d, a, b, c, 7, 10, 50];
        ii![c, d, a, b, 14, 15, 51];
        ii![b, c, d, a, 5, 21, 52];
        ii![a, b, c, d, 12, 6, 53];
        ii![d, a, b, c, 3, 10, 54];
        ii![c, d, a, b, 10, 15, 55];
        ii![b, c, d, a, 1, 21, 56];
        ii![a, b, c, d, 8, 6, 57];
        ii![d, a, b, c, 15, 10, 58];
        ii![c, d, a, b, 6, 15, 59];
        ii![b, c, d, a, 13, 21, 60];
        ii![a, b, c, d, 4, 6, 61];
        ii![d, a, b, c, 11, 10, 62];
        ii![c, d, a, b, 2, 15, 63];
        ii![b, c, d, a, 9, 21, 64];

        a = a.wrapping_add(aa);
        b = b.wrapping_add(bb);
        c = c.wrapping_add(cc);
        d = d.wrapping_add(dd);
    }

    /* Step 5. Output */
    let mut buf = [0; 16];
    LE::write_u32_into(&[a, b, c, d], &mut buf);
    BE::read_u128(&buf)
}

[cfg(test)]
mod tests {
    #[test]
    fn it_works() {
        use super::*;

        assert_eq!(md5(""), 0xd41d8cd98f00b204e9800998ecf8427e);
        assert_eq!(md5("a"), 0x0cc175b9c0f1b6a831c399e269772661);
        assert_eq!(md5("abc"), 0x900150983cd24fb0d6963f7d28e17f72);
        assert_eq!(md5("message digest"), 0xf96b697d7cb7938d525a2f31aaf161d0);
        assert_eq!(
            md5("abcdefghijklmnopqrstuvwxyz"),
            0xc3fcd3d76192e4007dfb496cca67e13b
        );
        assert_eq!(
            md5("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789"),
            0xd174ab98d277d9f5a5611c2c9f419d9f
        );
        assert_eq!(
            md5("12345678901234567890123456789012345678901234567890123456789012345678901234567890"),
            0x57edf4a22be3c955ac49da2e2107b67a
        );
    }
}

diy_md5/lib.rs at main · roiban1344/diy_md5

なお、LE, BEはそれぞれLittleEndian, BigEndianのエイリアスです。また、型推論が可能な場所の注釈は省略しました。

テストケースはASCIIの範囲しか扱っていませんが、漢字や絵文字を含んでいてもハッシュ化できます。

例えば、ゴジラSPにも登場した"解けばわかる"のハッシュ値が0x14980c8b8a96fd9e279796a61cf82c9cであることも確かめられます。

diy_md5/main.rs

fn main(){
    //ひらがなと漢字を含む文字列
    let hash = diy_md5::md5("解けばわかる");
    println!("{:032x}", hash);
    assert_eq!(0x14980c8b8a96fd9e279796a61cf82c9c, hash);

    //表記揺れで全く異なる値になる
    let hash = diy_md5::md5("解けば分かる");
    println!("{:032x}", hash); //606461eb515ea6d825117fbe76965899

    //emoji
    let hash = diy_md5::md5("🐶");
    println!("{:032x}", hash); //be0f7766d0c41a4386d47e18e8b91e15
}

diy_md5/main.rs at main · roiban1344/diy_md5

付録　sin(i)を計算

おまけです。Step 4.に登場する定数たち $T_i = \lfloor 4294967296\,|\sin(i)|\rfloor (1\leq i \leq 64)$ を計算しましょう。

テイラーの定理

テイラーの定理が基本です：

区間 $[0, a]\subset \mathbb{R}$ で $n+1$ 回微分可能な関数 $f:\mathbb{R}\rightarrow \mathbb{R}$ について、任意の $x\in [0,a]$ に対して

\begin{align*} f(x) = \sum_{k=0}^n \frac{f^{(k)}(0)}{k!}x^k + \frac{f^{(n+1)}(c)\,x^{n+1}}{(n+1)!} \end{align*}

となる $c\in (0,a)$ が存在する。

馴染み深い定理ですが、誤差を正しく評価して計算するとなるとなかなか大変です。

第 $n$ 部分和を $S_{n}:=\sum_{k=0}^n f^{(k)}(0)x^k/k!$ 、その誤差項を $R_{n}:=f^{(n+1)}(c)x^{n+1}/(n+1)!$ とします。

$f(x) = \sin x$ の場合、 $|f^{(n)}(c)|\leq 1$ であるため、全ての $x\in \mathbb{R}$ で $\lim_{n\rightarrow \infty} R_{n}=0$ と誤差項が $0$ に収束し、部分和は真の値に収束します。さらに、 $W_n:=|x^{n+1}/(n+1)!|$ とすると、

\begin{align} S_n-W_n \leq \sin x \leq S_n+W_n \end{align}

は厳密な評価になります。 $\lfloor 4294967296| \sin x|\rfloor$ の値が確定するためには、

$S_n-W_n$ と $S_n+W_n$ の符号が一致
$\lfloor 4294967296 |S_n-W_n|\rfloor = \lfloor 4294967296 |S_n+W_n|\rfloor$

となればよいです。

方法1. 力尽く

まずは愚直に全ての $\sin i$ に対して上の評価式を適用します。

$x$ が有理数（特に整数）なら上下限ともに有理数の範囲で計算できるため、巨大な整数さえ扱うことができればあとは計算するだけです。

巨大な整数やそれを成分に持つ有理数は num クレートによって扱えます。

num - crates.io: Rust Package Registry

実装

sin_brute_force/src/main.rs

use num::bigint::{BigInt, Sign};
use num::rational::BigRational;
use num::Signed;
use num::{One, Zero};
use std::result::Result;

const T: [u32; 65] = [
    0, 0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, 0xa8304613,
    0xfd469501, 0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e,
    0x49b40821, 0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x2441453, 0xd8a1e681,
    0xe7d3fbc8, 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8, 0x676f02d9,
    0x8d2a4c8a, 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60,
    0xbebfbc70, 0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x4881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8,
    0xc4ac5665, 0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, 0xffeff47d,
    0x85845dd1, 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1, 0xf7537e82, 0xbd3af235, 0x2ad7d2bb,
    0xeb86d391,
];

fn main() {
    let c = BigRational::from_integer(BigInt::new(Sign::Plus, vec![0, 1])); //=4294967296
    for x in 1..=64 {
        for i in 0.. {
            let (l, r) = approx_range(x, i);
            if l.signum() != r.signum() {
                continue;
            }
            let a = num::abs(&c * &l).floor().to_integer();
            let b = num::abs(&c * &r).floor().to_integer();
            if a == b {
                let t = to_u32(&a).unwrap();
                assert_eq!(T[x as usize], t);
                println!("{} {:x} {}", x, t, i);
                break;
            }
        }
    }
}

fn to_u32(n: &BigInt) -> Result<u32, ()> {
    let (sign, v) = n.to_u32_digits();
    if sign != Sign::Plus || v.len() != 1 {
        Result::Err(())
    } else {
        Result::Ok(v[0])
    }
}

fn approx_range(x: u32, n: u32) -> (BigRational, BigRational) {
    let mut sum = BigRational::zero();
    for i in 0..=n {
        sum += match i % 4 {
            1 => BigRational::new(BigInt::new(Sign::Plus, vec![x]).pow(i), factorial(i)),
            3 => BigRational::new(BigInt::new(Sign::Minus, vec![x]).pow(i), factorial(i)),
            _ => BigRational::zero(),
        };
    }
    let err = num::abs(BigRational::new(
        BigInt::new(Sign::Plus, vec![x]).pow(n + 1),
        factorial(n + 1),
    ));
    (&sum - &err, &sum + &err)
}

fn factorial(n: u32) -> BigInt {
    if n == 0 || n == 1 {
        BigInt::one()
    } else if n > 1 {
        BigInt::new(Sign::Plus, vec![n]) * factorial(n - 1)
    } else {
        panic!();
    }
}

diy_md5/main.rs at main · roiban1344/diy_md5

出力

1 d76aa478 13
2 e8c7b756 17
3 242070db 21
4 c1bdceee 25
5 f57c0faf 28
6 4787c62a 31
7 a8304613 34
8 fd469501 37
9 698098d8 41
10 8b44f7af 43
11 ffff5bb1 46
12 895cd7be 49
13 6b901122 53
14 fd987193 55
15 a679438e 57
16 49b40821 60
17 f61e2562 64
18 c040b340 66
19 265e5a51 68
20 e9b6c7aa 72
21 d62f105d 75
22 2441453 78
23 d8a1e681 80
24 e7d3fbc8 82
25 21e1cde6 86
26 c33707d6 89
27 f4d50d87 92
28 455a14ed 94
29 a9e3e905 96
30 fcefa3f8 101
31 676f02d9 104
32 8d2a4c8a 106
33 fffa3942 109
34 8771f681 111
35 6d9d6122 114
36 fde5380c 115
37 a4beea44 119
38 4bdecfa9 122
39 f6bb4b60 125
40 bebfbc70 126
41 289b7ec6 130
42 eaa127fa 132
43 d4ef3085 135
44 4881d05 137
45 d9d4d039 140
46 e6db99e5 145
47 1fa27cf8 147
48 c4ac5665 149
49 f4292244 152
50 432aff97 154
51 ab9423a7 158
52 fc93a039 161
53 655b59c3 163
54 8f0ccc92 166
55 ffeff47d 168
56 85845dd1 171
57 6fa87e4f 173
58 fe2ce6e0 176
59 a3014314 179
60 4e0811a1 183
61 f7537e82 184
62 bd3af235 186
63 2ad7d2bb 190
64 eb86d391 192

出力の各行は

i T_iの計算値 値が確定するのに必要な最小項数（上の式のn）

となっています。

assert_eq! でRFC1321に使われている値 T との一致を検証して、正しく計算できていることが確かめられました。

しかしこの方法は著しくエレガントさにかけます。結果に示した通り、 $i=64$ を計算するために、実際に $192!$ という巨大な値を計算しています。このとき、 $\sin 64$ の近似に用いられた区間の上下限値は

\begin{align*} S_n-W_n &= \frac{3651657157709632429621773107660790167468774167256674373270738569599674653579034218585481015992170214588323933855116913299410112610029504742220078937750554441083303523402000431070075574889112022950255030494439411312602987258791078297807139173049341047292712976990157835013751547718338437807893946114368}{3969080228748870734062672344953331518424707253458286072182412738360339508651264331715423002947493777172551554950501275192830478121245439671569386475282932022860772435425373575751535445841740829756037427703853970095093624941269207059267250157681021055299110269120835245786338418838568031787872314453125},\\ S_n+W_n &= \frac{8033645747959111499935260742681924724949722450998409251413353623870557859092481226822540268075253454840630505025928432680798495364051129427080970493420805669525683461644403784829709103014815311212742560679093383259998533966005777482736141375152977938267990628131907690468131425801523810792538686415552}{8731976503247515614937879158897329340534355957608229358801308024392746919032781529773930606484486309779613420891102805424227051866739967277452650245622450450293699357935821866653377980851829825463282340948478734209205974870792255530387950346898246321658042592065837540729944521444849669933319091796875}. \end{align*}

という過剰に巨大な分子分母を持つ有理数です。

方法2. 加法定理による再帰

求めたい $\sin$ の引数は $1$ ラジアンの整数倍であるため、 $\cos 1$ と組み合わせて加法定理を繰り返し使うことで計算できるはずです。

ただし、誤差を厳密に取り扱うため、和差積を取る際には簡易かつ限定的な区間演算^[11]を使います。

誤差は演算ごとに拡大するため、初期値となる $\sin 1, \cos 1$ はテイラー展開による評価で

\begin{align} \sin 1 &\in \left[ \frac{7761199951101802512}{2^{63}}, \frac{7761199951101802513}{2^{63}} \right],\\ \cos 1 &\in \left[ \frac{4983409179392355912}{2^{63}}, \frac{4983409179392355913}{2^{63}} \right] \end{align}

と誤差範囲を $1/2^{63}$ にまで狭めておきます。この値だけはテイラー展開を使って（第20項までで得られる）計算しておきます。

実装

sin_recursion/src/main.rs

use num::integer;
use std::ops::{Add, Mul, Sub};

//RFC1321記載の値
const T: [u32; 65] = [
    0, 0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, 0xa8304613,
    0xfd469501, 0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e,
    0x49b40821, 0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x2441453, 0xd8a1e681,
    0xe7d3fbc8, 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8, 0x676f02d9,
    0x8d2a4c8a, 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60,
    0xbebfbc70, 0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x4881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8,
    0xc4ac5665, 0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, 0xffeff47d,
    0x85845dd1, 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1, 0xf7537e82, 0xbd3af235, 0x2ad7d2bb,
    0xeb86d391,
];

const RANGE_MAX: usize = 64;

macro_rules! frac_from_interval {
    [$min:expr, $max:expr] => [
        Frac::new(Interval::new($min, $max))
    ]
}

fn main() {
    let mut sin = vec![Option::<Frac>::None; RANGE_MAX + 1];
    let mut cos = vec![Option::<Frac>::None; RANGE_MAX + 1];
    sin[1] = Some(frac_from_interval![
        7761199951101802512,
        7761199951101802513
    ]);
    cos[1] = Some(frac_from_interval![
        4983409179392355912,
        4983409179392355913
    ]);

    let print_data = |i: i32, frac: &Frac| -> () {
        let t = frac.times_4294967296_to_u32();
        assert_eq!(t, T[i as usize]);
        println!("{0} {1:#x} {2}", i, t, frac.num.max - frac.num.min)
    };

    print_data(1, &sin[1].unwrap());

    for i in 2..=RANGE_MAX {
        let s1 = sin[i / 2].unwrap();
        let c1 = cos[i / 2].unwrap();
        let s2 = sin[(i + 1) / 2].unwrap();
        let c2 = cos[(i + 1) / 2].unwrap();
        sin[i] = Some(s1 * c2 + c1 * s2);
        cos[i] = Some(c1 * c2 - s1 * s2);
        print_data(i as i32, &sin[i].unwrap());
    }
}

//区間演算に用いる「区間」。
[derive(Copy, Clone, Debug)]
struct Interval {
    min: i128,
    max: i128,
}

impl Interval {
    fn new(min: i128, max: i128) -> Interval {
        if min < 0 && 0 < max {
            panic!("Cannot contain 0.");
        } else if min > max {
            panic!("min = {} < max = {} is not satisfied.", min, max);
        }
        Interval { min, max }
    }

    fn signum(&self) -> i32 {
        if self.min > 0 && self.max > 0 {
            1
        } else if self.min < 0 && self.max < 0 {
            -1
        } else {
            panic!("Contains 0 {:?}", self);
        }
    }
}

//和差積の未実装。
impl Add for Interval {
    type Output = Self;
    fn add(self, rhs: Self) -> Self {
        Interval::new(self.min + rhs.min, self.max + rhs.max)
    }
}

impl Sub for Interval {
    type Output = Self;
    fn sub(self, rhs: Self) -> Self {
        Interval::new(self.min - rhs.max, self.max - rhs.min)
    }
}

impl Mul for Interval {
    type Output = Self;
    fn mul(self, rhs: Self) -> Self {
        let a = self.min;
        let b = self.max;
        let c = rhs.min;
        let d = rhs.max;
        let (min, max) = match (self.signum(), rhs.signum()) {
            (1, 1) => (a * c, b * d),
            (1, -1) => (b * c, a * d),
            (-1, 1) => (a * d, b * c),
            (-1, -1) => (b * d, a * c),
            _ => unreachable!(),
        };
        Interval::new(min, max)
    }
}

//2^63を分母に持ち、[-2^63,2^63]に含まれる区間を分母に持つ有理数。
[derive(Debug, Copy, Clone)]
struct Frac {
    num: Interval, //denom=2^63
}

impl Frac {
    fn new(num: Interval) -> Self {
        if num.min < -(1 << 63) || 1 << 63 < num.max {
            panic!("Out of range")
        }
        Self { num }
    }

    //2^32倍の整数部分を返す
    fn times_4294967296_to_u32(&self) -> u32 {
        let min = integer::div_floor(self.num.min.abs(), 1 << 31);
        let max = integer::div_floor(self.num.max.abs(), 1 << 31);
        if min == max {
            min as u32
        } else {
            println!("{:?}", self.num);
            panic!("min={}, max={}", min, max);
        }
    }
}

//和差積のみ実装。
impl Add for Frac {
    type Output = Self;
    fn add(self, rhs: Self) -> Self {
        Frac::new(self.num + rhs.num)
    }
}

impl Sub for Frac {
    type Output = Self;
    fn sub(self, rhs: Self) -> Self {
        Frac::new(self.num - rhs.num)
    }
}

impl Mul for Frac {
    type Output = Self;
    fn mul(self, rhs: Self) -> Self {
        let num_prod = self.num * rhs.num;
        let floor = integer::div_floor(num_prod.min, 1 << 63);
        let ceil = integer::div_ceil(num_prod.max, 1 << 63);
        frac_from_interval![floor, ceil]
    }
}

diy_md5/main.rs at main · roiban1344/diy_md5

出力

1 0xd76aa478 1
2 0xe8c7b756 4
3 0x242070db 9
4 0xc1bdceee 14
5 0xf57c0faf 20
6 0x4787c62a 22
7 0xa8304613 32
8 0xfd469501 40
9 0x698098d8 48
10 0x8b44f7af 54
11 0xffff5bb1 55
12 0x895cd7be 56
13 0x6b901122 73
14 0xfd987193 90
15 0xa679438e 96
16 0x49b40821 98
17 0xf61e2562 110
18 0xc040b340 130
19 0x265e5a51 139
20 0xe9b6c7aa 150
21 0xd62f105d 132
22 0x2441453 114
23 0xd8a1e681 137
24 0xe7d3fbc8 160
25 0x21e1cde6 179
26 0xc33707d6 198
27 0xf4d50d87 205
28 0x455a14ed 208
29 0xa9e3e905 237
30 0xfcefa3f8 272
31 0x676f02d9 257
32 0x8d2a4c8a 244
33 0xfffa3942 257
34 0x8771f681 276
35 0x6d9d6122 319
36 0xfde5380c 368
37 0xa4beea44 346
38 0x4bdecfa9 318
39 0xf6bb4b60 357
40 0xbebfbc70 398
41 0x289b7ec6 386
42 0xeaa127fa 368
43 0xd4ef3085 293
44 0x4881d05 232
45 0xd9d4d039 297
46 0xe6db99e5 378
47 0x1fa27cf8 405
48 0xc4ac5665 430
49 0xf4292244 423
50 0x432aff97 406
51 0xab9423a7 479
52 0xfc93a039 560
53 0x655b59c3 539
54 0x8f0ccc92 514
55 0xffeff47d 511
56 0x85845dd1 514
57 0x6fa87e4f 586
58 0xfe2ce6e0 672
59 0xa3014314 658
60 0x4e0811a1 624
61 0xf7537e82 655
62 0xbd3af235 680
63 0x2ad7d2bb 677
64 0xeb86d391 672

出力の各行は

i T_iの計算値 (区間の上限-下限)*2^63

となっています。

計算の各所に前提が崩れた場合に備えた意図的なpanic!を含んでいますが、必要な範囲内で無事に計算を終えることができました。assert_eq!でRFC1321記載の値と一致することも確かめられています。

ちなみにもっと先まで計算を進めた場合、i=8195で初めてpanic!します。

方法3. 浮動小数点数

浮動小数点数でも計算してみます。f32 で計算した場合は精度が足りず結果が一致しませんが、f64 は全てのケースでRFC1321記載の値と一致します。

f32 と f64 を統一的に扱うため、num-traits を利用しています。

num-traits - crates.io: Rust Package Registry

実装

sin_float/src/main.rs

use num_traits::{Float, FromPrimitive};

//RFC1321記載の値
const T: [u32; 65] = [
    0, 0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, 0xa8304613,
    0xfd469501, 0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e,
    0x49b40821, 0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x2441453, 0xd8a1e681,
    0xe7d3fbc8, 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8, 0x676f02d9,
    0x8d2a4c8a, 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60,
    0xbebfbc70, 0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x4881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8,
    0xc4ac5665, 0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, 0xffeff47d,
    0x85845dd1, 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1, 0xf7537e82, 0xbd3af235, 0x2ad7d2bb,
    0xeb86d391,
];

fn main() {
    for i in 1..=64 {
        let t_f32 = t_float::<f32>(i);
        let t_f64 = t_float::<f64>(i);
        let t = T[i as usize] as i64;
        let diff_f32 = (t_f32 - t).abs();
        let diff_f64 = (t_f64 - t).abs();
        assert_eq!(t, t_f64);
        println!("{:>2} {:08x} {:>3} {}", i, t, diff_f32, diff_f64);
    }
}

//floor(4294967296*abs(sin(i)))をf32, f64で計算
fn t_float<T: Float + FromPrimitive>(i: i64) -> i64 {
    let sin = T::from_i64(i).unwrap().sin();
    let x = sin.abs() * T::from_i64(4294967296_i64).unwrap();
    T::to_i64(&x.floor()).unwrap()
}

diy_md5/main.rs at main · roiban1344/diy_md5

出力

 1 d76aa478 120 0
 2 e8c7b756  86 0
 3 242070db  27 0
 4 c1bdceee  18 0
 5 f57c0faf  81 0
 6 4787c62a  42 0
 7 a8304613  19 0
 8 fd469501   1 0
 9 698098d8  40 0
10 8b44f7af  81 0
11 ffff5bb1  79 0
12 895cd7be  66 0
13 6b901122  34 0
14 fd987193 109 0
15 a679438e 114 0
16 49b40821  33 0
17 f61e2562  98 0
18 c040b340  64 0
19 265e5a51  17 0
20 e9b6c7aa  86 0
21 d62f105d  93 0
22 02441453   1 0
23 d8a1e681 127 0
24 e7d3fbc8  56 0
25 21e1cde6  26 0
26 c33707d6  42 0
27 f4d50d87 121 0
28 455a14ed  19 0
29 a9e3e905   5 0
30 fcefa3f8   8 0
31 676f02d9  39 0
32 8d2a4c8a 118 0
33 fffa3942  66 0
34 8771f681 127 0
35 6d9d6122  34 0
36 fde5380c  12 0
37 a4beea44  68 0
38 4bdecfa9  41 0
39 f6bb4b60  96 0
40 bebfbc70 112 0
41 289b7ec6   6 0
42 eaa127fa   6 0
43 d4ef3085 123 0
44 04881d05   3 0
45 d9d4d039  57 0
46 e6db99e5  27 0
47 1fa27cf8   8 0
48 c4ac5665 101 0
49 f4292244  68 0
50 432aff97  23 0
51 ab9423a7  89 0
52 fc93a039  57 0
53 655b59c3  61 0
54 8f0ccc92 110 0
55 ffeff47d 125 0
56 85845dd1  47 0
57 6fa87e4f  49 0
58 fe2ce6e0  32 0
59 a3014314  20 0
60 4e0811a1  33 0
61 f7537e82 126 0
62 bd3af235  53 0
63 2ad7d2bb   5 0
64 eb86d391 111 0

各行の値は、

i T[i] (f32で計算した場合の誤差) (f64で計算した場合の誤差)

です。再帰で実装した場合に計算可能な $1 \leq i\leq 8194$ の全ての場合で一致することも確かめられます。

とはいえ、これは「運が良かった」と捉えるべきものでしょう。 $\lfloor 4294967296 |\sin(i)| \rfloor$ という厳密な整数値を必要とする場合には信頼できない値です。

まとめ

本記事ではRFC1321に沿ってRustでMD5ハッシュアルゴリズムを実装しました。

なかなか慣れない所有権システムに代表される「Rust特有の躓きどころ」はあまりなかった一方で、バイト順を扱うための外部クレートの導入が実に快適だったり、文字列とバイトの列が厳格に区別されていたりする点にRustらしさを味わうことができるように思います。

本筋からは逸れた内容ですが、 $T_i$ の値という数学的に厳密な結果をこうしたプログラムで得られるのは何度経験しても良いものです。

SHA-1やSHA-2もほぼ同様の手順に基づくハッシュアルゴリズムです。いずれ実装してみましょう。

また、今回は含められなかったのがパフォーマンス測定です。ベンチマークテストの方法などを知ったうえで取り組むことは今後の課題です。

脚注

完全新作TVアニメシリーズ「ゴジラシンギュラポイント Godzilla Singular Point」公式サイト。特に第4・10話。試行錯誤では決して合成方法を発見できない物質「アーキタイプ」の喩えとして初めて登場。 ↩︎
RSA暗号のRの暗号学者。 ↩︎
command line - How to get the MD5 hash of a string directly in the terminal? - Ask Ubuntu ↩︎
Storing UTF-8 Encoded Text with Strings - The Rust Programming Language
String in std::string - Rust ↩︎
unicode - Isn’t on big endian machines UTF-8's byte order different than on little endian machines? So why then doesn’t UTF-8 require a BOM? - Stack Overflow ↩︎
1を追加するのはこの操作を可逆にして、自明な衝突を防ぐためでしょう。 ↩︎
rust - Is there some way to not show a warning for non snake case identifiers? - Stack Overflow ↩︎
ただ"independent and unbiased"という表現が気になります。"Unbiased"は数え上げればいいとして、"independent"という性質はこの結果から確かめられているのでしょうか？ ↩︎
u32 wrapping_add - Rust
How can integer overflow protection be turned off? - Stack Overflow ↩︎
u32 rotate_left - Rust ↩︎
区間演算の実装について(1) - kashiの日記などを参照。Rust用の精度保証付き数値計算用ライブラリを使いたいところですが、今回は実装してしまいました。 ↩︎

GitHubで編集を提案

概要

背景

MD5とは

バイト順の問題

Rustで実装

準備

Step 1. Append Padding Bits

Step 2. Append Length

Step 3. Initialize MD buffer

Step 4.Process Message in 16-word Blocks

4.1 補助関数

4.2 512ビットのブロックへの切り出し

4.3 バッファへ作用

Step 5. Output

付録 sin(i)を計算

テイラーの定理

方法1. 力尽く

方法2. 加法定理による再帰

方法3. 浮動小数点数

まとめ

Discussion

付録　sin(i)を計算