今年は実家に帰れなかった. 正月実家に戻らないのは初めてかも.
で, 年越しそばのための鰊(にしん)を買いにスーパーに行く. しかし見当たらない.
「にしん蕎麦のにしんはどこですか?」と定員に聞くとカップラーメンのコーナーに案内される.
なぜにカップラーメンとしばらく悩んで気が付いた. いや日清そばとちゃうから.
ああ, 京都, 奈良に戻りたい.
Excel列名変換問題で第２回社内プログラミングコンテストを開催してみた
『「Perlが苦手なんで解けませんでした」はNGとする』と書いてあったのでやってみる. ちなみに私はperlをtakesakoさんがやってるマインスイーパを解ける言語という程度しか知らない.
まず引数のとり方を調べる. $#ARGVとか$ARGV[0]を使うのね. 文字列の比較はeqと. shellっぽい. 配列は@とあったので変数$xを@xにすれば配列になるのかと思ったけどそうじゃなかった. 文字列のn番目の文字を取る方法が分からなかったのでsubstrを使うことにする. 文字コードの取得はordなのね. ここまで15分ぐらい.
一問目が動いた. ところがABCDEFGHIJKLMNを食わせると結果が小数になってしまう. 多倍長演算を使わないといかんのか. 調べてBigIntというのがあるらしい. 動いた. ここまで20分ぐらい.
次2問目. 関数の作り方を調べる時間もないし, そこまでのものでないのでフラットに書くことにする. else ifはelse ifでもelifでもなくelsifなのか. ついつい行末のセミコロンを忘れてエラーになってしまう. 文字の連結は+じゃなくてピリオドなのか. +だと0になるからはまったぞ.
BigIntにおかしな値を食わせたときどうなるのか. "NaN"ってなってるのでよく分かってないけどそれと比較してエラーにすることにする. 1問目の逆関数にすればよいと思ってとりあえず書いた. ところが26を食わすとおかしな値になる.
しばらく悩んでああそうか1ずらさないといけないんだと修正. 一通り動いてるっぽい. ここまで40分ぐらい.
use Math::BigInt;
if ($#ARGV != 1) {
print "ERR\n";
exit(1);
}
$mode = $ARGV[0];
$x = $ARGV[1];
if ($mode eq '0') {
$ret = Math::BigInt->new("0");
for ($i = 0; $i < length($x); $i++) {
$c = substr($x, $i, 1);
if ('A' le $c && $c le 'Z') {
$ret = $ret * 26 + ((ord $c) - (ord 'A') + 1);
} else {
print "bad input $x\n";
exit(1);
}
}
$str = substr($ret, 1);
print "$str\n";
} elsif ($mode eq '1') {
$v = Math::BigInt->new($x);
if ($v == 0 || $v eq "NaN") {
print "bad input $x";
exit(1);
}
$ret = "";
while ($v > 0) {
$q = ($v - 1) / 26;
$r = ($v - 1) % 26;
$ret = chr($r + (ord 'A')) . $ret;
$v = $q;
}
print "$ret\n";
} else {
print "input 0 or 1";
exit(1);
}
exit(0);
で, 答えを見る. 変数を配列にするのはsplit(//)を使うようだ. これはすぐには分からんな. でも, 引数の「上限なし」って言ってるわりに1位の人のプログラムでもABCDEFGHIJKL与えると小数になってるけどいいのかな. あと, 「TDDでやってる」のに例外テストが入ってない….
ついでにHaskellでも解いてみた. printfは文字を"ABX"ではなくABXで表示させたかったので使った. 簡単にできる方法あるのかな.
追記. melponnさんに教えていただいた. putStrLnを使えばよかっただけでした. ありがとうございます.
import System
import Data.Char
import Text.Printf
aton::String->Integer
aton str = foldl f 0 str
where f x c | 'A' <= c && c <= 'Z' = x * 26 + fromIntegral (ord c - ord 'A' + 1)
| otherwise = error $ "only [A-Z]:" ++ [c]
ntoa::Integer->String
ntoa 0 = error "bad input 0"
ntoa num = fst . until ((/=0).snd) f $ ([], num)
where f (os, num) = f (chr (fromIntegral r + ord 'A'):os, q)
where (q, r) = quotRem (num - 1) 26
main = do
args <- getArgs
case args of
["0", input] -> print . aton $ input
["1", input] -> printf "%s\n" . ntoa $ (read input::Integer)
_ -> error "usage: [0|1] input"
ああ, 負の値もはじくべきだった.
返り値の型でディスパッチ
これはC++03 Advent Calendar 2011の25日目の記事…ではありません(そんなの無いし). これは昨日の記事の番外編です.
昨日の記事を見て文字列と数値の変換関数で, std::to_string()はよいとしても, atoiやatol, atollがかっこ悪いと思われた方, 手を挙げて.
C++では引数の型が異なる関数を多重定義はできますが, 返り値の型が異なるだけではできません. したがって
int a = atoi("123");
long b = atol("123");
unsigned long c = atoul("123");
...
のように返り値に応じて関数名を変える必要があるのです. でもこれを
int a = my_atoi("123");
long b = my_atoi("123");
unsigned long c = my_atoi("123");
...
のようにできたらうれしいですね. 実はできたりします. しかもC++11の新機能やtemplateの黒魔術は不要です.
たとえばshort, unsigned short, intに対応したmy_atoiを作ってみましょう. my_atoiを関数ではなくクラスにするのがミソです. それとoperator XX()を組み合わせます.
#include <limits.h>
#include <stdlib.h>
#include <stdexcept>
struct my_atoi {
my_atoi(const char *str) : str_(str) {}
operator short() const
{
int x = atoi(str_);
if (x < SHRT_MIN || x > SHRT_MAX) throw std::invalid_argument("not in short");
return static_cast<short>(x);
}
operator unsigned short() const
{
int x = atoi(str_);
if (x < 0 || x > USHRT_MAX) throw std::invalid_argument("not in unsigned short");
return static_cast<unsigned short>(x);
}
operator int() const
{
return atoi(str_);
}
private:
const char *str_;
};
このようにすると次のように使うことができます.
#include <stdio.h>
int main() try
{
{
short x = my_atoi("-123");
printf("1. x=%d\n", x);
}
{
unsigned short x = my_atoi("123");
printf("2. x=%d\n", x);
}
{
int x = my_atoi("12345678");
printf("3. x=%d\n", x);
}
{
short x = my_atoi("12345678");
printf("4. x=%d\n", x); /* shortに入らないので例外が飛ぶ */
}
} catch (std::exception& e) {
printf("%s\n", e.what());
}
>1. x=-123
>2. x=123
>3. x=12345678
>not in short
意外と簡単ですね. ただ二つほど注意点があります.
一つはコンストラクタと評価のタイミングが分離されていることです. コンストラクタでstr_を初期化してからoperator()で評価するまでにstr_の先が消えるとうまく動きません. まあこの場合はコンストラクタでlong long val_;にいれてoperator()では範囲チェックだけをするというやりかたで回避できると思います.
もう一つは返り値の推論ができないとうまく動きません. たとえば
std::string s = std::to_string(my_atoi("123"));
はmy_atoiを何に変換すればよいのか分からないのでコンパイルエラーになります. 当然といえば当然ですが. こういうときは何に変換したいのかを明記する必要があります.
std::string s = std::to_string(static_cast<int>(my_atoi("123")));
おまけ: 今気が付いたけどC++11の規格書のp.1225のTable 150のSHRT_MAX, SHRT_MINなどの値の一覧, USHRT_MAXじゃなくてUSRT_MAXになってる. これってtypoやよね. 既に指摘されてましたっけ.
整数型long long
これはC++11 Advent Calendar 2011の24日目の記事です. もうあらかた皆さんが解説されて食傷ぎみですので簡単なところを紹介します.
C++11にはついに64bit(以上)の整数を表せるlong longが入りました.
# 規格書自体に64bitという文字は見当たらないのですが, climitsで定義されるLLONG_MAXを表せないといけなくて, C99ではその値は9223372036854775807(=263-1)以上とあるので大丈夫だと思います. 多分. 違ってたらすんません.
というわけで大きな整数値の型が入ったことでいくつか変更があります. まず4294967295(=232 - 1)の型はなんでしょうか.
今までのC++(C++03)でint, longが32bitの環境だった場合, 4294967295はintやlongでは表せないため未定義(!)でした. unsigned intやunsigned longとはならなかったのです. それがC++11ではlong longとして扱われるようになりました.
なお0xffffffffと書くとC++03でもC++11でもunsigned intとなります(cf. C/C++における定数クイズ).
long longを扱えるので文字列を整数に変換する関数atoiやatolの仲間atollも入りました. ついでにstd::stringからの変換関数stoi, stol, stollというのも追加されています.
std::string s = "9876543210932"; long long b = stoll(s);
stoi系は変換できないと例外が飛びます.
更に逆方法, つまり数値を文字列にする関数も追加されています.
long long a = 123456789012345LL; std::string t = std::to_string(a);
便利ですね.
大きなファイルを扱う.
gccとVCの巨大なファイルで書いたようにC++03では4GBを超えるファイルをうまく扱えませんでした. しかもfpos_t, off_type, streamoff, OFF_Tなど様々な型が入り乱れてわかりにくくなっていました.
C++11ではそれらの型の関係をできるだけすっきりさせました. std::streamoffはファイルサイズを表せる十分大きい整数型(普通long long)と明記されています. std::streamsizeは殆どsize_tなのですが-1を扱いたいのでPosixで定義されているssize_tを想定します. seekgなども正しく操作できます.
int main()
{
std::ifstream ifs("test.data", std::ios::binary);
ifs.seekg(4299999997LL, std::ios::beg);
std::cout << ifs.tellg() << std::endl;
ifs.seekg(0, std::ios::end);
std::cout << ifs.tellg() << std::endl;
}
test.dataが巨大なファイルの場合, 32bit環境においてもVC2010では上記リンク先にあるように独自拡張のifs.tellg().seekpos()なら正しく値を取れますが, ifs.tellg()ではおかしな値になっていました. VC2011ではifs.tellg()のままで正しい値を取れるようになっています. なおgcc-4.6.0ではまだうまく実装されていないようです.
以上long longにまつわるC++11の改良部分でした.
boost::mathで数値積分の式を作ってみる
これはBoost Advent Calendar 2011の23日目の記事です.
boost::mathには数学や物理などで使われる特殊関数と呼ばれるものがいくつか用意されています.
ここではそれらを使って数値積分の公式の一つであるGauss Legendreの公式を作ってみましょう.
以下, しばらくboost::mathのboostではなくmathよりの解説が続きます. すいません.
関数f(x)が与えられたときに区間[a, b]の積分を求めましょう.
数値積分は区間[a, b]の間の適当な点ziにおけるfの値f(zi)に適当な重みwiをかけて足しあわせて求めるのが基本です.
∫ab f(x)dx ≒ Σi f(zi) wi
たとえば区間[a, b]をN等分してh = (b - a) / N, zi = a + h i, wi = hとすると本来の積分の定義に近いものとなります.
台形公式やSimpsonの公式では積分区間を等分して求めるのに対し, Gauss Legendre公式はそうではありません. ziやwiをうまくとってよりよい近似を与えようという戦略です.
さて, その決め方ですが[-1, 1]区間上定義されたn次のLegendre多項式Pn(x)というものがあります.
Pn(x) = 1/(2nn!) (d/dx)n(x2 - 1)n.
それの零点(解)をziとして選びます. wiは
wi = 2 / (n Pn'(zi) Pn-1(zi))
で求められます. なぜこうすればうまくいくのかについては触れません. 「数値微分と数値積分[pdf]」(ソフトウェアとしての数値計算)などの資料を参考にしてください.
ではPn(x)の零点を求めてみましょう. 数式処理ソフトがあれば楽勝ですが一から求めるのはなかなか面倒です. C++11にいくつか数学関数が追加されましたがLegendre多項式は入っていません. そこでboost::mathを使ってみましょう. ここからはboostよりの話です. もちろんLegendre多項式も用意されています.
Pn(x)はboost::math::legendre_p(n, x)というので計算できます. また与えられた関数の零点を求める関数boost::math::tools::bisectが用意されています. 名前の通り与えられた区間を2分法で求めるものだと思われます. 関数f(x)とf(a)f(b) < 0となるような区間[a, b]を与えるとその間の解を求めてくれるラッパー関数を用意してみました.
#include <boost/math/tools/roots.hpp>
template<class F>
bool findRoot(double *x, F f, double a, double b)
try {
boost::math::tools::eps_tolerance<double>
tol(std::numeric_limits<double>::digits);
boost::uintmax_t max = 100;
std::pair<double, double> r = boost::math::tools::bisect(f, a, b, tol, max);
*x = fabs(f(r.first)) <= fabs(f(r.second)) ? r.first : r.second;
return true;
} catch (std::exception&) {
return false;
}
eps_tolerance<double> tolは区間を段々と縮小してa == bとみなす許容範囲(tolerance)を与えます. maxはiterationの最大回数です.
返り値は十分小さい区間r = [a, b]なのでそのうちよりfにおける値が小さい方を解とすることにしました. 求められないと例外が飛ぶのでfalseを返すことにします.
今回はPn(x)の全ての零点が必要です. 偶数次のLegendre多項式は偶関数なので[0, 1]区間をn等分して適当に解を求めさせることにしました.
ziが求まればあとはwiです. これにはPn(x)の微分が必要です. 差分で求めてもよいですが, Legendreの同伴(associate)関数との間に
legendre_p(n, m, x) = (-1)m (1 - x2)m/2 (d/dx)m Pn(x)
という関係があるのでそれを使います. legendre_pもboost::mathにあります.
/*
w_i = 2/(nP_n'(z_i)P_{n-1}(z_i))
legendre_p(n, m, x) = (-1)^m (1 - x^2)^(m/2) (d/dx)^m P_n(x)
*/
static double calc_w(double x)
{
const int m = 1;
double diff_legendre = boost::math::legendre_p(n, m, x)
/ (pow(-1.0, double(m)) * pow(1 - x * x, m / 2.0));
return 2 / (n * diff_legendre * boost::math::legendre_p(n - 1, x));
}
ziとwiが求まれば与えられた区間[a, b]に対して[-1, 1]に変数変換して重みつき和を求めることで積分を行います. もちろん, 求まったziとwiをハードコーディングしておけばboostは不要になります.
template<class F>
double integrate(F f, double a, double b) const
{
if (a == b) return 0;
bool minus = false;
if (a > b) {
minus = true;
std::swap(a, b);
}
double h = (b - a) / 2;
double ave = (b + a) / 2;
double ret = 0;
for (int i = 0; i < N; i++) {
double tmp = h * z[i];
ret += w[i] * (f(ave + tmp) + f(ave - tmp));
}
ret *= h;
if (minus) return -ret;
return ret;
}
では例としてexp(-x2)を[0, 0.9]で積分してみましょう. 12次のLegendre多項式を使ってみました. 厳密解は誤差関数erf()と呼ばれるものを2πで割ったものになります. erf()もboost::mathにあります.
double g(double x) { return std::exp(-x * x); }
int main()
{
const double t = 2.3;
const double pi = boost::math::constants::pi<double>();
GaussLegendre<> gl;
printf("integral(0, t) = %.18f\n", gl.integrate(g, 0, t));
printf("by erf = %.18f\n", boost::math::erf(t) / 2 * sqrt(pi));
}
>integral(0, t) = 0.885213811571632810
>by erf = 0.885213811571633040
なかなかよい感じです. Simpson法と比べてみましょう.
template<class F>
double simpson(F f, double a ,double b ,int N)
{
double s1 = 0;
const double h = (b - a) / double(N);
for (int i = 1; i < N/2; i++) {
double x = a + h * 2 * i;
s1 += f(x);
}
double s2 = 0;
for (int i = 1; i<= N/2; i++) {
double x = a + h * (2 * i - 1);
s2 += f(x);
}
return (b-a) / (3.0 * N) * (f(a) + f(b) + 2 * s1 + 4 * s2);
}
printf("sin(0, pi) = %.18f\n", gl.integrate((double (*)(double))sin, 0, pi));
printf("sin(0, pi) = %.18f\n", simpson((double (*)(double))sin, 0, pi, 100));
>sin(0, pi) = 2.000000000000000400 // Gauss Legendre
>sin(0, pi) = 2.000000010824504000 // Simpson
128点Simpsonより12点しかとらないGauss Legendreがずっと精度よいです. 初めてこの方式をみたときは感動したものです.
なお, 零点を求める方法として2分法の他にTOMS Algorithm 748という方法も実装されていました. 使い方は上記のboost::math::tools::bisectをboost::math::tools::toms748_solveに変えるだけです. 実験してみると2分法では値の確定までに50回程度呼ばれる(doubleの精度が52bitだから)のに対して, toms748_solveは10回程度で済みました. 関数の値の計算自体が重たいときはなかなか効率のよい手法のようです.
今回のソース全体はgauss_legendre.cppです. では.
やさしいCの書き方
これはEsolang Advent Calendar 2011の13日目の記事です.
といっても難解な言語ではなくやさしいCについての話をします. 題材は去年5分ネタでしゃべったカッコつけないプログラム[pdf]のp.19のソース. 今更ではありますが, 今でも時折「なんで動くの?」という感想を見かけるのでそのからくりについて.
ですので以下を読む前にとりあえず上記pdfをごらんください. で, ご存じの方は鼻で笑ってくださいまし. 以下は「え, なんで?」という方向けの説明です.
まずCでは変数の型を省略するとintになります(C++ではエラー). そしてポインタとintのsizeofが等しければ相互に変換できます(C++ではエラー).
main = 0x55535756;
mainは普通は関数ですが, intで宣言するとスタートアップ内から(&main)();と呼ばれます. すると 0x55535756;というバイト列を関数とみなしてプログラムが開始されることになります(そのためその領域に実行属性が与えられていないとエラーになります). ちなみに0xc3は関数から戻ってくるretを表すのですが
char main[] = "\xc3";
と
char *main = "\xc3";
とでは挙動が異なります. 前者はmainに突入するとすぐさま関数を抜けて戻って来るのに対して後者はsegmentation faultします. その違いはなんでしょうか. 考えてみるとよいでしょう.
さて0x55535756;をディスアセンブルすると
push esi push edi push ebx ...
と続きます. ですので後は基本的にはx86で階乗プログラム(X)を作れば終わりです. ただその関数の中からatoiやprintfなどのCの関数を呼ぶ必要があり, そのために少しだけテクニックを使ってます.
さて関数ポインタはWindowsとLinuxで異なるので両方で動かそうとするとXの中に埋め込むわけにはいきません. 間接的に呼び出すために, 位置をコントロールできる場所に関数ポインタを保持する必要があります. それが最初の2行です.
a = atoi; /* intのaをatoiの関数ポインタで初期化 */ b = printf; /* 同様 */
実はその次にはprintfに渡すための文字列を配置しています. つまりコードのレイアウトはこんな感じ.
(+ 0) pointer to atoi
(+ 4) pointer to printf
(+ 8) "%d! = %u\n"
(+20) push esi ; mainの先頭
(+21) push edi
...
main関数の中に突入した後, mainの中からatoiなどのポインタを保持した変数のアドレスを取得する必要があります. 相対アドレスはプログラム作成時に分かるのでmainのアドレスが分かればaやbのアドレスが分かります. ただmainのアドレスはOS毎に異なるので実行時にとる必要があります.
ところがx86自身には現在の実行アドレスを取得する命令がありません(x64にはあります). 代わりに自分の直後をcallで呼び出すことで現在アドレスをstackにpushさせ, それをpopでとるというテクニックを使いました.
call @f ; @@をcallして @@: pop esi ; popでesiにcall時のeipを設定
テクニック的にはそのぐらいなのですが, 実際に作るとなると多少の試行錯誤は必要です.
修正の度にバイナリコードサイズも変わり得るので関数ポインタを保持する変数への相対アドレスも変わったりします.
すると埋め込む即値も変わるのでなかなかやっかいです.
#asmで書いてアセンブルしたobjコードを見てサイズから手作業で書き直してhex化してコンパイルしてテスト? いやいやいや….
ですので私は「実行時に同じレイアウトとなるようなコードを生成してfac.cを生成するコード」を作成していました.
実行時コード生成にはXbyakを使います. 以下がそのコードです.
// fac-impl.cpp
#include <stdlib.h>
#include <stdio.h>
#include "xbyak/xbyak.h"
const size_t codeSize = 256;
struct {
int patoi; // -20
int pprintf; // -16
char str[12]; // -12
char code[codeSize];
} mem = {
(int)atoi,
(int)printf,
"%d! = %u\n",
};
#define GET_OFFSET(x) ((int)&(mem.code) - (int)&(mem.x))
struct Code : public Xbyak::CodeGenerator {
Code()
: Xbyak::CodeGenerator(codeSize, mem.code)
{
Xbyak::CodeArray::protect(mem.code, codeSize, true);
push(esi);
push(edi);
push(ebx);
push(ebp);
const size_t P = 4 * 4;
call("@f");
L("@@");
pop(esi);
const size_t adj = getSize() - 1;
mov(ebp, ptr [esp + P + 4]); // argc
cmp(ebp, 1);
je(".skip");
mov(ebp, ptr [esp + P + 8]); // argv
mov(ebp, ptr [ebp + 4]); // argv[1]
push(ebp);
mov(ebp, ptr [esi - adj - GET_OFFSET(patoi)]);
call(ebp);
add(esp, 4);
mov(ebp, eax);
L(".skip");
// ebp = n
sub(esp, 4 * 3);
lea(eax, ptr [esi - adj - GET_OFFSET(str)]);
mov(ptr [esp], eax);
mov(esi, ptr [esi - adj - GET_OFFSET(pprintf)]);
xor(ebx, ebx);
inc(ebx); // i = 1
mov(edi, ebx); // fac
L(".lp");
imul(edi, ebx);
mov(ptr [esp + 4], ebx);
mov(ptr [esp + 8], edi);
call(esi);
inc(ebx);
cmp(ebx, ebp);
jbe(".lp");
add(esp, 4 * 3);
pop(ebp);
pop(ebx);
pop(edi);
pop(esi);
ret();
}
};
int main(int argc, char *argv[])
{
// test code
Code code;
((void(*)(int,char*[]))code.getCode())(argc,argv);
// display code
puts("#include <stdio.h>");
puts("#include <stdlib.h>");
puts("a = atoi;");
puts("b = printf;");
char c = 'c';
size_t remain = sizeof(mem.str) + code.getSize();
const unsigned char *p = (unsigned char*)&mem.str;
do {
if (p == code.getCode()) {
printf("main = ");
} else {
printf("%c = ", c++);
}
union ic {
unsigned int i;
unsigned char c[4];
} ic;
size_t n = std::min(remain, 4U);
for (size_t i = 0; i < n; i++) {
ic.c[i] = *p++;
}
printf("0x%08x;\n", ic.i);
remain -= n;
} while (remain > 0);
}
GET_OFFSETは構造体のメンバの相対オフセットを得るためのマクロ. Codeの本体は与えられた引数から階乗を求めるプログラム.
mov(ebp, ptr [esi - adj - GET_OFFSET(patoi)]); call(ebp);
の部分でatoiのポインタを取得して呼び出しています. mainから後ろは生成した実行バイナリをa = ...; b = ...; の形で出力する部分です. byteでなくint単位で出力するためにパックしてます. これを実行すると問題のfac.cが生成されるのでした.
// cl fac-impl.cpp && ./fac-impl #include <stdio.h> #include <stdlib.h> a = atoi; b = printf; c = 0x20216425; d = 0x7525203d; e = 0x0000000a; main = 0x55535756; f = 0x000000e8; g = 0x6c8b5e00; h = 0xfd831424; i = 0x8b127401; j = 0x8b18246c; k = 0x8b55046d; l = 0xd5ffe36e; m = 0x8904c483; n = 0x0cec83c5; o = 0x89eb468d; p = 0x768b2404; q = 0x43db31e7; r = 0xaf0fdf89; s = 0x245c89fb; t = 0x247c8904; u = 0x43d6ff08; v = 0xee76eb39; w = 0x5d0cc483; x = 0xc35e5f5b;
それではみなさんもやさしいCを書いてみましょう.
#include <stdio.h>
FILE *openX(const char *fname)
{
FILE *fp = fopen(fname, "rb");
fprintf(stderr, "open %s %p\n", fname, fp);
return fp;
}
void closeX(FILE **fp)
{
fprintf(stderr, "close %p\n", *fp);
fclose(*fp);
}
int main()
{
puts("in");
{
__attribute__((cleanup(closeX)))
FILE *fp = openX("t.c");
}
puts("out");
return 0;
}
__attribute__のcleanupでmain内のfpのスコープが外れたときに呼ぶ関数を指定する.
その関数にはfpのポインタが渡される.
in open t.c 0x804a008 close 0x804a008 out# 釈迦に説法なんでしょうが….