うだうだWeblog2

2018-01-082018-01-09

内田春菊「がんまんが　～私たちは大病している～」ぶんか社

[ Amazon ] ISBN 978-4-8211-3566-0, \1000 + TAX

　内田春菊のエッセイ漫画に外れはない。流麗なペンタッチが白い画面に艶めかしい柔肌を描き出す。それでいて，内容はシビアで容赦ない。ダメンズを次々に恋人にして子供を作り，ダメと分かった時点で放り出す。内田本人の生命力・経済力が図抜けている分，割り切り方はスッパリしていて小気味良い。ワシはその力強い生き方を描いたエッセイ漫画である「私たちは繁殖している」は読んではいたものの，途中で脱落していたが，続刊が次々に出ているのを見るにつけ，このまま佐藤愛子か瀬戸内寂聴のように「完成」していくのであろうなと思っていたら，大腸がんになってしまったという。その顛末を，手術直後まで描いているのが本書である。その後については連載中なのでまた続刊が出るとのことである。

　20年以上前，オストメイトの方と仕事をしていたことがある。当時は人工肛門をサポートする技術が未熟だったらしく，アンモニアの匂いが漏れてはた目から見てても気の毒であった。さすがに最近は匂い消しフィルター付きの排便袋になっているらしく，昔よりずっと活動しやすいようである。その辺の事情については，内田のインタビュー記事（（１），（２），（３），（４））を読むとよい。本書を買って読む気のない人にはこちらを勧める。

　しかし，ワシとしては，本書まるごと読まれることをお勧めしたい。ガン発覚のいきさつから，治療方針が決まるまでの医療機関とのゴタゴタ，そして最終的に人工肛門を形成するに至るまでの一部始終が流麗なタッチで描かれていることが，がん治療というものの真実を知るいい資料になっていると感じるからである。

　内田春菊は我慢しない。言いたいことは言うし，やりたいことはやる。それ故に医者との軋轢もあったりするが，そのおかげで信頼できる医師とはしっかりしたコミュニケーションが取れており，家族（子息）からサポートも篤い。科学を無視した耳障りがいいだけの暴論がSNSに蔓延る昨今，客観的な科学の知識に基づき，自分の意志の表明を躊躇なく行うことの重要性を知らしめる本書は，がんを含めた病とともに生きていくワシら中高年にとっては，良い資料となるに違いないのである。

2018-01-052018-01-26

メルセンヌ素数で”FLOPS vs. Memory Bandwidth”を体験する

　50番目のメルセンヌ素数が見つかった，つーか，素数であることの検証が済んだ，という記事に呼応して下記のようなTweetをした。

見間違えたので再TW。古いXeonマシンでもGNU MP(GMP)のmpz_ui_pow_ui & mpz_sub_ui関数で計算して30秒足らず。ファイルDLするより計算した方が速い。素数検証の方がずっと大変。＞「史上最大の素数」約2年ぶりに更新最長だったものより約100万桁大 #ldnews https://t.co/lO9BwReyng pic.twitter.com/jEzQbfQhZu

— 幸谷智紀＠駿府(+髭) (@tkouya) 2018年1月5日

　これって，”FLOPS vs. Memory Bandwidth”に通じる話でもあるなと思ったので，ここでもう少し詳細にメモっておく。上記のTweetより少し手直ししたプログラムは下記の通り。

// Mersenne Prime Number Calculator by GNU MP
// Copyright (c) 2018 Tomonori Kouya
//
// Compile: cc mpz_mersennne.c -lgmp
//
#include <stdio.h>
#include <stdlib.h>

#include "gmp.h"

/* Mersenne Prime Number */
/* 2^p - 1

/* cf.) http://primes.utm.edu/mersenne/ */

int max_index = 50; // 2018-01-05
unsigned long known_mersenne_prime_p[] = {
	       2,        3,        5,        7,       13,       17,       19,       31,       61,       89, \
	     107,      127,      521,      607,     1279,     2203,     2281,     3217,     4253,     4423, \
	    9689,     9941,    11213,    19937,    21701,    23209,    44497,    86243,   110503,   132049, \
	  216091,   756839,   859433,  1257787,  1398269,  2976221,  3021377,  6972593, 13466917, 20996011, \
	24036583, 25964951, 30402457, 32582657, 37156667, 42643801, 43112609, 57885161, 74207281, 77232917 \
};

