CPU/GPU高速化セミナー 浮動小数点から文字列への高速変換の論文を読んでみた / cpugpu acceleration seminar 20230201

2023年2月1日(水)に実施した「CPU/GPU高速化セミナー ～浮動小数点から文字列への高速変換の論文を読んでみた～」の当日資料です。

Transcript

1. Copyright© Fixstars Group CPU/GPU高速化セミナー 浮動小数点から文字列への 高速変換の論文を読んでみた
2. Copyright© Fixstars Group 本日のAgenda • はじめに • フィックスターズのご紹介 • Ryū のアルゴリズム解説 • Ryū の GPU 移植 • 最新アルゴリズムの紹介 • Q&A / 告知 2
3. Copyright© Fixstars Group はじめに 3
4. Copyright© Fixstars Group 本講演の位置づけ 4 • 弊社がサービス展開している、ソフトウェアによる処理速度向上に関連した 様々な技術情報を発信しています • 性能モデルの理論と実践 (発表資料) ◦ ソフトウェア性能の理論的根拠となる 性能モデルの使用方法を実践例 • 今回の内容 ◦ Ryū のアルゴリズムについて、理論と実装を紹介 ◦ GPU移植の説明と、CPU との処理速度を比較 • こんな方に向いています ◦ アルゴリズム・高速計算の実装に興味のある方 ◦ システム性能向上に向けた、技術アプローチに興味のある方
5. Copyright© Fixstars Group 発表者紹介 冨田 明彦 ソリューションカンパニー 営業企画 2008年に入社。金融、医療業界において、 ソフトウェア高速化業務に携わる。その後、 新規事業企画、半導体業界の事業を担当し、 現職。 5 大野 真暉 ソリューション第一事業部 シニアエンジニア 2019年に新卒で入社。製造業における スケジューリング問題のソルバー開発に従事。 数理最適化、量子コンピュータ関連のプロジェ クトを手掛ける。
6. Copyright© Fixstars Group フィックスターズの ご紹介 6
7. Copyright© Fixstars Group フィックスターズの強み コンピュータの性能を最大限に引き出す、ソフトウェア高速化のエキスパート集団 ハードウェアの知見 アルゴリズム実装力 各産業・研究分野の知見 7 目的の製品に最適なハードウェアを見抜き、 その性能をフル活用するソフトウェアを開 発します。 ハードウェアの特徴と製品要求仕様に合わ せて、アルゴリズムを改良して高速化を実 現します。 開発したい製品に使える技術を見抜き、実 際に動作する実装までトータルにサポート します。
8. Copyright© Fixstars Group サービス提供分野 8 半導体 自動車 産業機器 生命科学 金融 • NAND型フラッシュメモリ向け ファームウェア開発 • 次世代AIチップの開発環境基盤 • 自動運転の高性能化、実用化 • 次世代パーソナルモビリティの 研究開発 • Smart Factory実現への支援 • マシンビジョンシステムの高速化 • ゲノム解析の高速化 • 医用画像処理の高速化 • AI画像診断システムの研究開発 • デリバティブシステムの高速化 • HFT(アルゴリズムトレード)の高速化
9. Copyright© Fixstars Group サービス領域 様々な領域でソフトウェア高速化サービスを提供しています。大量データの高速処理は、 お客様の製品競争力の源泉となっています。 9 組込み高速化 画像処理・ アルゴリズム開発 分散並列システム開発 GPU向け高速化 FPGAを活用した システム開発 量子コンピューティング AI・深層学習 自動車向け ソフトウェア開発 フラッシュメモリ向け ファームウェア開発
10. Copyright© Fixstars Group コンサルティング セットアップ支援、処理の分割や変換等のコンサル クラウド実行環境のご提供 クラウド経由での量子コンピュータ利用サービスを提供 ソフトウェア高速化・開発支援サービス 量子コンピュータを組み合わせてシステムの高速化を実現 セミナー・トレーニング 量子コンピュータの研究動向や活用事例、実際の利用方法等 量子コンピューティング活用支援 多様なハードウェアでのソフトウェア高速化サービスに加え、 量子コンピュータ活用支援とシステム開発を提供しています。 お客様の課題 ご支援内容 量子コンピューティングが課題の解決に役立 つか、確信が持てない 量子コンピューティングの検討をどう進めて いったら良いかわからない 作りたいアプリケーションがあるが、開発が 難しい 10
11. Copyright© Fixstars Group 物流倉庫における人員配置最適化の社会実装開発 内閣府主催SIP公開シンポジウム2021 11 1 通販向け物流倉庫で誰がどの仕事を担当す るか割り振る量子コンピューティングシス テムを開発 2 30分以上かかっていた割り当て作業が1秒 未満の計算で実現 3 業務負荷の偏りが解消し、全体の作業効率 も30％改善する割り当てパターンが計算で きた 作業者の負荷のばらつきを抑えつつ、全体の作業効率は30%向上 物流 分野 量子コンピューティング活用支援 サービス領域 ベルメゾンロジスコ様のご協力 を得て実施
12. Copyright© Fixstars Group Ryū 12
13. Copyright© Fixstars Group Ryū 2018年のPLDIで発表された論文 “Ryū: fast float-to-string conversion” を 紹介します。 浮動小数点を高速に文字列に変換するアルゴリズムを提案・実装した論文です。 固定幅整数の演算のみで文字列に変換し、正確かつ高速なことが特徴です。 Adams, Ulf. "Ryū: fast float-to-string conversion." Proceedings of the 39th ACM SIGPLAN Conference on Programming Language Design and Implementation. 2018. 13
14. Copyright© Fixstars Group Ryū - 目次 1. はじめに 2. 浮動小数点の復習：IEEE754 3. 文字列化のルール 4. 文字列化アルゴリズム • 理論面 • 実装面 5. 性能評価 6. GPU移植 7. 最新アルゴリズムの紹介 8. 参考資料 14
15. Copyright© Fixstars Group はじめに 浮動小数点の文字列化アルゴリズムはいたるところで用いられる、重要なアル ゴリズムである。 例えば、Javascriptでは基本的な数値型のNumber型をプリントするたびに浮動 小数点を文字列化している。 他にも、JSONフォーマットへ浮動小数点をシリアライズする際など、様々な箇 所で使用される。 15
16. Copyright© Fixstars Group 浮動小数点についての復習：IEEE754（１） IEEE754は浮動小数点数のビット表現（交換形式）を定義している。 ビット表現は三つの要素から構成される。 浮動小数点の形式ごとに、それぞれの要素のビット幅が決まっている。 1. 符号 sign 2. 指数 exponent 3. 仮数 fraction（mantissaとも） 16 Float64のビット表現 [5]より引用 符号のビット幅=1、指数のビット幅=11、仮数のビット幅=52
17. Copyright© Fixstars Group 浮動小数点についての復習：IEEE754（２） 浮動小数点のビット表現と数値の対応関係は下式で定義される。 • 𝑠：符号部の値 • 𝑠 = 0なら正の値、𝑠 = 1なら負の値 • 𝑒：指数部の値 • biasは浮動小数点の形式ごとに決まっている値 • Float32なら126で、Float64なら1022 • 𝑚：仮数部のビット列 • 1.mと0.mは2進数で記載した小数 17 浮動小数点の数値 [1]より引用
18. Copyright© Fixstars Group 文字列化のルール 浮動小数点をどのような文字列に変換するか？ 本発表では、次のルール1, 2に従う、正確かつ短い文字列への変換を考える。 a. 情報の保存：パーサ3は文字列化された値から元々の浮動小数点を復元できる b. 文字列は最短：文字列は可能な限り最短にする c. 正確な丸め：aとbのルールを満たす範囲で文字列は浮動小数点の数値に最も 近くする 1：How to print floating-point numbers accurately Guy L. Steele, Jr. and Jon L. White. 1990. 2：SteeleとWhiteは他のルールも定めているが、この三つのルールのみ用いられているケースが多い。 3：パーサは文字列化の逆向きの変換つまり、文字列から浮動小数点に変換するアルゴリズムのこと。 18
19. Copyright© Fixstars Group アルゴリズム概要 ~1.2を具体例に~（１） 1.2を例に文字列化アルゴリズムの概要を紹介する。 1.2のFloat64のビット表現は次のようになる。 0011111111110011001100110011001100110011001100110011001100110011 このビット表現を定義通り数値に変換すると次の値になる。 𝑓 =1.1999999999999999555910790149937383830547332763671875 この文字列は正確だがあまりに長い。 先ほどのルールの「2. 文字列は最短」に違反する。 19 符号 指数 仮数 浮動小数点の数値 [1]より引用
20. Copyright© Fixstars Group アルゴリズム概要 ~1.2を具体例に~（２） 数値𝒇の「近く」で、桁数の少ない数値を探す。 「近く」を正確に表現するために、下の数直線のように隣の浮動小数点数との境界を考える。 𝑓− ： 𝑓より一つ小さなFloat64 𝑓+ ： 𝑓より一つ大きなFloat64 𝑓の区間の中で、十進法の小数で表現したときに最も桁数の少ない実数が1.2なので、”1.2”と 文字列化する。 この文字列化アルゴリズムは先ほど定めた文字列化のルールに従う。 20 𝑓−の区間 𝑓− 𝑓+の区間 𝑓+ 𝑓の区間 𝑓 1.2 𝑓−と𝑓の中点 1.19999999999999984456877655247808434069156646728515625 𝑓+と𝑓の中点 1.20000000000000006661338147750939242541790008544921875
21. Copyright© Fixstars Group 文字列化アルゴリズム 文字列化アルゴリズムは複雑なので二つのステップに分けて解説する。 STEP1：理論的な話 任意精度の計算を駆使して浮動小数点を文字列に変換する。 STEP2：実装的な話 STEP1のアルゴリズムで任意精度が必要な計算を、全て固定幅整数の計算に置 き換える。 21
22. Copyright© Fixstars Group おことわり 正の浮動小数点の文字列化のみを考える。 つまり、以下の入力は除外する。 • NaN • Infinity • 負の浮動小数点 • 負の値を文字列化するときは、その正の値の文字列の先頭に「-」を付けると良い • +0.0 • アルゴリズムの都合上除外する また、指数部の値が変わる特殊なケースや丸めなどの話はスキップする。 （大事ではあるがアルゴリズムの概要を把握する分には不要） 22
23. Copyright© Fixstars Group 文字列化アルゴリズム：理論面 1. 浮動小数点を数値にデコードする 2. 隣の浮動小数点数との中点を計算する ◦ 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三つ組を 𝑢, 𝑣, 𝑤 × 2𝑒2とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2 は整数） 3. 2の冪を10の冪となるよう計算する ◦ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10 は整数） 4. 主要部の桁数が最も小さく正しく丸められた整数を計算する ◦ 主要部とは、整数をm × 10𝑒と書いた時のmのこととする ◦ 𝑎より大きく𝑐より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜 , 𝑒𝑜 を計算する 5. 整数を文字列としてプリントする ◦ 𝑑𝑜 × 10𝑒𝑜+𝑒10を文字列にする 23
24. Copyright© Fixstars Group STEP1：浮動小数点を数値にデコードする • 概要 • IEEE754の通りにデコードする • 小数部をはらって整数にする • 入力 • 𝑒 指数部 • 𝑚 仮数部 • 出力 • 𝑒𝑓 , 𝑚𝑓 24 浮動小数点の数値 [1]より引用 𝒆𝒇 , 𝒎𝒇 の計算式 [1]より引用 𝑙𝑒𝑛(𝑚)：仮数のビット数 𝑓 = 𝑚𝑓 × 2𝑒𝑓が成立する。
25. Copyright© Fixstars Group STEP1：具体例 1.2 25 0 s 1023 e 900719925474099 m 𝑓 = −1 0 × 252 + 900719925474099 × 21023 −1023 −52 𝑚𝑓 = 5404319552844595 𝑒𝑓 = −52 1. 浮動小数点を数値にデコードする
26. Copyright© Fixstars Group STEP2：隣の浮動小数点数との中点を計算する • 概要 • 一つ小さな浮動小数点数𝑓−との中点を計算する • 一つ大きな浮動小数点数𝑓+との中点を計算する • 隣の浮動小数点数はIEEE754の式とにらめっこすれば計算できる • 入力 • 𝑒, 𝑚 • 𝑒𝑓 , 𝑚𝑓 • 出力 • 𝑒2 • 𝑢, 𝑣, 𝑤 26 𝑢 ⋅ 2𝑒2 𝑓+ 𝑓 𝑤 ⋅ 2𝑒2 𝑣 ⋅ 2𝑒2 𝑓− 三つ組の計算式 𝑒2 = 𝑒𝑓 − 2 𝑢 = 4𝑚𝑓 − ቊ 1 if 𝑚 = 0 … 0 2 and 𝑒 > 0 … 01 2 2 otherwise 𝑣 = 4𝑚𝑓 𝑤 = 4𝑚𝑓 + 2 𝑚𝑓 × 2𝑒𝑓 = 𝑣 × 2𝑒2 が成立する。
27. Copyright© Fixstars Group STEP2：具体例 1.2 27 𝑢 ⋅ 2𝑒2 𝑓+ 𝑓 = 𝑚𝑓 × 2𝑒𝑓 𝑤 ⋅ 2𝑒2 𝑣 ⋅ 2𝑒2 𝑓− 𝑒2 = 𝑒𝑓 − 2 = −54 𝑣 = 4 × 𝑚𝑓 = 21617278211378380 𝑢 = 𝑣 − 2 = 21617278211378378 𝑤 = 𝑣 + 2 = 21617278211378382 2. 隣の浮動小数点との中点を計算する
28. Copyright© Fixstars Group STEP3：2の冪を10の冪となるよう計算する • 概要 • 10進法の数値として文字列化する際に2𝑛は邪魔なので展開する • 入力 • 𝑒2 • 𝑢, 𝑣, 𝑤 • 出力 • 𝑒10 • 𝑎, 𝑏, 𝑐 28 𝑢 × 2𝑒2 𝑓 = 𝑣 × 2𝑒2 𝑤 × 2𝑒2 𝑎 × 10𝑒10 𝑏 × 10𝑒10 𝑐 × 10𝑒10 三つ組（10進数）の計算式 𝑒2 ≥ 0 ⇒ 𝑒10 = 0, 𝑎, 𝑏, 𝑐 = 𝑢, 𝑣, 𝑤 ⋅ 2𝑒2 𝑒2 < 0 ⇒ 𝑒10 = 𝑒2 , 𝑎, 𝑏, 𝑐 = 𝑢, 𝑣, 𝑤 ⋅ 5−𝑒2 𝑎, 𝑏, 𝑐 × 10𝑒10 = 𝑢, 𝑣, 𝑤 × 2𝑒2 が成立する。
29. Copyright© Fixstars Group STEP3：具体例 1.2 29 𝑢 × 2𝑒2 𝑣 × 2𝑒2 𝑤 × 2𝑒2 𝑎 × 10𝑒10 𝑏 × 10𝑒10 𝑐 × 10𝑒10 𝑒2 ≥ 0 ⇒ 𝑒10 = 0, 𝑎, 𝑏, 𝑐 = 𝑢, 𝑣, 𝑤 ⋅ 2𝑒2 𝑒2 < 0 ⇒ 𝑒10 = 𝑒2 , 𝑎, 𝑏, 𝑐 = 𝑢, 𝑣, 𝑤 ⋅ 5−𝑒2 𝑒10 = 𝑒2 = −54 𝑎 = 𝑢 × 554 = 1199999999999999844568776552478084340691566467285156250 𝑏 = 𝑣 × 554 = 1199999999999999955591079014993738383054733276367187500 𝑐 = 𝑤 × 554 = 1200000000000000066613381477509392425417900085449218750 3. 2の冪を10の冪となるよう計算する
30. Copyright© Fixstars Group STEP4：主要部の桁数が最も小さく正しく丸められた整数を計算する • 概要 • compute_shortestで、𝑎より大きく𝑐より小さい整数の中で、主要部の桁数が最も小さい整数を求 める • 数学的に定式化すると、次の不等式を満たす整数𝑒𝑜 , 𝑑𝑜 を求めることに対応する • 𝒂 < 𝒅𝒐 × 𝟏𝟎𝒆𝒐 < 𝒄 • ただし𝑑𝑜 の桁数は最小とする • つまり𝑒𝑜 の値は最大 • 入力 • 𝑎, 𝑏, 𝑐 • 丸めのオプション • 出力 • 𝑒𝑜 , 𝑑𝑜 30 𝑎と𝑐が等しくなるまで 10で除算する compute_shortest [1]より引用
31. Copyright© Fixstars Group STEP4：（資料）正しいcompute_shortest 右の関数が真のcompute_shortestである。 前ページのcompute_shortestには一部嘘 があった。 嘘とは、文字列化のルールの一つに違反し ていること。 (3) 正確な丸め：1と2のルールを満たす文 字列で、文字列は浮動小数点の数値に最も 近くする 隣の浮動小数点数との中点である𝑎と𝑐のみ ならず、𝑏も一緒に計算することでこのルー ルも満たしたアルゴリズムになる。 31 def compute_shortest(a, b, c, accept_smaller, accept_larger, break_tie_down): i=0 a[0]=a, b[0]=b, c[0]=accept_larger ? c : c-1 all_a_zero[0]=True, all_b_zero[0]=True digit[0]=0 while floor(a[i]/10) < floor(c[i]/10): all_a_zero[i+1]=all_a_zero[i] and a[i]%10==0 all_b_zero[i+1]=all_b_zero[i] and digit[i]==0 a[i+1]=floor(a[i]/10), b[i+1]=floor(b[i]/10), c[i+1]=floor(c[i]/10) digit[i+1]=b[i]%10 i=i+1 if accept_smaller and all_a_zero[i]: while a[i]%10==0: all_b_zero[i+1]=all_b_zero[i] and digit[i]==0 a[i+1]=a[i]/10, b[i+1]=floor(b[i]/10), c[i+1]=floor(c[i]/10) digit[i+1]=b[i]%10 i=i+1 is_tie=(digit[i]==5) and all_b_zero[i] wand_round_down=(digit[i]<5) or (is_tie and break_tie_down) round_down=(want_round_down and (a[i]!=b[i] or all_a_zero[i])) or (b[i]+1>c[i]) return (d_o, e_o) = round_down ? (b_i, i) : (b_i+1, i)
32. Copyright© Fixstars Group STEP4：具体例 1.2 32 ⋮ 𝑎 / 1052 = 119 < 𝑐 / 1052 = 120 𝑎 / 1053 = 11 < 𝑐 / 1053 = 12 𝑎 / 1054 = 1 == 𝑐 / 1054 = 1 ⇒ 𝑑𝑜 = 12, 𝑒𝑜 = 53 4. 主要部の桁数が最も小さく正しく丸められた整数を計算する
33. Copyright© Fixstars Group STEP5：整数を文字列としてプリントする • 概要 • 𝑑𝑜 × 10𝑒𝑜+𝑒10を文字列にする • 入力 • 𝑒10 • 𝑒𝑜 , 𝑑𝑜 • 文字列化のオプション • 出力 • 文字列 33
34. Copyright© Fixstars Group STEP5：具体例 1.2 34 𝑑𝑜 × 10𝑒𝑜+𝑒10を文字列にする 𝑑𝑜 = 12 𝑒𝑜 + 𝑒10 = 53 − 54 = −1 ⇒”1.2” 5. 整数を文字列としてプリントする
35. Copyright© Fixstars Group 文字列化アルゴリズム：理論面（再掲） 1. 浮動小数点を数値にデコードする 2. 隣の浮動小数点数との中点を計算する ◦ 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三つ組を 𝑢, 𝑣, 𝑤 × 2𝑒2とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2 は整数） 3. 2の冪を10の冪となるよう計算する ◦ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10 は整数） 4. 主要部の桁数が最も小さく正しく丸められた整数を計算する ◦ 主要部とは、整数をm × 10𝑒と書いた時のmのこととする ◦ 𝑎より大きく𝑐より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜 , 𝑒𝑜 を計算する 5. 整数を文字列としてプリントする ◦ 𝑑𝑜 × 10𝑒𝑜+𝑒10を文字列にする 35
36. Copyright© Fixstars Group 文字列化アルゴリズム：具体例（まとめ） 36 S T E P 1 0 s 1023 e 900719925474099 m 𝑓 = −1 0 × 252 + 900719925474099 × 21023 −1023 −52 𝑚𝑓 = 5404319552844595 𝑒𝑓 = −52 𝑢 ⋅ 2𝑒2 𝑓+ 𝑓 𝑤 ⋅ 2𝑒2 𝑣 ⋅ 2𝑒2 𝑓− 𝑒2 = 𝑒𝑓 − 2 = −54 𝑣 = 4 × 𝑚𝑓 = 21617278211378380 𝑢 = 𝑣 − 2 = 21617278211378378 𝑤 = 𝑣 + 2 = 21617278211378382 S T E P 2
37. Copyright© Fixstars Group 文字列化アルゴリズム：具体例（まとめ） 37 𝑢 × 2𝑒2 𝑣 × 2𝑒2 𝑤 × 2𝑒2 𝑎 × 10𝑒10 𝑏 × 10𝑒10 𝑐 × 10𝑒10 𝑒2 ≥ 0 ⇒ 𝑒10 = 0, 𝑎, 𝑏, 𝑐 = 𝑢, 𝑣, 𝑤 ⋅ 2𝑒2 𝑒2 < 0 ⇒ 𝑒10 = 𝑒2 , 𝑎, 𝑏, 𝑐 = 𝑢, 𝑣, 𝑤 ⋅ 5−𝑒2 𝑒10 = 𝑒2 = −54 𝑎 = 𝑢 × 554 = 1199999999999999844568776552478084340691566467285156250 𝑏 = 𝑣 × 554 = 1199999999999999955591079014993738383054733276367187500 𝑐 = 𝑤 × 554 = 1200000000000000066613381477509392425417900085449218750 S T E P 3 S T E P 4 S T E P 5 ⋮ 𝑎 / 1052 = 119 < 𝑐 / 1052 = 120 𝑎 / 1053 = 11 < 𝑐 / 1053 = 12 𝑎 / 1054 = 1 == 𝑐 / 1054 = 1 ⇒ 𝑑𝑜 = 12, 𝑒𝑜 = 53 𝑑𝑜 × 10𝑒𝑜+𝑒10を文字列にする 𝑑𝑜 = 12 𝑒𝑜 + 𝑒10 = −1 ⇒”1.2”
38. Copyright© Fixstars Group 文字列化アルゴリズム：実装面 先ほどのアルゴリズムのボトルネックは任意精度の演算を行う部分である。 任意精度の演算を行うのはSTEP3とSTEP4である。 • STEP3：2の冪を処理して浮動小数点を「整数×10の冪」と表す • この整数が大きな整数となる • STEP4：二つの大きな整数が等しくなるまで10で除算する 38 1. 浮動小数点を数値にデコードする 2. 隣の浮動小数点数との中点を計算する ◦ 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三 つ組を 𝑢, 𝑣, 𝑤 × 2𝑒2とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2 は整数） 3. 2の冪を10の冪となるよう計算する ◦ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10 は整数） 4. 主要部の桁数が最も小さく正しく丸められた整数を計算する ◦ 主要部とは、整数をm × 10𝑒と書いた時のmのこととする ◦ 𝑎より大きく𝑐より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜 , 𝑒𝑜 を計算する 5. 整数を文字列としてプリントする ◦ 𝑑𝑜 × 10𝑒𝑜+𝑒10を文字列にする 𝑎, 𝑏, 𝑐 がとても 大きな整数
39. Copyright© Fixstars Group 具体例：Float32の𝒂, 𝒃, 𝒄の値 Float32 符号ビット：1、指数ビット：8、仮数ビット：23 bias：127 STEP1 𝑒𝑓 : 0 ~ 224 𝑚𝑓 : −149~104 STEP2 𝑢, 𝑣, 𝑤: 0~ 226 𝑒2 : −151~102 STEP3 𝑎, 𝑏, 𝑐: max 226 × 2102, 226 × 5151 ≒ 10131 39 STEP1 STEP2 STEP3 非常に大きい！
40. Copyright© Fixstars Group 具体例：Float64の𝒂, 𝒃, 𝒄の値 Float32 符号ビット：1、指数ビット：11、仮数ビット：52 bias：1023 STEP1 𝑒𝑓 : 0 ~ 253 𝑚𝑓 : −1074~971 STEP2 𝑢, 𝑣, 𝑤: 0~ 255 𝑒2 : −1076~969 STEP3 𝑎, 𝑏, 𝑐: max 255 × 2969, 255 × 51076 ≒ 10768 40 STEP1 STEP2 STEP3 非常に大きい！
41. Copyright© Fixstars Group 大きな整数の回避（１） STEP3とSTEP4を工夫できないか？ • STEP3：2の冪を処理して浮動小数点数を「整数×10の冪」とあらわす • この整数が大きな整数となる • STEP4：二つの大きな整数が等しくなるまで10で除算する STEP3で「大きな整数」があらわれ、STEP4で「大きな整数を10で除算した小 さな整数」にする。 STEP3とSTEP4を繋げることで、「大きな整数」を計算することなく、直接 「大きな整数を10で除算した小さな整数」を計算することを考える。 41
42. Copyright© Fixstars Group 大きな整数の回避（２） 今後、𝒆𝟐 ≥ 𝟎の場合に限定する。 STEP3とSTEP4で行う計算の概略を右に示す。 数学的な考察から、STEP4のwhileループは少なくとも次の𝑖 = 𝑞まで はまわることが分かる。 𝑞 = max(0, 𝑒2 × log10 2 − 1) 𝑎を10で𝑞回除算した値𝑎𝑞 は𝑢の100倍未満のため、小さな整数であ る： 𝑎𝑞 = 𝑎 / 10𝑞 ≒ 𝑢 × 2𝑒2 / 10𝑞 < 100𝑢 よって、whileループが𝑞回まわった時点の変数の値は小さい。 ⇒STEP3の段階で、STEP4のwhileループが𝒊 = 𝒒までまわった時の 変数の値を計算できれば、大きな整数の出現を回避できる！ 以降、 𝑎を計算せずに、直接𝑎𝑞 を計算する方法について考察する。 42 STEP3 𝑎 = 𝑢 × 2𝑒2, 𝑐 = 𝑤 × 2𝑒2 STEP4 while 𝑎 / 10𝑖 < 𝑐 − 1 / 10𝑖: 𝑖 = 𝑖 + 1 2𝑒2 ≒ 10𝑞 となる𝑞を選んでいる
43. Copyright© Fixstars Group 𝒂を10で𝒒回除算した値の計算 compute_shortestに登場する変数は色々あるが、 𝑎を10で𝑞回除算した値𝑎𝑞 の計算が最も困難 なので、 𝑎𝑞 の計算方法のみ解説する。 𝑎𝑞 を計算しやすい形に式変形する。 大きな値に赤枠を付ける。 非整数になりうる値に青枠を付ける。 𝑎𝑞 = 𝑢 × 2𝑒2 / 10𝑞 = 𝑢 × 2𝑒2−𝑞 / 5𝑞 （アイディア：大きな整数2𝑘で5𝑞の除算をなんとかする） = 𝑢 × 2𝑒2−𝑞−𝑘 × 2𝑘 / 5𝑞 （アイディア： 𝑘が十分大きいと次の式が成立する） = 𝑢 × 2𝑒2−𝑞−𝑘 × 2𝑘 / 5𝑞 + 1 43 式変形☆ 次スライド説明する
44. Copyright© Fixstars Group 𝒌の下界の算出：式変形☆の説明 𝑘の値が左式の不等式を満たすと前スライ ドの式変形☆が成立する。 証明は省略する。 青枠のことを「𝑘の下界」と呼ぶ。 赤枠の計算は難しいが、 𝑒2 ごとに予め計 算しておく量で、実行時には必要ない。 赤枠の計算方法の説明は省略する。 44 𝒌の下界の算出 [1]より引用
45. Copyright© Fixstars Group 変形した数式について ☆で式変形した下式は以下の点で優れている。 実行時変数は𝑢, 𝑒2 , 𝑞, 𝑘である。 • 𝑎𝑞 は三つの項の積である • 赤枠（ 2𝑒2−𝑞−𝑘 ）部分は𝑒2 − 𝑞 − 𝑘 < 0のため除算になるが、シフト演算で効率的に実装できる • 困難な5𝑞による除算が 2𝑘 / 5𝑞 に吸収されている • Float32やFloat64の場合、 𝑞の値はそれほど大きくない • 𝑞ごとに 2𝑘 / 5𝑞 の値をあらかじめ計算してテーブル化すれば、実行時には5𝑞による除算がなくなる 赤枠と緑枠の項が大きな値とならないような𝑘が存在すれば、効率的に𝑎𝑞 を計算できる。 ⇒2𝑘 ≒ 5𝑞となる𝑘 = log2 5𝑞付近で、下界より大きな𝑘 は存在するか？ 45 𝑎𝑞 = 𝑢 × 2𝑒2−𝑞−𝑘 × 2𝑘 / 5𝑞 + 1
46. Copyright© Fixstars Group 丁度良い𝒌の探索 • Float64で𝑘の下界を計算した（右図） • 横軸は𝑒2 • 𝑘の下界を青色でプロット • 2𝑘 / 5𝑞 の値をテーブルとして保持するのに必要な要素数 の最小値を橙色でプロット • おおよそ 𝑘 − log2 5𝑞 に対応する • 橙色はプラトーになっている！ • 最大値を𝐵0 とする • 𝑘 = 𝐵0 + log2 5𝑞 とすれば、𝑘は下界よりも大きな値と なる • さらに、実際に計算すると分かるが、この𝑘は丁度良い𝑘で ある！（次の二ページで計算する） • 2𝑘と5𝑞は大きな整数だが、 𝑞ごとに緑枠の値をあらかじめ 計算しておけば、実行時はテーブルを参照するだけで良い 46 𝑘の下界 [1]より引用
47. Copyright© Fixstars Group Float32の場合 • 変数の整理 • 𝑒2 : 0 ∼ 102, 𝑢: 0 ∼ 226 • 𝑞 = max(0, 𝑒2 × log10 2 − 1) • 𝑘 = 𝐵0 + log2 5𝑞 • 𝐵0 = 60 • 𝑞 + 𝑘 − 𝑒2 ≒ 𝐵0 + 𝑞 × 1 + log2 5 − 𝑒2 ≒ 𝐵0 + 𝑒2 × log10 2 − 1 × log2 10 − 𝑒2 = 𝐵0 + 𝑒2 − log2 10 − 𝑒2 ≒ 𝐵0 = 60 • 2𝑘 / 5𝑞 ≒ 2𝐵0+log2 5𝑞 / 5𝑞 ≒ 2𝐵0 = 260 • 赤枠：60ビット右シフト、緑枠：uint64_tの範囲内 • u（uint32_t）×緑枠（uint64_t）の結果を60ビット右シフトすると𝒂𝒒 （uint32_t） が得られる • uint128_tの範囲で計算できる 47
48. Copyright© Fixstars Group Float64の場合 • 変数の整理 • 𝑒2 : 0 ∼ 969, 𝑢: 0 ∼ 255 • 𝑞 = max(0, 𝑒2 × log10 2 − 1) • 𝑘 = 𝐵0 + log2 5𝑞 • 𝐵0 = 124 • 𝑞 + 𝑘 − 𝑒2 ≒ 𝐵0 + 𝑞 × 1 + log2 5 − 𝑒2 ≒ 𝐵0 + 𝑒2 × log10 2 − 1 × log2 10 − 𝑒2 = 𝐵0 + 𝑒2 − log2 10 − 𝑒2 ≒ 𝐵0 = 124 • 2𝑘 / 5𝑞 ≒ 2𝐵0+log2 5𝑞 / 5𝑞 ≒ 2𝐵0 = 2124 • 赤枠：124ビット右シフト、緑枠：uint128_tの範囲内 • u（uint64_t）×緑枠（uint128_t）の結果を124ビット右シフトすると𝒂𝒒 （uint64_t） が得られる • uint256_tの範囲で計算できる 48
49. Copyright© Fixstars Group 大きな整数の回避：まとめ • STEP3の計算を工夫することで、STEP4のwhileループで𝑖 = 𝑞の時点までス キップし、大きな整数の出現を回避できる • 工夫した計算はFloat32の場合、uint128_tで行える • Float64の場合、uint256_tで行える • 𝑞ごとに下式の緑枠の値を予め計算しておくことで、文字列化の際はテーブル 参照するだけで良い 49
50. Copyright© Fixstars Group 文字列化アルゴリズム：実装面（まとめ） 1. 浮動小数点数を数値にデコードする 2. 隣の浮動小数点数との中点を計算する ◦ 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三つ組を 𝑢, 𝑣, 𝑤 × 2𝑒2とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2 は整数） 3. 2の冪を10の冪となるよう計算する ◦ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10 は整数） ◦ 三つ組を10でq回除算した値(𝑎𝑞 , 𝑏𝑞 , 𝑐𝑞 )を計算する（𝑒10 も計算する） 4. 主要部の桁数が最も小さく正しく丸められた整数を計算する ◦ 主要部とは、整数をm × 10𝑒と書いた時のmのこととする ◦ 𝑎𝑞 より大きく𝑐𝑞 より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜 , 𝑒𝑜 を計算する 5. 整数を文字列としてプリントする ◦ 𝑑𝑜 × 10𝑒𝑜+𝑞+𝑒10を文字列にする 50
51. Copyright© Fixstars Group 性能評価： Ryū 浮動小数点を乱数で生成し、文字列化に要する時 間を比較した。 • 横軸：浮動小数点のビット表現を符号なし整 数とみなした時の値 • 縦軸：文字列化に要する時間[ns] 横軸の値がFloat32の場合は231の時、Float64の 場合は 263 の時、浮動小数点の符号が入れ替わる。 c, dのグラフが線対称なのはこのためである。 いずれの場合もRyūが最速だった。 51 浮動小数点を文字列化するのに要する時間の比較 [1]より引用 上段：C、下段：Java 左側：Float32、右側：Float64 実験環境 CPU：Intel(R) Core(TM) i7-4770K at 3.50GHz OS：Ubuntu Linux 17.10. C：clang 3.9 with the -O2 flag Java：OpenJDK 1.8.0_131.
52. Copyright© Fixstars Group GPU移植 RyūをGPUに移植し、大量の浮動小数点を文字列化するタスクの高速化を行った。 浮動小数点の配列と区切り文字列を受け取り、浮動小数点の文字列を区切り文 字でつなげた文字列を出力する関数を実装した。 void toString(const std::vector<T>& v, const std::string& sep, std::string* s); T = float or double 例えばv={1.2, 3.14, -0.7}, sep=”, ”のとき、”1.2, 3.14, -0.7”と文字列化する。 52
53. Copyright© Fixstars Group GPU移植：実装（１） GPUでのtoStringを３つの関数に分けて実装した。 1. f2dKernel, d2dKernel ◦ 浮動小数点を正しく丸めた10進数表現を計算する ▪ 文字列化アルゴリズムのSTEP1～STEP4に対応する ◦ 文字列化したときの文字数も計算する 2. thrust::exclusive_scan ◦ prefix sumで、文字列化した浮動小数点を書き込むインデックスを計算する 3. d2sKernel ◦ 10進数表現を文字列化し、前ステップで計算したインデックスに文字列を書き込む ▪ 文字列化アルゴリズムのSTEP5に対応する 53 float to decimal double to decimal という意味
54. Copyright© Fixstars Group GPU移植：実装（２） 54 1.2 3.14 -0.7 0.006 … v 12 314 -7 6 … f2dKernel, d2dKernel GPUの各スレッドが配列の各要素を処理する -1 -2 -1 -3 … d e 3 4 4 5 … size 0 3 7 11 16 … index thrust::exclusive_scan size配列のprefix sumを計算する 1 . 2 3 . 1 4 - 0 . 7 0 . 0 0 6 … s d2sKernel GPUの各スレッドが配列の各要素を文字列化し、 indexで指定されたインデックスに書き込む 区切り文字列が空文字列の場合
55. Copyright© Fixstars Group GPU移植：性能評価（１） 3つのアルゴリズムを比較した。 1. Ryū（CPU） 2. Ryū（GPU） 3. yyjson（CPU） ◦ Cで実装された、高速なJSONライブラリ ◦ 浮動小数点の文字列化はSchubfachベースのアルゴリズムで行う ◦ ※浮動小数点配列を文字列化した結果はRyūと異なる ◦ ※フェアな比較対象ではないが、参考として比較した Ryū（CPU）とyyjson（CPU）はCPUのシングルスレッドで動作するプログラムであ る。 55
56. Copyright© Fixstars Group GPU移植：性能評価（２） 実験環境 ◦ CPU：Intel(R) Xeon(R) CPU E5-2698 v4 @ 2.20GHz ◦ GPU：Tesla V100 計測区間 1. Ryū（CPU）：配列を文字列化する時間 2. Ryū（GPU）：GPU上の配列をGPU上で文字列化する時間 3. yyjson（CPU）：文字列のマロック＋JSONオブジェクトを文字列化する時 間 ◦ yyjson_mut_writeの時間 56
57. Copyright© Fixstars Group GPU移植：性能評価（３） • Ryū（CPU）は浮動小数点を一要素当たり、floatは43[ns]で、doubleは60[ns]で文字列化した • おおよそ論文通りのオーダーだった • Ryū（GPU）はCPU版よりfloatで200倍ほど、doubleで100倍ほど高速となった • 実行時間の8割はd2sKernelで、これを高速化すればもっと速くなる • yyjson（CPU）はRyū（CPU）よりも高速だった • 本実験でも、SchubfachはRyūよりも高速だった 57 配列の要素数 Ryū（CPU） Ryū（GPU） 1,000,000 46.2 0.513 10,000,000 433 2.42 100,000,000 4,335 22.2 配列の要素数 Ryū（CPU） Ryū（GPU） yyjson（CPU） 1,000,000 60.3 0.795 53.7 10,000,000 604 5.60 533 100,000,000 6,037 54.1 5,280 floatの配列を文字列化 するのに要する時間 [ms] yyjson（CPU）は浮動小 数点を全てdoubleで扱う ため除外 doubleの配列を文字列化 するのに要する時間 [ms]
58. Copyright© Fixstars Group 最新アルゴリズムの紹介：概要 Dragonboxを最新のアルゴリズムとして、その周辺の次の4つのアルゴリズムを 簡単に紹介する。 • Schubfach • Ryū • Grisu-Exact • Dragonbox 58
59. Copyright© Fixstars Group 最新アルゴリズムの紹介：比較 59 初出 時期 著者 アルゴリズム 名前の由来 Schubfach OpenJDKのcore- libs-devのパッチ 2018年4月, 11 月 Raffaello Giulietti Ryūのほぼ上位互換 Ryūのcompute_shortest をループ無しで行う 鳩ノ巣原理のドイツ語の Schubfachprinzipに由来する Ryū PLDI 2018年6月 Ulf Adams 本発表で解説した Grisuの別名のDragonを日本 語に訳した Grisu-Exact GitHub 2020年 jk-jeon Ryū + Grisu Dragonbox GitHub 2020年 jk-jeon Schubfach + Grisu Grisu⇒Dragon Schubfach⇒box (鳩ノ巣原理＝箱入れ原理) ※Grisu-ExactとDragonboxのアルゴリズムの「+ Grisu」について RyūとSchubfachは三つ組で計算する。 Float64の場合、64bit×128bitを三回も計算する必要があり、高コストである。 Grisuのアイディアを利用して、三つ組の計算を省略し、一度の64bit×128bitの計算で済ませるよ うにした。
60. Copyright© Fixstars Group 最新アルゴリズムの紹介：性能評価 • 文字列化に要する時間を4つのアルゴリズムで比較した ([6]より引用) • 上段：Float32、下段：Float64 • 左側：最短文字列がN桁の浮動小数点を文字列化するのに要する時間（Nが横軸） • 右側：浮動小数点をビットとして一様に生成し、文字列化に要する時間をプロット • 結果 • Dragonboxが最速だった • 最短文字列の桁数が小さいほど、Ryūは他のアルゴリズムよりも遅かった • compute_shortest でのナイーブなループが遅い 60 実験環境 CPU：Intel (R) Core™ i7- 7700HQ CPU @2.80GHz C++：clang-cl compiler shipped with Visual Studio 2022 17.0.4 凡例を拡大
61. Copyright© Fixstars Group 参考資料 • [1] https://dl.acm.org/doi/10.1145/3192366.3192369 • Ryūの論文 • [2] https://github.com/ulfjack/ryu • Ryū著者による実装 • [3] https://github.com/abolz/Drachennest • 浮動小数点をminimal length decimalに変換するアルゴリズムの実行時間を比較する • 文字列としてプリントするのに要する時間を除いて計測し、比較している • [4] https://kubo39.hatenablog.com/entry/2022/06/09/浮動小数点数の10進表記についてまとめてみ た • 各種文字列化アルゴリズムがどのライブラリで使用されているか調べている • [5] https://ja.wikipedia.org/wiki/IEEE_754 • Wikipedia IEEE754 • [6] https://github.com/jk-jeon/dragonbox/blob/master/other_files/Dragonbox.pdf • Dragonboxのアルゴリズムを説明したpdf • [7] https://github.com/ibireme/yyjson • yyjson 61
62. Copyright© Fixstars Group Thank you! お問い合わせ窓口 : hr-seminar@fixstars.com 62