8ビットの整数で53ビットの倍精度を組み立てる。中国剰余定理をGPUに載せたOzaki Scheme II

数学を食わず嫌いして後悔したかも

안녕하신게라！パナソニック コネクト株式会社クラウドソリューション部の加賀です。

前回の記事では、AI向けGPUで倍精度の回路が減るなか、Ozaki Scheme IIがCUDA 13.4のcuBLASに入ったことを紹介しました。そのうえで、GPUを選ぶときの要件は「倍精度の回路があること」ではなく、本当に欲しい「倍精度と同じ精度が出ること」で書ける、という話をしました。その中で、Ozaki Scheme IIについては「中国剰余定理を使う方式」とだけ書き、中身は別の記事で扱うとしていました。

今回はその中身です。前半では、中国剰余定理がどんな定理で、どのように広がってきたかを説明します。後半では、Ozaki Scheme IIがこの定理をどう使い、GPUで動かすために何を工夫したかを説明します。最後に、エミュレーションで足りるかを判断する観点をまとめます。

正直に言うと、1500年以上前の算術書にあった問題から始まった定理が、いまGPUで倍精度の行列積を支えているとまでは思っていませんでした。調べるほど、数学をもっとちゃんとやっておけば良かったかなー？と、今さら少し後悔しています。そんなわけで、この記事は難しい数式を使わず、小学校で習う割り算の余りの範囲で書きます（数式で書けるほど理解できていない、というのが本音です）。

記事の中では、次の言葉を使います。

言葉	意味
法	割る数のことです。「7で割った余り」なら、7が法です
行列積	行列どうしの掛け算です。中身は、掛け算と足し算だけの繰り返しです
テンソルコア	GPUの中にある、行列積専用の演算器です。AI向けのGPUでは、8ビットの整数（INT8）の行列積を特に速く計算できます

中国剰余定理の歴史と仕組み

孫子算経からガウスまで

3〜5世紀ごろの中国の算術書「孫子算経」に、次のような問題があります。

数の分からない物がある。3つずつ数えると2余り、5つずつ数えると3余り、7つずつ数えると2余る。物はいくつあるか。

割った余りだけが分かっているときに、元の数を求める問題です。13世紀には、秦九韶が「数書九章」（1247年）で、この種の問題の一般的な解き方（大衍総数術）をまとめ、暦の計算のような問題に当てはめています。その中の計算の一部が、後で紹介する「大衍求一術」です。その後、ガウスが「整数論（Disquisitiones Arithmeticae）」（1801年）で、整数の性質についての定理として整理しました。

「中国剰余定理」という名前を誰が付けたのかが気になって調べてみましたが、突き止められませんでした。確認できたいちばん古い例は、アメリカの数学者ディクソンの教科書「Introduction to the Theory of Numbers」（1929年）に出てくる英語の「Chinese remainder theorem」です。孫子算経や秦九韶の解き方は、1850年代にイギリスの宣教師アレクサンダー・ワイリーが英語で紹介し、ヨーロッパに知られるようになっていました。確認できる範囲では、名前が登場するのはそれから70年以上後です¹。

時計の針で見る定理の中身

割り算の余りは、時計の針と同じ動きをします。3で割った余りは、目盛りが0、1、2の3つしかない時計の針だと考えられます。数が1増えるたびに針が1目盛り進み、1周すると0に戻ります。

目盛りが3つ、5つ、7つの時計を並べて、0から同時に針を進めていくとします。3つの針が最初にそろって0に戻るのは、105（3×5×7）を数えたときです。それまでの0から104では、3つの針の位置の組が重なることはありません。つまり、針の位置の組（余りの組）を見れば、元の数がひとつに決まります。これが中国剰余定理の中身です。

3、5、7は例にすぎません。どの2つを選んでも1以外に共通の約数がない数の組（互いに素な数の組）なら、法はどんな数でも、何個使ってもかまいません。区別できる数の個数は、法をすべて掛けた数になります。反対に、互いに素でない組では成り立ちません。目盛りが2つと4つの時計では、4を数えた時点で両方の針が0に戻ってしまうので、区別できるのは4個の数だけです。

