2023年2月1日(水)に実施した「CPU/GPU高速化セミナー ～浮動小数点から文字列への高速変換の論文を読んでみた～」の当日資料です。
Copyright© Fixstars GroupCPU/GPU高速化セミナー浮動小数点から文字列への高速変換の論文を読んでみた
View Slide
Copyright© Fixstars Group本日のAgenda● はじめに● フィックスターズのご紹介● Ryū のアルゴリズム解説● Ryū の GPU 移植● 最新アルゴリズムの紹介● Q&A / 告知2
Copyright© Fixstars Groupはじめに3
Copyright© Fixstars Group本講演の位置づけ4● 弊社がサービス展開している、ソフトウェアによる処理速度向上に関連した様々な技術情報を発信しています● 性能モデルの理論と実践 (発表資料)○ ソフトウェア性能の理論的根拠となる性能モデルの使用方法を実践例● 今回の内容○ Ryū のアルゴリズムについて、理論と実装を紹介○ GPU移植の説明と、CPU との処理速度を比較● こんな方に向いています○ アルゴリズム・高速計算の実装に興味のある方○ システム性能向上に向けた、技術アプローチに興味のある方
Copyright© Fixstars Group発表者紹介冨田 明彦ソリューションカンパニー 営業企画2008年に入社。金融、医療業界において、ソフトウェア高速化業務に携わる。その後、新規事業企画、半導体業界の事業を担当し、現職。5大野 真暉ソリューション第一事業部 シニアエンジニア2019年に新卒で入社。製造業におけるスケジューリング問題のソルバー開発に従事。数理最適化、量子コンピュータ関連のプロジェクトを手掛ける。
Copyright© Fixstars Groupフィックスターズのご紹介6
Copyright© Fixstars Groupフィックスターズの強みコンピュータの性能を最大限に引き出す、ソフトウェア高速化のエキスパート集団ハードウェアの知見 アルゴリズム実装力 各産業・研究分野の知見7目的の製品に最適なハードウェアを見抜き、その性能をフル活用するソフトウェアを開発します。ハードウェアの特徴と製品要求仕様に合わせて、アルゴリズムを改良して高速化を実現します。開発したい製品に使える技術を見抜き、実際に動作する実装までトータルにサポートします。
Copyright© Fixstars Groupサービス提供分野8半導体自動車産業機器生命科学金融● NAND型フラッシュメモリ向けファームウェア開発● 次世代AIチップの開発環境基盤● 自動運転の高性能化、実用化● 次世代パーソナルモビリティの研究開発● Smart Factory実現への支援● マシンビジョンシステムの高速化● ゲノム解析の高速化● 医用画像処理の高速化● AI画像診断システムの研究開発● デリバティブシステムの高速化● HFT(アルゴリズムトレード)の高速化
Copyright© Fixstars Groupサービス領域様々な領域でソフトウェア高速化サービスを提供しています。大量データの高速処理は、お客様の製品競争力の源泉となっています。9組込み高速化画像処理・アルゴリズム開発分散並列システム開発GPU向け高速化FPGAを活用したシステム開発量子コンピューティングAI・深層学習自動車向けソフトウェア開発フラッシュメモリ向けファームウェア開発
Copyright© Fixstars Groupコンサルティングセットアップ支援、処理の分割や変換等のコンサルクラウド実行環境のご提供クラウド経由での量子コンピュータ利用サービスを提供ソフトウェア高速化・開発支援サービス量子コンピュータを組み合わせてシステムの高速化を実現セミナー・トレーニング量子コンピュータの研究動向や活用事例、実際の利用方法等量子コンピューティング活用支援多様なハードウェアでのソフトウェア高速化サービスに加え、量子コンピュータ活用支援とシステム開発を提供しています。お客様の課題 ご支援内容量子コンピューティングが課題の解決に役立つか、確信が持てない量子コンピューティングの検討をどう進めていったら良いかわからない作りたいアプリケーションがあるが、開発が難しい10
Copyright© Fixstars Group物流倉庫における人員配置最適化の社会実装開発内閣府主催SIP公開シンポジウム2021111 通販向け物流倉庫で誰がどの仕事を担当するか割り振る量子コンピューティングシステムを開発2 30分以上かかっていた割り当て作業が1秒未満の計算で実現3 業務負荷の偏りが解消し、全体の作業効率も30％改善する割り当てパターンが計算できた作業者の負荷のばらつきを抑えつつ、全体の作業効率は30%向上物流分野 量子コンピューティング活用支援サービス領域ベルメゾンロジスコ様のご協力を得て実施
Copyright© Fixstars GroupRyū12
Copyright© Fixstars GroupRyū2018年のPLDIで発表された論文 “Ryū: fast float-to-string conversion” を紹介します。浮動小数点を高速に文字列に変換するアルゴリズムを提案・実装した論文です。固定幅整数の演算のみで文字列に変換し、正確かつ高速なことが特徴です。Adams, Ulf. "Ryū: fast float-to-string conversion." Proceedings of the 39th ACM SIGPLANConference on Programming Language Design and Implementation. 2018.13
Copyright© Fixstars GroupRyū - 目次1. はじめに2. 浮動小数点の復習：IEEE7543. 文字列化のルール4. 文字列化アルゴリズム• 理論面• 実装面5. 性能評価6. GPU移植7. 最新アルゴリズムの紹介8. 参考資料14
Copyright© Fixstars Groupはじめに浮動小数点の文字列化アルゴリズムはいたるところで用いられる、重要なアルゴリズムである。例えば、Javascriptでは基本的な数値型のNumber型をプリントするたびに浮動小数点を文字列化している。他にも、JSONフォーマットへ浮動小数点をシリアライズする際など、様々な箇所で使用される。15
Copyright© Fixstars Group浮動小数点についての復習：IEEE754（１）IEEE754は浮動小数点数のビット表現（交換形式）を定義している。ビット表現は三つの要素から構成される。浮動小数点の形式ごとに、それぞれの要素のビット幅が決まっている。1. 符号 sign2. 指数 exponent3. 仮数 fraction（mantissaとも）16Float64のビット表現 [5]より引用符号のビット幅=1、指数のビット幅=11、仮数のビット幅=52
Copyright© Fixstars Group浮動小数点についての復習：IEEE754（２）浮動小数点のビット表現と数値の対応関係は下式で定義される。• 𝑠：符号部の値• 𝑠 = 0なら正の値、𝑠 = 1なら負の値• 𝑒：指数部の値• biasは浮動小数点の形式ごとに決まっている値• Float32なら126で、Float64なら1022• 𝑚：仮数部のビット列• 1.mと0.mは2進数で記載した小数17浮動小数点の数値 [1]より引用
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
Copyright© Fixstars Groupアルゴリズム概要 ~1.2を具体例に~（１）1.2を例に文字列化アルゴリズムの概要を紹介する。1.2のFloat64のビット表現は次のようになる。0011111111110011001100110011001100110011001100110011001100110011このビット表現を定義通り数値に変換すると次の値になる。𝑓 =1.1999999999999999555910790149937383830547332763671875この文字列は正確だがあまりに長い。先ほどのルールの「2. 文字列は最短」に違反する。19符号 指数 仮数浮動小数点の数値 [1]より引用
Copyright© Fixstars Groupアルゴリズム概要 ~1.2を具体例に~（２）数値𝒇の「近く」で、桁数の少ない数値を探す。「近く」を正確に表現するために、下の数直線のように隣の浮動小数点数との境界を考える。𝑓− ： 𝑓より一つ小さなFloat64𝑓+ ： 𝑓より一つ大きなFloat64𝑓の区間の中で、十進法の小数で表現したときに最も桁数の少ない実数が1.2なので、”1.2”と文字列化する。この文字列化アルゴリズムは先ほど定めた文字列化のルールに従う。20𝑓−の区間𝑓−𝑓+の区間𝑓+𝑓の区間𝑓1.2𝑓−と𝑓の中点1.19999999999999984456877655247808434069156646728515625𝑓+と𝑓の中点1.20000000000000006661338147750939242541790008544921875
Copyright© Fixstars Group文字列化アルゴリズム文字列化アルゴリズムは複雑なので二つのステップに分けて解説する。STEP1：理論的な話任意精度の計算を駆使して浮動小数点を文字列に変換する。STEP2：実装的な話STEP1のアルゴリズムで任意精度が必要な計算を、全て固定幅整数の計算に置き換える。21
Copyright© Fixstars Groupおことわり正の浮動小数点の文字列化のみを考える。つまり、以下の入力は除外する。• NaN• Infinity• 負の浮動小数点• 負の値を文字列化するときは、その正の値の文字列の先頭に「-」を付けると良い• +0.0• アルゴリズムの都合上除外するまた、指数部の値が変わる特殊なケースや丸めなどの話はスキップする。（大事ではあるがアルゴリズムの概要を把握する分には不要）22
Copyright© Fixstars Group文字列化アルゴリズム：理論面1. 浮動小数点を数値にデコードする2. 隣の浮動小数点数との中点を計算する○ 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三つ組を𝑢, 𝑣, 𝑤 × 2𝑒2とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2は整数）3. 2の冪を10の冪となるよう計算する○ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10は整数）4. 主要部の桁数が最も小さく正しく丸められた整数を計算する○ 主要部とは、整数をm × 10𝑒と書いた時のmのこととする○ 𝑎より大きく𝑐より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜, 𝑒𝑜を計算する5. 整数を文字列としてプリントする○ 𝑑𝑜× 10𝑒𝑜+𝑒10を文字列にする23
Copyright© Fixstars GroupSTEP1：浮動小数点を数値にデコードする• 概要• IEEE754の通りにデコードする• 小数部をはらって整数にする• 入力• 𝑒 指数部• 𝑚 仮数部• 出力• 𝑒𝑓, 𝑚𝑓24浮動小数点の数値 [1]より引用𝒆𝒇, 𝒎𝒇の計算式 [1]より引用𝑙𝑒𝑛(𝑚)：仮数のビット数𝑓 = 𝑚𝑓× 2𝑒𝑓が成立する。
Copyright© Fixstars GroupSTEP1：具体例 1.2250s1023e900719925474099m𝑓 = −1 0 × 252 + 900719925474099 × 21023 −1023 −52𝑚𝑓= 5404319552844595 𝑒𝑓= −521. 浮動小数点を数値にデコードする
Copyright© Fixstars GroupSTEP2：隣の浮動小数点数との中点を計算する• 概要• 一つ小さな浮動小数点数𝑓−との中点を計算する• 一つ大きな浮動小数点数𝑓+との中点を計算する• 隣の浮動小数点数はIEEE754の式とにらめっこすれば計算できる• 入力• 𝑒, 𝑚• 𝑒𝑓, 𝑚𝑓• 出力• 𝑒2• 𝑢, 𝑣, 𝑤26𝑢 ⋅ 2𝑒2𝑓+𝑓𝑤 ⋅ 2𝑒2𝑣 ⋅ 2𝑒2𝑓−三つ組の計算式𝑒2= 𝑒𝑓− 2𝑢 = 4𝑚𝑓− ቊ1 if 𝑚 = 0 … 0 2and 𝑒 > 0 … 01 22 otherwise𝑣 = 4𝑚𝑓𝑤 = 4𝑚𝑓+ 2𝑚𝑓× 2𝑒𝑓 = 𝑣 × 2𝑒2 が成立する。
Copyright© Fixstars GroupSTEP2：具体例 1.227𝑢 ⋅ 2𝑒2𝑓+𝑓 = 𝑚𝑓× 2𝑒𝑓𝑤 ⋅ 2𝑒2𝑣 ⋅ 2𝑒2𝑓− 𝑒2= 𝑒𝑓− 2 = −54𝑣 = 4 × 𝑚𝑓= 21617278211378380𝑢 = 𝑣 − 2 = 21617278211378378𝑤 = 𝑣 + 2 = 216172782113783822. 隣の浮動小数点との中点を計算する
Copyright© Fixstars GroupSTEP3：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 が成立する。
Copyright© Fixstars GroupSTEP3：具体例 1.229𝑢 × 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= 12000000000000000666133814775093924254179000854492187503. 2の冪を10の冪となるよう計算する
Copyright© Fixstars GroupSTEP4：主要部の桁数が最も小さく正しく丸められた整数を計算する• 概要• compute_shortestで、𝑎より大きく𝑐より小さい整数の中で、主要部の桁数が最も小さい整数を求める• 数学的に定式化すると、次の不等式を満たす整数𝑒𝑜, 𝑑𝑜を求めることに対応する• 𝒂 < 𝒅𝒐× 𝟏𝟎𝒆𝒐 < 𝒄• ただし𝑑𝑜の桁数は最小とする• つまり𝑒𝑜の値は最大• 入力• 𝑎, 𝑏, 𝑐• 丸めのオプション• 出力• 𝑒𝑜, 𝑑𝑜30𝑎と𝑐が等しくなるまで10で除算するcompute_shortest [1]より引用
Copyright© Fixstars GroupSTEP4：（資料）正しいcompute_shortest右の関数が真のcompute_shortestである。前ページのcompute_shortestには一部嘘があった。嘘とは、文字列化のルールの一つに違反していること。(3) 正確な丸め：1と2のルールを満たす文字列で、文字列は浮動小数点の数値に最も近くする隣の浮動小数点数との中点である𝑎と𝑐のみならず、𝑏も一緒に計算することでこのルールも満たしたアルゴリズムになる。31def compute_shortest(a, b, c, accept_smaller, accept_larger, break_tie_down):i=0a[0]=a, b[0]=b, c[0]=accept_larger ? c : c-1all_a_zero[0]=True, all_b_zero[0]=Truedigit[0]=0while floor(a[i]/10) < floor(c[i]/10):all_a_zero[i+1]=all_a_zero[i] and a[i]%10==0all_b_zero[i+1]=all_b_zero[i] and digit[i]==0a[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]%10i=i+1if accept_smaller and all_a_zero[i]:while a[i]%10==0:all_b_zero[i+1]=all_b_zero[i] and digit[i]==0a[i+1]=a[i]/10, b[i+1]=floor(b[i]/10), c[i+1]=floor(c[i]/10)digit[i+1]=b[i]%10i=i+1is_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 GroupSTEP4：具体例 1.232⋮𝑎 / 1052 = 119 < 𝑐 / 1052 = 120𝑎 / 1053 = 11 < 𝑐 / 1053 = 12𝑎 / 1054 = 1 == 𝑐 / 1054 = 1⇒ 𝑑𝑜= 12, 𝑒𝑜= 534. 主要部の桁数が最も小さく正しく丸められた整数を計算する
Copyright© Fixstars GroupSTEP5：整数を文字列としてプリントする• 概要• 𝑑𝑜× 10𝑒𝑜+𝑒10を文字列にする• 入力• 𝑒10• 𝑒𝑜, 𝑑𝑜• 文字列化のオプション• 出力• 文字列33
Copyright© Fixstars GroupSTEP5：具体例 1.234𝑑𝑜× 10𝑒𝑜+𝑒10を文字列にする𝑑𝑜= 12𝑒𝑜+ 𝑒10= 53 − 54 = −1⇒”1.2”5. 整数を文字列としてプリントする
Copyright© Fixstars Group文字列化アルゴリズム：理論面（再掲）1. 浮動小数点を数値にデコードする2. 隣の浮動小数点数との中点を計算する○ 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三つ組を𝑢, 𝑣, 𝑤 × 2𝑒2とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2は整数）3. 2の冪を10の冪となるよう計算する○ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10は整数）4. 主要部の桁数が最も小さく正しく丸められた整数を計算する○ 主要部とは、整数をm × 10𝑒と書いた時のmのこととする○ 𝑎より大きく𝑐より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜, 𝑒𝑜を計算する5. 整数を文字列としてプリントする○ 𝑑𝑜× 10𝑒𝑜+𝑒10を文字列にする35
Copyright© Fixstars Group文字列化アルゴリズム：具体例（まとめ）36STEP10s1023e900719925474099m𝑓 = −1 0 × 252 + 900719925474099 × 21023 −1023 −52𝑚𝑓= 5404319552844595 𝑒𝑓= −52𝑢 ⋅ 2𝑒2𝑓+𝑓𝑤 ⋅ 2𝑒2𝑣 ⋅ 2𝑒2𝑓− 𝑒2= 𝑒𝑓− 2 = −54𝑣 = 4 × 𝑚𝑓= 21617278211378380𝑢 = 𝑣 − 2 = 21617278211378378𝑤 = 𝑣 + 2 = 21617278211378382STEP2
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= 1200000000000000066613381477509392425417900085449218750STEP3STEP4STEP5⋮𝑎 / 1052 = 119 < 𝑐 / 1052 = 120𝑎 / 1053 = 11 < 𝑐 / 1053 = 12𝑎 / 1054 = 1 == 𝑐 / 1054 = 1⇒ 𝑑𝑜= 12, 𝑒𝑜= 53𝑑𝑜× 10𝑒𝑜+𝑒10を文字列にする𝑑𝑜= 12𝑒𝑜+ 𝑒10= −1⇒”1.2”
Copyright© Fixstars Group文字列化アルゴリズム：実装面先ほどのアルゴリズムのボトルネックは任意精度の演算を行う部分である。任意精度の演算を行うのはSTEP3とSTEP4である。• STEP3：2の冪を処理して浮動小数点を「整数×10の冪」と表す• この整数が大きな整数となる• STEP4：二つの大きな整数が等しくなるまで10で除算する381. 浮動小数点を数値にデコードする2. 隣の浮動小数点数との中点を計算する○ 小さな浮動小数点数との中点、元の浮動小数点数、大きな浮動小数点数との中点の三つ組を𝑢, 𝑣, 𝑤 × 2𝑒2とする（𝑢, 𝑣, 𝑤は正の整数、 𝑒2は整数）3. 2の冪を10の冪となるよう計算する○ 三つ組を 𝑎, 𝑏, 𝑐 × 10𝑒10の形に変形する（𝑎, 𝑏, 𝑐は正の整数、 𝑒10は整数）4. 主要部の桁数が最も小さく正しく丸められた整数を計算する○ 主要部とは、整数をm × 10𝑒と書いた時のmのこととする○ 𝑎より大きく𝑐より小さい整数の中で主要部の桁数が小さい整数𝑑𝑜, 𝑒𝑜を計算する5. 整数を文字列としてプリントする○ 𝑑𝑜× 10𝑒𝑜+𝑒10を文字列にする𝑎, 𝑏, 𝑐 がとても大きな整数
Copyright© Fixstars Group具体例：Float32の𝒂, 𝒃, 𝒄の値Float32符号ビット：1、指数ビット：8、仮数ビット：23bias：127STEP1𝑒𝑓: 0 ~ 224𝑚𝑓: −149~104STEP2𝑢, 𝑣, 𝑤: 0~ 226𝑒2: −151~102STEP3𝑎, 𝑏, 𝑐: max 226 × 2102, 226 × 5151 ≒ 1013139STEP1STEP2STEP3非常に大きい！
Copyright© Fixstars Group具体例：Float64の𝒂, 𝒃, 𝒄の値Float32符号ビット：1、指数ビット：11、仮数ビット：52bias：1023STEP1𝑒𝑓: 0 ~ 253𝑚𝑓: −1074~971STEP2𝑢, 𝑣, 𝑤: 0~ 255𝑒2: −1076~969STEP3𝑎, 𝑏, 𝑐: max 255 × 2969, 255 × 51076 ≒ 1076840STEP1STEP2STEP3非常に大きい！
Copyright© Fixstars Group大きな整数の回避（１）STEP3とSTEP4を工夫できないか？• STEP3：2の冪を処理して浮動小数点数を「整数×10の冪」とあらわす• この整数が大きな整数となる• STEP4：二つの大きな整数が等しくなるまで10で除算するSTEP3で「大きな整数」があらわれ、STEP4で「大きな整数を10で除算した小さな整数」にする。STEP3とSTEP4を繋げることで、「大きな整数」を計算することなく、直接「大きな整数を10で除算した小さな整数」を計算することを考える。41
Copyright© Fixstars Group大きな整数の回避（２）今後、𝒆𝟐≥ 𝟎の場合に限定する。STEP3とSTEP4で行う計算の概略を右に示す。数学的な考察から、STEP4のwhileループは少なくとも次の𝑖 = 𝑞まではまわることが分かる。𝑞 = max(0, 𝑒2× log102 − 1)𝑎を10で𝑞回除算した値𝑎𝑞は𝑢の100倍未満のため、小さな整数である：𝑎𝑞= 𝑎 / 10𝑞 ≒ 𝑢 × 2𝑒2 / 10𝑞 < 100𝑢よって、whileループが𝑞回まわった時点の変数の値は小さい。⇒STEP3の段階で、STEP4のwhileループが𝒊 = 𝒒までまわった時の変数の値を計算できれば、大きな整数の出現を回避できる！以降、 𝑎を計算せずに、直接𝑎𝑞を計算する方法について考察する。42STEP3𝑎 = 𝑢 × 2𝑒2, 𝑐 = 𝑤 × 2𝑒2STEP4while 𝑎 / 10𝑖 < 𝑐 − 1 / 10𝑖:𝑖 = 𝑖 + 12𝑒2 ≒ 10𝑞となる𝑞を選んでいる
Copyright© Fixstars Group𝒂を10で𝒒回除算した値の計算compute_shortestに登場する変数は色々あるが、 𝑎を10で𝑞回除算した値𝑎𝑞の計算が最も困難なので、 𝑎𝑞の計算方法のみ解説する。𝑎𝑞を計算しやすい形に式変形する。大きな値に赤枠を付ける。非整数になりうる値に青枠を付ける。𝑎𝑞= 𝑢 × 2𝑒2 / 10𝑞 = 𝑢 × 2𝑒2−𝑞 / 5𝑞（アイディア：大きな整数2𝑘で5𝑞の除算をなんとかする）= 𝑢 × 2𝑒2−𝑞−𝑘 × 2𝑘 / 5𝑞（アイディア： 𝑘が十分大きいと次の式が成立する）= 𝑢 × 2𝑒2−𝑞−𝑘 × 2𝑘 / 5𝑞 + 143式変形☆次スライド説明する
Copyright© Fixstars Group𝒌の下界の算出：式変形☆の説明𝑘の値が左式の不等式を満たすと前スライドの式変形☆が成立する。証明は省略する。青枠のことを「𝑘の下界」と呼ぶ。赤枠の計算は難しいが、 𝑒2ごとに予め計算しておく量で、実行時には必要ない。赤枠の計算方法の説明は省略する。44𝒌の下界の算出 [1]より引用
Copyright© Fixstars Group変形した数式について☆で式変形した下式は以下の点で優れている。実行時変数は𝑢, 𝑒2, 𝑞, 𝑘である。• 𝑎𝑞は三つの項の積である• 赤枠（ 2𝑒2−𝑞−𝑘 ）部分は𝑒2− 𝑞 − 𝑘 < 0のため除算になるが、シフト演算で効率的に実装できる• 困難な5𝑞による除算が 2𝑘 / 5𝑞 に吸収されている• Float32やFloat64の場合、 𝑞の値はそれほど大きくない• 𝑞ごとに 2𝑘 / 5𝑞 の値をあらかじめ計算してテーブル化すれば、実行時には5𝑞による除算がなくなる赤枠と緑枠の項が大きな値とならないような𝑘が存在すれば、効率的に𝑎𝑞を計算できる。⇒2𝑘 ≒ 5𝑞となる𝑘 = log25𝑞付近で、下界より大きな𝑘 は存在するか？45𝑎𝑞= 𝑢 × 2𝑒2−𝑞−𝑘 × 2𝑘 / 5𝑞 + 1
Copyright© Fixstars Group丁度良い𝒌の探索• Float64で𝑘の下界を計算した（右図）• 横軸は𝑒2• 𝑘の下界を青色でプロット• 2𝑘 / 5𝑞 の値をテーブルとして保持するのに必要な要素数の最小値を橙色でプロット• おおよそ 𝑘 − log25𝑞 に対応する• 橙色はプラトーになっている！• 最大値を𝐵0とする• 𝑘 = 𝐵0+ log25𝑞 とすれば、𝑘は下界よりも大きな値となる• さらに、実際に計算すると分かるが、この𝑘は丁度良い𝑘である！（次の二ページで計算する）• 2𝑘と5𝑞は大きな整数だが、 𝑞ごとに緑枠の値をあらかじめ計算しておけば、実行時はテーブルを参照するだけで良い46𝑘の下界 [1]より引用
Copyright© Fixstars GroupFloat32の場合• 変数の整理• 𝑒2: 0 ∼ 102, 𝑢: 0 ∼ 226• 𝑞 = max(0, 𝑒2× log102 − 1)• 𝑘 = 𝐵0+ log25𝑞• 𝐵0= 60• 𝑞 + 𝑘 − 𝑒2≒ 𝐵0+ 𝑞 × 1 + log25 − 𝑒2≒ 𝐵0+ 𝑒2× log102 − 1 × log210 − 𝑒2= 𝐵0+ 𝑒2− log210 − 𝑒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
Copyright© Fixstars GroupFloat64の場合• 変数の整理• 𝑒2: 0 ∼ 969, 𝑢: 0 ∼ 255• 𝑞 = max(0, 𝑒2× log102 − 1)• 𝑘 = 𝐵0+ log25𝑞• 𝐵0= 124• 𝑞 + 𝑘 − 𝑒2≒ 𝐵0+ 𝑞 × 1 + log25 − 𝑒2≒ 𝐵0+ 𝑒2× log102 − 1 × log210 − 𝑒2= 𝐵0+ 𝑒2− log210 − 𝑒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
Copyright© Fixstars Group大きな整数の回避：まとめ• STEP3の計算を工夫することで、STEP4のwhileループで𝑖 = 𝑞の時点までスキップし、大きな整数の出現を回避できる• 工夫した計算はFloat32の場合、uint128_tで行える• Float64の場合、uint256_tで行える• 𝑞ごとに下式の緑枠の値を予め計算しておくことで、文字列化の際はテーブル参照するだけで良い49
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
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.50GHzOS：Ubuntu Linux 17.10.C：clang 3.9 with the -O2 flagJava：OpenJDK 1.8.0_131.
Copyright© Fixstars GroupGPU移植RyūをGPUに移植し、大量の浮動小数点を文字列化するタスクの高速化を行った。浮動小数点の配列と区切り文字列を受け取り、浮動小数点の文字列を区切り文字でつなげた文字列を出力する関数を実装した。void toString(const std::vector& 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
Copyright© Fixstars GroupGPU移植：実装（１）GPUでのtoStringを３つの関数に分けて実装した。1. f2dKernel, d2dKernel○ 浮動小数点を正しく丸めた10進数表現を計算する■ 文字列化アルゴリズムのSTEP1～STEP4に対応する○ 文字列化したときの文字数も計算する2. thrust::exclusive_scan○ prefix sumで、文字列化した浮動小数点を書き込むインデックスを計算する3. d2sKernel○ 10進数表現を文字列化し、前ステップで計算したインデックスに文字列を書き込む■ 文字列化アルゴリズムのSTEP5に対応する53float to decimaldouble to decimal という意味
Copyright© Fixstars GroupGPU移植：実装（２）541.2 3.14 -0.7 0.006 …v12 314 -7 6 …f2dKernel, d2dKernelGPUの各スレッドが配列の各要素を処理する-1 -2 -1 -3 …de3 4 4 5 …size0 3 7 11 16 …indexthrust::exclusive_scansize配列のprefix sumを計算する1 . 2 3 . 1 4 - 0 . 7 0 . 0 0 6 …sd2sKernelGPUの各スレッドが配列の各要素を文字列化し、indexで指定されたインデックスに書き込む区切り文字列が空文字列の場合
Copyright© Fixstars GroupGPU移植：性能評価（１）3つのアルゴリズムを比較した。1. Ryū（CPU）2. Ryū（GPU）3. yyjson（CPU）○ Cで実装された、高速なJSONライブラリ○ 浮動小数点の文字列化はSchubfachベースのアルゴリズムで行う○ ※浮動小数点配列を文字列化した結果はRyūと異なる○ ※フェアな比較対象ではないが、参考として比較したRyū（CPU）とyyjson（CPU）はCPUのシングルスレッドで動作するプログラムである。55
Copyright© Fixstars GroupGPU移植：性能評価（２）実験環境○ 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
Copyright© Fixstars GroupGPU移植：性能評価（３）• 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.51310,000,000 433 2.42100,000,000 4,335 22.2配列の要素数 Ryū（CPU） Ryū（GPU） yyjson（CPU）1,000,000 60.3 0.795 53.710,000,000 604 5.60 533100,000,000 6,037 54.1 5,280floatの配列を文字列化するのに要する時間 [ms]yyjson（CPU）は浮動小数点を全てdoubleで扱うため除外doubleの配列を文字列化するのに要する時間 [ms]
Copyright© Fixstars Group最新アルゴリズムの紹介：概要Dragonboxを最新のアルゴリズムとして、その周辺の次の4つのアルゴリズムを簡単に紹介する。• Schubfach• Ryū• Grisu-Exact• Dragonbox58
Copyright© Fixstars Group最新アルゴリズムの紹介：比較59初出 時期 著者 アルゴリズム 名前の由来Schubfach OpenJDKのcore-libs-devのパッチ2018年4月, 11月RaffaelloGiuliettiRyūのほぼ上位互換Ryūのcompute_shortestをループ無しで行う鳩ノ巣原理のドイツ語のSchubfachprinzipに由来するRyū PLDI 2018年6月 Ulf Adams 本発表で解説した Grisuの別名のDragonを日本語に訳したGrisu-Exact GitHub 2020年 jk-jeon Ryū + GrisuDragonbox GitHub 2020年 jk-jeon Schubfach + Grisu Grisu⇒DragonSchubfach⇒box(鳩ノ巣原理＝箱入れ原理)※Grisu-ExactとDragonboxのアルゴリズムの「+ Grisu」についてRyūとSchubfachは三つ組で計算する。Float64の場合、64bit×128bitを三回も計算する必要があり、高コストである。Grisuのアイディアを利用して、三つ組の計算を省略し、一度の64bit×128bitの計算で済ませるようにした。
Copyright© Fixstars Group最新アルゴリズムの紹介：性能評価• 文字列化に要する時間を4つのアルゴリズムで比較した ([6]より引用)• 上段：Float32、下段：Float64• 左側：最短文字列がN桁の浮動小数点を文字列化するのに要する時間（Nが横軸）• 右側：浮動小数点をビットとして一様に生成し、文字列化に要する時間をプロット• 結果• Dragonboxが最速だった• 最短文字列の桁数が小さいほど、Ryūは他のアルゴリズムよりも遅かった• compute_shortest でのナイーブなループが遅い60実験環境CPU：Intel (R) Core™ i7-7700HQ CPU @2.80GHzC++：clang-cl compiler shippedwith Visual Studio 2022 17.0.4凡例を拡大
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• yyjson61
Copyright© Fixstars GroupThank you!お問い合わせ窓口 : [email protected]62
fixstars
poteboy
heyitsmohit
_tanako
cocoeyes02
yosuke_furukawa
kazue
dalinaum
skanehira
tapompom
workspace93
iuchimasatoshi
shohata
moore
bcantrill
eileencodes
malarkey
ammeep
sachag
revolveconf
frogandcode
tmm1
rasmusluckow
scottboms
lemiorhan