CPU/GPU高速化セミナー ～高速文字列変換アルゴリズムの論文を読んでみた～（2023/02/01）

CPU/GPU高速化セミナー 浮動小数点から文字列への 高速変換の論文を読んでみた Copyright© Fixstars Group

本日のAgenda ● はじめに ● フィックスターズのご紹介 ● Ryū のアルゴリズム解説 ● Ryū の GPU 移植 ● 最新アルゴリズムの紹介 ● Q&A / 告知 Copyright© Fixstars Group 2

はじめに Copyright© Fixstars Group 3

本講演の位置づけ ● 弊社がサービス展開している、ソフトウェアによる処理速度向上に関連した 様々な技術情報を発信しています ● 性能モデルの理論と実践 (発表資料) ○ ソフトウェア性能の理論的根拠となる 性能モデルの使用方法を実践例 ● 今回の内容 ○ ○ Ryū のアルゴリズムについて、理論と実装を紹介 GPU移植の説明と、CPU との処理速度を比較 ● こんな方に向いています ○ ○ アルゴリズム・高速計算の実装に興味のある方 システム性能向上に向けた、技術アプローチに興味のある方 Copyright© Fixstars Group 4

• https://speakerdeck.com/fixstars/gpugao-su-hua-semina-xing-neng-moderufalseli-lun-toshi-jian-li-lun-bian

発表者紹介 冨田 明彦 大野 真暉 ソリューションカンパニー 営業企画 ソリューション第一事業部 シニアエンジニア 2008年に入社。金融、医療業界において、 ソフトウェア高速化業務に携わる。その後、 新規事業企画、半導体業界の事業を担当し、 現職。 2019年に新卒で入社。製造業における スケジューリング問題のソルバー開発に従事。 数理最適化、量子コンピュータ関連のプロジェ クトを手掛ける。 Copyright© Fixstars Group 5

フィックスターズの ご紹介 Copyright© Fixstars Group 6

フィックスターズの強み コンピュータの性能を最大限に引き出す、ソフトウェア高速化のエキスパート集団 ハードウェアの知見 アルゴリズム実装力 各産業・研究分野の知見 目的の製品に最適なハードウェアを見抜き、 その性能をフル活用するソフトウェアを開 発します。 ハードウェアの特徴と製品要求仕様に合わ せて、アルゴリズムを改良して高速化を実 現します。 開発したい製品に使える技術を見抜き、実 際に動作する実装までトータルにサポート します。 Copyright© Fixstars Group 7

サービス提供分野 半導体 産業機器 金融 自動車 ● NAND型フラッシュメモリ向け ファームウェア開発 ● 次世代AIチップの開発環境基盤 生命科学 ● Smart Factory実現への支援 ● マシンビジョンシステムの高速化 ● 自動運転の高性能化、実用化 ● ゲノム解析の高速化 ● 次世代パーソナルモビリティの 研究開発 ● 医用画像処理の高速化 Copyright© Fixstars Group ● デリバティブシステムの高速化 ● HFT(アルゴリズムトレード)の高速化 ● AI画像診断システムの研究開発 8

サービス領域 様々な領域でソフトウェア高速化サービスを提供しています。大量データの高速処理は、 お客様の製品競争力の源泉となっています。 組込み高速化 GPU向け高速化 AI・深層学習 画像処理・ アルゴリズム開発 FPGAを活用した システム開発 分散並列システム開発 量子コンピューティング 自動車向け フラッシュメモリ向け ソフトウェア開発 ファームウェア開発 Copyright© Fixstars Group 9

量子コンピューティング活用支援 多様なハードウェアでのソフトウェア高速化サービスに加え、 量子コンピュータ活用支援とシステム開発を提供しています。 お客様の課題 ご支援内容 量子コンピューティングが課題の解決に役立 つか、確信が持てない セミナー・トレーニング 量子コンピュータの研究動向や活用事例、実際の利用方法等 量子コンピューティングの検討をどう進めて いったら良いかわからない コンサルティング セットアップ支援、処理の分割や変換等のコンサル 作りたいアプリケーションがあるが、開発が 難しい クラウド実行環境のご提供 クラウド経由での量子コンピュータ利用サービスを提供 ソフトウェア高速化・開発支援サービス 量子コンピュータを組み合わせてシステムの高速化を実現 Copyright© Fixstars Group 10

