13

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

17

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).

enter image description here

\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.


binomials B

Grouping three by three the digits and allowing line breaks with \threebythree macro:

three by three


First answer:

binomial coeffs

\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
  • Just to make sure that I get the latest/'best' code snippers above, can I make you create a complete example with them collected? :) – Svend Tveskæg Feb 24 '14 at 9:14
  • Somewhat embarrassing but how do I make \xintFor* #1 in {\xintSeq{0}{10}} \do {\xintFor* #2 in {\xintSeq{0}{#1}} \do {\binomialb{#1}{#2}\xintifForLast{\par}{, \hskip 0pt plus 1pt}} } into a macro with the single argument being the size of Pascal's triangle (in this case 10)? My problem is that I already have #1 and #2 but they are not the arguments I want for the macro. – Svend Tveskæg Feb 24 '14 at 16:45
  • just double the #: ##1 and ##2. And then {10} in \xintSeq {0}{10} may be replaced by {#1} if your #1 holds the 10 (possibly something like \value{mycounter} or \arabic{mycounter}), or something like {\x} where \x expands to 10. May even be 5+5. – user4686 Feb 24 '14 at 17:10
  • I mean... whoa! – Manuel Mar 15 '16 at 19:32
13

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...

enter image description here

% !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}
  • @jfbu - Could you be more specific as to what a "multiprecision" library is? Lua(La)TeX does load the entire Lua math library. – Mico Feb 23 '14 at 17:59
  • @jfbu - I've amended the example to allow the printing of arbitrary-length integers. I must confess I don't know enough about the innards of Lua's math library to judge if arbitrary (user-specified) precision is possible. – Mico Feb 23 '14 at 18:40
  • ah I see, thus, the computation was done exactly all along but the previous version printed as a float, right? I thought the computation itself had been done with floats. – user4686 Feb 23 '14 at 18:44
  • @jfbu - Indeed, the earlier (and now deleted) example relied on the default settings of the C sprintf function (which underlies the tex.sprint function) for "truncating" a long integer to scientific notation. The string "%.0f" tells string.format to print the function's second argument as a floating-point (not "scientific") number with 0 decimals; I suppose this is a slightly round-about approach to printing an arbitrarily long integer. – Mico Feb 23 '14 at 18:59
  • @jfbu - Just as in your answer, no input-sanity checking is performed in this MWE. Moreover, it is assumed (implicitly!) that the user will know to write \mchoose{1000}{100} rather than \mchoose{1000}{900}. :-) – Mico Feb 23 '14 at 19:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.