DouKaku? Language detail: Matlab

Matlabらしくベクトル計算で。

a = kron(1:9, repmat(1, 1, 9));
b = repmat(1:9, 1, 9);
c = a.*b;
fprintf('%d * %d = %2d\n', [a; b; c])

正整数のゲーデル数化？ (Nested Flatten)

所詮は浮動小数点数なので大きなnについての精度は無考慮。

primes(100)は100までの素数25個を返す。 25と決めうちしているのは気持ち悪いが、double型が正確に表わせるのは高々十数桁であり、その範囲で使うぶんにはとりあえずはOKとする。桁数に応じた数の素数をとりだしたい場合は、第4行を次のように置き換えればよい:

if n < 6
  p = primes(11);  % Five smallest primes
else
  u = n*(log(n) + log(log(n)));  % An upper bound of the value of the nth prime
  p = primes(u);
end

上界uはピエール・デザルトの定理による。

function g = goedel(n)
  s = num2str(n) - '0';
  n = length(s);
  p = primes(100);
  g = prod(p(1:n).^s);

２進数の記述 (Nested Flatten)

ビルトイン関数 bin2dec() があります．

see: MATLAB Function Reference: bin2dec

1	bin2dec('01101001')

格子点の列挙 (Nested Flatten)

範囲の計算が変なまま出してしまったので再投稿。

同じく最後三行のコメントを外すと求めた格子点図示する。等距離の場合の順序をx軸を起点とした偏角の順にソートするという条件も満たしている。1000個目は [-8, 16]。

function r = latticepointsaroundorigin(num)
% r = latticepointssorted(n) retruns a list of lattice points,
% in the order of  distance from the origin and polar angle.
% Uncomment the last three lines to show the locations of the points.
% (ja.doukaku.org Q65)
rg = ceil(sqrt(num/pi));
[x y] = meshgrid(-rg:rg);
pts = [x(:) y(:) sqrt(x(:).^2+y(:).^2) atan2(-y(:),-x(:))];
r = sortrows(pts, [3,4]);
r = r(1:num,1:2);
% figure
% axis([-rg rg -rg rg])
% text(r(:,1),r(:,2),num2cell(1:num))

求める個数を引数として与えると、格子点を行列を返す。最後の三行のコメントを外すと求めた格子点をグラフ上に表示する。

function r = latticepointsaroundorigin(num)
% r = latticepointssorted(n) retruns a list of lattice points,
% in the order of  distance from the origin and polar angle.
% Uncomment the last three lines to show the locations of the points.
% (ja.doukaku.org Q65)
rg = ceil(sqrt(num)/2);
[x y] = meshgrid(-rg:rg);
pts = [x(:) y(:) x(:).^2+y(:).^2 atan2(-y(:),-x(:))];
r = sortrows(pts, [3,4]);
r = r(1:num,1:2);
% figure
% axis([-rg rg -rg rg])
% text(r(:,1),r(:,2),num2cell(1:num))

仲間はずれの判定 (Nested Flatten)

本質の部分はuniqueで種類を取り出してhistcで数えてsortで数の多い順にするところまでだが、その後各場合によって条件分岐して出力を整えるところはあまりすっきりしてない。

function [type majority singular] = uniformity(x)
%[type majority singular] = uniformity(x)
%  checks uniformity of vector x and returns
%  1 when x is uniform (or all elements are the same),
%  2 when x is uniform except for a single element with a different value,
%  or 3 when x is multiform (contains three or more values).
%  In the cases of 1 and 2, the value of the majority is returned. In the case
%  of 2, the value of the singular element is also returned.
% (ja.doukaku.org Q53)
majority = -1; singular = -1;
r = unique(x);
h = histc(x,r);
[hs k] = sort(h,'descend');
type = length(r);
if type==1
    majority = r(k(1));
elseif type==2 && hs(2)==1
    majority = r(k(1));
    singular = r(k(2));
else
    type = 3;
end

重複無し乱数 (Nested Flatten)

