Compile integer functions

Question

I need to speed up a function that uses Binomial and QuotientRemainder. I used Compile, but unfortunately MA resorts to uncompiled evaluation

inxC = Compile[{{n, _Integer}, {k, _Integer}}, 
  Module[{q, ix, ic, iv},
   q = k;
   ix = 0;
   While[q > 0, q = q - Binomial[n, ix]^2; ix++];
   ix = ix - 1;
   q = q + Binomial[n, ix]^2;
   {iv, ic} = QuotientRemainder[q - 1, Binomial[n, ix]];
   {ix, iv + 1, ic + 1}
   ]]
inxC[6, 887]
(*{4, 15, 15}*)

How such a function can be compiled? It seems the optimizer does not know Binomial and QuotientRemainder.

Answer 1

inxC = With[{Binomial = Evaluate@Round@FunctionExpand@Binomial[#, #2] &}, 
  Compile[{{n, _Integer}, {k, _Integer}},
   Module[{q, ix, ic, iv}, q = k;
    ix = 0;
    While[q > 0, q = q - Binomial[n, ix]^2; ix++];
    ix = ix - 1;
    q = q + Binomial[n, ix]^2;
    iv = Floor[(q - 1)/Binomial[n, ix]];
    ic = Mod[q - 1, Binomial[n, ix]];
    {ix, iv + 1, ic + 1}]]];
inxC[6, 887]

{4, 15, 15}

Answer 2

Very similar approach to @chyanog 's. One could use equivalent, compiled functions for QuotientRemainder and Binomial and express them in terms of Gamma, Quotient and Mod.

bin = Compile[{{a, _Integer}, {b, _Integer}}, 
               Floor[Gamma[a + 1]/(Gamma[b + 1] Gamma[a - b + 1])]]

inxC3 = Compile[{{n, _Integer}, {k, _Integer}}, 
   Module[{q, ix, ic, iv}, q = k; 
          ix = 0;
      While[q > 0, q = q - bin[n, ix]^2; ix++]; ix = ix - 1; 
            q = q + bin[n, ix]^2; {iv, ic} = {Quotient[q - 1, bin[n, ix]], Mod[q - 1, bin[n, ix]]};
            {ix, iv + 1, ic + 1}], {{bin[_, _], _Integer}}, 
  CompilationOptions -> {"InlineExternalDefinitions" -> True, "InlineCompiledFunctions" -> True}]

inxC3[6, 887]
(* {4, 15, 15} *)

Answer 3

It seems that manual expressions for Binomial and QuotientRemainder are unavoidable. I follow @chyanog and @b.gates.you.know.what but exploit a recursive relation of binomial coefficients

$$\binom{n}{k} = \frac{n-k+1}{k} \binom{n}{k-1}, \quad \binom{n}{0}=1$$

for constant $n$ to avoid the manipulation of real numbers.

inxC = Compile[{{n, _Integer}, {k, _Integer}}, 
  Module[{q, ix, ic, iv, b},
   q = k;
   ix = 0;
   b = 1;
   While[q > 0, b = If[ix == 0, 1, Quotient[b*(n - ix + 1), ix] ]; 
    q = q - b^2; ix++];
   ix = ix - 1;
   q = q + b^2;
   iv = Quotient[q - 1, b];
   ic = Mod[q - 1, b];
   {ix, iv + 1, ic + 1}]]

Now some timing results:

ns = 12
{tm1, a} = Timing[Table[inxC[ns, j], {j, Binomial[2 *ns, ns]}]];
{tm2, a} = Timing[Table[inxCch[ns, j], {j, Binomial[2 *ns, ns]}]];
{tm3, a} = Timing[Table[inxCbg[ns, j], {j, Binomial[2 *ns, ns]}]];
{tm1, tm2, tm3}
(*{2.72751, 6.13963, 8.32903}*)

Here inxCch is defined in the chyanog's post and inxCbg is defined in the b.gates.you.know.what's post.