/*                                                                                                                           
** binomod.gp ver. 1.1 (c) 2007 by Max Alekseyev                                                                           
**                                                                                                                           
** binomod(n,k,p)	computes binomial(n,k) modulo prime p using Lucas theorem
** binomod1(n,k,p)	tests if binomial(n,k) is co-prime to prime p
** binomod2(n,k,m)	tests if binomial(n,k) is co-prime to integer m
** binoval(n,k,p)	valuation of n! with respect to prime p
*/

{  \\ binomial(n,k) mod PRIME p
   binomod(n,k,p) = local(r,v,f,np,kp);
   if(k<0,return(Mod(0,p)));
   if(n<0,return( (-1)^k*binomod(-n+k-1,k,p) ));
   if(k>n,return(Mod(0,p)));
   
   v = vector(p);
   while(k,
     np = n%p;
     kp = k%p;
     if(kp>np,return(Mod(0,p)));
     if(kp && np-kp, v[np]++; v[kp]--; v[np-kp]--);
     n\=p; k\=p;
   );
   r=Mod(1,p); 
   f=r; for(i=2,p-1,f*=i;if(v[i],r*=f^v[i]));
   r
}

{  \\ is binomial(n,k)!=0 mod PRIME p
   binomod1(n,k,p) = local(np,kp);
   if(k<0,return(0));
   if(n<0,n=-n+k-1);
   if(k>n,return(0));
   
   while(k,
     np = n%p;
     kp = k%p;
     if(kp>np,return(0));
     n\=p; k\=p;
   );
   1
}

{  \\ is binomial(n,k) COPRIME to m
   binomod2(n,k,m) = local(f);
   if(k<0,return(0));
   if(n<0,n=-n+k-1);
   if(k>n,return(0));

   f=factorint(m)[1,];
   for(i=1,#f,if(!binomod1(n,k,f[i]),return(0)));
   1
}

{  \\ valuation( binomial(n,k), p ), p is prime
   binoval(n,k,p) = local(m=[n,k,n-k],r=0);
   while(m,
     m \= p;
     r += m[1] - m[2] - m[3];
   );
   r
}