横軸が無限大でγに収束します。

実際に計算するため、またMPFRを使ってみます。

MPFR Library

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <mpfr.h>

int main()
{
  mpfr_t r, one, work, nn, sum;

  mpfr_init2(r, 100);
  mpfr_init2(one, 100);
  mpfr_init2(nn, 100);
  mpfr_init2(work, 100);
  mpfr_init2(sum, 100);

  unsigned int x = 3000000000;
  mpfr_set_ui(one, 1, MPFR_RNDD);
  mpfr_set_ui(sum, 0, MPFR_RNDD);
  for(unsigned int i=1;i<=x;i++){
    mpfr_set_ui(nn, i, MPFR_RNDD);
    mpfr_div(work, one, nn, MPFR_RNDD);
    mpfr_add(sum, sum, work, MPFR_RNDD);
  }
  mpfr_log_ui(work, x, MPFR_RNDD);
  mpfr_sub(r, sum, work, MPFR_RNDD);
  mpfr_out_str(stdout, 10, 0, r, MPFR_RNDD);
  mpfr_clear(r);
}

Wikipedia
https://ja.wikipedia.org/wiki/%E3%82%AA%E3%82%A4%E3%83%A9%E3%83%BC%E3%81%AE%E5%AE%9A%E6%95%B0
では、

およそ0.57721 56649 01532 86060 65120 90082 40243 10421 59335 93992 35988 05767 23488 48677 26777 66467 09369 47063 29174 67495…

とありますが、このブログラムの結果は

5.7721566506819952722608849107166e-1

real 3m8.076s
user 3m7.427s
sys 0m0.585s

となり、そこそこ標準的なPCでこれだけ時間をかけて小数点第8位までしか一致しませんでした。
なかなか収束が遅いようです。
オイラーが6位まで計算したとのことですが、PCのない時代にすごいです。

現代数学のホットな研究対象であるゼータ関数論との関係も深いこの定数ですが、この分野に何か惹きつけられるものがあり勉強中です。こういうときプログラム言語というのは理解を深めるためにとても役に立つと実感します。
最近のPCの性能のおかげで、数学のミステリーに迫る楽しみを味わえるのは、いい時代です。

MPFR Library

通常、複素数演算は実部と虚部を分けて計算すればいいのですが、累乗に複素数を使うようなイメージがしにくい計算では、このようなライブラリが有効なのではと思います。
そこで、下記はこれを使う有名なオイラーの公式ですが、これをライブラリで計算してみることにします。

xにπを代入すると、世界でもっとも美しいといわれる公式になります。

両辺とも直接的な関係はないのですが、それぞれテーラー展開という無限級数にして計算すると一致するという不思議な関係です。
虚数乗というなんと解釈していいかわからないようなものでも、数学の世界では表現できるという、ここが面白いところです。

環境: Ubuntu 18.04

インストール・ビルド

sudo apt install libmpc-dev
cc -o euler euler.c -lmpc -lmpfr -lgmp

euler.c

#include <stdio.h>
#include <mpc.h>

int main() {

  mpfr_t pi, sn;
  mpc_t work,  e, re1, im1, ipi;
  
  mpfr_init2(pi, 128);
  mpfr_init2(sn, 128);
  
  mpc_init2(work, 128);
  mpc_init2(e, 128);
  mpc_init2(ipi, 128);
  
   mpc_init2(re1, 128);
   mpc_init2(im1, 128);
   
  mpc_set_ui_ui(re1,  1,  0, MPC_RNDNN);
  mpc_set_ui_ui(im1,  0,  1, MPC_RNDNN);
  
  mpc_exp(e, re1, MPC_RNDDD);
  mpfr_const_pi(pi, MPFR_RNDD);

  mpc_out_str(stdout, 10, 0, e, MPC_RNDDD);
  printf("\n");
  mpfr_out_str(stdout, 10, 0, pi, MPFR_RNDD);
  printf("\n");
   
  mpc_mul_fr(ipi, im1, pi, MPC_RNDDD);
  mpc_pow(work, e, ipi, MPC_RNDNN);
  mpc_out_str(stdout, 10, 0, work, MPC_RNDDD);
  printf("\n");
  
  mpfr_sin(sn, pi, MPFR_RNDD);
  mpfr_out_str(stdout, 10, 0, sn, MPFR_RNDD);
  printf("\n");
  
  mpc_clear(work);
  ;
  return 0;
}

結果

(2.718281828459045235360287471352662497747 0)
3.141592653589793238462643383279502884195
(-1.000000000000000000000000000000000000000 1.338432504205046150309804279711632522358e-38)
1.883041077660785116745909548456034940273e-39

まず演算に使う、虚数単位(im1)、ネイピア数e、と円周率piを用意します。(eにmpfr_const_eみたいなものがない?)
re1はeの値を導き出すために使います。
結果の3行目が、求めたい値です。実部は-1で問題ないですが、虚部は0にならずかなり小さい数となっています。sin(π)も単独で演算してみましたが、一致はしませんでした。
コンピュータの計算では当然誤差がでてくるため、その範囲としますが、一桁違うのはちょっと気になります。

複素数型なのかどうかヘッダーファイルにある関数の型とにらめっこして慎重に型変換をしていきました。
なかなか手間ですが、やっていることがわかりやすく応用を考えたとき、各種ラッパーライブラリよりも使いやすいと思っています。

その前に現在の標準的な64bitコンピュータで扱える演算精度を確認したいと思います。

#include<stdio.h>

int main()
{
	long n1 = 0xffffffffffffffff;
	printf("%lu : %ld\n", n1, n1);

	double d1 = 1.0, d2 = 3.0, d3;
	d3 = d1 / d2;
	printf("%.20lf\n", d3);
}

18446744073709551615 : -1
0.33333333333333331483

64bit unsignedの整数で20桁、浮動小数仮数部52bitで16桁の有効桁を確認できます。
実用的にはこれで十分な桁数がありますが、数学の世界では無限に必要となるようです。

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <mpfr.h>

int main()
{
  mpfr_t n1, n2;

  mpfr_init2(n1, 200);
  mpfr_init2(n2, 200);
	                                               
  mpfr_set_str(n1, "1111111111111111111111111", 10, MPFR_RNDD); 

  mpfr_mul(n2, n1, n1, MPFR_RNDD);
  mpfr_out_str(stdout, 10, 0, n2, MPFR_RNDD);                
  printf("\n");

  mpfr_clear(n1); 
  mpfr_clear(n2); 

  return 0;
}

1.2345679012345679012345676543209876543209876543210000000000000e48

200bitの精度の整数を掛け算した結果です。桁数が多いとC言語の数字として表現できないので文字列にして初期化しています。

次に円周率です。

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <mpfr.h>

int main()
{
  mpfr_t pi;

  mpfr_init2(pi, 2000);                                                  

  mpfr_const_pi(pi, MPFR_RNDD);
  mpfr_out_str(stdout, 10, 0, pi, MPFR_RNDD);
  printf("\n");

  mpfr_clear(pi);                                    
}

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405131995

2000bitの精度で初期化し、πを求める関数があるのでそれを利用しました。
このくらいだとすぐに結果を出してくれます。(もっといけますが、結果が長すぎるで・・、素数判定は64bitでかなり遅いので、これを使うレベルでない)

たまにはPCもこういう使い方すると、スペックの進化も実感することができます。

環境: Ubuntu 18.04

https://www.mpfr.org/