物流倉庫における人員配置最適化の社会実装開発 内閣府主催SIP公開シンポジウム2021 分野 1 物流 サービス領域 量子コンピューティング活用支援 通販向け物流倉庫で誰がどの仕事を担当す るか割り振る量子コンピューティングシス ベルメゾンロジスコ様のご協力 を得て実施 テムを開発 2 30分以上かかっていた割り当て作業が1秒 3 業務負荷の偏りが解消し、全体の作業効率 未満の計算で実現 も30％改善する割り当てパターンが計算で きた Copyright© Fixstars Group 作業者の負荷のばらつきを抑えつつ、全体の作業効率は30%向上 11

Ryū Copyright© Fixstars Group 12

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. Copyright© Fixstars Group 13

Ryū - 目次 1. はじめに 8. 参考資料 2. 浮動小数点の復習：IEEE754 3. 文字列化のルール 4. 文字列化アルゴリズム • 理論面 • 実装面 5. 性能評価 6. GPU移植 7. 最新アルゴリズムの紹介 Copyright© Fixstars Group 14

はじめに 浮動小数点の文字列化アルゴリズムはいたるところで用いられる、重要なアル ゴリズムである。 例えば、Javascriptでは基本的な数値型のNumber型をプリントするたびに浮動 小数点を文字列化している。 他にも、JSONフォーマットへ浮動小数点をシリアライズする際など、様々な箇 所で使用される。 Copyright© Fixstars Group 15

浮動小数点についての復習：IEEE754（１） IEEE754は浮動小数点数のビット表現（交換形式）を定義している。 ビット表現は三つの要素から構成される。 浮動小数点の形式ごとに、それぞれの要素のビット幅が決まっている。 1. 符号 sign 2. 指数 exponent 3. 仮数 fraction（mantissaとも） Float64のビット表現 [5]より引用 符号のビット幅=1、指数のビット幅=11、仮数のビット幅=52 Copyright© Fixstars Group 16

浮動小数点についての復習：IEEE754（２） 浮動小数点のビット表現と数値の対応関係は下式で定義される。 浮動小数点の数値 [1]より引用 • 𝑠：符号部の値 • 𝑠 = 0なら正の値、𝑠 = 1なら負の値 • • biasは浮動小数点の形式ごとに決まっている値 Float32なら126で、Float64なら1022 • 1.mと0.mは2進数で記載した小数 • 𝑒：指数部の値 • 𝑚：仮数部のビット列 Copyright© Fixstars Group 17

文字列化のルール 浮動小数点をどのような文字列に変換するか？ 本発表では、次のルール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：パーサは文字列化の逆向きの変換つまり、文字列から浮動小数点に変換するアルゴリズムのこと。 Copyright© Fixstars Group 18

アルゴリズム概要 ~1.2を具体例に~（１） 1.2を例に文字列化アルゴリズムの概要を紹介する。 1.2のFloat64のビット表現は次のようになる。 符号 指数 仮数 0011111111110011001100110011001100110011001100110011001100110011 このビット表現を定義通り数値に変換すると次の値になる。 𝑓 =1.1999999999999999555910790149937383830547332763671875 この文字列は正確だがあまりに長い。 先ほどのルールの「2. 文字列は最短」に違反する。 Copyright© Fixstars Group 浮動小数点の数値 [1]より引用 19

アルゴリズム概要 ~1.2を具体例に~（２） 数値𝒇の「近く」で、桁数の少ない数値を探す。 「近く」を正確に表現するために、下の数直線のように隣の浮動小数点数との境界を考える。 𝑓 − ： 𝑓より一つ小さなFloat64 𝑓 + ： 𝑓より一つ大きなFloat64 𝑓の区間の中で、十進法の小数で表現したときに最も桁数の少ない実数が1.2なので、”1.2”と 文字列化する。 この文字列化アルゴリズムは先ほど定めた文字列化のルールに従う。 𝑓 − の区間 𝑓− 𝑓 + の区間 𝑓+ 𝑓の区間 𝑓 1.2 𝑓 − と𝑓の中点 1.19999999999999984456877655247808434069156646728515625 𝑓 + と𝑓の中点 1.20000000000000006661338147750939242541790008544921875 Copyright© Fixstars Group 20