#2264の指摘をみてなるほどと思ったので簡単なテスト。順位相関を使ってみた。

#2285のコード(MATLAB標準のrandpermも同じ方法)は最初に乱数を要求された数列の個数分だけ作り、その順序を返す。MATLABのsortは値が同じ場合は順序を保存するので、たとえば1番目と2番目の乱数が同じ数字だと必ず1が2の前になってしまう。

投稿のコードは乱数の有効桁数(最初の行の最後の引数)を制限し、生成した数列とその位置の順位相関の分布と平均を表示する。有効桁数を1にした場合はかなり偏りが生じるが、3桁くらいでほとんど偏りが無い。順位相関を見る限りでは、普通の精度の乱数を使えばこのアルゴリズムでの偏りを気にする必要は無さそう。

x = bingo(25,1000,1);
rho = corr((1:size(x,2))',x(1:size(x,1),:)','type','Kendall');
figure(1); hist(rho,-1:0.1:1); axis([-1 1 0 300])
fprintf('Mean = %f, skewness = %f\n', mean(rho), skewness(rho));

function r = bingo(n,k,p)
% r = bingo(n,k) returns a k-by-n array, each row of which is a random
% permutation of 1:n, using random numbers with p significant digits .
[m r] = sort(randr([k n],p),2);

function r = randr(M,p)
r = floor(rand(M).*10^p)./10^p;

最初の投稿でえらそうに別解です、なんて書いたけどMATLABのrandpermは実は.mファイル(つまりMATLABで書かれたスクリプト)でその中身は投稿したコードと全く同じだということに気付いた。情けないのであらためて別解。前と同様第二引数で何個生成するか指定する。sortを使う方法よりずいぶん遅い。次のようにすると計算時間の測定と順序が均等かどうかのおおまかな確認ができる。

tic;r=bingo(100,10000);bar(sum(r,1));toc

function r = bingo(n,k)
% r = bingo(n) retruns a random permutation of 1:n.
% r = bingo(n,k) returns a k-by-n array, each row of which is a random
% permutation of 1:n.
% (ja.doukaku.org Q46)
if nargin<2
    k = 1;
end
[s r] = meshgrid(1:k,1:n);
% q = zeros(k,2);
for p = n:-1:2
    q(:,1) = ceil(rand(k,1)*p);
    q(:,2) = p;
    q = q+s(1:2,:)'*n-n;
    r(q) = r(q(:,[2 1]));
end
r = r';

標準で必要な関数があるのにあえて別解を用いるメリットが最初の例では特に無かったので、少し書き換えてk個の列を一度に返せるようにした。標準のrandpermでk個生成するにはループでk回呼ばないといけない。余談だがMATLAB7.4ではデフォルトの乱数生成アルゴリズムとして#2269さんのリンク先でも紹介されているメルセンヌ・ツイスタが採用されている。

function r = bingo(n,k)
% r = bingo(n) retruns a random permutation of 1:n.
% r = bingo(n,k) returns a k-by-n array, each row of which is a random
% permutation of 1:n.
% (ja.doukaku.org Q46)
if nargin<2
    k = 1;
end
[m r] = sort(rand([k n]),2);

標準のrandpermという関数が問のbingoと全く同じ動作をします。投稿のコードはsortがソート前の順序を返すことを使った別解です。

1 2	function r = bingo(n) [m r] = sort(rand(1,n));

ダブル完全数 (Nested Flatten)

ループを使わないで書いてみた。手元のマシンで2分くらいかかった。メモリを食うし、ループを使う方がずっと速い。

function r = doubleperfectnumber(rmax)
% Return every integer, the double of which is equal to the sum of its divisors
% (ja.doukaku.org Q25).
[m n] = meshgrid(uint16(1:rmax-1),uint16(1:rmax));
m(find(m(:)>=n(:))) = 0;
m(find(mod(n,m)~=0)) = 0;
r = find(sum(m,2)==n(:,1)*2);
% plot(1:rmax,sum(m,2),1:rmax,n(:,1));

倍数になる13進数 (Nested Flatten)

