Summation mod$ (2^{16} +1)$ How can we calculate $\displaystyle\sum\limits_{k=1}^{2^{16}} \binom{2k}{k}(3\times 2^{14} +1)^k (k-1)^{2^{16}-1}$ mod $(2^{16} +1)$? I am aware that $2^{16} +1$ is a prime. 
 A: It is congruent to 28673 modulo $2^{16}+1$.
The main problem with evaluating this equation is efficiently calculating binomial coefficients modulo p.  Fortunately Lucas' Theorem allows us to do this reasonably efficiently.  I implemented this calculation in GAP:
BinomialModP:=function(m,n,p)
  local M,N;
  M:=[];
  N:=[];
  while(m>0 or n>0) do
    M:=Concatenation([m mod p],M);
    N:=Concatenation([n mod p],N);
    m:=Int(m/p);
    n:=Int(n/p);
  od;
  return Product([1..Size(M)],i->Binomial(M[i],N[i]) mod p) mod p;
end;;

h:=0;
p:=2^16+1;
for k in [1..p-1] do
  b:=BinomialModP(2*k,k,p);
  c:=PowerModInt(49153,k,p);
  d:=PowerModInt(k-1,65535,p);
  inc:=(b*c*d) mod p;
  h:=(h+inc) mod p;
  # Print("k=",k," inc=",inc," (mod ",p,") h=",h," (mod ",p,")\n");
od;
Print(h,"\n");

which was quite slow (but got there in the end).  I also implemented it in C using the GMP "bignum" library (which was much faster):
#include <stdio.h>
#include <gmp.h>
#include <limits.h>

const unsigned long long int p=65537; // 2^16+1
mpz_t p_GMP;
const unsigned long long int e=65535; // 2^16-1
const unsigned long long int z=49153; // 3*2^14+1

const unsigned long long int second_term[16]={1, 49153, 61441, 64513, 65281, 65473, 65521, 65533, 65536, 16384, 4096, 1024, 256, 64, 16, 4};

// An implementation of Lucas' Theorem; see: http://en.wikipedia.org/wiki/Lucas'_theorem
unsigned long long int binomial_mod_p_lucas(unsigned long long int m, unsigned long long int n) {
  mpz_t bin,inc;
  mpz_init(bin);
  mpz_set_ui(bin,1);
  mpz_init(inc);
  mpz_set_ui(inc,1);

  while(m>0 || n>0) {
    mpz_bin_uiui(inc,m%p,n%p);
    mpz_mod_ui(inc,inc,p);
    mpz_mul(bin,bin,inc);
    mpz_mod_ui(bin,bin,p);
    m=m/p;
    n=n/p;
  }
  unsigned long long int result=mpz_get_ui(bin);
  mpz_clear(bin);
  mpz_clear(inc);
  return result;
}

unsigned long long int k_minus_one_to_the_e(unsigned long long int k) {
  mpz_t pow;
  mpz_init(pow);
  mpz_set_ui(pow,k-1);
  mpz_powm_ui(pow,pow,e,p_GMP);
  unsigned long long int result=mpz_get_ui(pow);
  mpz_clear(pow);
  return result;
}

unsigned long long int main() {
  mpz_init(p_GMP);
  mpz_set_ui(p_GMP,p);

  printf("max unsigned long long int: %llu  max possible u-int: %llu\n",ULLONG_MAX,(p-1)*(p-1));

  unsigned long long int h=0,inc;
  int k;

  for(k=1;k<=32768;k++) {
    inc=binomial_mod_p_lucas(2*k,k);
    inc=(inc*second_term[k%16])%p;
    inc=(inc*k_minus_one_to_the_e(k))%p;
    h=(h+inc)%p;
    // printf("k=%llu inc=%llu (mod %llu) h=%llu (mod %llu)\n",k,inc,p,h,p);
  }

  printf("found h=%llu (mod %llu).\n",h,p);
  mpz_clear(p_GMP);
  return 0;
}

Both implementations returned $\equiv 28673 \pmod {65537}$.  Two additional simplifications were used in the C code:


*

*$z:=3 \times 2^{14}+1$ has multiplicative order $16$ modulo $p$ (so I store the 16 values of $z^i \pmod p$ in memory and recall them from there when needed), and

*${2k \choose k} \equiv 0 \pmod p$ if $2k \geq p$ and $k<p$ (so we only need to sum over $k=1,\ldots,2^{15}$).


(This code is valid only on a 64-bit machine in -std=c99 mode.)