負の数も扱えます。針を逆に戻すことを負の数と考え、区別できる数の後ろ半分を負の数として読めば、0を挟んだ範囲の数を表せます。

繰り上がりのない計算

この定理が計算に役立つのは、足し算と掛け算を余りのまま進められるからです。

数を足すことは、それぞれの時計の針を、足す数の分だけ進めることにあたります。各時計の針がどこに止まるかは、その時計の今の位置と、進める量だけで決まり、ほかの時計は関係しません。掛け算は足し算の繰り返しなので、同じことが言えます。

そのため、足し算と掛け算だけでできた計算は、大きな元の数を使わずに、時計ごとに余りのまま進められます。筆算では下の桁の繰り上がりを待たないと上の桁が決まりませんが、時計と時計の間にはこの繰り上がりがありません。時計ごとの計算は、お互いを待たずに別々に進められます。

最後に元の数に戻すときは、「1つの時計だけを1目盛り進め、ほかの時計は動かさない数」を時計ごとに用意しておきます。3と5と7なら、70、21、15がそれにあたります。たとえば70は、5と7の時計を動かさないように5×7＝35の倍数から探し、3の時計が1を指す最初の数（35の2倍）として見つかります。21と15も、同じ考え方で見つかります。

この数を、それぞれの針の位置の回数だけ足し合わせます。105以上になった分は、時計が1周した分なので105ずつ引けば、元の数に戻ります。孫子算経の問題なら、針の位置は3の時計が2、5の時計が3、7の時計が2なので、70を2回、21を3回、15を2回足して233です。ここから105を2回引いた23が答えです。

秦九韶の大衍求一術

70のような「1つの時計だけを1目盛り進める数」を、法が大きくても手順どおりに求められるようにしたのが、秦九韶の方法です。数書九章では、次の名前を使って手順を説明しています²。

数書九章での名前	意味	孫子算経の問題での値
定母	法	3、5、7
衍母	法をすべて掛けた数	105
衍数	衍母を自分の法で割った数（ほかの時計を動かさない数）	35、21、15
奇数	衍数を自分の法で割った余り（奇数・偶数の「奇数」ではありません）	2、1、1
乗率	奇数に掛けると、自分の法で割った余りが1になる数	2、1、1
用数	衍数×乗率（1つの時計だけを1目盛り進める数）	70、21、15

大衍求一術は、このうち乗率を求める部分の名前です。田の字の4つのマスに、右上に奇数、右下に定母、左上に1、左下に0を置きます。右の上下で大きいほうを小さいほうで割り、余りに置き換えます。そのたびに、商と反対側の左の数を掛けて、割られた側の左の数に足します。右上が1になったら止め、そのときの左上が乗率です。

70を求める場合で追ってみます。奇数は2、定母は3です。

1. 右下の3を右上の2で割ると、商は1、余りは1です。右下を1にし、左下に「商1×左上の1」を足して1にします。
2. 右上の2を右下の1で割ります。右上に1を残したいので、商を1、余りを1とします。右上を1にし、左上に「商1×左下の1」を足して2にします。
3. 右上が1になったので止めます。左上の2が乗率で、衍数35×2＝70が用数です。

右の列でしているのは、2つの数を、割った余りで置き換え続ける計算です。これはユークリッドの互除法そのものです。左上の数は、奇数をその数倍すると、定母で割った余りが右上の数と同じになるように更新されていきます。右上が1になった時点で、左上は「奇数に掛けると余りが1になる数」、つまり乗率になっています。今の言い方では、拡張ユークリッドの互除法で、法を割った余りの世界での逆数（逆元）を求める計算にあたります。2の手順のように、右上に1が残るように最後の商を調整するところが、「求一」（1を求める）という名前に表れています。

定理が保証すること

中国剰余定理は、数学の定理として証明されています。互いに素な法を何個使っても、法をすべて掛けた数より小さい範囲では、どの余りの組にも元の数がちょうど1つ対応します。そのため、足し算と掛け算を余りのまま計算して元に戻せば、答えの整数を誤差なく求められます。