文字列化アルゴリズム 文字列化アルゴリズムは複雑なので二つのステップに分けて解説する。 STEP1：理論的な話 任意精度の計算を駆使して浮動小数点を文字列に変換する。 STEP2：実装的な話 STEP1のアルゴリズムで任意精度が必要な計算を、全て固定幅整数の計算に置 き換える。 Copyright© Fixstars Group 21

おことわり 正の浮動小数点の文字列化のみを考える。 つまり、以下の入力は除外する。 • NaN • Infinity • 負の浮動小数点 • 負の値を文字列化するときは、その正の値の文字列の先頭に「-」を付けると良い • アルゴリズムの都合上除外する • +0.0 また、指数部の値が変わる特殊なケースや丸めなどの話はスキップする。 （大事ではあるがアルゴリズムの概要を把握する分には不要） Copyright© Fixstars Group 22

文字列化アルゴリズム：理論面 1. 浮動小数点を数値にデコードする 2. 隣の浮動小数点数との中点を計算する ○ 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三つ組を 𝑢, 𝑣, 𝑤 × 2𝑒2 とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2 は整数） 3. 2の冪を10の冪となるよう計算する ○ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10 の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10 は整数） 4. 主要部の桁数が最も小さく正しく丸められた整数を計算する ○ ○ 主要部とは、整数をm × 10𝑒 と書いた時のmのこととする 𝑎より大きく𝑐より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜 , 𝑒𝑜 を計算する 5. 整数を文字列としてプリントする ○ 𝑑𝑜 × 10𝑒𝑜+𝑒10 を文字列にする Copyright© Fixstars Group 23

STEP1：浮動小数点を数値にデコードする • 概要 浮動小数点の数値 [1]より引用 • IEEE754の通りにデコードする • 小数部をはらって整数にする • 𝑒 指数部 • 𝑚 仮数部 • 𝑒𝑓 , 𝑚𝑓 • 入力 • 出力 𝒆𝒇 , 𝒎𝒇 の計算式 [1]より引用 𝑙𝑒𝑛(𝑚)：仮数のビット数 𝑓 = 𝑚𝑓 × 2𝑒𝑓 が成立する。 Copyright© Fixstars Group 24

STEP1：具体例 1.2 1. 浮動小数点を数値にデコードする 0 1023 900719925474099 s e m 0 52 𝑓 = −1 × 2 + 900719925474099 × 21023 −1023 −52 𝑚𝑓 = 5404319552844595 Copyright© Fixstars Group 𝑒𝑓 = −52 25

STEP2：隣の浮動小数点数との中点を計算する • 概要 𝑓+ 𝑓 𝑓− • 一つ小さな浮動小数点数𝑓 − との中点を計算する • 一つ大きな浮動小数点数𝑓 + との中点を計算する • 隣の浮動小数点数はIEEE754の式とにらめっこすれば計算できる • 𝑒, 𝑚 • 𝑒𝑓 , 𝑚𝑓 𝑢 ⋅ 2𝑒2 𝑣 ⋅ 2𝑒2 𝑤 ⋅ 2𝑒2 • 入力 • 出力 • 𝑒2 • 𝑢, 𝑣, 𝑤 三つ組の計算式 𝑒2 = 𝑒𝑓 − 2 1 if 𝑚 = 0 … 0 2 and 𝑒 > 0 … 01 𝑢 = 4𝑚𝑓 − ቊ 2 otherwise 𝑣 = 4𝑚𝑓 𝑤 = 4𝑚𝑓 + 2 2 𝑚𝑓 × 2𝑒𝑓 = 𝑣 × 2𝑒2 が成立する。 Copyright© Fixstars Group 26

STEP2：具体例 1.2 2. 隣の浮動小数点との中点を計算する 𝑓− 𝑓 = 𝑚𝑓 × 2𝑒𝑓 𝑢 ⋅ 2𝑒2 𝑣 ⋅ 2𝑒2 𝑓+ 𝑤 ⋅ 2𝑒2 Copyright© Fixstars Group 𝑒2 = 𝑒𝑓 − 2 = −54 𝑣 = 4 × 𝑚𝑓 = 21617278211378380 𝑢 = 𝑣 − 2 = 21617278211378378 𝑤 = 𝑣 + 2 = 21617278211378382 27

STEP3：2の冪を10の冪となるよう計算する • 概要 • 10進法の数値として文字列化する際に2𝑛 は邪魔なので展開する • 𝑒2 • 𝑢, 𝑣, 𝑤 • 𝑒10 • 𝑎, 𝑏, 𝑐 • 入力 • 出力 𝑢 × 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 が成立する。 Copyright© Fixstars Group 28

