憂鬱な午後のひととき
最近, 靴の裏にローラーが入ったものを履いてる子供をちらほら見かける.
それも広い遊び場でじゃなくて, スーパーや駅とかの場所で. 転んだりぶつかったりしてるの見るとはらはらする.
ちょっと前は, 柔らかいサンダルを履いててエスカレーターで挟まる事故があったけど, こっちの方が遥かに危険だろう.
なんであんなもの履かせるのか理解できない. 重大な事故にあわないと分からんのか.
『臓器農場』(帚木蓬生)
立て続けに同じ作者の. 無脳症の胎児を意図的に発生させて, それを移植のために使う(脳が無いからそもそも人ではないという論理). かなり不気味.
ただ, 無脳症の胎児の描写が怖い印象が強くて, 「その胎児も生きてる」というヒロインの主張がやや弱い感じ.
論理なんかではなく, 問答無用にそう思わせるような場面が欲しかった.
あと, 邪魔になった人を悪役がさくっと殺してしまうパターンはもうちょっとなんとかならんか.
四則演算ルーチンの再開. まあ世の中ええのはいっぱいあるので自分で作るのはどうかと思うが, やっぱり細かい制御は自分でやりたいしな.
fmathをblogで
紹介したよとのメール. 少しでも反応があるのはやっぱり嬉しいもんだ.
東京に来て初めてまともな本棚を注文した. 実家にある本も少しずつ持ってきたいな.
本社で社内技術発表会.
amachang,
TAKESAKOさん,
奥さんと営業の方たちとの座談会がウケていた.
立場が違うと視点も違うなと, なかなか興味深かった.
『千日紅の恋人』(帚木蓬生)
古いアパートの管理人である時子は, 読み始めた当初, 60歳ぐらいのイメージだったが途中で38歳と判明して戸惑った.
あれれ, 印象と大分違う. くみ取り式のトイレとか, 家賃が3万とかの時代に葬式代200万って高すぎちゃうのかな.
まあ, それはともかく管理人とアパートの癖のある住人たちとのトラブルや助け合いが淡々と丁寧に(ややくどい嫌いも)述べられる.
重大な事件も起こりそうで起こらない. 私はこういうの嫌いやないけど, じれったい人も多そうだ. 起承承承承...承承結(ハッピーエンド)みたいなの.
_s の戻り値としてそれが返ってきているのはなんかよろしくない気が
???
MSDNに
errno_t strcpy_s(
char *strDestination,
size_t numberOfElements,
const char *strSource
);
>正常終了した場合は 0 を返し、それ以外の場合はエラーを返します。
>実行の継続が許可された場合、これらの関数は EINVAL を返し、errno を EINVAL に設定します。
と明記されている以上, 問題は全くないように思いますが, 何にひっかかってらっしゃるのでしょう.
str*_s系
>なんだこの 22 は
Invalid argumentを表すerrnoですね. *_s系は失敗するとerrnoが設定されるのでそれを見れば原因が分かる.
printf("err%s\n", strerror(22));
>Invalid argument
_snprintf_sでも同じです(cf.
例外ハンドラ).
『とめはねっ！鈴里高校書道部』(ドラマ)
原作は未読. でも配役がぴったりという印象だったし, 演出もくどすぎず, おもろかった.
グループでパフォーマンスしながら書道をするというのが新鮮. ラストの字もうまいことそれぞれの個性が出ていたように思う.
久しぶりに東京現代美術館に行って,
サイバーアーツジャパン―アルスエレクトロニカの30年というのを見た.
残念ながら自分の感性がずれていたのか全然来るものが無かった. 井上雄彦の大壁画はすごかったけど. よーあんな大きいの書いたな.
休日は混んでるだろうからとせっかく有休とって行ったのに, 乗ってみたかった
media vehicleなどいくつか体験できるものは平日はやってなかったのも悲しい.
『パズル・パレス』(ダン・ブラウン, 越前敏弥訳, 熊谷千寿訳)
実につまらんかった. 正直ダ・ヴィンチ・コードもそんなに面白いと思わなかったのだけど, こちらはストーリー以前に突っ込み所が多すぎてすらすらと読めなかった.
unordered_mapってkeyとvalueのペアが一致していたとしても入力の順序が異なるとoperator==()がfalseになることがあるんだな.
当たり前と言っちゃ当たり前だけどちょっと悩んだ.
下記コードはstd::mapだとsameだけどunordered_mapではdiffになる. operator==()はない方がいい気もする... ってgccだとoperator==()は無いみたい.
規格だとどうだっけ.
# Unsupported expressionsとなってた.
#include <stdio.h>
#include <map>
#include <string>
#include <unordered_map>
template<class T>
void put(const T& x)
{
for (typename T::const_iterator i = x.begin(), ie = x.end(); i != ie; ++i) {
printf("%s:%d\n", i->first.c_str(), i->second);
}
}
template<class T>
void test()
{
const char tbl[][10] = { "std", "map" };
const size_t n = sizeof(tbl) / sizeof(*tbl);
T h1, h2;
for (size_t i = 0; i < n; i++) {
h1[tbl[i]] = i;
h2[tbl[n - 1 - i]] = n - 1 - i;
}
if (h1 == h2) {
puts("same");
} else {
puts("diff");
puts("h1");
put(h1);
puts("h2");
put(h2);
}
}
int main(int argc, char *argv[])
{
puts("unordered_map");
test<std::tr1::unordered_map<std::string, int>>();
puts("map");
test<std::map<std::string, int> >();
}
『エンブリオ』(帚木蓬生)
こないだ読んだ『インターセックス』の前作. のっけから衝撃的な手術が行われている.
主人公の医師は倫理委員会などを無視して最先端の生殖医療を進めるが, 必ずしも金銭や名誉が第一の目的というわけではない.
患者が望むものが第一, そして医療だけでなく経営や看護士たちの教育にも力を入れてるあたりがよくあるパターンと違う. 多忙なのに患者の愚痴や不満もじっくり聞くというのがすごいな.
これだけの頭脳と技術があればわざわざリスクの多い犯罪なんかしなくても別の方法がとれたんじゃないかと思うのがちょっと残念.
という内部のアルゴリズムは分からん, という方でも
fmathはfmath.hppをincludeしてexpやlogをfmath::expやfmath::logにするだけで使えるから,
そういう計算を多用してるなら(特にgccで効果あり)は試してもらえるとうれしい(多少の誤差は注意してね).
syoyoさんも速くなったよって言って下さってるし.
二つ目.
なぜroundが欲しかったかというと x を離散的な値 n と残り t = x - n に分けたかったから.
そうすると
exp(x) = exp(t) * exp(n)
と分割できて後半をテーブル引きできる.
前半の t はテイラー展開による近似計算をしたいので, できるだけ abs(t) が小さいようにしたい.
n として round() を使うとせいぜい abs(t) <= 0.5 としか言えない.
もし n が t に最も近い半整数(m/2, m ∈ Z)ならば abs(t) <= 0.25 とできる.
もっと一般に n = m/2
s, m ∈ Z が t に最も近い数値ならば abs(t) <= 1/2
s + 1 となる.
そういう形をとるには x を a 倍してから round して a で割ればよい.
t2 = x - round(x * a) / a とおくと(a = 2
s).
abs((x * a) - round(x * a)) <= 0.5
だから t2 <= 0.5 / a となる.
floatの仮数部は 23bit だから s = 11 ならばテイラー展開 exp(t) = 1 + t + t
2 / 2 + ... の
2次の項が1/2
12 * 2で抑えられる(***). つまり無視できる.
exp(t) = 1 + t
三つ目.
上記のとき, テーブル部は n = 1/2
s = 1/2
11 で離散化することになり, テーブルサイズは(xmax - xmin) * 2
11 ≒ (8 - (-88)) * 2048 = 360448
と大きくなってしまう. ここでテーブルサイズを小さくすることを考える.
離散化の部分を見直してみる. a は t2 を 1/2
12 で抑えたかっただけだから a 自体が2の巾である必要はない.
テーブル引きでは exp(round(x * a) / a) = exp(n / a) の計算が欲しいだけである.
a' = a / log(2)とおいてみよう.
exp(n / a) = exp(n / (a' * log(2))) = pow(2, n / a')
とeの巾ではなく2の巾にすれば pow(2, round(n / a')) の部分はシフトだけで実現できて, 残りの端数のみをテーブル引きにすることができる(cf.
powの近似計算 または iccのexp()).
a / log(2)が(***)を満たすためには a = 2
11であればよい.
今回はちょっとけちってa = 2
10にした.
そのときa' = 1024 / log(2) = 1477.3197. これが
aの正体. cはそれの逆数.
クイズの答え一つ目.
b = 12582912 = (1 << 23) + (1 << 22)のとき, x + b - bはxを整数に丸める機能を持つ(ただしxの範囲は[0, (1 << 22) - 1]).
整数に丸めるだけなら1 << 23でよい(その場合xの範囲も1bit増やせる).
今回は丸めと同時にx + bの仮数部の値も同時に利用したい. そのときの今回のbを使うと負のときでも整合性を保てる.
union fi {
unsigned int i;
float f;
} fi;
// 正の場合はどちらでもOK
fi.f = 1234.567017 + (1 << 23);
printf("%d\n", (int)((fi.i << 9)) >> 9); // 1234
fi.f = 1234.567017 + (1 << 23) + (1 << 22);
printf("%d\n", (int)((fi.i << 10)) >> 10); // 1234
// 負の場合は(1 << 23)ではうまくいかない
fi.f = -1234.567017 + (1 << 23);
printf("%d\n", (int)((fi.i << 9)) >> 9); // -2469
fi.f = -1234.567017 + (1 << 23) + (1 << 22);
printf("%d\n", (int)((fi.i << 10)) >> 10); // -1235 // こちらはOK
ただこのコードはfloatとdoubleを区別しなかったり, x + b - bをxにしてしまうようなアグレッシグな最適化では動作しないだろう(asmで書くなら関係ないけど).
#include <stdio.h>
float round(float x)
{
const float b = (1 << 23) + (1 << 22);
x += 12582912;
return x - 12582912;
}
int main()
{
const double tbl[] = {
0, 0.4, 0.5, 0.6, 1.2, 1.5, 1.9, 2.5, 3.5, 123.456,
(1 << 22) - 1,
(1 << 22) - 0.5,
(1 << 22),
};
for (size_t i = 0; i < sizeof(tbl) / sizeof(*tbl); i++) {
float x = tbl[i];
printf("round( %f)= %f\n", x, round(x));
printf("round(%f)=%f\n", -x, round(-x));
}
}
// 結果
round( 0.000000)= 0.000000
round(-0.000000)=0.000000
round( 0.400000)= 0.000000
round(-0.400000)=0.000000
round( 0.500000)= 0.000000
round(-0.500000)=0.000000
round( 0.600000)= 1.000000
round(-0.600000)=-1.000000
round( 1.200000)= 1.000000
round(-1.200000)=-1.000000
round( 1.500000)= 2.000000
round(-1.500000)=-2.000000
round( 1.900000)= 2.000000
round(-1.900000)=-2.000000
round( 2.500000)= 2.000000
round(-2.500000)=-2.000000
round( 3.500000)= 4.000000
round(-3.500000)=-4.000000
round( 123.456001)= 123.000000
round(-123.456001)=-123.000000
round( 4194303.000000)= 4194303.000000
round(-4194303.000000)=-4194303.000000
round( 4194303.500000)= 4194304.000000
round(-4194303.500000)=-4194304.000000
round( 4194304.000000)= 4194304.000000
round(-4194304.000000)=-4194304.000000
・とりあえずC++でのみ書いて可変長で四則演算はできるようになった.
・クイズは, まあ考えてみるのがいいんじゃないかな. 情報落ちを利用した整数化(浮動小数点の位置を固定化)するテクニックは, 『ハッカーの楽しみ』にも載ってたし, 昔の
LAMEでも使われてたように思う.
JANTのため筑波へ.
遠かったので電車の中でひたすら簡易多倍長ライブラリ作る. 昔未踏で作ったやつをベースにしてもよかったのだけど,
今回は64bitにも対応したいし速度よりも汎用性を重視したいので一から作る.
・η
TペアリングのCインタフェースをこんな感じで作ってるけどどうだろう.
アプリを作るには逆数用の多倍長ライブラリもいるなあ. 簡易版でよいから作るか.
// 初期化(今のところF_3^193上のηTペアリングのみ. そのうち他も?)
MIE_System *sys;
MIE_Err err = MIE_CreateSystem(&sys, MIE_ETAT_M193);
// 基礎体の元(F_3^m)
MIE_Fp *x = sys->FpCreate();
// 拡大体の元((F_3^m)^6)
MIE_Fq *s = sys->FqCreate();
MIE_Fq *t = sys->FqCreate();
// 楕円曲線上(y^2 = x^3 - x + 1)の点
MIE_Ec *P = sys->EcCreate();
MIE_Ec *Q = sys->EcCreate();
MIE_Ec *R = sys->EcCreate();
// x = 0としてPをセット．yは自動的に決まる
sys->EcSetX(P, x);
// x = 21000222020202022202202221111としてQをセット．
sys->FpSetByStr(x, "21000222020202022202202221111");
sys->EcSetX(Q, x);
// Rもセット．
sys->FpSetByStr(x, "1111222000002220202020222022022211111");
sys->EcSetX(R, x);
sys->pairing(s, P, R); // s = e(P, R)
sys->pairing(t, Q, R); // t = e(Q, R)
// ペアリングとってから掛けたものをsとする
sys->FqMul(s, s, t); // s = e(P, R) * e(Q, R)
// 点どうして加算して
sys->EcAdd(P, P, Q); // P = P + Q
// ペアリングとったものをtとする
sys->pairing(t, P, R); // t = e(P + Q, R)
// それらが等しければOK
TEST_ASSERT(sys->FqEqual(s, t)); // e(P, R) * e(Q, R) = e(P + Q, R)
// 後始末
sys->FpDestroy(x);
sys->FqDestroy(s);
sys->FqDestroy(t);
sys->EcDestroy(P);
sys->EcDestroy(Q);
sys->EcDestroy(R);
MIE_DestroySystem(sys);
・
解答編
トライしてもらった. 面白かったと言われたのでよかった. 「exp(x)の導関数はx」ではなくてexp(x)ですね.
exp(x)の一次近似の項がx.
ヒントを少し.
一つ目: a = 1.4773197e3, b = 12582912, c = 6.7690155e-4とする.
t1 = x * a + b;
t2 = x - (t1 - b) * c = x - x * a * c = (1 - a * c) * x
と書いてみると, bを足して同じbを引いてるのはなんだよ, a * c = 0.9999999947(ほぼ1)なのでt2は"ほぼ0な値 * x"?と思った方, 数学的には正しいけど今回のプログラムでは間違い.
実際, 下記の例のようにx = 1のときでもt2 = 2.164107e-004とそれなりにある. それは何故か. またその意図は.
# これを自力で解けたらかなりなもの.
二つ目: y = (t1 - b) * cとおくと, (2)からx = t2 + y. つまりexp(x) = exp(t2) * exp(y).
t2が十分小さいから(5)でexp(t2) = 1 + t2としてる. とするとその精度はどう見積もってるのか.
三つ目: 残りのexp(y)はテーブル引き + αで求めてる. そのままなら結構大きなテーブルサイズになるはずなのに1024個しか用意していないのは何故.
指数関数に関わる浮動小数点数のクイズ
次の関数は
fmat.hppのexp()に使われているロジックをCで書いたものである.
そこに現れる各定数は何を意味するのか, どういうロジックなのか答えよ(ただし -88 <= x <= 88としてよい -- そうじゃなければ0か∞なので).
# 本当は
並カンの"SIMD パズルの森"で話してみたかったが, あまりにニッチ過ぎる感があってお蔵入りしたままになってしまった.
unsigned int tbl[1024]; // 実行時に一度だけ初期化される
union fi {
float f;
unsigned int i;
};
float expC(float x)
{
float t1 = x * 1.4773197e3 + 12582912; // (1)
fi fi;
fi.f = t1;
float t2 = x - (t1 - 12582912) * 6.7690155e-4; // (2)
unsigned int u = ((fi.i + 130048) >> 10) << 23; // (3)
fi.i = u | tbl[fi.i & 1023]; // (4)
return (1 + t2) * fi.f; // (5)
}
たとえばx = 1のときは1行目でt1 = 1.258439e+7, (2)でt2 = 2.164107e-004, u = 0x40000000, tbl = 0x2deeb1となって2.718282を出力する.
同様にx = -4.3ならt1 = 1.257656e+7, t2 = -3.215451e-4, u=0x3c000000, tbl = 0x5e60f5で0.013569.
x = 5.8ならt1 = 1.259148e+7, t2 = 3.077103e-4, u = 0x43800000, tbl=0x251958で330.299591となる(3/6 説明のため(1)のtと(2)のtの名前を変えました).
『インターセックス』(帚木蓬生)
染色体が男性で外見が女性, あるいはその逆, あるいは両方の場合をそういうらしい. 性同一性障害とはまた違う.
何か陰謀を企む医院長 vs そこに引き抜かれた性差医療の専門医師という話かと思ったら, わりとあっさり医院長が専門医師に説得されて肩すかし. 更に裏があるのかと邪推したけどそちらは主眼ではなかったようだ(前作(未読)の続き?らしい).
それはともかく当事者たちの悩みや, 医者自身が持つ偏見, コミュニティの話などルポタージュ的な部分が圧巻.
ユートピア的な病院や上流階級の描写が多少多い気もするけどお薦め.
VAIOのリモコンが動かなくなった. テレビ見るのがめんどくさい.
眼鏡のフレームが壊れた. 直しにいかなきゃ.
去年は
憩室炎とか
ようわからん検査とかなんやかんやで医療費が10万越えた. 確定申告しに行こう.
・
Coverityの静的解析ツールの体験版を見せてもらった.
社内で実際に開発されたものに対してデモしてもらう. かなり精度よくピンポイントにレポートがあがる感じだった.
ある関数の使われ方の統計も取っていて, たとえば「返り値のチェックをある箇所でだけやってない」ようなのも報告していた. これは単なるパターンマッチだけではできない.
自分が書いたところもcriticalではないものの若干見つかった. VCやboostのバグも出てたなあ. 大規模なプロジェクトだと導入を検討する価値は十分あると思った.
・
ペアリングフォーラムで
ηTペアリングの実装なんかもC++じゃなくてCのインタフェースがあるとうれしいという話があった.
そのうち用意する.
御意見は光成滋生<herumi@nifty.com>までお願いします