Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question

1 Answer 1

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.