ただし、証明されているのは整数の計算についてだけです。小数をどう整数にするか、丸めの誤差をどう扱うかについては、定理は何も言っていません。また、答えが区別できる範囲からはみ出すと、時計が1周して別の数と区別できなくなります。12時間の時計で、3時と15時の区別がつかないのと同じです。

暗号や計算機代数での応用

中国剰余定理は、計算機の世界でも広く使われてきました。

• 暗号の分野では、RSA暗号の復号で、大きな数の計算を2つの素数で割った余りに分けて進め、中国剰余定理で戻して速くする方法が使われています。RSAの秘密鍵の標準的な形式には、このための値があらかじめ含まれています³。
• 数を余りの組で表して計算する方式は、剰余数系と呼ばれます。繰り上がりがないので、足し算や掛け算を法ごとに並列に計算できます。1988年には、余りで表した数の行列積を専用の回路で計算する研究も発表されています⁴。
• 計算機代数の分野では、整数の行列を厳密に計算するために、いくつかの素数で割った余りごとに計算し、中国剰余定理で元に戻す方法が使われています。

Ozaki Scheme IIは、この流れを倍精度の行列積に当てはめたものです。

Ozaki Scheme IIのGPU実装

行列積への当てはめ方

行列積の中身は、掛け算と足し算だけです。そのため、中国剰余定理の「余りのまま計算できる」性質を使えます。行列の各要素を余りに置き換えれば、法ごとに、余りの行列どうしの行列積を計算するだけで済みます。法を256以下にしておけば余りは8ビットに収まるので、この行列積は、AI向けのGPUが得意なINT8の行列積になります。

ただし、中国剰余定理が扱えるのは整数だけで、答えが区別できる範囲に収まっている必要もあります。Ozaki Scheme IIは、定理の前後に手順を足して、この2つの条件を満たすようにしています。全体の流れは次のとおりです。

中国剰余定理そのものを使うのは手順4で、手順2と3は「余りのまま計算できる」性質を使う部分です。図のように、手順2と3は法の数だけ横に並び、それぞれがほかの法の結果を待たずに進みます。論文のアルゴリズムでも、各手順を倍精度、整数、高い精度のどれで計算するかが区別されています⁵。

2の累乗による整数化（手順1と5）

倍精度の行列には小数が入っているので、そのままでは余りを取れません。そこで、行ごと、列ごとに2の累乗を掛けて大きな数にし、小数点以下を切り捨てて整数にします。最後の手順5で、同じ2の累乗で割り戻します。

2の累乗を掛けたり割ったりしても、結果が倍精度で表せる範囲に収まっていれば、コンピュータの中では数の指数の部分がずれるだけなので、誤差は出ません。手順1と5で誤差が出るのは、小数点以下を切り捨てるところだけです。

掛ける2の累乗は、大きいほど切り捨てる部分が小さくなり、精度が上がります。ただし、大きくしすぎると、行列積の答えが中国剰余定理で区別できる範囲からはみ出します。論文は、この上限を見積もる方法として2つの版を用意しています。速さを優先する版（OS II-fast）は、行と列それぞれのベクトルの長さ（ノルム）から、コーシー・シュワルツの不等式で手早く上限を見積もります。精度を優先する版（OS II-accu）は、見積もりそのものにINT8の行列積を1回使い、より正確な上限を求めて、使える桁を増やします⁵。

256以下の法の選び方（手順2）

法には、2から256までの数のうち、互いに素で、できるだけ大きい数を選びます。法が大きいほど、1回の行列積で区別できる範囲が広がるからです。論文では、法を16個使う場合に、256、255、253、251、247、239、233、229、227、223、217、211、199、197、193、191を挙げています⁵⁶。

法は素数でなくてもかまいません。中国剰余定理が求めているのは、互いに素であることだけだからです。256（2の8乗）、255（3×5×17）、253（11×23）、247（13×19）、217（7×31）は素数ではありませんが、ほかの法と共通の約数を持たないので使えます。素数だけで16個選ぶと251から167までになり、区別できる数の個数は2の123乗ほどです。素数でない数も使うと2の125乗ほどになり、同じ16回の行列積で、区別できる範囲が約4倍（2ビット分）広がります。

