floor関数を実装する

2023-01-15
#algorithm #math

C言語の数学関数の floor を自前で実装するにはどうするか。

floor関数

floor 関数は日本語だと「床関数」、与えられた浮動小数点数を超えない最大の整数を返す関数で、math.hで宣言されている：

double floor(double x);

実装に入る前にdoubleがプロセッサ上でどのように表現されているかについて。

浮動小数点数の内部表現(double)

浮動小数点数をコンピュータ内部でどうやって扱っているかというと、大抵の場合IEEE 754という形式で扱っている。 double の場合、全体で64ビット長で、符号部1、指数部11、仮数部52となる。

  • 仮数部は0.5以上1.0未満に正規化され、最上位ビットが省略される
  • 指数部は1022(= 0x3fe)底上げされている
  • 0.0や±InfinityやNaNは特別扱い

例えば1.0は、0x3ff_0000000000000（0.5 x 2 ^ 1）となる。

floor関数の実装方法

double を扱う関数になるが、浮動小数点数演算ではなくビット演算で処理する。 でひとまず正の場合のみを考える。

仮数部に対する小数点の位置は ((52+1) - 指数) となる。 実際には指数は1022底上げされているので、指数部をeとすると(52+1) - (e - 1022)となる。

それで小数点の位置によって場合わけ：

  • 仮数部にかかる場合(1022 < e <= 1074)：小数部をビットマスクで落とす
  • 仮数部より上になる場合(e <= 1022)：値が小さくて元々整数部がないので、0.0を返す
  • 仮数部より下になる場合(e > 1074)：値が大きくて元々小数部がないので、元の値を返す

入力が0.0の場合、e = 0で条件２となるが、一応-0.0を考慮して元の値を返すことにする。 入力がNaNまたは±Infinityの場合、e = 0x7ffで条件３で元の値が返されるのでOK。

入力が負の場合

値が負の場合は、符号を無視して仮数部を見た場合に小数となる部分が0じゃない場合は切り上げてやる。

ソース

というわけで、実装したソースは floor.c

#include <stdint.h>

#define FRAC_BIT   (52)
#define EXPO_POS   FRAC_BIT
#define EXPO_BIT   (11)
#define EXPO_BIAS  (1022)
#define GET_BIASED_EXPO(hex)  (((int)((hex) >> EXPO_POS)) & ((1 << EXPO_BIT) - 1))

double myfloor(double x) {
  int64_t q = *(int64_t*)&x;
  int e = GET_BIASED_EXPO(q);
  if (e <= EXPO_BIAS + FRAC_BIT && e != 0) {
    if (e <= EXPO_BIAS)
      return q >= 0 ? 0.0 : -1.0;

    int64_t one = (int64_t)1 << ((EXPO_BIAS + FRAC_BIT + 1) - e);
    if (q > 0) {
      int64_t r = q & -one;
      return *(double*)&r;
    } else {
      int64_t frac = q & (one - 1);
      if (frac != 0) {
        int64_t r = (q & -one) + one;
        return *(double*)&r;
      }
    }
  }
  return x;
}

#include <stdio.h>
int main() {
  printf("%f\n", myfloor(1.23));
}

ceil関数

ceil関数も、floorの正負を逆にして同様に：ceil.c

雑感

doubleの値のビット表現を取り出すためにC言語ではint64_tを指すポインタ経由で取り出す必要があって、それがいいのかどうか。

それなりによく使いそうなもんだしそんなに回路複雑にならなそうだけど、プロセッサの命令として実装されてないだろうか？

関連記事

新着記事