Puzzler: Binomial Coefficients
Author Message
Puzzler: Binomial Coefficients

Julien Constantin writes on Tuesday, April 2:

Quote:
> In this thread on computing binomial coefficients, it was suggested to obtain
> them by means of the exponent E(N,K,P) of each prime P in the prime
> factorization of K!N.    A formula for this exponent is given in
>     Goetgheluck, P.  Computing binomial coefficients,
>     The Amer. Mathematical Monthly, April 1987, p. 360-365.
> This formula is embodied in the following APL function BINEXP, where the
> exponent is given by Z and NKK is the triple N,K,N-K.

> {del}Z{<-}P BINEXP NNK
> Z{<-}+/1{neg}1{neg}1+.{times}{floor}NNK{jot}.{divide}P*{iota}{floor}P{log}N
> {del}
> Note in particular that E(N,K,P)=1 if N-K <P <= N
>                     and         =0 if N/2 <P <= N-K

> The function BINOME will produce all the exponents in the factorization of N!K.

> {del}Z{<-}K BINOME N;T;PRE
> PRE{<-}PRIMESTO N        {lamp} any function giving all the primes to N
> T{<-}{neg}1+(PRE>N{divide}2){iota}1
> Z{<-}PRE>N-K{<-}K{min}N-K
> Z[{iota}T]{<-}(PRE[{iota}T]) BINEXP {each}{enclose}N,K,N-K
> {del}

This is the "factorial factors" algorithm that Jim Weigang discovered
in Knuth, Volume 1, page 46; an implementation of it in J was previously

improved here):

[ p=.ple 20
2 3 5 7 11 13 17 19

p pexp 11                        NB. exponents for !k
8 4 2 1 1 0 0 0
p pexp 20                        NB. exponents for !n
18 8 4 2 1 1 1 1
p pexp 20-11                     NB. exponents for !n-k
7 4 1 1 0 0 0 0
p pexp 11,20,20-11
8 4 2 1 1 0 0 0
18 8 4 2 1 1 1 1
7 4 1 1 0 0 0 0

-/ p pexp 11,20,20-11
_3 0 _1 0 0 _1 _1 _1
- -/ p pexp 11,20,20-11          NB. same as _1 1 _1 +/ .* matrix
3 0 1 0 0 1 1 1
11 bce 20
3 0 1 0 0 1 1 1

*/ p ^ 11 bce 20
167960
11 ! 20
167960

Quote:
> For example, the exponents of the primes in 1000!6000 are obtained in
> about 0.4 sec. with APL68000 on a PowerMac or with APL2 under OS/2
> on a Pentium 100.

On a P90 running WinNT and J3.02x,  1000 bce 6000  takes 0.15 seconds,
of which  ple 6000  takes 0.004 seconds.

Sun, 20 Sep 1998 03:00:00 GMT

 Page 1 of 1 [ 1 post ]

Relevant Pages