256は、8ビットで表せる値の個数そのものなので、INT8で扱える法の中でいちばん大きい数です。偶数は、約数の2が共通になるので、法の組の中に1つしか入れられません。その1つを2の累乗の256にすると、2を約数に持つ法の中でいちばん大きい数を使えます。254（2×127）や252を選ぶと、256より小さいうえに、ほかの偶数はやはり使えません。また、256で割った余りは、2進数で表したときの下8ビットそのものです。

余りは、0を挟んだ範囲で持ちます。たとえば7で割った余りなら、0から6ではなく-3から3で持ちます。こうすると、256以下のどの法でも、余りは-128から127に収まります。これはINT8でそのまま表せる範囲です。256で割った余りが128になる場合だけは、この範囲を1つ超えます。ただ、INT8では128が-128として扱われ、256の時計の上でも128と-128は同じ位置なので、問題は起きません。

法の組と、手順4で使う「1つの時計だけを進める数」は、法の個数ごとに前もって計算し、表にしておきます。計算のたびに求め直す必要はありません。

誤差のないINT8行列積（手順3）

法ごとの行列積には、cuBLASに以前からあるINT8の行列積の関数を使います。テンソルコアは、8ビットの整数を受け取り、32ビットの整数（INT32）で足し合わせます。

INT8どうしの掛け算の答えは、絶対値でいちばん大きくても128×128（16384）です。このため、足し合わせる数（行列の内側の大きさ）が2の17乗（131072）未満なら、合計がINT32からあふれることはなく、行列積は誤差なしで計算できます。それより大きい行列は、小さなブロックに分けて計算します⁵。

この手順が、GPUとの相性の良さを決めています。

1つ目に、計算の大部分が、ふつうの大きな行列積になります。法1つにつき、INT8の行列積は1回です。テンソルコアは、大きな行列積をまとめて計算するときに、いちばん速さを出せます。

2つ目に、法ごとの行列積はお互いに独立しています。時計と時計の間に繰り上がりがないので、ほかの法の結果を待たずに計算できます。

3つ目に、手順2、4、5は行列の1マスごとに独立した計算で、手順1も行や列ごとに並行して進められます。数万から数億のマスを、GPUのたくさんの演算器で同時に処理できます。

行列が大きくなるほど、行列積の計算量は、ほかの手順の計算量より速く増えます。そのため、行列が大きいほど、全体の時間のうちテンソルコアでの行列積が占める割合が大きくなります。反対に、小さな行列では、ネイティブの倍精度より遅くなる例が論文に示されています⁵。

高い精度での復元（手順4）

法を16個使うと、区別できる数の個数は2の125乗ほどになります。倍精度の数（53ビット）にも、64ビットの整数にも収まらない大きさです。

そこで論文では、元の数に戻す計算だけを、倍精度より高い精度で行います。実装では、倍精度の数を2つ組み合わせて、およそ2倍の桁を持たせる計算方法（疑似的なdouble-double演算）を使っています⁵。この計算も1マスごとに独立しているので、GPUで同時に処理できます。

行列積の回数と性能

法を1つ増やすと、区別できる範囲はその法の大きさの分（191倍〜256倍）、ビット数でいえば8ビット弱広がります。行列積では2つの数を掛けるので、そのうち各要素に回せるのは約4ビットです。

論文によると、法を16個使えば、各要素をおよそ53ビット（倍精度の桁数）の整数にでき、行列積16回で倍精度を模倣できる見込みです。以前の方式であるOzaki Scheme Iは、同じ精度を出すのに28〜35回の行列積が必要でした⁵。法の個数がほぼそのまま行列積の回数になるので、回数を半分ほどに減らせることが、Ozaki Scheme IIの利点です。

論文の精度の実験では、試した条件の入力に対して、法14〜15個でcuBLASのネイティブの倍精度と同程度の精度が出ています。性能の測定では、NVIDIA GH200で、16384行16384列の行列積を法14個で計算したときに80.2 TFLOPSが出ています。cuBLASのネイティブの倍精度（60.9 TFLOPS）を上回る値です⁵。

