Factorial (or, even better, binomial coefficient) function

Question

Is there a function in l3fp for calculating factorials or, even better, binomial coefficients directly?

Answer 1

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

Answer 2

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}