STEP3：具体例 1.2 3. 2の冪を10の冪となるよう計算する 𝑢 × 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 Copyright© Fixstars Group 29

STEP4：主要部の桁数が最も小さく正しく丸められた整数を計算する • 概要 • • • 入力 compute_shortestで、𝑎より大きく𝑐より小さい整数の中で、主要部の桁数が最も小さい整数を求 める 数学的に定式化すると、次の不等式を満たす整数𝑒𝑜 , 𝑑𝑜 を求めることに対応する • 𝒂 < 𝒅𝒐 × 𝟏𝟎𝒆𝒐 < 𝒄 • ただし𝑑𝑜 の桁数は最小とする compute_shortest [1]より引用 • つまり𝑒𝑜 の値は最大 • • 𝑎, 𝑏, 𝑐 丸めのオプション • 𝑒𝑜 , 𝑑𝑜 • 出力 𝑎と𝑐が等しくなるまで 10で除算する Copyright© Fixstars Group 30

STEP4：（資料）正しいcompute_shortest 右の関数が真のcompute_shortestである。 前ページのcompute_shortestには一部嘘 があった。 嘘とは、文字列化のルールの一つに違反し ていること。 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 (3) 正確な丸め：1と2のルールを満たす文 字列で、文字列は浮動小数点の数値に最も 近くする 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) Copyright© Fixstars Group 31

STEP4：具体例 1.2 4. 主要部の桁数が最も小さく正しく丸められた整数を計算する ⋮ 𝑎 / 1052 = 119 < 𝑐 / 1052 = 120 𝑎 / 1053 = 11 < 𝑐 / 1053 = 12 𝑎 / 1054 = 1 == 𝑐 / 1054 = 1 ⇒ 𝑑𝑜 = 12, 𝑒𝑜 = 53 Copyright© Fixstars Group 32

STEP5：整数を文字列としてプリントする • 概要 • 𝑑𝑜 × 10𝑒𝑜+𝑒10 を文字列にする • 𝑒10 • 𝑒𝑜 , 𝑑𝑜 • 文字列化のオプション • 文字列 • 入力 • 出力 Copyright© Fixstars Group 33

STEP5：具体例 1.2 5. 整数を文字列としてプリントする 𝑑𝑜 × 10𝑒𝑜+𝑒10 を文字列にする 𝑑𝑜 = 12 𝑒𝑜 + 𝑒10 = 53 − 54 = −1 ⇒”1.2” Copyright© Fixstars Group 34

文字列化アルゴリズム：理論面（再掲） 1. 浮動小数点を数値にデコードする 2. 隣の浮動小数点数との中点を計算する ○ 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三つ組を 𝑢, 𝑣, 𝑤 × 2𝑒2 とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2 は整数） 3. 2の冪を10の冪となるよう計算する ○ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10 の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10 は整数） 4. 主要部の桁数が最も小さく正しく丸められた整数を計算する ○ ○ 主要部とは、整数をm × 10𝑒 と書いた時のmのこととする 𝑎より大きく𝑐より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜 , 𝑒𝑜 を計算する 5. 整数を文字列としてプリントする ○ 𝑑𝑜 × 10𝑒𝑜+𝑒10 を文字列にする Copyright© Fixstars Group 35

文字列化アルゴリズム：具体例（まとめ） S T E P 1 0 1023 900719925474099 s e m 𝑓 = −1 0 × 252 + 900719925474099 × 21023 −1023 −52 𝑚𝑓 = 5404319552844595 S T E P 2 𝑓+ 𝑓 𝑓− 𝑢 ⋅ 2𝑒2 𝑣 ⋅ 2𝑒2 𝑤 ⋅ 2𝑒2 Copyright© Fixstars Group 𝑒𝑓 = −52 𝑒2 = 𝑒𝑓 − 2 = −54 𝑣 = 4 × 𝑚𝑓 = 21617278211378380 𝑢 = 𝑣 − 2 = 21617278211378378 𝑤 = 𝑣 + 2 = 21617278211378382 36