ただし、論文自身も書いているとおり、実際の精度は入力の行列によって変わります。2の累乗は行の中でいちばん大きい要素に合わせて決めるので、同じ行に大きさのかけ離れた要素が混ざると、小さい要素ほど多くの桁を切り捨てられます。前回の記事で紹介したとおり、cuBLASのエミュレーション（Ozaki Scheme IとIIの両方が対象）では、ADPという仕組みが入力を調べ、エミュレーションでは精度を保てないと判断すると、ネイティブの倍精度に切り替えるとされています⁷。

実装の公開とcuBLASへの搭載

論文の実装は、理研R-CCSのGEMMul8として公開されています。手元のGPUで試す方法（cuBLASのエミュレーションを有効にする環境変数や、GEMMul8の使い方）は、前回の記事の「試すには」の節にまとめています。

NVIDIAのcuBLASでは、CUDA 13.4から、倍精度のエミュレーションでOzaki Scheme Iより速くなる場合に、Ozaki Scheme IIが使われるようになりました。対象はAmpere世代以降のGPUです。B200では、倍精度の行列積で最大175 TFLOPSとされています。次の世代のRubinでは、新しいデータ型（TI16）を使う命令でOzaki Scheme IとIIを動かすと説明されています⁸。cuBLASの内部の実装は公開されていないので、論文の手順とどこまで同じかは分かりません。

エミュレーションで足りるかを判断する観点

ここまでの内容から、自分の計算をエミュレーションに任せてよいかを考えるときの観点をまとめます。

観点	確かめること
値の大きさのばらつき	左側の行列の同じ行、右側の行列の同じ列に、大きさが何桁も違う値が混ざっていないか。2の累乗は行や列ごとに共通なので、混ざっているほど小さい値の桁が切り捨てられ、必要な法の個数が増えます59
必要な精度	本当に倍精度と同じ精度が要るのか、もっと粗くてよいのか。cuBLASでは、使うビット数を固定したり、ADPが決めるビット数を増減したりできます9
行列の大きさ	行列が小さいと、行列積以外の手順の時間が目立ち、ネイティブの倍精度より遅くなることがあります。また、行列の内側の大きさが大きいほど、各要素に回せるビット数は少し減ります5
倍精度への切り替え	ADPがネイティブの倍精度に切り替えた計算は、GPUの倍精度の回路で動きます。倍精度の回路が少ないGPUで、切り替わる入力が多いと、期待した速さは出にくくなります7
InfとNaN	入力にInfやNaNが入る場合、ふつうの倍精度と扱いが異なることがあります。正しく扱う設定にすると、遅くなることがあります9
再現性とメモリ	毎回同じ結果が必要なら、作業用のメモリ（ワークスペース）を前もって渡す必要があります。必要な量は条件によって変わり、最大で8GBです9

迷ったときは、実際のデータでネイティブの倍精度とエミュレーションの両方を計算し、結果の差を比べるのが確実です。

まとめ

中国剰余定理は、「孫子算経」の3と5と7の問題に始まり、秦九韶やガウスを経て、互いに素な法なら何個でも使える定理として証明されました。足し算と掛け算を余りのまま、法ごとに別々に計算できる性質から、暗号や剰余数系、整数の行列計算などに使われてきました。

Ozaki Scheme IIは、この性質を倍精度の行列積に当てはめた方式です。定理の外側にある「小数を整数にする」「区別できる範囲に収める」「大きな数を戻す」といった課題を、2の累乗による誤差のない拡大、256以下の法の選び方、高い精度での復元で解決しています。そのうえで、計算の大部分を、お互いに独立したINT8の行列積にしたことで、GPUのテンソルコアの速さをそのまま生かせるようになりました。

一方で、何個の法で足りるかは入力の行列によって変わります。エミュレーションで足りるかどうかは、扱うデータを知っている側でなければ判断できません。前回の記事の「欲しいのは回路ではなく精度」という話にも、ここでつながります。

