プリクラ問題

covering designのt=2に相当、
La Jolla Covering Repository Tablesに全部載ってます…
リンク先クリックすればいいんだね。
Combinatorial covering designs に文献とかあるから、なんかもう・・・
もういいや

fortran90によるのコードはこのページの下の方です
だんだんと条件を付け加えて更新していくつもりです。n=9~10でm=3,4程度は厳しいです。



力技でやってみたよ!
主眼が回数になっているので、実際の組み合わせは?と気になってやったものです。

事の発端

コンピュータによる総当たりでは今のままではn=10位が限界です。
だんだんと改良して増やせていければ、と。
コードは後ほど僕に時間ができ次第、公開します。



より↓

念のための確認もかねて、計算は表の値-1から確かめて総当たりを行った結果です。なので具体的な組み合わせが載っているものは正しいはずです。
組み合わせの総数はおおよそ
\(\displaystyle \sum_{q=3}^{q^{Answer}} {}_{{}_nC_m}C_{q}\)
という膨大な数の組み合わせとなります。ここで\(q^{Answer}\)はプリクラ問題の最小回数です。
\(q=3\)はあり得ない、全てを検証する必要はない、など確実にわかる要素は多いですが、例えばn=9,m=3,解q=12を考えると
\(\displaystyle {}_{{}_9C_3}C_{12}=112,992,892,764,570 \ \ \sim 10^{14}(=100T)\)回
の計算が必要です。理想的な状況を考え、cpu10GHz,コア数10で並列処理で行う場合、計算時間は
\(\displaystyle \frac{10^{14}回}{10\times 10^{9}Hz\cdot 10\mbox{コア}}=1000[秒]\)
なのでぎりぎり計算できそうかなーくらいになります。
しかし、このあたりが限界で、n=10,m=3,解q=17になると、
\(\displaystyle {}_{{}_{10}C_3}C_{17}=189,916,591,435,396,829,640\ \ \sim 10^{20}(100E)\)回
の計算となります。先ほどの設定で行うならば、\(10^{9}\)秒~31.7年かかります。
スパコンの京(10PHz)をフルに使えるならば1000[s]で終わります。しかしそこで終わりです。n=11は…もう無理です。

効率のいい方法を考えない限り、解は得られません。

n\m 2 3 4 5 6 7 8 9 10 11
2 1 1 1 1 1 1 1 1 1 1
3 3 1 1 1 1 1 1 1 1 1
4 6 3 1 1 1 1 1 1 1 1
5 10 4 3 1 1 1 1 1 1 1
6 15 6 3 3 1 1 1 1 1 1
7 21 7 5 3 3 1 1 1 1 1
8 28 11 6 4 3 3 1 1 1 1
9 36 12 8 5 3 3 3 1 1 1
10 45 17 9 6 4 3 3 3 1 1
11 55 19 11 7 6 4 3 3 3 1

(3,2)=3

1 2
1 3
2 3

(4,2)=6

(4,3)=3

1 2 3
1 2 4
1 3 4

(5,2)=10

(5,3)=4

1 2 3
1 2 4
1 2 5
3 4 5

(5,4)=3

1 2 3 4
1 2 3 5
1 2 4 5

(6,2)=15

(6,3)=6

1 2 3
1 2 4
1 5 6
2 5 6
3 4 5
3 4 6

(6,4)=3

1 2 3 4
1 2 5 6
3 4 5 6

(6,5)=3

1 2 3 4 5
1 2 3 4 6
1 2 3 5 6

(7,2)=21

(7,3)=7

1 2 3
1 4 5
1 6 7
2 4 6
2 5 7
3 4 7
3 5 6

(7,4)=5

1 2 3 4
1 2 3 5
1 2 3 6
1 2 3 7
4 5 6 7
or
1 2 3 5
1 4 6 7
2 3 4 5
2 4 6 7
3 5 6 7

(7,5)=3

1 2 3 4 5
1 2 3 6 7
1 4 5 6 7

(7,6)=3

1 2 3 4 5 6
1 2 3 4 5 7
1 2 3 4 6 7

(8,2)=28

(8,3)=11

計算中…

(8,4)=6

1 2 3 4
1 2 5 6
1 2 7 8
3 4 5 6
3 4 7 8
5 6 7 8

(8,5)=4

1 2 3 4 5
1 2 3 4 6
1 2 3 7 8
4 5 6 7 8

(8,6)=3

1 2 3 4 5 6
1 2 3 4 7 8
1 2 5 6 7 8

(8,7)=3

1 2 3 4 5 6 7
1 2 3 4 5 6 8
1 2 3 4 5 7 8

(9,2)=36

(9,3)=12

計算中…

(9,4)=8

1 2 3 5
1 2 3 6
1 2 4 7
1 2 8 9
1 3 4 8
1 3 7 9
4 5 6 9
5 6 7 8


info)cputime=3.0415[s], Nc=10523209, com=1 2 35 115 126

(9,5)=5

1 2 3 4 5
1 2 3 4 6
1 2 7 8 9
3 4 7 8 9
5 6 7 8 9

(9,6)=3

1 2 3 4 5 6
1 2 3 7 8 9
4 5 6 7 8 9

(9,7)=3

1 2 3 4 5 6 7
1 2 3 4 5 8 9
1 2 3 6 7 8 9

(9,8)=3

1 2 3 4 5 6 7 8
1 2 3 4 5 6 7 9
1 2 3 4 5 6 8 9

(10,2)=45

(10,3)=17

計算厳しい

(10,4)=9

計算中(4日後くらい?)

(10,5)=6

(10,6)=4

1 2 3 4 5 6
1 2 3 4 7 8
1 2 3 4 9 10
5 6 7 8 9 10

(10,7)=3

1 2 3 4 5 6 7
1 2 3 4 8 9 10
1 5 6 7 8 9 10

(10,8)=3

1 2 3 4 5 6 7 8
1 2 3 4 5 6 9 10
1 2 3 4 7 8 9 10

(10,9)=3

1 2 3 4 5 6 7 8 9
1 2 3 4 5 6 7 8 10
1 2 3 4 5 6 7 9 10

(11,2)=55

(11,3)=19

(11,4)=11

(11,5)=7

(11,6)=

(11,7)=4

1 2 3 4 5 6 7
1 2 3 4 5 6 8
1 2 3 4 9 10 11
5 6 7 8 9 10 11

(11,8)=3

1 2 3 4 5 6 7 8
1 2 3 4 5 9 10 11
1 2 6 7 8 9 10 11

(11,9)=3

1 2 3 4 5 6 7 8 9
1 2 3 4 5 6 7 10 11
1 2 3 4 5 8 9 10 11

(11,10)=3

1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 11
1 2 3 4 5 6 7 8 10 11



fortran90のコード。総当たりで調べます。
別にinputファイルが必要です。n:m=3:2以上の時は考慮していません。3回で終わるからね。
下の例はn=6, m=3に相当しています。(inputのrはmに相当します。)
qiは調べ始める最低回数、qeは調べ終わる最高回数です。
出力はこうなります。

 M          35