ループを使う方法は既に出ているので、ループを使わない方法で書いてみた。順番にチェックするのではなく全部計算するので、調べる範囲の上限を与えなければいけないしメモリも食う。条件を満たす最小の数を返せというこの問題に対しては効率のいい方法とは言えない。例えばmからnまでの整数で条件を満たすものを全て見つけろ、というような問題だとするとループより速いかも? 最後2行のコメントを外すと13進数とみなす場合と10進数とみなす場合との比をグラフに表示する。この比が整数となる数字列が条件を満たすことになる。最後から3行目のminを外すと範囲内の全部の答を出力する。

function n = finddoubleinb13(rmax)
% Returns a sequence of digits that is twice bigger in base 13 than in base 10.
% The numbers from 10 to rmax in base 10 will be checked (ja.doukaku.org Q14).
n10 = (10:rmax)';
n13 = base2dec(num2str(n10),13);
idx = find(~mod(n13(:),n10(:)));
n = min(n10(idx));
% plot(n10,n13(:)./n10(:))
% grid on

全ての組み合わせ (Nested Flatten)

とりあえず御言葉に甘えて数字のみ対応。

function c = multisetcomb(varargin)
% return all combinations of items in given column vectors by picking up one
% item from each at once.
s1 = varargin{1};
s2 = varargin{2};
c = [repmat(s1,size(s2,1),1) sortrows(repmat(s2,size(s1,1),1))];
if length(varargin) > 2
    c = multisetcomb(c,varargin{3:end});
end

長方形の交差判定 (Nested Flatten)

標準でrectintという関数があるので引数を合わせただけ。コメントアウトしてある最後の三行は交差の状態を図示。

function a = checkrectoverlap(r1, r2)
%  returns positive value if two rectangles are overlapped, zero if not.
%  r1,r2: [left top right bottom]
rr1=r1-[0 0 r1(1:2)];
rr2=r2-[0 0 r2(1:2)];
a = rectint(rr1,rr2);
% rectangle('Position',rr1);
% rectangle('Position',rr2);
% axis equal

アレイのuniq (Nested Flatten)

MATLABの標準関数uniqueはソートした結果を返すが、オプションを与えることによって各要素のソート前の位置を得ることができる。それを用いて元の順序を復元する。

function y = unsortuniq(x)
[b i] = unique(x,'first');
bi = sortrows([b;i]',2)';
y = bi(1,:);

呼んだのは誰？ (Nested Flatten)

dbstackでファイル名やそのファイル名の中の何行目で呼ばれたかも取得可能。デバッグ目的の関数なので。

1
2
3

function foo
st = dbstack;
st(2).name

リストを逆順に表示 (Nested Flatten)

MATLABは多くの関数が入力としてスカラーだけでなくベクトルを受け入れる。fprintfもそういう関数の一つなのでベクトルを与えると個々の要素に対してフォーマット付き出力をしてくれる。これを利用。ただし要素の型に応じてフォーマット指定文字列の部分を変更しないといけないのが難点。投稿のコードは整数を想定。ベクトルを逆順にする部分はfliplrという関数が用意されてるけどあえて違う書き方で。

1 2	function revdisp(v) fprintf('%d\n',v(length(v):-1:1));

「組合せ型の最小完全ハッシュ関数」の逆関数 (Nested Flatten)

とりあえず簡潔に。

function data = invoptcombhash(x, n, m)
combpos = nchoosek(1:n,m);
data = zeros(1,n);
data(combpos(x,:)) = 1;

与えられた数字のケタ数 (Nested Flatten)

常用対数でやってみました。

使用例：
>> [d,m] = doukaku40(321321)
d =
     6
m =
      100000

1
2
3

function [d,m] = doukaku40(n)
  d = floor(log10(n)) + 1;
  m = 10^(d-1);

Hello, world! (Nested Flatten)

Matlabで単純に

1	disp('Hello, world!');

DouKaku?

Language detail: Matlab

Unsolved challenges

codes

Index

Feed

Other

Link

Pathtraq