CPP-QUIZ from Unreal 2017 ( part-2/2; 答え編 )

（この記事は前の投稿 CPP-QUIZ from Unreal 2017 ( part-2/2; 出題編 ) に対応する答え編です。）

2. こたえ

Q1. float 型の [ 0.0 .. 1.0 ] の値を uint8 型の [ 0 .. 256 ) へ写像したかった

出題のソースコードでは、 "出力が 255 となる元の値の範囲" が 0 から 254 のそれぞれに対応する範囲に対して均等とならない問題が発生する。

例えば、値がグレイスケールカラーの輝度値を表す場合に出題のソースコードを用いると、入力が 1.0 未満の場合に本来よりも僅かに明るくシフトしてしまう結果となる。グレイスケールカラーで人間が目視するだけであれば、大きな問題にはならない事も多いが、変換を繰り返したり、その値での制御を一定時間続ける装置に適用したりすると問題と成り得る事もある。

出題のソースコードと in と out の関係の一部:

in	out
[ 0 / 255.f .. 1 / 255.f )	0
[ 1 / 255.f .. 2 / 255.f )	1
[ 2 / 255.f .. 3 / 255.f )	2
[ 254 / 255.f .. 255 / 255.f )	254
[ 255 / 255.f .. 255 / 255.f ]	255

修正方法:

256 未満で最大の float の数を掛けると [ 0 .. 256 ) に可能な限り均等に写像できる。

int out = in * 255.9999847412109375f;

おまけ: 255.9999847412109375f はどこから来たか？

const int32 x = 0b0'10000110'11111111111111111111111;
const float y = *(float*)&x;
std::cout << std::fixed << std::setprecision( 32 ) << y << '\n';

Q2. 意図しない整数のオーバーフローと問題の遮蔽

出題のコードでは TV で 32-Bit 以上の整数型、例えば int32 を扱うと問題が起こりうる。例えば、次のような意図しない結果が得られる:

auto x = 
  Lerp
  ( std::numeric_limits< int32 >::max()
  , std::numeric_limits< int32 >::min()
  , 1.0f
  );

x には +2,147,483,647 相当の値が期待されるが、実際には -2,147,483,648 が得られる。 b - a は -2,147,483,648 - 2,147,483,647 = -1 となり、 ratio = 1.0f と乗算され -1.0f 、 これが -2,147,483,648 と乗算されると float 型の +2,147,483,648 となるが、 return 直前の (TV) により int32 型へキャストされる。 +2,147,483,648 を int32 へキャストすると -2,147,483,648 になる。