q  -->    3
q  -->    4
q  -->    5
q  -->    6
q  -->    7
 number of calculation     1276464
 cpu_time  0.905050993     [sec]
+------------+
  1  2  3
  1  4  5
  1  6  7
  2  4  6
  2  5  7
  3  4  7
  3  5  6
+------------+

ここからfortran90のコード。gfortranで確かめています。

program main
  implicit none
  integer::n,r,qi,qe,qc
  integer::i,j
  real::t1,t2
  integer,allocatable::groups(:,:)
 
  n=7
  r=3
  qi=3
  qe=10
  allocate(groups(1:qe,1:r))
  groups=0

  call cpu_time(t1)  
  call pkprob(n,r,qi,qe,qc,groups)
  call cpu_time(t2)
  write(6,*)"cpu_time",t2-t1,"[sec]"

  write(6,'(A)')"+------------+"
  do i=1,qc
     do j=1,r
        write(6,'(i3)',advance='no')groups(i,j)
     enddo
     write(6,*)          
  enddo
  write(6,'(A)')"+------------+"
 
  stop
end program main
!==========================
subroutine pkprob(n,r,qi,qe,qc,groups)
  !Covering Design t=2
  implicit none
  integer,intent(in)::n,r,qi,qe
  integer,intent(out)::qc,groups(1:qe,1:r)
 
  integer::q,M,ic,ncalc
  integer::i,j,q2,k
  integer,allocatable::ini(:),tmp(:),cgroup(:,:),pk(:,:)
  integer,allocatable::choose(:)
 
  integer::check1,nCr
  integer::csign
  !  logical::csign,gsign
  external::nCr

  !Combination M=nCr
  M=nCr(n,r)  
  write(6,*)"M",M

  !csign      : Sign of critical combination of array  
  !ic         : Critical value of i which number of combination of combination
  !qc         : Critical value of q which number of combination of combination
  !pk(1:n,1:n): Checking array if every number fill(1) combination or not(0)
  !ncount(1:n): To check calculation number of times
  allocate(ini(1:r),tmp(1:r),cgroup(1:M,1:r))
  ini(1:r)=0
  cgroup(1:M,1:r)=0
  tmp(1:r)=0
 
  do j=1,r
     ini(j)=j
     tmp(j)=ini(j)
  enddo
  !Substitute combination to cgroup.
  do i=1,M
     cgroup(i,1:r)=tmp(1:r)
     call combination(n,r,ini,tmp)
  enddo
  deallocate(ini,tmp)
  !csign      : Sign of critical combination of array  
  !ic         : Critical value of i which number of combination of combination
  !qc         : Critical value of q which number of combination of combination
  !pk(1:n,1:n): Checking array if every number fill(1) combination or not(0)
  !ncount(1:n): To check calculation number of times
  csign=0
  ic=0
  qc=0
  allocate(pk(1:n,1:n))
  pk=0

  ncalc=0

  !q : We need to consider below combination.
  !  Example ) n=5, r=3 case.
  !
  !cgroup(i,1:r)
  !i\r| 1 2 3
  !----------
  ! 1 | 1 2 3
  ! 2 | 1 2 4
  ! 3 | 1 2 5
  ! 4 | 1 3 4
  ! 5 | 1 3 5
  ! 6 | 1 4 5
  ! 7 | 2 3 4
  ! 8 | 2 3 5
  ! 9 | 2 4 5
  !10 | 3 4 5
  !
  !We have to choose q= 3 or 4 or 5 ...combination of "cgroup".
  !Sum total of these combination can write D=_{nCr}C_q.
  !Specifically,
  !---choose q groups from cgroup
  ! +----q=3----+ +-----q=4-----+ +------q=5------+
  ! |choose(1:q)| | choose(1:q) | |  choose(1:q)  |
  ! | l |       | | l |         | | l |           |
  ! |===|=======| |===|=========| |===|===========|
  ! | 1 | 1 2 3 | | 1 | 1 2 3 4 | | 1 | 1 2 3 4 5 |
  ! | 2 | 1 2 4 | | 2 | 1 2 3 5 | | 2 | 1 2 3 4 6 |
  ! | 3 | 1 2 5 | | 3 | 1 2 3 6 | | 3 | 1 2 3 4 7 |
  ! | 4 | 1 2 6 | | 4 | 1 2 4 7 | | 4 | 1 2 3 4 8 |
  ! | 5 | 1 2 7 | | 5 | 1 2 4 8 | | 5 | 1 2 3 4 9 |
  ! | 6 | 1 2 8 | | 6 | 1 2 5 9 | | 6 | 1 2 3 4 10|
  ! | 7 | 1 2 9 | | 7 | 1 2 5 10| | 7 | 1 2 3 5 6 |
  ! |...| . . . | |...| . . . . | |...| . . . .   |
  ! |120| 4 5 6 | |210| 7 8 9 10| |252| 6 7 8 9 10|
  ! +-----------+ +-------------+ +---------------+
  !
  ! If q=4, l=3
  !  We choose below combination, and check satisfy condition or not by "pk".
  ! 1 | 1 2 3
  ! 2 | 1 2 4
  ! 3 | 1 2 5
  ! 6 | 1 4 5
 
  do q=qi,qe
     write(6,'(A,i5)')"q  -->",q
     allocate(choose(1:q),ini(1:q))
     choose(1:q)=0
     ini(1:q)=0
     
     do j=1,q
        ini(j)=j
        choose(j)=ini(j)
     enddo
         
     csign=0
     do while(csign.eq.0)
        !condition 1 to reduce calculation
        check1=cgroup(choose(1),1)
        if(check1.ne.1)exit

        !condition 2 to reduce calculation
        check1=cgroup(choose(2),1)
        if(check1.ge.2)exit

        !condition 3 to reduce calculation for less than n:r=3:2
        check1=cgroup(choose(q),1)
        if(check1.ge.3)then
           pk=0
           do q2=1,q
              do j=1,r
                 do k=1,r
                    pk(cgroup(choose(q2),j),cgroup(choose(q2),k))=1
                 enddo
              enddo
           enddo
           !memory calculation number of times
           ncalc=ncalc+1
           if(minval(pk).eq.1)then
              csign=1
              qc=q
           endif

           if(csign.eq.1)exit
        endif

        !go forward choose(1:q)
        call combination(M,q,ini,choose)
     enddo
     
     if(csign.eq.1)exit
     deallocate(choose,ini)
  enddo
  write(6,*)"number of calculation",ncalc

  groups(1:qe,1:r)=0
  do i=1,qc
     do j=1,r
        groups(i,j)=cgroup(choose(i),j)
     enddo
  enddo
 
  return
end subroutine pkprob
!---------------------------------
subroutine combination(n,r,ini,arr)
  implicit none  
  integer(8)::i,n,r,ini(1:r),bef(1:r),arr(1:r)
  integer(8)::numx
  logical::key(1:r)
 
  bef(1:r)=arr(1:r)-ini(1:r)
  arr(1:r)=0
  key(1:r)=.false.

  numx=n-r
  do i=1,r
     if(bef(i).eq.numx)key(i)=.true.
  enddo

  do i=1,r-1
     if(key(i+1))then
        if(key(i))then
           if(i.ne.1)arr(i)=arr(i-1)      
        else
           arr(i)=bef(i)+1
        endif
     else
        arr(i)=bef(i)
     endif
  enddo
  if(key(r))then
     arr(r)=arr(r-1)
  else
     arr(r)=bef(r)+1
  endif
 
  arr(1:r)=arr(1:r)+ini(1:r)

  return