int main(void)
{
	mpz_t a;
	int index;

	mpz_init(a);

	// input
	fprintf(stderr, "2^p[index] - 1 = ? ");
	fprintf(stderr, "index(1 - %d) = ", max_index);
	scanf("%d", &index);
	if((index > max_index) || (index <= 0))
	{
		fprintf(stderr, "ERROR: index = %d is illegal!\n", index);
		return EXIT_FAILURE;
	}

	index--;
	printf("p[%d] = %ld\n", index + 1, known_mersenne_prime_p[index]);

	// a := 2^p - 1
	mpz_ui_pow_ui(a, 2UL, known_mersenne_prime_p[index]);
	//gmp_printf("2^%ld     = %Zd\n", known_mersenne_prime_p[index], a);
	mpz_sub_ui(a, a, 1UL);
	gmp_printf("2^%ld - 1 = %Zd\n", known_mersenne_prime_p[index], a);

	mpz_clear(a);

	return 0;
}

// Mersenne Prime Number Calculator by GNU MP

// Compile: cc mpz_mersennne.c -lgmp

#include <stdio.h>

#include <stdlib.h>

#include "gmp.h"

/* Mersenne Prime Number */

/* 2^p - 1

/* cf.) http://primes.utm.edu/mersenne/ */

int max_index = 50; // 2018-01-05

unsigned long known_mersenne_prime_p[] = {

2, 3, 5, 7, 13, 17, 19, 31, 61, 89, \

107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, \

9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, \

216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917, 20996011, \

24036583, 25964951, 30402457, 32582657, 37156667, 42643801, 43112609, 57885161, 74207281, 77232917 \

};

int main(void)

{

mpz_t a;

int index;

mpz_init(a);

// input

fprintf(stderr, "2^p[index] - 1 = ? ");

fprintf(stderr, "index(1 - %d) = ", max_index);

scanf("%d", &index);

if((index > max_index) || (index <= 0))

{

fprintf(stderr, "ERROR: index = %d is illegal!\n", index);

return EXIT_FAILURE;

}

index--;

printf("p[%d] = %ld\n", index + 1, known_mersenne_prime_p[index]);

// a := 2^p - 1

mpz_ui_pow_ui(a, 2UL, known_mersenne_prime_p[index]);

//gmp_printf("2^%ld = %Zd\n", known_mersenne_prime_p[index], a);

mpz_sub_ui(a, a, 1UL);

gmp_printf("2^%ld - 1 = %Zd\n", known_mersenne_prime_p[index], a);

mpz_clear(a);

return 0;

}

　Tweetに挙げたのは$2^p$(mpz_ui_pow_ui関数で計算)と$2^{p – 1}$(mpz_sub_ui関数で計算)の両方をファイルに書き出す奴だったので，後者のみ書き出すようにするとほぼ半分の時間で済む。

　10進数テキストで書き出すと大体45MBぐらいのファイルになる訳だが，ネット越しだとダウンロードする方が時間がかかる。100BASE-TXの環境でもイーブン。大体，計算そのものにかかる時間よりファイルに書き出す時間の方が問題で，更に言うと，メインメモリ内部のデータ移動の法に時間が費やされている。128MBぐらいのL3キャッシュがあれば，FLOPS(この場合は整数演算だからMIPSだけど)上げるよりもずっと高速に計算終了している筈で，さらに言えば，ファイルに書き出しせずにオンメモリで済むならもっと高速。つまりCPU内部の演算効率よりメモリ帯域(memory bandwidth)の方がよっぽど問題，というお話になるわけで。

　という自分用メモ。

[2018-01-26(Fri) 追記] 「多倍長数値計算入門（仮）」執筆開始につき，上記のCプログラムをJuliaスクリプトで書き直してみた。なーるほど，これは使いやすいや。

#  Mersenne Prime Number Calculator by GNU MP
#  Copyright (c) 2018 Tomonori Kouya
# 
#  Run: julia mersennne.jl 
# 

# Mersenne Prime Number
# 2^p - 1

# cf.) http:# primes.utm.edu/mersenne/

max_index = 50; #  2018-01-05
known_mersenne_prime_p  = [
	       2,        3,        5,        7,       13,       17,       19,       31,       61,       89,
	     107,      127,      521,      607,     1279,     2203,     2281,     3217,     4253,     4423,
	    9689,     9941,    11213,    19937,    21701,    23209,    44497,    86243,   110503,   132049,
	  216091,   756839,   859433,  1257787,  1398269,  2976221,  3021377,  6972593, 13466917, 20996011,
	24036583, 25964951, 30402457, 32582657, 37156667, 42643801, 43112609, 57885161, 74207281, 77232917
];