文字列化アルゴリズム：具体例（まとめ） S T E P 3 S T E P 4 S T E P 5 𝑢 × 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 ⋮ 𝑎 / 10 = 119 < 𝑐 / 1052 = 120 𝑎 / 1053 = 11 < 𝑐 / 1053 = 12 𝑎 / 1054 = 1 == 𝑐 / 1054 = 1 ⇒ 𝑑𝑜 = 12, 𝑒𝑜 = 53 52 𝑑𝑜 × 10𝑒𝑜+𝑒10 を文字列にする 𝑑𝑜 = 12 𝑒𝑜 + 𝑒10 = −1 ⇒”1.2” Copyright© Fixstars Group 37

文字列化アルゴリズム：実装面 先ほどのアルゴリズムのボトルネックは任意精度の演算を行う部分である。 任意精度の演算を行うのはSTEP3とSTEP4である。 • STEP3：2の冪を処理して浮動小数点を「整数×10の冪」と表す • • この整数が大きな整数となる STEP4：二つの大きな整数が等しくなるまで10で除算する 1. 2. 浮動小数点を数値にデコードする 隣の浮動小数点数との中点を計算する ○ 3. 𝑎, 𝑏, 𝑐 がとても 大きな整数 4. 5. 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三 つ組を 𝑢, 𝑣, 𝑤 × 2𝑒2 とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2 は整数） 2の冪を10の冪となるよう計算する ○ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10 の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10 は整数） ○ ○ 主要部とは、整数をm × 10𝑒 と書いた時のmのこととする 𝑎より大きく𝑐より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜 , 𝑒𝑜 を計算する ○ 𝑑𝑜 × 10𝑒𝑜+𝑒10 を文字列にする 主要部の桁数が最も小さく正しく丸められた整数を計算する 整数を文字列としてプリントする Copyright© Fixstars Group 38

具体例：Float32の𝒂, 𝒃, 𝒄の値 Float32 符号ビット：1、指数ビット：8、仮数ビット：23 bias：127 STEP1 STEP1 STEP2 STEP3 𝑒𝑓 : 0 ~ 224 𝑚𝑓 : −149~104 𝑢, 𝑣, 𝑤: 0~ 226 𝑒2 : −151~102 STEP2 非常に大きい！ 𝑎, 𝑏, 𝑐: max 226 × 2102 , 226 × 5151 ≒ 10131 STEP3 Copyright© Fixstars Group 39

具体例：Float64の𝒂, 𝒃, 𝒄の値 Float32 符号ビット：1、指数ビット：11、仮数ビット：52 bias：1023 STEP1 STEP1 STEP2 STEP3 𝑒𝑓 : 0 ~ 253 𝑚𝑓 : −1074~971 𝑢, 𝑣, 𝑤: 0~ 255 𝑒2 : −1076~969 STEP2 非常に大きい！ 𝑎, 𝑏, 𝑐: max 255 × 2969 , 255 × 51076 ≒ 10768 Copyright© Fixstars Group STEP3 40

大きな整数の回避（１） STEP3とSTEP4を工夫できないか？ • • STEP3：2の冪を処理して浮動小数点数を「整数×10の冪」とあらわす • この整数が大きな整数となる STEP4：二つの大きな整数が等しくなるまで10で除算する STEP3で「大きな整数」があらわれ、STEP4で「大きな整数を10で除算した小 さな整数」にする。 STEP3とSTEP4を繋げることで、「大きな整数」を計算することなく、直接 「大きな整数を10で除算した小さな整数」を計算することを考える。 Copyright© Fixstars Group 41

大きな整数の回避（２） 今後、𝒆𝟐 ≥ 𝟎の場合に限定する。 STEP3とSTEP4で行う計算の概略を右に示す。 数学的な考察から、STEP4のwhileループは少なくとも次の𝑖 = 𝑞まで はまわることが分かる。 𝑞 = max(0, 𝑒2 × log10 2 − 1) 𝑎を10で𝑞回除算した値𝑎𝑞 は𝑢の100倍未満のため、小さな整数であ る： 𝑎𝑞 = 𝑎 / 10𝑞 ≒ 𝑢 × 2𝑒2 / 10𝑞 < 100𝑢 よって、whileループが𝑞回まわった時点の変数の値は小さい。 STEP3 𝑎 = 𝑢 × 2𝑒2 , 𝑐 = 𝑤 × 2𝑒2 STEP4 while 𝑎 / 10𝑖 < 𝑐 − 1 / 10𝑖 : 𝑖 =𝑖+1 2𝑒2 ≒ 10𝑞 となる𝑞を選んでいる ⇒STEP3の段階で、STEP4のwhileループが𝒊 = 𝒒までまわった時の 変数の値を計算できれば、大きな整数の出現を回避できる！ 以降、 𝑎を計算せずに、直接𝑎𝑞 を計算する方法について考察する。 Copyright© Fixstars Group 42