end subroutine combination
!----------------------------------
function nCr(n,r)
  implicit none
  integer(8)::n,r,i,r0,nCr

  r0=n-r
  if(r.le.r0)r0=r

  nCr=1
  do i=1,r0
     nCr=nCr*(n-r0+i)
     nCr=nCr/i
  enddo
 
  return
end function nCr

LaTeXテンプレート

メモ書き用のLaTeXテンプレートです。

\documentclass[11pt, a4j]{jarticle}
%用紙のサイズ
%setlength{paperwidth}{210mm}
\setlength{\paperwidth}{200mm}
\setlength{\paperheight}{297mm
}
%上側余白
\setlength{\voffset}{-18mm}
%\setlength{\voffset}{0mm}
\setlength{\topmargin}{0mm}
%ヘッダー
%\setlength{\headsep}{0mm}
%\setlength{\headheight}{0mm}
%傍注欄
%\setlength{\marginparsep}{0mm}
%\setlength{\marginparwidth}{0mm}
%\setlength{\marginparpush}{0mm}
%本文スペースの大きさ
\setlength{\textheight}{257mm}
\setlength{\textwidth}{170mm
}
%%ページ左側の余白設定
\setlength{\hoffset}{0mm}
\setlength{\oddsidemargin}{-4mm}
\setlength{\evensidemargin}{-4mm
}
%フッター
%\setlength{\footskip}{10mm}