#  input
print("2^p[index] - 1 = ? ");
print("index(1 - $max_index) = ");
index = parse(Int, chomp(readline(STDIN)));
if (index > max_index) || (index <= 0)
	print("ERROR: index = $index is illegal!\n");
	return EXIT_FAILURE;
end

println("p[$index] = ", known_mersenne_prime_p[index]);

#  a := 2^p - 1
a = BigInt(2)^known_mersenne_prime_p[index] - 1;
println("2^", known_mersenne_prime_p[index], " - 1 = $a");

# Mersenne Prime Number Calculator by GNU MP

# Run: julia mersennne.jl

# Mersenne Prime Number

# 2^p - 1

# cf.) http:# primes.utm.edu/mersenne/

max_index = 50; # 2018-01-05

known_mersenne_prime_p = [

2, 3, 5, 7, 13, 17, 19, 31, 61, 89,

107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423,

9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049,

216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917, 20996011,

24036583, 25964951, 30402457, 32582657, 37156667, 42643801, 43112609, 57885161, 74207281, 77232917

];

# input

print("2^p[index] - 1 = ? ");

print("index(1 - $max_index) = ");

index = parse(Int, chomp(readline(STDIN)));

if (index > max_index) || (index <= 0)

print("ERROR: index = $index is illegal!\n");

return EXIT_FAILURE;

end

println("p[$index] = ", known_mersenne_prime_p[index]);

# a := 2^p - 1

a = BigInt(2)^known_mersenne_prime_p[index] - 1;

println("2^", known_mersenne_prime_p[index], " - 1 = $a");

2018-01-01

2018年1月1日（月・祝）？

　明けましておめでとうございます。本年もよろしくお願い致します。

　昨年は殆ど更新しなかったが，今年はせめて週報とか月報程度のことは書いていきたいものである。遅ればせながら英語やり始めたり，「LAPACK/BLAS入門」の出張講義したり(PDF)，出るあてがあるんだかないんだか分らん本のリライトを始めたりして，年なりにできることをチマチマやってきたので，本年はそれを継続していきたい。

　ま，守れないのはいつものことだが，年初ぐらいは抱負を語るのも悪くないかな。ともかく，今後ともよろしくお願い致しまする。

2017-08-12

8/12(土) 駿府・？

　今週は月曜日午後から職場がロックアウトされているので，自宅でシコシコダラダラ過ごす。積ん読していた漫画類をほぼクリアした以外は，神さんとWindowsマシンとワシの自宅MacBook AirをWindows 10 1703にアップデートする作業に費やしていた。

　神さんのマシンはWin10化していたのだが，何故だか7月のupdateしてたらブート途中でブルースクリーンになる症状が出てオジャンになっていたのを，Windows 7に戻してから(Office 2010 Personal入れなきゃならんので）改めてWindows 10にするという，途中のUpdate時間だけがやたらにかかる作業。ハードウェアのトラブルではないらしく，HDDは正常だったので原因は不明だが，一応元に戻せたのでよしとする。一つ感心したのは，Windows 10のライセンスは，以前，無料アップデート期間中にWin7 -> Win10にした時にアクティベーションしたものがそのまま使えたこと。当たり前のことなんだけど，改めてライセンス取り直しで金払うとなればアホらしいからねぇ。

　ワシのMacBook Airは，Boot Camp環境のWin 8.1を10にしたこと。そのままアップデートできればよかったのだが，あいにくDreamSparkのWin8.1 Enterprise版を入れてしまっていたので，Win10のEnterprise版にするには改めてクリーンインストールせねばならず，それならばPro版にしておこうということで改めてインストールし直ししたのであった。

　結局２台ともクリーンインストールではあったが，古い環境を一掃できたのはよかったかなと。Cloud時代になって，必要なデータファイルはiCloud，Google Drive(Backup & Sync)，OneDriveに入れてあるのも幸いした。おかげでmacOS SerraでもWin10でもアプリさえ入れなおせばすぐに作業が続けられるのも時代であるなぁとシミジミ。

　そういえば，今日はVPSへ移動させたワシのWebアプリをPerl CGiからPHPに変換したのであった。今更Perl CGI環境を維持したところで良いことないし，セキュリティ的にもイヤ。時代的にもPHPは必要最小限，サーバ側のデータ入出力だけにとどめ，殆どのUIはJavaScriptで実装という流れだし，Webアプリもどんどん分業化が進む昨今である。

　ボチボチ寝ます。

2017-07-132017-07-13

破綻する区間演算の事例

　精度保証計算（という言葉はあまり好きではないのだけれど）の事例も紹介しようかと，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 写像

// 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;
}

// 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;

}

2018年4月
日	月	火	水	木	金	土
« 3月
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30