破綻する区間演算の事例

　精度保証計算（という言葉はあまり好きではないのだけれど）の事例も紹介しようかと，ARB（半径）とMPFI（区間）を使ってロジスティック写像\(x_{i+1} := 4x_i(1-x_i)\)の計算をさせてみたところ，桁数が足りないと破綻することが分かった（今更）。つーことで，128bits計算の事例（上）と256bitsの事例（下）をば。

　有効桁数がなくなった途端に区間や半径が急激にデカくなって破綻していることが一目でわかるのはいいとして，精度はともかく計算は最後まで実行したい時には使えん。とりあえず，こういう桁落ちが酷い事例には区間演算だろうとそうでなかろうと桁数増やしてやんなさいという結論になるな。

　あと，実係数2次方程式の事例も作ってみたが，判別式で実数解と複素数解の区別をするif文をどうしようとか，実用に当たってはやはりメンドクサイことを考えないといけないのでその辺をどうかして工夫すると研究になる訳で，こっから先は精度保証の方々にお任せしようっと。まぁ「最初から複素係数でやれ」という結論になるんだろうけど，それはARBの演習問題にしておけばいいか。

　つーことで，使用したプログラム(logistic_arb.c, logisitc_mpfi.c)を最後に掲載しておく。使用時にはARB, MPFIの他，MPFR/GMPが必要。

// logistic 写像
// ARB版
#include <stdio.h>
#include "arb.h"

int main()
{
 int i;
 slong prec;
 arb_t x[102];

 // 初期化
 for(i = 0; i < 102; i++)
  arb_init(x[i]);

 printf("prec(bits) = "); scanf("%ld", &prec);

 // 初期値
 //x[0] = 0.7501;
 arb_set_str(x[0], "0.7501", prec);

 for(i = 0; i <= 90; i++)
 {
  if((i % 10) == 0)
  {
   printf("%5d, ", i);
   arb_printd(x[i], 17);
   printf("\n");
  }

  //x[i + 1] = 4 * x[i] * (1 - x[i]);
  arb_sub_ui(x[i + 1], x[i], 1UL, prec);
  arb_neg(x[i + 1], x[i + 1]);
  arb_mul(x[i + 1], x[i + 1], x[i], prec);
  arb_mul_ui(x[i + 1], x[i + 1], 4UL, prec);

 }

 // 消去
 for(i = 0; i < 102; i++)
  arb_clear(x[i]);

 return 0;
}

// logistic 写像
// MPFR & MPFI版
#include <stdio.h>
#include "mpfr.h"
#include "mpfi.h"

int main()
{
 int i;
 unsigned long prec;
 mpfi_t x[102];
 mpfr_t relerr;

 printf("prec(bits) = "); scanf("%ld", &prec);
 mpfr_set_default_prec(prec);

 // 初期化
 mpfr_init(relerr);
 for(i = 0; i < 102; i++)
  mpfi_init(x[i]);

 // 初期値
 //x[0] = 0.7501;
 mpfi_set_str(x[0], "0.7501", 10);

 for(i = 0; i <= 90; i++)
 {
  if((i % 10) == 0)
  {
   printf("%5d, ", i);
   mpfi_out_str(stdout, 10, 17, x[i]);
   mpfi_diam(relerr, x[i]);
   mpfr_printf("%10.3RNe\n", relerr);
  }

  //x[i + 1] = 4 * x[i] * (1 - x[i]);
  mpfi_ui_sub(x[i + 1], 1UL, x[i]);
  mpfi_mul(x[i + 1], x[i + 1], x[i]);
  mpfi_mul_ui(x[i + 1], x[i + 1], 4UL);

 }

 // 消去
 mpfr_clear(relerr);
 for(i = 0; i < 102; i++)
  mpfi_clear(x[i]);

 return 0;
}

T.Kouya