#include<stdio.h>
#include<stdlib.h>
#include<stdint.h>
#include<gmp.h>
#include<math.h>

int main(int ac, char **av){
  uint64_t v, w;
  int q, k, j, b, i;
  double sum = 0, eps;
  mpz_t *pow, s, tmp, tmp2, tmp3, p;

  if(ac != 3){
    printf("Usage: E_q_k q k\n");
    exit(0);
  }
  q = atoi(av[1]);
  k = atoi(av[2]);
  pow = (mpz_t*)malloc(2*k*sizeof(mpz_t));
  for(i = 0; i < 2*k; i++) mpz_init(pow[i]);
  mpz_set_ui(*pow, 1);
  mpz_init_set_ui(s, 1); mpz_init(tmp); mpz_init(tmp2); mpz_init(p);
  for(i = 1; i < 2*k; i++) mpz_mul_ui(pow[i], pow[i-1], q);
  for(i = 1; i < k; i++) mpz_add(s, s, pow[i]);
  for(b = 2; b <= k; b++) if((b&3) != 3 && q%b == 1 && (b+1)%q) for(v = 2; v < q; v++){
    w = v;
    for(i = 0; i < b; i++) w = (w*v)%q;
    if(w != 1) continue;
    mpz_set_ui(p, v);
    for(j = 1; j <= k; j++){
      mpz_pow_ui(tmp3, p, b);
      mpz_mul(tmp, tmp3, p);
      mpz_sub_ui(tmp, tmp, 1);
      mpz_divexact(tmp, tmp, pow[j]);
      mpz_neg(tmp, tmp);
      mpz_mod(tmp, tmp, pow[1]);
      if(!mpz_sgn(tmp)) continue;
      mpz_powm_ui(tmp2, p, b, pow[1]);
      mpz_mul_ui(tmp2, tmp2, b+1);
      mpz_invert(tmp2, tmp2, pow[1]);
      mpz_mul(tmp, tmp, tmp2);
      mpz_mod(tmp, tmp, pow[1]);
      mpz_mul(tmp, tmp, pow[j]);
      if((!(b&1) || mpz_congruent_ui_p(p, 1, 4)) && mpz_cmp(tmp3, pow[2*j]) < 0 && mpz_probab_prime_p(p, 10)){
        eps = 2*j-b*log(mpz_get_d(p))/log(q);
        sum += eps;
        printf("p=%s j=%d b=%d eps=%f sum=%f\n", mpz_get_str(NULL, 10, p), j, b, eps, sum);
      }
      mpz_add(p, tmp, p);
    }
  }
  printf("sum=%f\n", sum);
}
/*

--------------        BCR results (in order to check)      --------------

serveur6:~/code/misc$ ./E_q_k 7 172
p=2 j=1 b=2 eps=1.287586
p=10744682090246617 j=19 b=2 eps=0.060748
p=12927056705861584875687630399167656943998305343900319630427909162956962804028786608263511417448696587 j=119 b=2 eps=1.077194
p=3376853 j=8 b=2 eps=0.549693
p=936579478224094047977 j=25 b=2 eps=0.368966
p=677822686658425215407071060694728215733634249273881 j=61 b=2 eps=1.703623
p=4020563766827570048511998229655093775332634764300364251622805949471752457774955475739650336310509313359098137765181615640612833 j=150 b=2 eps=0.379646
sum=5.427456

NB: factor p=2 should not be considered...

serveur6:~/code/misc$ ./E_q_k 3221 42
p=11 j=1 b=4 eps=0.812548
sum=0.812548

serveur6:~/code/misc$ ./E_q_k 612067 22
p=63601 j=1 b=2 eps=0.339855
p=52915283185303600217046163660913353511950087333134106429889256916343584191707324489815681786937851 j=17 b=2 eps=0.225341
sum=0.565195

serveur6:~/code/misc$ ./E_q_k 3169 36
p=97 j=1 b=2 eps=0.865001
p=5875516237 j=3 b=2 eps=0.419159
p=266602399893630549086712579594816998201 j=11 b=2 eps=0.048148
sum=1.332307

--------------        William's roadblocks       --------------

serveur6:~/code/misc$ ./E_q_k 2801 82
p=7 j=1 b=4 eps=1.019412
sum=1.019412

serveur6:~/code/misc$ ./E_q_k 547 106
p=17883403903191919 j=6 b=2 eps=0.128177 /// 3*229*1468723*11832589
sum=0.128177

serveur6:~/code/misc$ ./E_q_k 13 268
p=3 j=1 b=2 eps=1.143365
p=14600137274456049185171 j=20 b=2 eps=0.205594 /// 3121*3593793042070502653
p=6779605828531585362945329435261372873457736192243 j=44 b=2 eps=0.327311 /// 7,..D
p=38837129956841413196435465514153003370303860787375480621584612861879637753841042549799899 j=80 b=2 eps=0.944790 /// 19,..D
p=22384066606615859143319242320402727 j=31 b=2 eps=0.327315 /// 3*13597*2739724132861*131632421278843711
sum=2.948375

serveur6:~/code/misc$ ./E_q_k 191 102
sum=0.000000

*/