割り算の余りという小学校の算数から始まった定理が、1500年以上たって、GPUで倍精度の行列積を支えています。数学は、正しく使えばこれほど強い道具になるのだと、書きながら改めて感じました。

---

お断り
記事内容は個人の見解であり、所属組織の立場や戦略・意見を代表するものではありません。
あくまでエンジニアとしての経験や考えを発信していますので、ご了承ください。

1. 名前の初出は、J. Miller「Earliest Known Uses of Some of the Words of Mathematics (C)」の「CHINESE REMAINDER THEOREM」の項によります。ワイリーによる紹介は、英語版Wikipediaの「Chinese remainder theorem」と「Alexander Wylie (missionary)」によります。ワイリーの「Jottings on the Science of the Chinese」は、1853年に上海で出版されています。 ↩

2. 用語と手順は、数書九章 巻一「大衍類」の記述をまとめた中国語版Wikipediaの「大衍求一术」によります。同ページでは、孫子算経の問題に当てはめた例として、衍数35、21、15、奇数2、1、1、乗率2、1、1、用数70、21、15が示されています。秦九韶の方法全体については、MacTutorの「Qin Jiushao」も参考になります。 ↩

3. RSA暗号の標準であるPKCS #1（RFC 8017）の秘密鍵の形式には、中国剰余定理を使った計算のための値（exponent1、exponent2、coefficient）が含まれています。 ↩

4. C. A. Papachristou, S. Hwang: A systolic array structure for matrix multiplication in the residue number system, Supercomputing, Lecture Notes in Computer Science, vol. 297, Springer, 1988, pp. 716-731. DOI: 10.1007/3-540-18991-2_40。Ozaki Scheme IIの論文の参考文献に挙げられています。 ↩

5. K. Ozaki, Y. Uchino, T. Imamura: Ozaki scheme II: A GEMM-oriented emulation of floating-point matrix multiplication using an integer modular technique, The International Journal of High Performance Computing Applications, 2026. DOI: 10.1177/10943420261467787 / arXiv:2504.08009。法の選び方とINT32での誤差のない計算の条件は3.1節、アルゴリズムはAlgorithm 1、OS II-fastとOS II-accuの違い、法14〜15個での精度の結果、性能の測定は4.1節によります。法16個で約53ビットという見込みは、要素の大きさに大きな差がない場合の数字です。 ↩ ↩² ↩³ ↩⁴ ↩⁵ ↩⁶ ↩⁷ ↩⁸ ↩⁹ ↩¹⁰

6. 論文の一覧には素数の241が入っていません。GEMMul8の実装の表（src/oz2/common/table.hpp）では、INT8を使う場合の法は256、255、253、251、247、241、239、233、229、227、223、217、211、199、197、193の順に並び、17個目以降に191、181、179、173が続きます。この16個なら、区別できる数の個数は2の125.4乗ほどです。 ↩

7. NVIDIA「cuBLAS」ドキュメントの Fixed-Point の節。固定小数点のエミュレーションはOzaki Scheme IとIIに基づいており、既定の動作（Dynamic Mantissa Control）ではADPが必要なビット数を求め、上限を超える場合はネイティブの倍精度に切り替えると説明されています。ADPのアルゴリズムは、A. Schwarz ほか: Guaranteed DGEMM Accuracy While Using Reduced Precision Tensor Cores Through Extensions of the Ozaki Scheme, SCA/HPCAsia '26, pp. 91-101. DOI: 10.1145/3773656.3773670 で説明されています。同論文が評価しているのはOzaki Scheme Iです。 ↩ ↩²

8. NVIDIA「CUDA Toolkit Release Notes」の cuBLAS: Release 13.4。 ↩

9. NVIDIA「cuBLAS」ドキュメントの Fixed-Point の節（各行・各列で2の累乗を共有する表現、Fixed Mantissa Control、Fixed-Point Workspace Requirements、Fixed-Point Performance Guide）と、Support For Floating Point Special Values の節によります。固定小数点のエミュレーションでは、InfとNaNはそのままでは扱えず、扱う設定にすると追加の処理が入るとされています。 ↩ ↩² ↩³ ↩⁴