\\% q-factorial, q-binomial, multinomial, q-multinomial, rising/falling factorial power.
\\ Author: Joerg Arndt
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2011-August-04 (10:05)


qpoch(n,x=x,q=q)=
\\ q-Pochhammer symbol:
\\  qpoch(n,a,q) = (a;q)_n = (1-a)*(1-a*q)*...*(1-a*q^(n-1))
\\ (taken from Michael Somos' tricks.gp)
{
    if ( n>=0,
        prod(i=0,n-1,1-x*q^i)
        ,
        1/qpoch(-n,x*q^n,q)
    );
} /* ----- */


qpochv(n,v)=
{
    return(  prod(k=1, #v-1, qpoch(n, v[k], v[#v]) ) );
} /* ----- */

qpochve(n,v,e)=
{
    return(  prod(k=1, #v-1, qpoch(n, v[k], v[#v])^e[k] ) );
} /* ----- */


\\ q-factorial:
qfact(n,q=q)= { qpoch(n,q,q) / (1-q)^n; }
\\qfact(n,q=q)= { prod(k=1,n, (1-q^k) ) / (1-q)^n; }


\\ q-binomial coefficients
\\ (taken from Michael Somos' tricks.gp)
\\ (note:  uses one possible convention for m<0 or m>n)
qbinomial(n, m, q='q)=
{
    if ( q==1,  return(binomial(n,m)) );
    if ( n-m>m,  m=n-m );
    if ( (n>=0) & (n<m),  return(0) );
    if ( n<0, return( (-1)^m*q^(m*n-(m^2-m)/2) * qbinomial(m-1-n,m,q)) );
    return ( prod(i=1,n, polcyclo(i,q)^(n\i-m\i-(n-m)\i)) );
} /* ----- */

multinomial(v)=
\\ Multinomial coefficient s!/prod(j=1,n, v[j])
\\ where n = #v and s = sum(j=1,n,v[j]).
\\ Same as: number of permutations of the multiset
\\ with v[1] elements of type 1, v[2] elements of type 2, etc.
{
    local(n, s, ret);
    n = #v;  \\ Number of types of elements
    s = sum(j=1,n,v[j]); \\ total number of elements
\\    return(  s! / prod(j=1,n, v[j]!) );
    ret = 1;
    for (k=1, n-1,
        ret *= binomial(s, v[k]);
        s -= v[k];
    );
    return( ret );
} /* ----- */

qmultinomial(v, q='q)=
\\ q-multinomial coefficient.
{
    local(n, s, ret);
    n = #v;
    s = sum(j=1,n,v[j]);
    ret = 1;
    for (k=1, n-1,
        ret *= qbinomial(s, v[k], q);
        s -= v[k];
    );
    return( ret );
} /* ----- */


ffactpow(n,e)=
\\ Falling factorial power.
{
    local(f=1);
    for (k=0, e-1, f*=(n-k) );
    return( f );
} /* ----- */

rfactpow(n,e)=
\\ Rising factorial power.
{
    local(f=1);
    for (k=0, e-1, f*=(n+k) );
    return( f );
} /* ----- */


\\ ==== end of file ====