std::cout
  << std::bitset< 32 >( (int32) +2'147'483'647 ) << ' ' << std::showpos << (int32) +2'147'483'647 << '\n'
  << std::bitset< 32 >( (int32) +2'147'483'648 ) << ' '                 << (int32) +2'147'483'648 << '\n'
  << std::bitset< 32 >( (int32) +2'147'483'649 ) << ' '                 << (int32) +2'147'483'649 << '\n'
  ;

01111111111111111111111111111111 +2147483647
10000000000000000000000000000000 -2147483648
10000000000000000000000000000001 -2147483647

この問題は TV が整数型でも int8, uint8, int16, uint16 の場合には表面化しない。 a, b が int よりも狭いビット幅の整数型の場合に b - a は signed int ないし unsigned int で計算され、結果もその計算に使われた型のままとなる。

int16 a = +32'767;
int16 b = -32'768;
auto c = b - a;
std::cout << typeid( decltype( c ) ).name() << ' ' << std::to_string( c );

s s i -65535

b - a が int16 ではなく int32 で計算され、結果も int32 となるので -65535 を余裕をもって扱えてしまう。結果、 int 未満のビット幅の整数型を Lerp の TV で扱う限りは問題が表面化せず、一見すると任意の整数型に対して意図通りに動作するように見えてしまう。

もし、 Lerp に単体試験を網羅的に行うよう書いていたとしても、計算時間と "整数型" の概念を漠然と考えて、 int8 と uint8、あるいは int16 と uint16 で試験を実装してしまっていたなら、この問題はその単体試験が成功してしまうために何かが起きた際に原因を探すのが少し難しくなる事もあるかもしれない。

この問題が起こらない Lerp の実装例は次の通り:

  return (TV)( a * ( 1 - ratio ) + b * ratio );

数学的には等価な ( a + ratio * ( b - a ) ) と ( a * ( 1 - ratio ) + b * ratio ) も現実の計算機では等価とならない事もある。

Q3. 浮動小数点数型

IsSameAltitude の return では宇宙船Aと宇宙船Bの人工惑星Oからのそれぞれの高度を計算し、最終的に == で float 値を比較している。この処理では float 型の値が完全に一致する必要があるが、ストーリーのニーズからは float 値としての完全な一致ではなく、ゲームのプレイヤーが認識する高度の数値として感覚的に妥当な一致が必要です。

出題の実装ではミニマップにはほとんど何も表示されないか、ほんの一瞬（たいていのゲームではそれは 16 ms や 33 ms くらいの視認できるかどうかわからないくらいに）だけ、極稀に何か表示され…た気がする…、そんなような現実が発生し、ミニマップにまともに同じ高度の宇宙船が映らないというバグチケットが（おそらく）上がる事になるでしょう。

この系は float で十分に表現可能という条件から、宇宙船の高度は10進数にして高々7桁に収まるよう扱われます。 73.1 km のように。

しかし、 float の値としては、 例えソースコード上で 73.1f と書いたとしても、 std::cout << 73.1f; は 73.1 と表示されるとしても、実際の値は 73.099'998'474'121'093'750 相当となります。一般的な x64 アーキテクチャーやそれと互換の処理系での float は IEEE754/binary32 なので。

これはもちろん FVector の x, y, z それぞれについても、また operator- の結果も、 Length の計算の乗算や加算や std::sqrt の結果も全てが計算機にとってはこのバイナリー表現の世界が現実となります。

よって、 例え宇宙船Aと宇宙船Bの高度がプリントエフデバッグやログ出力などで確認すると同じ 73.1 と出力される状況でも、 IsSameAltitude は false を返す意図せず多発、というよりほとんどの状況となる。

float x = 73.1f;
float y = 73.1f + 4.0e-6f;
std::cout
  << x << '\n'
  << y << '\n' 
  << std::fixed << std::setprecision(32)
  << std::bitset< 32 >( *(int32*)&x ) << ' ' << x << '\n'
  << std::bitset< 32 >( *(int32*)&y ) << ' ' << y << '\n'
  ;

1.41421
73.1
73.1
01000010100100100011001100110011 73.09999847412109375000000000000000
01000010100100100011001100110100 73.10000610351562500000000000000000

この点を意識せずにデバッグすると原因の特定に少しだけ手間取るかもしれない。あるいは知識として把握していても、うっかり == で float の比較を書いてしまい予期しないバグを発生させてしまう事は新人や極度に疲弊した状態のプログラマーにはしばしば生じる事があるし、コンパイラーもこの意図は汲んでエラーや警告を出してくれるほど今のところは親切でないし、一般的な処理系の float は IEEE754/binary32 のままです。

この問題の対応としては、 IsNearlyEqual を用意し、また、 IsSameAltitude には許容誤差を明示的に実引数として渡せるようにする事です:

constexpr float default_error_tolerance = 1.0e-3f;

template < typename T >
bool IsNearlyEqual( const T a, const T b, const T error_tolerance = (T)default_error_tolerance )
{
  static_assert( std::is_floating_point< T >::value, "" );
  return std::abs( a - b ) <= error_tolerance;
}

bool IsSameAltitude
( const FVector& a
, const FVector& b
, const FVector& o
, const float error_tolerance = default_error_tolerance
)
{
  auto o_to_a = a - o;
  auto o_to_b = b - o;
  auto altitude_of_a = o_to_a.Length();
  auto altitude_of_b = o_to_b.Length();
  return IsNearlyEqual( altitude_of_a, altitude_of_b, error_tolerance );
}

Q4. 遅すぎた 2 の指数 [ 1, 2, 4, 8, 16, 32, .. ] の判定

template < typename T >
bool IsPowerOfTwo( const T in )
{
  static_assert( std::is_integral< T >::value, "" );
  return in == 0 ? false : ( ( in & ( in - 1 ) ) == 0 );
}

この実装は執筆時の wandbox で出題時と同様に簡易的に計測したところ Optimization=OFF で 5,041,452 #/sec （出題実装比 1.89 倍高速）、 Optimization=ON で 8,007,129 #/sec （出題実装比 1.94 倍高速）で動作した。

Q5. 正弦と余弦も速くしたい

象限を判定し最も精度良く計算可能な低角度域へシフトしつつマクローリン展開を必要な有効数字が十分に得られる程度計算する。

#define PI 3.14159265358979f
std::tuple< float, float > SinCosB( const float angle_in_radians )
{
  float quotient = angle_in_radians * 0.5f / PI;
  if ( angle_in_radians >= 0.0f )
    quotient = (float)( (int)( quotient + 0.5f ) );
  else
    quotient = (float)( (int)( quotient - 0.5f ) );
  float a = angle_in_radians - quotient * 2.0 * PI;
  float s = 0.0f;
  if ( a > 0.5f * PI )
  {
    a = PI - a;
    s = -1.0f;
  }
  else if ( a < -0.5f * PI )
  {
    a = -PI - a;
    s = -1.0f;
  }
  else
    s = +1.0f;
  float p = a * a;
  // [ 1 / 2!, 1 / 3!, 1 / 4!, .., 1 / 11! ]
  constexpr float f2 = 1.0f / 2.0f;
  constexpr float f3 = f2 / 3.0f;
  constexpr float f4 = f3 / 4.0f;
  constexpr float f5 = f4 / 5.0f;
  constexpr float f6 = f5 / 6.0f;
  constexpr float f7 = f6 / 7.0f;
  constexpr float f8 = f7 / 8.0f;
  constexpr float f9 = f8 / 9.0f;
  constexpr float f10 = f9 / 10.0f;
  constexpr float f11 = f10 / 11.0f;
  return std::make_tuple
  ( a * ( +1.0f + p * ( -f3 + p * ( +f5 + p * ( -f7 + p * ( +f9 * p * ( -f11 ) ) ) ) ) )
  , s * ( 1.0f + p * ( -f2 + p * ( +f4 + p * ( -f6 + p * ( +f8 + p * ( -f10  ) ) ) ) ) )
  );
}

この実装は執筆時の wandbox で出題時と同様に簡易的に計測したところ Optimization=OFF で 2,898,145 #/sec （出題実装比 1.06 倍高速）、 Optimization=ON で 5,914,322 #/sec （出題実装比 1.45 倍高速）で動作した。

（追記: 2018-12-12）

↑めでたしめでたし、と思われたが、実は↑の「こたえ」のコードにも新たなバグが埋め込まれてしまっていました。

CPP-QUIZ Q5 のこたえ、マクローリン展開による sin/cos の高速化に潜んでいた大きな誤差を生じるバグ - C++ ときどき ごはん、わりとてぃーぶれいく☆

と、いうわけで、改めまして、正しくは、

#define PI 3.14159265358979f
std::tuple< float, float > SinCosB( const float angle_in_radians )
{
  float quotient = angle_in_radians * 0.5f / PI;
  if ( angle_in_radians >= 0.0f )
    quotient = (float)( (int)( quotient + 0.5f ) );
  else
    quotient = (float)( (int)( quotient - 0.5f ) );
  float a = angle_in_radians - quotient * 2.0 * PI;
  float s = 0.0f;
  if ( a > 0.5f * PI )
  {
    a = PI - a;
    s = -1.0f;
  }
  else if ( a < -0.5f * PI )
  {
    a = -PI - a;
    s = -1.0f;
  }
  else
    s = +1.0f;
  float p = a * a;
  // [ 1 / 2!, 1 / 3!, 1 / 4!, .., 1 / 11! ]
  constexpr float f2 = 1.0f / 2.0f;
  constexpr float f3 = f2 / 3.0f;
  constexpr float f4 = f3 / 4.0f;
  constexpr float f5 = f4 / 5.0f;
  constexpr float f6 = f5 / 6.0f;
  constexpr float f7 = f6 / 7.0f;
  constexpr float f8 = f7 / 8.0f;
  constexpr float f9 = f8 / 9.0f;
  constexpr float f10 = f9 / 10.0f;
  constexpr float f11 = f10 / 11.0f;
  return std::make_tuple
  ( a * ( +1.0f + p * ( -f3 + p * ( +f5 + p * ( -f7 + p * ( +f9 + p * ( -f11 ) ) ) ) ) )
  , s * ( 1.0f + p * ( -f2 + p * ( +f4 + p * ( -f6 + p * ( +f8 + p * ( -f10  ) ) ) ) ) )
  );
}

となります😃 こちらも @nekketsuuu さんからのご指摘が基でバグを発見できました。ありがとうございます😃

Q6. ニアリーイコール、再び

出題の実装には2つの問題が潜んでいる。

1. a=π/2（ =90° ）, b=2π+π/2 （ =360°+90°=450°） に対して return が false となるが、一般的には何周多く回っても結果的に示す角度は等価なため、本来は true が返ると期待される。
2. a=0（ =0° ）, b=2π - 1.0e-4f （ =360°-1.0e-4=359.99 ）は標準の許容誤差 e=1.0e-3f 範囲内の角度で周回もしていないが abs(a-b) は abs( 0 - (2π-1.0e-4f) ) となり、 -2π+1.0e-4f < e は false と評価されてしまう。本来は true が返ると期待される。

対応として、 a と b の角度の差を正しく評価する FindDeltaAngleRadians を定義する:

#define PI 3.14159265358979f
/// 弧度法単位の角度の差を計算
/// @param a angle of A in radians
/// @param b angle of B in radians
/// @return a と b の角度の差
float FindDeltaAngleInRadians( const float a, const float b )
{
  auto d = b - a;
  
  if ( ! std::isfinite( d ) )
    return d;
  
  while ( d > PI )
    d -= 2 * PI;
  while ( d < -PI )
    d += 2 * PI;
  
  return d;
}

それから、 IsNearlyEqualAngleInRadians の角度の差の計算に FindDeltaAngleInRadians を用いる:

bool IsNearlyEqualAngleInRadians( const float a, const float b, const float e = 1.0e-3f )
{
  return std::abs( FindDeltaAngleInRadians( a, b ) ) < e;
}

(2017-12-06 追記)

クイズの作成時に題意としては考慮していませんでしたが、後日「こんな答えも含まれるのでは？」と鋭いご指摘を頂き、確かにこのクイズの他の設問で扱っているような出題と答えからするとそれもクイズの答えの1つとして記述した方が自然と私も思い、頂いたご指摘の Tweet を引用する形で追記・紹介いたします😃

お、鋭いですねー😃 楽しんで頂けて嬉しいです。さて、ご指摘の件は今回クイズとしてQ6を用意した時点では出題者の意図としては考慮しませんでしたが、このクイズの他の設問からはこれも答え編で触れた方が自然ですね。Tweet引用で内挿しておきます。ありがとうございます。

— Usagi Ito (@USAGI_WRP) 2017年12月5日

ぴんと来ない方も Tweet に添付して頂いた wandbox で状況を確認できるコード を見ると、たいへん分かり易いと思います。ありがとうございます。

Q7. 人間が読みやすい "数値" にしたかった

出題の実装では次の3つの問題が発生する。

1. int32 型の値は 999,999,999 を超える値も取り得る。例えば 1234567890 をこの関数の実引数として与えると、出力は 1234,567,890 となり期待される桁区切りが足りずに、奇妙で読み難い、意図しない出力が得られる。
2. int32 型の値は負の値を取り得る。例えば -123456789 を与えると 出力は -123456789 となり桁区切りがまったく行われず、意図しない出力が得られる。
3. 世界は広い。日本やアメリカ合衆国の文化圏では数値の桁区切りが3桁ごとに , （カンマ）で行われる慣例が一般化しているが、例えばドイツやイタリアでは . （ドット）を桁区切りに使い、 , を小数点に使う。フランスやロシアでは （空白文字）を桁区切りに使い、 , を小数点に使う。スイスでは ' （シングルクォート）を桁区切りに使い、 , を小数点に使う。 int32 は整数型なので小数点の考慮は必要無いが、少なくとも製品が使用される文化圏を限定できない場合には桁区切り文字は , とは限らない事を考慮する必要がある。

(1) int32 型の値は [ -2,147,483,648 .. +2,147,483,647 ] を取り得るので、桁区切りは3つ目まで考慮した実装とする必要がある。

// 出題の実装に倣うなら、3回目の桁区切りを追加すれば
// (1) については意図通りの結果が得られるようになる。
  if ( in > 999'999'999 )
  {
    out = ',' + buffer.substr( buffer.size() - 3, 3 ) + out;
    buffer.resize( buffer.size() - 3 );
  }

(2) は負の値も桁区切りの判定を正常に行えるよう修正する必要がある。

// 出題の実装を最小限の改修で済ませるならば if の条件に負の値の判定も || で盛り込めば 
// (2) については意図通りの結果が得られるようになる。
  if ( in > 999 || in < -999 )

(3) は , を任意の文字へ変更可能に対応する。

std::string FormatIntToHumanReadable( const int32 in, const char splitter = ',' )
{
  // 中略 //
    out = splitter + buffer.substr( buffer.size() - 3, 3 );
  // 後略 //

(1), (2), (3) の修正を加え、ついでに if の3連続を while 1つに整理すると次の実装となる:

std::string FormatIntToHumanReadable( const int32 in, const char splitter = ',' )
{
  auto buffer = std::to_string( in );
  std::string out;

  while ( buffer.size() > 3 && buffer[ buffer.size() - 4 ] != '-' )
  {
    out = splitter + buffer.substr( buffer.size() - 3, 3 ) + out;
    buffer.resize( buffer.size() - 3 );
  }

  out = buffer + out;

  return out;
}

なお、このコードは "文字列処理" としては分かり易い側面もありますが、 C++ の実装としてはやや "暢気" な実装です。また、そもそも std::ostream へロケールを設定して << するだけでよしなに桁区切りしてくれるとか、そういった指摘もあるかと思います。それらは出題の「この実装ではどのような問題が起こり得るだろうか」の題意からは外れますので、興味に応じてエクストラにお楽しみになられたら良いと思います。

Q8. 整数型、再びトラブる

同じ原因に起因するバグが2箇所に潜んでいる。

1 つは、 float intensity = std::abs( in ) の結果は常に正の数とはならない。例えば、 T=int32, in=std::numeric_limits<int32>::min() が関数へ与えられた場合、 intensity は -2147483648.0f となる。 std::abs( in ) は int 型以上のビット幅の符号付きの整数型に対して、その最小値が与えられた場合、プログラマーが数学的に意図する絶対値は返せず、結果的には最小値がそのまま返ってくる。これを float 型へ変換したものが intensity へ保持され、結果、負の値を T 型の最大値である正の値で除算した値が return され、意図しない -1 が得られてしまう事になる。

この関数の仕様は return が [ 0 .. 1 ] なので、 -1 が返る想定をしない何らかの処理がとんでもない動作を引き起こす可能性がある。

もう 1 つは、 intensity / std::numeric_limits< T >::max() の結果が T 型の最大値と最小値の数学的な意味での絶対値が 1 だけ違うために 0 を中心として符号を除去した intensity を T 型の最大値で除算してしまうと、元の値が負の値だった場合に意図した意図した値域の写像とならない。例えば、 T=int8, in=-128 の場合、本来は 1.0f が期待されるが、この出題の実装では 1.007874011993408203125f となる。

この関数の実装上の問題は、以上の2点が起こり得る事。何れも符号付きの整数型の最小値に起因する。

また、この関数の仕様についても1点明確ではない問題が潜んでいる。 T が符号付き整数型の場合に、負の強度の最大値は T 型の 最小値 として T型の値の正、負それぞれの全域を強度とするのか、負の値については 最小値 + 1 までを正常値として正、負の値は符号が異なっても同じ強度とするのかが明確でない。

// 正は正、負は負でそれぞれの数学的な意味での絶対値が最大の値を最大の強度とする場合の修正例
template < typename T >
float GetIntensityToUNormFloat( const T in )
{
  return in >= 0
    ? in / (float)std::numeric_limits< T >::max()
    : in / (float)std::numeric_limits< T >::min()
    ;
}

// 正も負も数学的な意味での絶対値が同じ値は同じ強度とし、
// 符号付き整数型の最小値は不正値とする場合の修正例
template < typename T >
float GetIntensityToUNormFloat( const T in )
{
  check( std::is_unsigned< T >::value || in > std::numeric_limits< T >::min() );
  return std::abs( in / (float)std::numeric_limits< T >::max() );
}

Q9. G.C.D.

5 つの実装例と簡易的な速度評価結果を紹介します。

(1) はじめに Unreal Engine 4.18 の実装を出題に併せて調整した実装:

template < typename T >
T GreatestCommonDivisor( T a, T b )
{
  static_assert( std::is_integral< T >::value, "" );
  check( a >= 0 )
  check( b >= 0 )
  while ( b != 0 )
  {
    T t = b;
    b = a % b;
    a = t;
  }
  return a;
}

(2) while タイプの亜種:

template < typename T >
T GreatestCommonDivisor( T a, T b )
{
  static_assert( std::is_integral< T >::value, "" );
  check( a >= 0 )
  check( b >= 0 )
  if ( b )
    while ( ( a %= b ) && ( b %= a ) );
  return ( a + b );
}

(3) std::function recursive タイプ:

template < typename T >
T GreatestCommonDivisor( T a, T b )
{
  static_assert( std::is_integral< T >::value, "" );
  check( a >= 0 )
  check( b >= 0 )
  std::function< T( T, T ) > f;
  f = [&f]( auto a, auto b ){ return ( b != 0 ) ? f( b, a % b ) : a; };
  return f( a, b );
}

(4) function recursive タイプ:

template < typename T >
T GreatestCommonDivisor_Impl( T a, T b )
{ return ( b != 0 ) ? GreatestCommonDivisor_Impl( b, a % b ) : a; }

template < typename T >
T GreatestCommonDivisor( T a, T b )
{
  static_assert( std::is_integral< T >::value, "" );
  check( a >= 0 )
  check( b >= 0 )
  return GreatestCommonDivisor_Impl( a, b );
}

(5) C++17er:

template < typename T >
T GreatestCommonDivisor( T a, T b )
{
  static_assert( std::is_integral< T >::value, "" );
  check( a >= 0 )
  check( b >= 0 )
  // std::gcd がツールチェインごとにどのようなソースで実装されているかはさておき
  return std::gcd( a, b );
}

それぞれの簡易的な実行速度評価結果:

実装パターン	Optimization	#/sec
(1) UE4 while	OFF	1,137,116
(1) UE4 while	ON	1,495,536
(2) while 亜種	OFF	1,121,166
(2) while 亜種	ON	1,503,565
(3) std::function recursive	OFF	313,537
(3) std::function recursive	ON	1,236,848
(4) function recursive	OFF	954,043
(4) function recursive	ON	1,534,016
(5) C++17er	OFF	960,103
(5) C++17er	ON	1,481,497

Q10. L.C.M.

2 つの実装例と簡易的な速度評価結果を紹介します。

(1) はじめに Unreal Engine 4.18 の実装を出題に併せて調整した実装:

template < typename T >
T LeastCommonMultiplier( T a, T b )
{
  static_assert( std::is_integral< T >::value, "" );
  check( a >= 0 )
  check( b >= 0 )
  // Q9 G.C.D. の実装を使う
  T gcd = GreatestCommonDivisor( a, b );
  return gcd == 0 ? 0 : ( a / gcd ) * b;
}

(2) C++17er:

template < typename T >
T LeastCommonMultiplier( T a, T b )
{
  static_assert( std::is_integral< T >::value, "" );
  check( a >= 0 )
  check( b >= 0 )
  // std::gcd がツールチェインごとにどのようなソースで実装されているかはさておき
  return std::lcd( a, b );
}

実装パターン	Optimization	#/sec
(1) UE4 a / GCD * b	OFF	1,126,620
(1) UE4 a / GCD * b	ON	1,456,371
(2) C++17er	OFF	889,107
(2) C++17er	ON	1,418,852

おまけ: G.C.D. & L.C.M. の wandbox での簡易的な実行速度評価に使用したコード

G.C.D., L.C.M. 以外でも概ねは同様のコードで簡易的な実行速度評価を行っています。

#include <iostream>
#include <limits>
#include <type_traits>
#include <chrono>
#include <random>
#include <array>
#include <algorithm>
#include <cassert>
#define check(X) { assert(X); }
using int8 = std::int8_t;
using uint8 = std::uint8_t;
using int16 = std::int16_t;
using uint16 = std::uint16_t;
using int32 = std::int32_t;
using uint32 = std::uint32_t;
using int64 = std::int64_t;
using uint64 = std::uint64_t;

// Q9 のこたえを入れた header file
#include "GCD.h"
// Q10 のこたえを入れた header file
#include "LCM.h"

int main()
{
  std::array< uint64, 1024 > numbers;
  std::mt19937_64 rng( 0 );
  std::generate( numbers.begin(), numbers.end(), rng );
  double x = 0.0f;
  size_t count = 0;
  auto ta = std::chrono::steady_clock::now();
  while ( std::chrono::steady_clock::now() <= ta + std::chrono::seconds( 10 ) )
  {
    auto a = numbers[   count % numbers.size() ];
    auto b = numbers[ ++count % numbers.size() ];
    x += GreatestCommonDivisor( a, b );
    //x += LeastCommonMultiplier( a, b );
  }
  std::cout << ( count / 10 ) << ' ' << x << '\n';
}