𝒂を10で𝒒回除算した値の計算 compute_shortestに登場する変数は色々あるが、 𝑎を10で𝑞回除算した値𝑎𝑞 の計算が最も困難 なので、 𝑎𝑞 の計算方法のみ解説する。 𝑎𝑞 を計算しやすい形に式変形する。 大きな値に赤枠を付ける。 非整数になりうる値に青枠を付ける。 𝑎𝑞 = 𝑢 × 2𝑒2 / 10𝑞 = 𝑢 × 2𝑒2 −𝑞 / 5𝑞 （アイディア：大きな整数2𝑘 で5𝑞 の除算をなんとかする） = 𝑢 × 2𝑒2 −𝑞−𝑘 × 2𝑘 / 5𝑞 （アイディア： 𝑘が十分大きいと次の式が成立する） = 𝑢 × 2𝑒2 −𝑞−𝑘 × 式変形☆ 次スライド説明する 2𝑘 / 5𝑞 + 1 Copyright© Fixstars Group 43

𝒌の下界の算出：式変形☆の説明 𝒌の下界の算出 [1]より引用 𝑘の値が左式の不等式を満たすと前スライ ドの式変形☆が成立する。 証明は省略する。 青枠のことを「𝑘の下界」と呼ぶ。 赤枠の計算は難しいが、 𝑒2 ごとに予め計 算しておく量で、実行時には必要ない。 赤枠の計算方法の説明は省略する。 Copyright© Fixstars Group 44

変形した数式について ☆で式変形した下式は以下の点で優れている。 実行時変数は𝑢, 𝑒2 , 𝑞, 𝑘である。 𝑎𝑞 = 𝑢 × 2𝑒2−𝑞−𝑘 × 2𝑘 / 5𝑞 + 1 • 𝑎𝑞 は三つの項の積である • 赤枠（ 2𝑒2−𝑞−𝑘 ）部分は𝑒2 − 𝑞 − 𝑘 < 0のため除算になるが、シフト演算で効率的に実装できる • 困難な5𝑞 による除算が 2𝑘 / 5𝑞 に吸収されている • • Float32やFloat64の場合、 𝑞の値はそれほど大きくない 𝑞ごとに 2𝑘 / 5𝑞 の値をあらかじめ計算してテーブル化すれば、実行時には5𝑞 による除算がなくなる 赤枠と緑枠の項が大きな値とならないような𝑘が存在すれば、効率的に𝑎𝑞 を計算できる。 ⇒2𝑘 ≒ 5𝑞 となる𝑘 = log 2 5𝑞 付近で、下界より大きな𝑘 は存在するか？ Copyright© Fixstars Group 45

丁度良い𝒌の探索 • • • Float64で𝑘の下界を計算した（右図） • 横軸は𝑒2 • 𝑘の下界を青色でプロット • 2𝑘 / 5𝑞 の値をテーブルとして保持するのに必要な要素数 𝑘の下界 [1]より引用 の最小値を橙色でプロット • おおよそ 𝑘 − log 2 5𝑞 に対応する 橙色はプラトーになっている！ • 最大値を𝐵0 とする 𝑘 = 𝐵0 + log2 5𝑞 とすれば、𝑘は下界よりも大きな値と なる • さらに、実際に計算すると分かるが、この𝑘は丁度良い𝑘で • ある！（次の二ページで計算する） 2𝑘 と5𝑞 は大きな整数だが、 𝑞ごとに緑枠の値をあらかじめ 計算しておけば、実行時はテーブルを参照するだけで良い Copyright© Fixstars Group 46

Float32の場合 • 変数の整理 • • • 𝑒2 : 0 ∼ 102, 𝑢: 0 ∼ 226 𝑞 = max(0, 𝑒2 × log10 2 − 1) 𝑘 = 𝐵0 + log 2 5𝑞 • 𝐵0 = 60 • 𝑞 + 𝑘 − 𝑒2 ≒ 𝐵0 + 𝑞 × 1 + log 2 5 − 𝑒2 ≒ 𝐵0 + 𝑒2 × log10 2 − 1 × log 2 10 − 𝑒2 = 𝐵0 + 𝑒2 − log 2 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の範囲で計算できる Copyright© Fixstars Group 47

