8 gcc の __float128

(2016/2/25)

「__float128 の利用法」 がわかりやすい。

gcc の version 4.6 以降では __float128, __complex128 というのが使えるとか。 IEEE 754 2008 で導入された binary128 だとか。

/*
 * testfloat128.c -- __float128 のテスト
 *   gcc testfloat128.c -lquadmath
 */

#include <stdio.h>
#include <quadmath.h>
#define MAXN (100)

int main(void)
{
  int n;
  char message[128];
  __float128 e, newe, x;
  x = 1.0Q; e = x;
  for (n = 1; n <= MAXN; n++) {
    x /= n;
    newe = e + x;
    quadmath_snprintf(message, sizeof(message), "%.40Qf", e);
    puts(message);
    if (e == newe)
      break;
    else
      e = newe;
  }
  printf("2.7182818284590452353602874713526624977572470937000 --- N[E,50]\n");
  printf("10進法で34桁程度の精度があるというけれど、本当？\n");
  return 0;
}

[chronos:~/work] mk% gcc testfloat128.c -lquadmath
[chronos:~/work] mk% ./a.out
1.0000000000000000000000000000000000000000
2.0000000000000000000000000000000000000000
2.5000000000000000000000000000000000000000
2.6666666666666666666666666666666665382713
2.7083333333333333333333333333333330765427
2.7166666666666666666666666666666663841969
2.7180555555555555555555555555555554613990
(途中略)
2.7182818284590452353602874713492665146805
2.7182818284590452353602874713525463733758
2.7182818284590452353602874713526596180565
2.7182818284590452353602874713526634699164
2.7182818284590452353602874713526624977572470937000 --- N[E,50]
10進法で34桁程度の精度があるというけれど、本当？
[chronos:~/work] mk%

なるほど、なるほど。

(2016/12/1)

__float128 と double で、 1000次正方行列同士の乗算にどれくらい時間がかかるか。

/*
 * gcc -O3 -funroll-loops -o testfloat128 testfloat128.c -I/usr/local/include -L/usr/local/lib -lmatrix -lm
 */

#include <stdio.h>
#include <math.h>
#include <matrix.h>
#include <quadmath.h>

int main(void)
{
  int i, j, k, n;
  char message[128];
  __float128 norm;
  matrix128 a,b,c;
  n = 1000;
  a = new_matrix128(n,n);
  b = new_matrix128(n,n);
  c = new_matrix128(n,n);
  if (a == NULL || b == NULL || c == NULL) {
    fprintf(stderr, "cannot allocate\n");
    exit(1);
  }
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++) {
      a[i][j] = i + j;
      b[i][j] = i - j;
    }
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++) {
      c[i][j] = 0;
      for (k = 0; k < n; k++)
	c[i][j] += a[i][k] * b[k][j];
    }
  norm = 0;
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      norm += fabsq(c[i][j]);
  quadmath_snprintf(message, sizeof(message), "%.7Qf", norm);
  printf("%s\n", message);
  return 0;
}

[chronos:~/work] mk% ccmg testfloat128.c -lquadmath ; time ./testfloat128
257388424419000.0000000
122.528u 0.179s 2:02.73 99.9%	0+0k 0+0io 2pf+0w

double バージョンを作ってみたが、あまりスピード出ない (5〜10秒)。MATLABだと、掛け算 $1$ 秒もかからない。