Factorial (or, even better, binomial coefficient) function
Since release 1.2f
, xintexpr
has a binomial
function, both for exact and for float evaluations (there are corresponding macros in xint
and xintfrac
respectively).
\documentclass[a4paper]{article}
\usepackage[hscale=0.83]{geometry}
\usepackage{xintexpr}[2016/03/12]% release 1.2f or later
\usepackage[english]{babel}
\usepackage[np, autolanguage]{numprint}
\usepackage{amsmath}
\begin{document}
\begingroup
\hrule\medskip
\edef\N{\xinttheiiexpr 0..[30]..300\relax}% 0, 30, 60, ..., 270, 300
\xintFor #1 in\N\do{$\binom{300}{#1}=\xinttheiiexpr binomial(300,#1)\relax$\par\smallskip}%
\medskip\hrule\medskip
\edef\N{\xinttheiiexpr [\N][1:-1]\relax}% cut out first and last
\xintFor #1 in\N\do{$\binom{300}{#1}\approx\np{\xintthefloatexpr binomial(300,#1)\relax}$\par\smallskip}%
\medskip\hrule\medskip
\xintFor #1 in {500, 1000, 1500, 2000, 2500}\do
{$\binom{3000}{#1}\approx\np{\xintthefloatexpr binomial(3000,#1)\relax}$\par\smallskip}
\medskip\hrule\medskip
\endgroup
\end{document}
The rest of this answer is now obsoleted.
edit 2015-11-05 because recent versions of xint
do not load xinttools
anymore.
First, an implementation of binomial(n,k) = n choose k
which uses only
\numexpr
. Will fail if the actual value is at least 2^31
(the first
too big ones are 2203961430
= binomial(34,16)
and
2333606220
= binomial(34,17)
). The 2-arguments macro \binomialb
is expandable, hence it can be used inside \num
of siunitx
or \numprint
of numprint
, or an \edef
or \write
, or in an \ifnum
test.
Its arguments are given to a \numexpr
, thus they may be infix expressions like 2+3
or a counter like \value{mycount}
, etc.
The algorithm is described in the comments of the macro \binomialb
. No package needed.
Then, the same algorithm but with the long-integer macros of xint
. The arguments to this \binomialB
are treated exactly like those of \binomialb
. Again the macro is (f)-expandable and expands fully in two steps.
\binomialB
is able to handle n = 1000
, k = 500
.
Then a \threebythree
macro to print long numbers grouping digits by three and allowing line breaks. This macro needs only xinttools
for its \xintLength
macro, xint
(which loads xinttools
) is not needed.
Finally, a macro \binomialA
which was my first answer; it is somewhat slower for big results as it does n*(n-1)*..*(n-k+1)
and then divides by the factorial of k
. For example factorial(500)
has 1135 digits whereas binomial(1000,500)
turns out to have only 300 digits. Thus the computation via \binomialA
does one gigantic division, which turns out to be about 2.2x
slower than \binomialB
with its 499 divisions with small divisors. (observed with xint 1.2
).
\documentclass {article}
\usepackage{xint, xinttools}
% for testing of line breaks + digits grouping
% \usepackage{siunitx} % no line breaks
% \usepackage{numprint} % no line breaks in math, line breaks in text possible
% Expandably computing binomial(n,k)=n choose k
% after having replaced k by the smallest of k and n-k, and checked if
% k=0, either one of the following products produces integers at each
% mutiply/divide steps:
% n * (n-1)/2 * (n-2)/3 * .... * (n-k+1)/k
% or
% (n-k+1) * (n-k+2)/2 * (n-k+3)/3 * ... * n/k
% eTeX \numexpr does multiply/divide in one "double-precision" step,
% thus arithmetic overflow should not happen, as long as the result is <
% 2^31 (and naturally the initial n, binomial (2147483648,0) will not
% work
% For no special reason I chose the second product.
% (notice that as k<n-k+1 also the first product is increasing, no
% intermediate thing can cause overflow if the final thing does not)
% Each (n-k+j)/j step could be seen as (n-k)/j + 1, thus only j would need
% incrementing; up to the price of an extra addition, and I preferred to
% carry around both an u=n-k+j and a v=j
% ALGORITHM
% replace k by the smallest of k and n-k
% if k=0 return 1
% else set w=n-k+1
% u=n-k+2
% v=2
% endif
% if v>k return w
% else
% w<-w*u/v
% u<-u+1
% v<-v+1
% repeatif
% Constraint: expandability. Adding +1 has a cost and fetching a list of
% tokens also has one. To use one macro less, or not do twice u->u+1,
% u and v are shifted from the start by 1 to be usable directly in the
% updating of w.
% no check on validity of inputs
%-----------------------------------------------------------
% expandable macro \binomialb. No package needed.
\catcode`_ 11
\def\binomialb #1#2{\romannumeral0\expandafter
\binomialb_a\the\numexpr #1\expandafter.\the\numexpr #2.}
\def\binomialb_a #1.#2.{\expandafter\binomialb_b\the\numexpr #1-#2.#2.}
\def\binomialb_b #1.#2.{\ifnum #1<#2 \expandafter\binomialb_ca
\else \expandafter\binomialb_cb
\fi {#1}{#2}}
\def\binomialb_ca #1{\ifnum#1=0 \expandafter \binomialb_one\else
\expandafter \binomialb_d\fi {#1}}
\def\binomialb_cb #1#2{\ifnum #2=0 \expandafter\binomialb_one\else
\expandafter\binomialb_d\fi {#2}{#1}}
\def\binomialb_one #1#2{ 1}
\def\binomialb_d #1#2{\expandafter\binomialb_e \the\numexpr #2+1.#1!}
% n-k+1.k! -> u=n-k+2.v=2.w=n-k+1.k!
\def\binomialb_e #1.{\expandafter\binomialb_f \the\numexpr #1+1.2.#1.}
% u.v.w.k!
\def\binomialb_f #1.#2.#3.#4!%
{\ifnum #2>#4 \binomialb_end\fi
\expandafter\binomialb_f
\the\numexpr #1+1\expandafter.%
\the\numexpr #2+1\expandafter.%
\the\numexpr #1*#3/#2.#4!}
\def\binomialb_end #1*#2/#3!{\fi\space #2}
\catcode`_ 8
%-----------------------------------------------------------
% expandable macro \binomialB. Requires package xint
\catcode`_ 11
\def\binomialB #1#2{\romannumeral0\expandafter
\binomialB_a\the\numexpr #1\expandafter.\the\numexpr #2.}
\def\binomialB_a #1.#2.{\expandafter\binomialB_b\the\numexpr #1-#2.#2.}
\def\binomialB_b #1.#2.{\ifnum #1<#2 \expandafter\binomialB_ca
\else \expandafter\binomialB_cb
\fi {#1}{#2}}
\def\binomialB_ca #1{\ifnum#1=0 \expandafter \binomialB_one\else
\expandafter \binomialB_d\fi {#1}}
\def\binomialB_cb #1#2{\ifnum #2=0 \expandafter\binomialB_one\else
\expandafter\binomialB_d\fi {#2}{#1}}
\def\binomialB_one #1#2{ 1}
\def\binomialB_d #1#2{\expandafter\binomialB_e \the\numexpr #2+1.#1!}
% n-k+1.k! -> u=n-k+2.v=2.w=n-k+1.k!
\def\binomialB_e #1.{\expandafter\binomialB_f \the\numexpr #1+1.2.#1.}
% u.v.w.k!
\def\binomialB_f #1.#2.#3.#4!%
{\ifnum #2>#4 \binomialB_end\fi
\expandafter\binomialB_f
\the\numexpr #1+1\expandafter.%
\the\numexpr #2+1\expandafter.%
\romannumeral0\xintiiquo{\xintiiMul{#1}{#3}}{#2}.#4!}
\def\binomialB_end #1\romannumeral0\xintiiquo #2#3!%
{\fi\binomialB_enda #2}
\def\binomialB_enda\xintiiMul #1#2{ #2}
\catcode`_ 8
%-----------------------------------------------------------
% The three by three macro. Requires xinttools (not automatically loaded by xint)
% \threebythree for printing 3 by 3
% argument supposed to be f-expandable
% defines the thousand separator, for example:
\DeclareRobustCommand*{\threebythreesep}{\ifmmode
\allowbreak\mskip 6mu plus 1mu minus 1mu
\else
\hskip .3333em plus 0.0555em minus 0.0555em \fi }
\catcode`_ 11
% submitting the #1 to an \edef could be better
\newcommand*{\threebythree}[1]{\expandafter\threebythree_a
\romannumeral-`0#1.}
\def\threebythree_a #1{\if#1-\expandafter\xint_firstoftwo\else
\expandafter\xint_secondoftwo\fi
{-\threebythree_b}{\threebythree_b #1}}
% #1 is assumed to be not empty
\def\threebythree_b #1.{\expandafter\threebythree_c\expandafter
{\romannumeral0\xintlength {#1}}#1...}
\def\threebythree_c #1{\ifcase \numexpr #1+3-3*((#1+2)/3)\relax
\expandafter\threebythree_da
\or\expandafter\threebythree_db
\or\expandafter\threebythree_dc
\fi }
\def\threebythree_da #1#2#3{#1#2#3\threebythree_e}
\def\threebythree_db #1{#1\threebythree_e}
\def\threebythree_dc #1#2{#1#2\threebythree_e}
\def\threebythree_e #1#2#3{\if#1.\else \threebythreesep
#1#2#3\expandafter\threebythree_e\fi}
\catcode`_ 8
%-----------------------------------------------------------
% a less efficient algorithm doing n(n-1)...(n-k+1) divided by k!:
% expandable, expands in two steps, requires xint
\catcode`_ 11
\def\binomialA #1#2{\romannumeral0%
\expandafter\binomialA_a\the\numexpr #1\expandafter.\the\numexpr #2.}
\def\binomialA_a #1.#2.{\expandafter\binomialA_b\the\numexpr #1-#2.#2.{#1}}
\def\binomialA_b #1.#2.{\ifnum #1<#2 \expandafter\binomialA_c
\else \expandafter\binomialA_d
\fi {#1}{#2}}
\def\binomialA_c #1#2#3{\xintiiquo {\xintiiPrd {\xintSeq [-1]{#3}{#2+1}}}%
{\xintiFac {#1}}}
\def\binomialA_d #1#2#3{\xintiiquo {\xintiiPrd {\xintSeq [-1]{#3}{#1+1}}}%
{\xintiFac {#2}}}
\catcode`_ 8
%-----------------------------------------------------------
\begin{document}\thispagestyle{empty}
\fdef\test {\binomialB {10}{5}}
\noindent
\texttt{\string\binomialB \{10\}\{5\} expanded twice: \meaning\test!!!
(<-space check)}
\xintFor* #1 in {\xintSeq {0}{20}}
\do
{\xintFor* #2 in {\xintSeq {0}{#1}}
\do {\binomialb {#1}{#2}\xintifForLast{\par}{, \hskip 0pt plus 1pt}}%
% comparison for checking everything ok
% \xintFor* #2 in {\xintSeq {0}{#1}}
% \do {\binomialB {#1}{#2}\xintifForLast{\par}{, \hskip 0pt plus 1pt}}%
% \xintFor* #2 in {\xintSeq {0}{#1}}
% \do {\binomialA {#1}{#2}\xintifForLast{\par}{, \hskip 0pt plus 1pt}}%
}
% \noindent\pdfresettimer
% \fdef\test {\binomialB {1000}{500}}
% (algorithm B: \the\pdfelapsedtime)\par
% text mode
${1000 \choose 500}={}$\threebythree{\binomialB {1000}{500}}
% or \threebythree {\test} etc...
% math mode
${1000 \choose 575}=\threebythree{\binomialB {1000}{575}}$
% \medskip
% \noindent\pdfresettimer
% \fdef\test {\binomialA {1000}{500}}
% (algorithm A: \the\pdfelapsedtime)\par
% % ${1000 \choose 500}={}$\threebythree{\test}
% Testing printing with package siunitx or numprint
% the \num of siunitx does not allow line breaks
% the \numprint of numprint allows line breaks in text mode, but NOT in
% math mode, if its thousand separator is suitably configured.
% \npthousandsep {\hskip .1666em plus .5pt minus .5pt}
% \numprint {\binomialB {1000}{500}}
\end{document}
remark 2015/11/05: the timing of algorithm B vs algorithm A below is obsolete; with
xint 1.2
release, the \binomialA{1000}{500} wasn't anymore 12x
but only 2.2x
slower than \binomialB{1000}{500}.
remark 2016/03/15: it seems that \binomialB has gotten about 20% faster with 1.2f's but anyhow 1.2f's binomial(1000,500) is more than 2x
faster still.
Grouping three by three the digits and allowing line breaks with \threebythree
macro:
First answer:
\input xint.sty
\input xinttools.sty
\hsize 19cm
\vsize 27.7cm
\hoffset \dimexpr -1in+1cm\relax
\voffset \dimexpr -1in+1cm\relax
\nopagenumbers
\catcode`_ 11
\def\binomial #1#2{\expandafter\binomial_a\the\numexpr #1\expandafter.%
\the\numexpr #2.}
\def\binomial_a #1.#2.{\expandafter\binomial_b\the\numexpr #1-#2.#2.{#1}}
\def\binomial_b #1.#2.{\ifnum #1<#2 \expandafter\binomial_c
\else \expandafter\binomial_d
\fi {#1}{#2}}
\def\binomial_c #1#2#3{\xintiiQuo {\xintiiPrd {\xintSeq [-1]{#3}{#2+1}}}%
{\xintiFac {#1}}}
\def\binomial_d #1#2#3{\xintiiQuo {\xintiiPrd {\xintSeq [-1]{#3}{#1+1}}}%
{\xintiFac {#2}}}
\catcode`_ 8
\xintFor* #1 in {\xintSeq {0}{20}}
\do
{\xintFor* #2 in {\xintSeq {0}{#1}}
\do {\binomial {#1}{#2}\xintifForLast{\par}{, \hskip 0pt plus 1pt}}%
}
\def\allowsplits #1{\ifx #1\relax \else #1\hskip 0pt plus 1pt\relax
\expandafter\allowsplits\fi}%
\def\printnumber #1{\expandafter\allowsplits \romannumeral-`0#1\relax }%
\xintFor* #1 in {\xintSeq {0}{50}}\do
{\printnumber{\binomial {50}{#1}}\xintifForLast{\par}{\unskip, }}%
\xintFor* #1 in {\xintSeq {0}{100}}\do
{\printnumber{\binomial {100}{#1}}\xintifForLast{\par}{\unskip, }}%
\bye
Here's a solution that uses the math capabilities of Lua(La)TeX. A TeX-side macro called \mchoose{n}{k}
and a Lua-side function named mchoose(n,k)
do most of the work. An auxiliary Lua function called fwrite
allows the printing of arbritrary-length integers; no provision, though, is made in this example to introduce line breaks when printing overly long numbers...
% !TEX TS-program = lualatex
\documentclass{article}
\usepackage[a4paper,margin=2cm]{geometry} % just for this example
\usepackage{luacode}
%%% Lua-side code
\begin{luacode}
-- 'mchoose' is patterned after the posting
-- in http://stackoverflow.com/a/15302448.
-- Thanks, @egreg, for providing the pointer to this posting!
function mchoose ( n , k ) -- 'n' should be no smaller than 'k'
if ( k == 0 or k == n ) then
return 1
else
return ( n * mchoose ( n-1, k-1 ) / k )
end
end
-- Print an arbitrary-length integer as a string
function fwrite ( z )
tex.sprint ( string.format("%.0f" , z ) )
end
\end{luacode}
%%% LaTeX-side code
\newcommand{\mchoose}[2]{\directlua{fwrite(mchoose(#1,#2))}}
\begin{document}
\verb+\mchoose{6}{3}+: \mchoose{6}{3}.
\verb+\mchoose{1000}{70}+: \tiny \mchoose{1000}{70}.
\end{document}