Float64の場合 • 変数の整理 • • • 𝑒2 : 0 ∼ 969, 𝑢: 0 ∼ 255 𝑞 = max(0, 𝑒2 × log10 2 − 1) 𝑘 = 𝐵0 + log 2 5𝑞 • 𝐵0 = 124 • 𝑞 + 𝑘 − 𝑒2 ≒ 𝐵0 + 𝑞 × 1 + log 2 5 − 𝑒2 ≒ 𝐵0 + 𝑒2 × log10 2 − 1 × log 2 10 − 𝑒2 = 𝐵0 + 𝑒2 − log 2 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の範囲で計算できる Copyright© Fixstars Group 48

大きな整数の回避：まとめ • STEP3の計算を工夫することで、STEP4のwhileループで𝑖 = 𝑞の時点までス キップし、大きな整数の出現を回避できる • 工夫した計算はFloat32の場合、uint128_tで行える • Float64の場合、uint256_tで行える • 𝑞ごとに下式の緑枠の値を予め計算しておくことで、文字列化の際はテーブル 参照するだけで良い Copyright© Fixstars Group 49

文字列化アルゴリズム：実装面（まとめ） 1. 浮動小数点数を数値にデコードする 2. 隣の浮動小数点数との中点を計算する ○ 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三つ組を 𝑢, 𝑣, 𝑤 × 2𝑒2 とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2 は整数） 3. 2の冪を10の冪となるよう計算する ○ ○ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10 の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10 は整数） 三つ組を10でq回除算した値(𝑎𝑞 , 𝑏𝑞 , 𝑐𝑞 )を計算する（𝑒10 も計算する） ○ ○ 主要部とは、整数をm × 10𝑒 と書いた時のmのこととする 𝑎𝑞 より大きく𝑐𝑞 より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜 , 𝑒𝑜 を計算する ○ 𝑑𝑜 × 10𝑒𝑜+𝑞+𝑒10 を文字列にする 4. 主要部の桁数が最も小さく正しく丸められた整数を計算する 5. 整数を文字列としてプリントする Copyright© Fixstars Group 50

性能評価： Ryū 浮動小数点を乱数で生成し、文字列化に要する時 間を比較した。 • 横軸：浮動小数点のビット表現を符号なし整 数とみなした時の値 • 縦軸：文字列化に要する時間[ns] 横軸の値がFloat32の場合は231 の時、Float64の 場合は 263 の時、浮動小数点の符号が入れ替わる。 c, dのグラフが線対称なのはこのためである。 浮動小数点を文字列化するのに要する時間の比較 [1]より引用 上段：C、下段：Java 左側：Float32、右側：Float64 いずれの場合もRyūが最速だった。 実験環境 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. Copyright© Fixstars Group 51

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”と文字列化する。

Copyright© Fixstars Group

52

GPU移植：実装（１） GPUでのtoStringを３つの関数に分けて実装した。 float to decimal double to decimal という意味 1. f2dKernel, d2dKernel ○ ○ 浮動小数点を正しく丸めた10進数表現を計算する ■ 文字列化アルゴリズムのSTEP1～STEP4に対応する 文字列化したときの文字数も計算する 2. thrust::exclusive_scan ○ prefix sumで、文字列化した浮動小数点を書き込むインデックスを計算する 3. d2sKernel ○ 10進数表現を文字列化し、前ステップで計算したインデックスに文字列を書き込む ■ 文字列化アルゴリズムのSTEP5に対応する Copyright© Fixstars Group 53

GPU移植：実装（２） v 1.2 3.14 -0.7 区切り文字列が空文字列の場合 … 0.006 f2dKernel, d2dKernel GPUの各スレッドが配列の各要素を処理する d 12 314 -7 6 … e -1 -2 -1 -3 … size 3 4 4 5 … 0 index 3 7 11 16 thrust::exclusive_scan size配列のprefix sumを計算する … d2sKernel GPUの各スレッドが配列の各要素を文字列化し、 indexで指定されたインデックスに書き込む s 1 . 2 3 . 1 4 - 0 . Copyright© Fixstars Group 7 0 . 0 0 6 … 54

