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.

-

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.)

-