%画像表示のために必要
\usepackage[dvipdfmx]{graphicx}
%数式内で太字を使えるようにするbm{A}。http://www.latex-cmd.com/equation/vector.html
\usepackage{bm}
%華文字フォントを使えるようにする。http://www.biwako.shiga-u.ac.jp/sensei/kumazawa/tex/mathrsfs.html
\usepackage{mathrsfs}
%複雑な数式記号、数式フォントに対応する。http://hooktail.sub.jp/tex/symbol/amsmathsymbols.html
\usepackage{amsmath}
%ディラックのブラケット記号を使えるように。http://dayinthelife.at.webry.info/200705/article_2.html
\usepackage{braket}
%図のキャプションの文字を小さくする。http://ozuma.txt-nifty.com/blog/2004/01/latexcaption.html、http://www.ctan.org/tex-archive/macros/latex209/contrib/misc
\usepackage{hangcaption}
%for use hangcaption.sty
\captionwidth=0.9\textwidth %本文サイズの90%
\newcommand{\captionfonts}{\footnotesize} %キャプションをsmallフォントで
\makeatletter % Allow the use of @ in command names
\long\def\@makecaption#1#2{%
\vskip\abovecaptionskip
\sbox\@tempboxa{{\captionfonts #1: #2}}%
\ifdim \wd\@tempboxa >\hsize
{\captionfonts #1: #2\par}
\else
\hbox to\hsize{\hfil\box\@tempboxa\hfil
}%
\fi
\vskip\belowcaptionskip}
\makeatother % Cancel the effect of \makeatletter
%for use hangcaption.sty end.
%数式がページを跨いでもよいようにする。http://tex.stackexchange.com/questions/57338/how-can-i-allow-gather-to-break-between-pagesより。
\allowdisplaybreaks[4]
%式や図を参照したい時の新たな定義式を設定。
\newcommand{\eq}[1]{式~(\ref{#1})}
\newcommand{\fig}[1]{図~\ref{#1}
}
%図とキャプションの間を狭める。
\setlength \abovecaptionskip{-5pt}
%jreportだと"関連図書"になってしまうのを"参考文献"に変える。
%|-- この時、
%|-- jreportの場合は\bibname
% -- jarticleの場合は\refnameにすること。
\renewcommand{\refname}{参考文献}
%画像の場所を.tex以外の場所に置きたい時、.texからのパスを記述。
\graphicspath{{./fig/pic1/}{./fig/pic2/}{./fig/pic1/aa1/}}
   
\begin{document}
 \section{a}
    ここに本文を書いていきます。
    \subsetion{a_sub}
    ...
\section{b
}
    ...
%   画像を載せる場合   
%   \begin{figure}[th]
%      \centering
%      \includegraphics[width=0.8\textwidth]{test.eps}
%      %キャプション内の文章の間隔を変更
%      \renewcommand{\baselinestretch}{1}
%      \caption{ここにキャプションを}
%      \label{fig:TBA}
%    \end{figure}

    ...
\end{document}

latexコマンド



alignの省略はできない!?


\newcommandによる省略はできません。
LATEX community, newcommand with \end{align*}
によると一応できます。しかし、推奨はされてません。
そのやり方は、

\makeatletter
\newenvironment{ma}
  {\start@align\@ne\st@rredtrue\m@ne}
  {\endalign}
\makeatothe

をlatexのプリアンブルに記述し、

\begin{ma}
  a &=b\\
  &=c
\end{ma}

とすれば、確かに省略はできます。ただし、\beginは書かなくてはならないので、若干省略とは違う感じです。

figure環境について


figure環境テンプレ


プリアンブルに記述すること

\usepackage[dvipdfmx]{graphicx}
    %画像挿入のために必要なパッケージ
\setlength \abovecaptionskip{-5pt}
    %図とキャプションの間を狭める
\graphicspath{{./fig/}{./fig/abc/}}
    %パス(.tex以外の場所の画像を参照したい時)

figure環境内に記述すること

\begin{figure}
    \centering
    \includegraphics[width=\hsize]{data.eps}
    \label{data1234}
    \renewcommand{\baselinestretch}{1}
    \caption{この図は(-ω-;)を表す。}
\end{figure}

別のフォルダにある画像を参照する


dir/
├ test.tex
├ test.aux
├ test.dvi
├ …
data.eps
fig/
      ├ data2.eps
      └ abc/
            └ rrr.eps
通常のままでは.texと同じ階層に画像ファイルを置いておかなくては参照されません。
すなわち上の構造である場合、参照されるのはdata.epsのみであり、data2.eps,rrr.epsは参照されません。
参照されるようにするためにはlatexのプリアンブルに
\graphicspath{{./fig/}{./fig/abc/}}
と記述します。そのうえで
\includegraphics[width=\hsize]{data2.eps}
と書くと参照されるようになります。

図とキャプションの間を調節する


図とキャプションの間は通常
+—fig通常—+
位の広さを持ちます。この図とキャプションの間を狭くするためにはプリアンブルに
\setlength \abovecaptionskip{-5pt}
と記述します。ここでは-5ptとしていますが、各々調節して最適な値を探してください。
+—fig(狭くした例)—+

図のキャプションの間隔を狭くする


通常の図は、キャプションの間隔が広いと感じます。
+—fig—+
これの調節は、figure環境の中に
\renewcommand{\baselinestretch}{1}
という言葉を入れましょう。{1}は行間の倍率を表すそうです。行間の調整より。
+—fig—+

LXDE

LinuxMint 17.1 “Rebecca” においてデスクトップ環境”lxde”をインストールするための手順です。

LXDEの導入

lxdeとは?
Arch Linux LXDE (日本語)より引用

LXDE (“Lightweight X11 Desktop Environment”) は極めて高いパフォーマンスを持ち省エネルギーなデスクトップ環境です。開発者たちの国際的なコミュニティにより管理され、美しいインターフェース、多言語対応、標準的なキーボードショートカットとタブファイラーのような追加機能を含みます。LXDE は CPU も RAM も他のデスクトップ環境に比べ消費が少ないです。特に、ネットブック、モバイル端末や、古い計算機などによるクラウドコンピューティングにあわせてデザインされています。

ぶっちゃければ軽いよ!ということです。
windowsだったらこの変更に相当する操作はありません。
これからやろうとしているのは
OS → バージョン → デスクトップ環境
の順で書けば、
LinuxMint → 17.1(Rebecca) → cinnamon
から
LinuxMint → 17.1(Rebecca) → LXDE
に変更する、という操作です。cinnamon, LXDEの他にMate, KDE, Xfceなどがあります。
詳しくは調べてください。

設定→Synapticパッケージマネージャ
と進みます。
Synapticパッケージマネージャの検索ボックスでlxdeを検索してチェック、チェックして出てくる関連のものも全てチェックします。
さらに、lxdmもチェックします。
どうやらこのlxdmがないとログイン画面でどのデスクトップ環境でスタートするか?の選択ができません。

あとはSynapticパッケージマネージャの”適用”を押せば終了です。ログアウトして、
ログイン画面で(多分デフォルトでは)λになっているところを押してlxdeで起動を選択、その後ログインしましょう。


lxterminal のカラーを変えたい!

  …できません。

Thread: Colour Scheme on lxterminal in lubuntu

LXTerminal Colors?
より。

At this time we do not offer any more customization than foreground/background. You are encouraged to run any terminal emulator that meets your needs.

訳すと、”前景色と背景色以上のカスタマイズは提供しないよ。カスタムしたかったら別のターミナルをインストールしてねミ☆”
だそうです。

→guake terminalというのがよさそうです。

sudo apt-get install guake

で、インストールが終わったら端末で

guake

と入力し、F12キーを押せば立ち上がります。guake terminalの設定は
スタートボタン→設定→Guake Preferences
で設定できます。


LXDEのスタートアップ(起動時に自動的に実行させるプログラム)の設定がないよ!

Arch LXDEによると、
ディレクトリの場所:~/.config/lxsession/LXDE/autostart
にLXDEを起動させたときに実行させるプログラムが書いてあります。その中身に何が書いてあるかというと、

@xscreensaver -no-splash
@lxpanel --profile LXDE
@pcmanfm --desktop --profile LXDE
@/usr/lib/policykit-1-gnome/polkit-gnome-authentication-agent-1

のような文が記述されていると思います。意味はArch LXDEによると@以下の文が起動時のプログラムとして実行されます。
なので、ここを書き換えて、たとえばguake端末を起動時に立ち上げたいのならば

@xscreensaver -no-splash
@lxpanel --profile LXDE
@pcmanfm --desktop --profile LXDE
@/usr/lib/policykit-1-gnome/polkit-gnome-authentication-agent-1
@guake

という風に@guakeを付け加えてあげましょう。


guakeの横幅を変えたいよ!

How To Adjust Guake Terminal Width を見てみるとUbuntu12.10以降はguakeの設定ファイルは

/usr/bin/guake

にあります。LinuxMint17.1でも同じで、上のディレクトリにあります。この中身を編集していきます。

sudo emacs /usr/bin/guake

で開き、820~830行目あたりにある

def get_final_window_rect(self):
    """Gets the final size of the main window of guake. The height
    is the window_height property, width is window_width and the
    horizontal alignment is given by window_alignment.
    """

    screen = self.window.get_screen()
    height = self.client.get_int(KEY('/general/window_height'))
    width = 100
    halignment = self.client.get_int(KEY('/general/window_halignment'))

のwidth=100になっているところ書き換えてwidth=50とか80とかにします。
ちなみにこのwidthはウインドウの横幅に対するパーセンテージです。

※サイトで検索した時、下の場所にあるからねーっていう記述がかなりあったけれどこれはUbuntu12.04以前の場合です。
/usr/lib/guake/guake.py
昔のだね。

フォントサイズ変更

emacsでフォントサイズを変更したい場合は?環境はemacs24。

(調べてもよくわからないし出てこない…。)
一度だけ変更したい場合は、Ctrl-x-± で可能です。
ではずっと変更したい場合は?
emacsの設定ファイルに記述します。emacsの設定ファイルの場所は
emacswiki – InitFileより

  1. For GnuEmacs, it is ~/.emacs or _emacs or ~/.emacs.d/init.el.
  2. For XEmacs, it is ~/.xemacs or ~/.xemacs/init.el.
  3. For AquamacsEmacs, it is ~/.emacs or ~/Library/Preferences/Aquamacs Emacs/Preferences.el

にあります。環境によって場所が異なります。
Linuxmint17.1,”sudo apt-get install emacs24″で行った場合、おそらく~/.emacsでしょう。

21.8 Fontsを見るとデフォルトでは
フォント(の種類) : monospace font
フォントサイズ : 10pt
として設定されているそうです。

テキストエディタで~/.emacsを開き、
(custom-set-variables


)
の括弧の

(add-to-list 'default-frame-alist
                       '(font . "DejaVu Sans Mono-12"))

という文を記述します。
現在、フォントのタイプの変更は僕は分かりません。多分DejaVu…を変えればよさそう。
12のところを10とか14とかに変えればデフォルトのフォントサイズを変更することができます。

見出しに使っている画像は8,9,10,11,12,20の場合の比較です。参考にどうぞ。

僕の使ってる~./emacsはこちら

LinuxMint導入時にやること

とはいってもVMwarePlayerを使った環境を想定しています。

環境は
LinuxMint 17.1 Quina

  • アップデート
  • sudo apt-get update
    sudo apt-get dist-upgrade
  • VMware tools
  • ※導入されているかの確認は ($ service vmware-tools status で動いていたらvmware-tools start/running という文字列が出ます。)
    20160123-151006_c
    を押してからの手順となります。

    sudo su
        パスワード入力
    mkdir /mnt/cdrom
    mount /dev/cdrom /mnt/cdrom
            mount:ブロックデバイス /dev/sr0 は書き込み禁止です、読み込み専用でマウントします. という文が出てくる。
    cd /tmp
    rm -rf vmware-tools-distrib
    tar zxpf /mnt/cdrom/VMwareTools-9.6.2-1688356.tar.gz
        /mnt/cdrom/VMw くらいまで打ってTabキーで補完してください。もしかしたらバージョンが上のではないかもしれない。(2014/12/28 確認)
    cd vmware-tools-distrib
    ./vmware-install.pl
  • emacs24
  • sudo apt-get install emacs24

    ホームディレクトリに『.emacs』という名前のファイルを作り、以下を記述する。

    (custom-set-variables
     ;; custom-set-variables was added by Custom.
     ;; If you edit it by hand, you could mess it up, so be careful.
     ;; Your init file should contain only one such instance.
     ;; If there is more than one, they won't work right.
     ;;'(font-lock-maximum-size 512000)
     '(font-use-system-font t)
     '(inhibit-startup-screen t)
    )
    (add-to-list 'default-frame-alist
                           '(font . "DejaVu Sans Mono-12"))

    (custom-set-faces
     ;; custom-set-faces was added by Custom.
     ;; If you edit it by hand, you could mess it up, so be careful.
     ;; Your init file should contain only one such instance.
     ;; If there is more than one, they won't work right.
     )
  • 日本語環境
  • ibus+Anthy
    Software Manager(日本語だったら “ソフトウェアの管理” )

    “ibus”と”ibus-anthy”を検索してそれらをインストール

    再起動したのち、右下のところにのキーボードのアイコンがあるのでそれをクリック、日本語-Anthy を選択すれば日本語入力ができる。

  • Linuxのディレクトリ名を英語に変更
  • LANG=C xdg-user-dirs-gtk-update
  • htop
  • sudo apt-get install htop
  • gnuplot
  • sudo apt-get install gnuplot

    ホームディレクトリに初期設定用のファイル『.gnuplot』を作り、以下を記述する。

    set terminal wxt noraise
    set palette defined(0"#ffffff",0.8"#00008b",1.8"#2ca9e1",3"#008000",4.2"#ffff00",5"#eb6101",5.5"#8b0000")
  • gfortran
  • sudo apt-get install gfortran
  • ファイルオープン数の変更とスタックサイズの変更(あまり推奨はされてません。)
  • sudo emacs /etc/security/limits.conf

    で、

        <domain>        <type>      <item>      <value>
        sikino         soft    stack          131072
        sikino         hard    stack          131072
        sikino         soft    nofile         16384
        sikino         hard    nofile         16384

    ※sikinoは僕のユーザー名であり、あなた自身のユーザー名にそれぞれ変更してください。

  • imagemagick(convertコマンドとか使えるようにするあれ)
  • sudo apt-get install imagemagick

    (http://blog.mirakui.com/entry/20110123/1295795409)

  • exfat形式の認識
  • →VMwareの場合、|Player(P)▼|から、管理→仮想マシンの設定を選択、
    USBコントローラ 
    接続| USBの互換性 [USB3.0 ▼] に変更
    が良かったのか、もしくは

    sudo apt-get install subversion scons libfuse-dev gcc
    cd ~
    svn co http://exfat.googlecode.com/svn/trunk/ exfat-read-only
    cd exfat-read-only
    scons
    sudo scons install
    cd ..
    rm -rf exfat-read-only

    が良かったのか…
    多分前者だけでいい気がしますが、よく分かりません。
    どちらか、もしくは両方行えば認識してくれます。
    Mount ExFAT filesystems under Linux Mint / Ubuntu

  • ffmpeg
  • ※アップデートは割と早いです。以下のコマンドがいつまで使えるかは分かりません(以下のは2014/12/28確認したものです)。
    ここのサイト(Compile FFmpeg on Ubuntu, Debian, or Mint)を見て、更新されていればそちらを使ってください。

    sudo apt-get update
    sudo apt-get -y install autoconf automake build-essential libass-dev libfreetype6-dev libgpac-dev \
      libsdl1.2-dev libtheora-dev libtool libva-dev libvdpau-dev libvorbis-dev libx11-dev \
      libxext-dev libxfixes-dev pkg-config texi2html zlib1g-dev
    mkdir ~/ffmpeg_sources
    sudo apt-get install yasm
    sudo apt-get install libx264-dev
    sudo apt-get install unzip
    cd ~/ffmpeg_sources
    wget -O fdk-aac.zip https://github.com/mstorsjo/fdk-aac/zipball/master
    unzip fdk-aac.zip
    cd mstorsjo-fdk-aac*
    autoreconf -fiv
    ./configure --prefix="$HOME/ffmpeg_build" --disable-shared
    make
    make install
    make distclean
    sudo apt-get install libmp3lame-dev
    sudo apt-get install libopus-dev
    cd ~/ffmpeg_sources
    wget http://webm.googlecode.com/files/libvpx-v1.3.0.tar.bz2
    tar xjvf libvpx-v1.3.0.tar.bz2
    cd libvpx-v1.3.0
    PATH="$HOME/bin:$PATH" ./configure --prefix="$HOME/ffmpeg_build" --disable-examples
    PATH="$HOME/bin:$PATH" make
    make install
    make clean
    cd ~/ffmpeg_sources
    wget http://ffmpeg.org/releases/ffmpeg-snapshot.tar.bz2
    tar xjvf ffmpeg-snapshot.tar.bz2
    cd ffmpeg
    PATH="$HOME/bin:$PATH" PKG_CONFIG_PATH="$HOME/ffmpeg_build/lib/pkgconfig" ./configure \
      --prefix="$HOME/ffmpeg_build" \
      --extra-cflags="-I$HOME/ffmpeg_build/include" \
      --extra-ldflags="-L$HOME/ffmpeg_build/lib" \
      --bindir="$HOME/bin" \
      --enable-gpl \
      --enable-libass \
      --enable-libfdk-aac \
      --enable-libfreetype \
      --enable-libmp3lame \
      --enable-libopus \
      --enable-libtheora \
      --enable-libvorbis \
      --enable-libvpx \
      --enable-libx264 \
      --enable-nonfree \
      --enable-x11grab
    PATH="$HOME/bin:$PATH" make
    make install
    make distclean
    hash -r
  • FireFoxにFlashを導入する
  • Flash プラグインをインストールする

    1. Adobe.com の Flash Player ダウンロードページ を開く
    2. install_flash_player_(バージョン)_linux.(プロセッサ).tar.gz ファイルを保存する
    3.   ( 2014/12/28時点では64bitコンピュータの場合、install_flash_player_11_linux.x86_64.tar.gz という名前でした。)

    4. ファイルをダウンロードしたディレクトリへ移動します (例: cd /home/user/Downloads)
    5. ダウンロードしたファイルをtarコマンドで解凍
    6. libflashplayer.so ファイルが出てくる。

      tar -zxvf install_flash_player_(バー ジョン)_linux.(プロセッサ).tar.gz
    7. スーパーユーザ 権限で、抽出した libflashplayer.so ファイルを Firefox をインストールしたディレクトリ内の plugins サブディレクトリにコピーしてください。
    8. 例えば、Firefox が /usr/lib/mozilla にインストールされている場合は、

      sudo cp libflashplayer.so /usr/lib/mozilla/plugins

      コマンドを実行し、スーパーユーザのパスワードを入力します。(LinuxMint17.1,Rebecca,cinnamonで特に変更していない場合、デフォルトでは上の場所に置けばいいです。)

    終わったら、firefoxを再起動してください。

カラーマップの上に線を描く

gnuplotでカラーマップの上に2次元で書かれたグラフを書きたいとします。

この場合は

set pm3d map
splot "fort.11" u 1:2:3 with pm3d, "fort.10" u 1:2:($2-$2) with point

とすればokです。
カラーマップの情報はfort.10, 書きたい線のデータはfort.11に書かれているとします。
考えは、2次元のデータをあたかも3次元のデータとして扱うことで解決します。

coloronplot
上のデータを得るためのfortranコードはこんな感じです。

program main
  implicit none
  double precision::x,y,h
 
  h=5.d-2
  y=-3.d0
  do while(y.le.3.d0)
     x=-3.d0
     do while(x.le.3.d0)
        write(11,*)x,y,x+y
        x=x+h
     enddo
     write(11,*)
     
     write(10,*)y,sin(y)
     y=y+h
  enddo

  return
end program main

818 re(1):pm3dによる等高線図(カラーマップ)に2次元グラフを重ねる方法

組み合わせ

組み合わせを求めるfortranコードです。
順番を気にしないで、1~nの番号が振られたボールの中からr個取りだす組み合わせは具体的に何か?のプログラムを作ります。
2つのやり方を考えます。

  1. 再帰サブルーチン”recursion”を用いる方法
  2. 漸化式を用いる方法

再帰サブルーチンを用いる方法は、組み合わせを求める手法として一般的な手法です。ここでは公開されているサブルーチンをもとに紹介します。
漸化式を用いる方法は、僕がある問題を解く上で大きな数の具体的な組み合わせ(具体的には\({}_{nCr}C_{10}\)とか)が必要になったわけです。再帰サブルーチンの方法は確かにいいのですが、一度に配列のメモリを確保して一度に書き出すことしかできないため(僕の知識がないだけかと思うけど…。)、メモリ不足になり太刀打ちできなくなりました。そこで漸化式が記述できれば、ということで考えたものです。

再帰サブルーチン”recursion”を用いる方法


再帰サブルーチンというものを使います。コードのソースはCombinations -Rosettacode.orgです。
ただし、上のリンクでは内部サブルーチンを用いて、画面出力するだけだったのでこれらを変更しました。

直したものが↓になります。

subroutine cc(n,r,nCr,cgroup)
  !sikinote
  implicit none
  integer::n,m,r,nCr,cgroup(1:nCr,1:r),comb(1:r)
  integer::k,nx,rx
 
  cgroup=0
  nx=n
  rx=r
  m=1
  k=1
  call combination(m,k,nCr,nx,rx,cgroup,comb)
 
  return
end subroutine cc

recursive subroutine combination(m,k,nCm,n_max,m_max,cgroup,comb)
  !sikinote
  implicit none
  integer,intent(in)::m
  integer::k,n,nCm,n_max,m_max,cgroup(1:nCm,1:m_max),comb(1:m_max)
  logical::check
 
  if (m > m_max) then
     cgroup(k,1:m_max)=comb(1:m_max)
     k=k+1
  else
     do n = 1, n_max
        check=.false.
        if(m.eq.1)then
           check=.true.
        else
           if((m == 1) .or. (n > comb (m - 1)))check=.true.
        endif
        if (check) then
           comb (m) = n
           call combination(m+1,k,nCm,n_max,m_max,cgroup,comb)
        end if
     end do
  end if
end subroutine combination

再帰サブルーチンは自分自身を呼び出すもので、if文を繰り返し使いたい時によく使われます。例えば、ある変数iに応じて
i=1の時

if(条件式)then
endif

i=2の時

if(条件式)then
    if(条件式)then
    endif
endif

i=3の時

if(条件式)then
    if(条件式)then
        if(条件式)then
        endif
    endif
endif

を実装する際に使えます。詳しくは調べていないのでこれ以上説明するのは避けます。
とにかく上で、call cc(n,r,nCr,cgroup)で使えます。

nCrを求めるプログラム(関数)は

function nCr(n,r)
  !sikinote
  implicit none
  integer::n,r,i,r0,nCr

  r0=n-r
  if(r.le.r0)r0=r

  nCr=1
  do i=1,r0
     nCr=nCr*(n-r0+i)
     nCr=nCr/i
  enddo
 
  return
end function nCr

でokです。\({}_{29}C_{14}\)まで計算可能です。詳しくは大きい数を扱うで書いたので参考にしてください。

メインのプログラムを書くと、

program main
  !sikinote
  implicit none
  integer::i,j,n,r,nCr,C
  integer,allocatable::cgroup(:,:)
  external::nCr
 
  n=6
  r=4
  C=nCr(n,r)

  allocate(cgroup(1:C,1:r))
  cgroup(1:C,1:r)=0
  call cc(n,r,C,cgroup)

  write(6,'(A)')"+------------+"
  do i=1,C
     do j=1,r
        write(6,'(i3)',advance='no')cgroup(i,j)
     enddo
     write(6,*)
  enddo
  write(6,'(A)')"+------------+"

  stop
end program main

と書いて実行すると、以下の出力が得られます。

+------------+
  1  2  3  4
  1  2  3  5
  1  2  3  6
  1  2  4  5
  1  2  4  6
  1  2  5  6
  1  3  4  5
  1  3  4  6
  1  3  5  6
  1  4  5  6
  2  3  4  5
  2  3  4  6
  2  3  5  6
  2  4  5  6
  3  4  5  6
+------------+

漸化式を用いる方法


+------------+
  1  2  3  4
  1  2  3  5
  1  2  3  6
  1  2  4  5
  1  2  4  6
  1  2  5  6
  1  3  4  5
  1  3  4  6
  1  3  5  6
  1  4  5  6
  2  3  4  5
  2  3  4  6
  2  3  5  6
  2  4  5  6
  3  4  5  6
+------------+

をrecursionを使わないで得るにはどうしよう…と考えます。
何度か考えてやったことは全ての行を初めの値で引く、です。実際引いてみると

+------------+
  0  0  0  0
  0  0  0  1
  0  0  0  2
  0  0  1  1
  0  0  1  2
  0  0  2  2
  0  1  1  1
  0  1  1  2
  0  1  2  2
  0  2  2  2
  1  1  1  1
  1  1  1  2
  1  1  2  2
  1  2  2  2
  2  2  2  2
+------------+

となります。なんとなく3進数っぽく見えたのでこれの漸化式さえ分かれば後で復元して…でいけそうです。
繰り上がりの時に一番大きな数を超えない、ということに注意していけば、

subroutine combination(n,r,ini,arr)
  !Update array "arr"
  !Example{n=6(maximun number)
  !        r=4(column of "ini" and "arr" array)
  !        ini(1)~ini(4)=1,2,3,4 (initial number of combination)
  !        arr(1)~arr(4)  (array what you want to update)  }
  ! - If you input   arr(1)~arr(4) = 1,2,3,4
  !      will update arr(1)~arr(4) = 1,2,3,5
  !
  ! - If you input   arr(1)~arr(4) = 1,3,5,6
  !      will update arr(1)~arr(4) = 1,4,5,6
  !
  ! - If you input   arr(1)~arr(4) = 1,4,5,6
  !      will update arr(1)~arr(4) = 2,3,4,5
  !
  !sikinote
  implicit none  
  integer::i,n,r,ini(1:r),bef(1:r),arr(1:r)
  integer::numx
  logical::key(1:r)
 
  !Memory array to bef(1~r) before update in "(n-r+1)ary" picture.
  bef(1:r)=arr(1:r)-ini(1:r)
  !Initialize "arr"
  arr(1:r)=0
  !Memory moves up or not to logical array "key"
  key(1:r)=.false.

  !numx means moves up number (n-r).
  numx=n-r
  do i=1,r
     if(bef(i).eq.numx)key(i)=.true.
  enddo
   
  do i=1,r-1
     if(key(i+1))then
        if(key(i))then
           if(i.ne.1)arr(i)=arr(i-1)      
        else
           arr(i)=bef(i)+1
        endif
     else
        arr(i)=bef(i)
     endif
  enddo
  if(key(r))then
     arr(r)=arr(r-1)
  else
     arr(r)=bef(r)+1
  endif

  !Restore from "(n-r+1)ary" picture to "10ary" picture
  arr(1:r)=arr(1:r)+ini(1:r)

  return
end subroutine combination

というサブルーチンでいけそうです。”key”は繰り上がりのある桁の場所をマークするためのものです。
先ほど紹介したnCrを求めるプログラムnCrと共にメインのプログラムを書くと、

program main
  !sikinote
  implicit none
  integer::i,j,n,r,C
  integer,allocatable::ini(:),cgroup(:)
  integer::nCr
  external::nCr
 
  n=6
  r=4
  C=nCr(n,r)
  allocate(ini(1:r),cgroup(1:r))
  ini(1:r)=0
  cgroup(1:r)=0
 
  do j=1,r
     ini(j)=j
     cgroup(j)=ini(j)
  enddo
 
  write(6,'(A,i0,A,i0)')"Combination @ n -> ",n,", r -> ",r
  do i=1,C
     do j=1,r
        write(6,'(i3)',advance='no')cgroup(j)
     enddo
     write(6,*)
     call combination(n,r,ini,cgroup)
  enddo
 
  stop
end program main

となります。出力は望んでいた通り、

Combination @ n -> 6, r -> 4
  1  2  3  4
  1  2  3  5
  1  2  3  6
  1  2  4  5
  1  2  4  6
  1  2  5  6
  1  3  4  5
  1  3  4  6
  1  3  5  6
  1  4  5  6
  2  3  4  5
  2  3  4  6
  2  3  5  6
  2  4  5  6
  3  4  5  6

めでたしめでたし。

時間計測

fortran90において、
このプログラムは何秒かかるか?
早い?遅い?
を調べるためにはもちろん計算時間を計るのが一番です。

  1. 実時間とCPU時間の違い
  2. 実時間計測
  3. CPU時間計測

実時間とCPU時間の違い


『計算時間』には大きく2種類あります。実時間CPU時間です。

  • 実時間 → 現実の世界での経過時間
  • CPU時間 → プログラムの実行でCPUを使った時間

※[1]より。参考文献とは程遠いですが、僕の認識とあっていて、綺麗な説明と感じたので載せておきます。

例えば、計算時間はそんなにかかってないのに詳細な画像を得たいがために書き出すデータ量を増やしている場合、実時間は長く、CPU時間は短くなります。
また、あるデータ点1つを得るために計算時間は膨大にかかる場合、実時間≈CPU時間になります。

fortranでは時刻を得るためのサブルーチンが既に用意されています。
ここで紹介するのは、実時間を得るための
date_and_time

system_clock
の2種類。
CPU時間を得るためのサブルーチンは
cpu_time
の1種類です。

実時間計測(date_and_time)


その場所で何時何分何秒を出力させたい時、その場所に

  write(6,'(3A,i0,A,i0,A,i0,A,i0,A,i0,A,i0,A)')"  == ", &
       c," ",ti(1),"/",ti(2),"/",ti(3), &
       "  ",ti(5),":",ti(6),":",ti(7),", (yyyy/mm/dd  hh:mm:ss)"

を入れましょう[2]。character(10)::b(1:3), integer::ti(1:8)です。

使う際はこうすると良いでしょう。

program main
  implicit none

  call current_time("program start ")  

  !program here
  write(6,*)"press Enter"
  read *
   
  call current_time("program finish")  
 
  stop
end program main

subroutine current_time(c)
  implicit none
  character(*),intent(in)::c
  integer::ti(1:8)
  character(10)::b(1:3)
 
  call date_and_time(b(1),b(2),b(3),ti)
  write(6,*)
  write(6,'(3A,i0,A,i0,A,i0,A,i0,A,i0,A,i0,A)')"  == ", &
       c," ",ti(1),"/",ti(2),"/",ti(3), &
       "  ",ti(5),":",ti(6),":",ti(7),", (yyyy/mm/dd  hh:mm:ss)"
  write(6,*)
 
  return
end subroutine current_time

サブルーチンcurrent_timeを定義して、それを呼び出すだけのほうが分かりやすいかと思います。
引数は何という文言をその場所に指定したいかを示すものです。
実行して適当なときにEnterキーを押してプログラムを進めると、

$ gfortran main2.f90
$ ./a.out
   
  == program start  2016/2/15  9:16:38, (yyyy/mm/dd  hh:mm:ss)

 press Enter

  == program finish 2016/2/15  9:16:40, (yyyy/mm/dd  hh:mm:ss)

$

という結果が得られるでしょう。もちろん、プログラムを動かした時点での時間が表示されます。
一瞬で終わるプログラムの場合は差は見られ無いと思います。

実時間計測(system_clock)


system_clockはfortran標準実装なので汎用性は高いです[4]。system_clockは組み込まれているサブルーチンなので下のプログラムをそのままコピペして使う事が出来ます。

program main
  implicit none
  integer::ti,tf,tr ! ti: initial time, tf: finish time, tr: time rate
  integer::i

  call system_clock(ti)

  !program here
  write(6,*)"press Enter"
  read *

  call system_clock(tf,tr)
 
  write(6,'(f10.3,A)')(tf-ti)/dble(tr),"[s]"
 
  return
end program main

上の文を実行して適当なときにEnterキーを押してプログラムを進めると、

$ gfortran main2.f90
$ ./a.out
    press Enter
      1.283[s]
$

という出力が得られます。

CPU時間の計測(cpu_time)


CPU時間の計測するには

    call cpu_time(t0)
    ...
    call cpu_time(t1)
    write(6,'(f10.3,A)')(t1-t0),"[CPU sec]"

を入れましょう[3]。real::t1,t2です。

program main
  implicit none
  real::t0,t1
  integer::i

  call cpu_time(t0)

  !program here
  write(6,*)"press Enter"
  read *
 
  call cpu_time(t1)
 
  write(6,'(f10.3,A)')(t1-t0),"[CPU sec]"
 
  stop
end program main

でokです。適当な時間待って、Enterキーを押すと

$ gfortran main2.f90
$ ./a.out
    press Enter
     0.000[s]

という出力が得られるでしょう。
これはcpuを動かしていないためであり、cpuが動いていた時間というのは限りなく0だ、ということです。

参考文献

[1]ユーザ時間とシステム時間の違いを教えてください。
[2]1.4.7.1 date_and_time: 日付と時刻の取得 -ORACLE®
[3]8 移植性のある時間計測の方法 -nag
[4]9.254 SYSTEM_CLOCK — Time function

大きな数を扱う

大きい数=工夫or対数!

大きい数を扱う…例えば階乗を計算しろ、という問題や組み合わせ\(_nC_r\)を計算しろ、という問題に時々あたります。

まず、典型的な例として組み合わせ
\(\displaystyle _nC_r=\frac{n!}{(n-r)!r!}\)
を計算しましょう。
この計算はそのまま定義のまま計算しようとすると階乗が出てくるので、途中の値は物凄くでかい数字になるのですが、最後の答えは小さい値になるはずです。それなのに計算できないのはおかしくない!?工夫次第で何とかなるものです。
下の表は3つの計算方法についてまとめたものです。

integer(4バイト整数型±2 147 483 647)を用いた計算
定義のまま計算 工夫して計算 対数を利用して計算
最大のn,r \(_{12}C_6\) \(_{29}C_{14}\) \(_{33}C_{16}\)
理由 階乗の計算によるオーバーフロー “割り切れる”ための制約 (限界値)

integer(8バイト整数型±9 223 372 036 854 775 807)を用いた計算
計算方法 定義のまま計算 工夫して計算 対数を利用して計算
最大のn,r \(_{20}C_{10}\) \(_{61}C_{30}\) \(_{48}C_{24}\)
理由 階乗の計算によるオーバーフロー “割り切れる”ための制約 倍精度型変数の有効桁数が足りない
定義のまま計算

例えばfortranでintegerの宣言により定義した場合、階乗が正しく計算できる範囲(オーバーフローしない範囲)は、\(12!=479\ 001\ 600\)
まで計算できます。
変数の型がいくつまで計算できるかは、組み込み関数hugeを使うことで確認できます。

program main
  implicit none
  integer::n

  write(6,*)"huge",huge(n)

  stop
end program main

これで確認すると、整数型は
\(\pm 2\ 147\ 483\ 647\)
までの値なら代入することができるため、
\(13!=6\ 227\ 020\ 800\)は計算できないことがわかります。
一応関数はこんな感じ。

function nCr_fact(n,r)
  !sikinote
  implicit none
  integer::nCr_fact,fact,n,r
  external::fact
 
  nCr_fact=fact(n)/(fact(n-r)*fact(r))
 
  return
end function nCr_fact
function fact(n)
  implicit none
  integer::fact,i,n
 
  if(n.le.-1)write(6,*)"####warning#### parameter of factorial has negative value"
  fact=1
  do i=2,n
     fact=fact*i
  enddo

  return
end function fact

( ※余談ですが、8バイト型整数で宣言、すなわちinteger(8)で宣言すると
\(9\ 223\ 372\ 036\ 854\ 775\ 807\)
まで計算できます(19桁、\(20!=2\ 432\ 902\ 008\ 176\ 640\ 000\)まで)。)

なので定義通りの計算はn=12までが限界です。実用的じゃありません。

工夫して計算

\(\displaystyle _nC_r=\frac{n!}{(n-r)!r!}\)を工夫します。

\(
\begin{align}
\displaystyle _nC_r&=\frac{n!}{(n-r)!r!} \\
&=\frac{n(n-1)(n-2)\cdots(n-r+2)(n-r+1)}{r(r-1)(r-2)\cdots2\cdot 1} \\
&=\frac{n}{r}\cdot \frac{n-1}{r-1}\cdot \frac{n-2}{r-2}\cdots \frac{n-r+2}{2}\cdot \frac{n-r+1}{1} \\
&=\prod_{i=1}^r\frac{n-r+i}{i}
\end{align}
\)
として計算します。こうすれば階乗みたいに大きな数を計算する必要がなくなります。
関数を書けばこうなります。

function nCr(n,r)
  !sikinote
  implicit none
  integer::n,r,i,r0,nCr

  r0=n-r
  if(r.le.r0)r0=r

  nCr=1
  do i=1,r0
     nCr=nCr*(n-r0+i)
     nCr=nCr/i
  enddo
 
  return
end function nCr

で\(_nC_r\)を求めることができます。
上のプログラムでは、最大\(_{29}C_{14}=77\ 558\ 760\)まで求めることができます。
これは絶対に割り切れることを前提にしているためにおこるもので、

nCr=nCr*(n-r0+i)

の文のところで最大値が決まってしまうためです。
定式化すれば、
\(_nC_r\)に対し、\(_{n-1}C_{r-1}\cdot n\)が扱える整数値を超えないこと
と表現できます。

\(_{29}C_{14}=77\ 558\ 760\)が計算できるのは、
\(_{29-1}C_{14-1}*29=1\ 085\ 822\ 640\)
であるので計算可能です。次の\(_{30}C_{15}\)を計算するためには
\(_{30-1}C_{15-1}*30=2\ 326\ 762\ 800\)が
integerの最大の整数値\(\pm 2\ 147\ 483\ 647\)
を超えてしまうのでこれ以上計算はできません。
※もしもinteger(8)を用いると\(_{61}C_{30}\)まで計算できます。

対数を利用して計算

もう一つ対数を利用した方法を紹介します。
対数を利用すると、
\(\displaystyle _nC_r=\frac{n!}{(n-r)!r!}\)は、\(S= _nC_r\)とおくと
\(
\begin{align}
\displaystyle \ln(S)&=\ln\left[\frac{n!}{(n-r)!r!}\right] \\
&=\ln(n!)-\ln\{(n-r)!\}-\ln(r!) \\
&=\sum_{i=n-r+1}^n\ln(i)-\sum_{i=2}^{r}\ln(i)
\end{align}
\)
と書けます。よって、
\(\displaystyle S=\exp\left[\sum_{i=n-r+1}^n\ln(i)-\sum_{i=2}^{r}\ln(i)\right]\)
です。プログラムを書けば

function nCr_log(n,r)
  !sikinote
  implicit none
  integer::n,r,i,r0,nCr2
  double precision::tmp
 
  r0=n-r
  if(r.le.r0)r0=r

  tmp=0.d0
  do i=n-r0+1,n
     tmp=tmp+log(dble(i))
  enddo
  do i=2,r0
     tmp=tmp-log(dble(i))
  enddo

  nCr2=nint(exp(tmp))  

  return
end function nCr_log

となります。この場合は\(_{n}C_{r}\)が整数型で扱える範囲を越えなければいいという条件だけになるので、
integerでは\(_{33}C_{16}\)までokです。

integer(8)では\(_{66}C_{33}\)までできます…が、上のプログラムそのままではできません。
おそらく、組み込み関数nint(実数型を整数型に四捨五入して代入する)がinteger(8)に対応していない(かも)ということと、倍精度の有効桁数がだんだんと足らなくなっていくので厳しいです。後者による問題によって、確認した範囲では\(_{48}C_{24}\)位が限界です。