GPU移植：性能評価（１） 3つのアルゴリズムを比較した。 1. Ryū（CPU） 2. Ryū（GPU） 3. yyjson（CPU） ○ ○ ○ ○ Cで実装された、高速なJSONライブラリ 浮動小数点の文字列化はSchubfachベースのアルゴリズムで行う ※浮動小数点配列を文字列化した結果はRyūと異なる ※フェアな比較対象ではないが、参考として比較した Ryū（CPU）とyyjson（CPU）はCPUのシングルスレッドで動作するプログラムであ る。 Copyright© Fixstars Group 55

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の時間 Copyright© Fixstars Group 56

GPU移植：性能評価（３） • • • Ryū（CPU）は浮動小数点を一要素当たり、floatは43[ns]で、doubleは60[ns]で文字列化した • おおよそ論文通りのオーダーだった Ryū（GPU）はCPU版よりfloatで200倍ほど、doubleで100倍ほど高速となった • 実行時間の8割はd2sKernelで、これを高速化すればもっと速くなる yyjson（CPU）はRyū（CPU）よりも高速だった • 本実験でも、SchubfachはRyūよりも高速だった floatの配列を文字列化 するのに要する時間 [ms] yyjson（CPU）は浮動小 数点を全てdoubleで扱う ため除外 配列の要素数 Ryū（GPU） 1,000,000 46.2 0.513 10,000,000 433 2.42 100,000,000 4,335 22.2 配列の要素数 doubleの配列を文字列化 するのに要する時間 [ms] Ryū（CPU） 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 Copyright© Fixstars Group 57

最新アルゴリズムの紹介：概要 Dragonboxを最新のアルゴリズムとして、その周辺の次の4つのアルゴリズムを 簡単に紹介する。 • Schubfach • Ryū • Grisu-Exact • Dragonbox Copyright© Fixstars Group 58

最新アルゴリズムの紹介：比較 初出 時期 著者 アルゴリズム 名前の由来 Schubfach OpenJDKのcorelibs-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の計算で済ませるよ うにした。 Copyright© Fixstars Group 59

最新アルゴリズムの紹介：性能評価 • • 凡例を拡大 文字列化に要する時間を4つのアルゴリズムで比較した ([6]より引用) 結果 • • • • • 上段：Float32、下段：Float64 左側：最短文字列がN桁の浮動小数点を文字列化するのに要する時間（Nが横軸） 右側：浮動小数点をビットとして一様に生成し、文字列化に要する時間をプロット Dragonboxが最速だった 最短文字列の桁数が小さいほど、Ryūは他のアルゴリズムよりも遅かった compute_shortest でのナイーブなループが遅い • 実験環境 CPU：Intel (R) Core™ i77700HQ CPU @2.80GHz C++：clang-cl compiler shipped with Visual Studio 2022 17.0.4 Copyright© Fixstars Group 60

参考資料 • • • • • • • [1] https://dl.acm.org/doi/10.1145/3192366.3192369 • Ryūの論文 • Ryū著者による実装 • • 浮動小数点をminimal length decimalに変換するアルゴリズムの実行時間を比較する 文字列としてプリントするのに要する時間を除いて計測し、比較している • 各種文字列化アルゴリズムがどのライブラリで使用されているか調べている • Wikipedia IEEE754 • Dragonboxのアルゴリズムを説明したpdf • yyjson [2] https://github.com/ulfjack/ryu [3] https://github.com/abolz/Drachennest [4] https://kubo39.hatenablog.com/entry/2022/06/09/浮動小数点数の10進表記についてまとめてみ た [5] https://ja.wikipedia.org/wiki/IEEE_754 [6] https://github.com/jk-jeon/dragonbox/blob/master/other_files/Dragonbox.pdf [7] https://github.com/ibireme/yyjson Copyright© Fixstars Group 61

• https://dl.acm.org/doi/10.1145/3192366.3192369
• https://github.com/ulfjack/ryu
• https://github.com/abolz/Drachennest
• https://kubo39.hatenablog.com/entry/2022/06/09/浮動小数点数の10進表記についてまとめてみた
• https://ja.wikipedia.org/wiki/IEEE_754
• https://github.com/jk-jeon/dragonbox/blob/master/other_files/Dragonbox.pdf
• https://github.com/ibireme/yyjson

Thank you! お問い合わせ窓口 : [email protected] Copyright© Fixstars Group 62