Euler-Kronecker Constants

euler.cc

/*  Josip Smolcic 2019 summer research project in Number Theory, with Amir Akbary.
 *  Computation of Euler-Kronecker constants of cyclotomic fields.
 *  Adapted from section 3.2 of https://arxiv.org/abs/1108.3805
 *
 *  Compile with:
 *  g++ -std=c++11 -lfftw3l fftw.cc -o fft-gq
 */

#include <bits/stdc++.h>

// uncomment this and use boost digamma for faster runtime (O(n log(n)) instead of O(n^2)).
// #include <boost/math/special_functions/digamma.hpp>

// fftw is used for multi-dim DFT in O(n log(n)).
// fftwl is used for long double precision.
#include <fftw3.h> 
using namespace std;

// constants and short hands

#define SQR(X) (X) * (X)
#define fst first
#define snd second
#define all(X) begin(X), end(X)

typedef long double ld;
const int maxq = 1001001, maxd = 100, terms = 1000;
const ld pi = 3.14159265358979323846264338327950288419;
const ld eg = 0.57721566490153286060651209008240243104;
const ld eps = 1e-4;

// ---------------------------------------------------------------------------
// helper functions

int gcd(int a, int b) {
  if (b == 0) return a;
  else return gcd(b, a%b);
}

// Assumes b <= a.
pair<int, int> egcd(int a, int b) {
  assert(b <= a);
  int s = 0, t = 1, u = 1, v = 0, q = 0;
  while (b != 0) {
    q = a/b; swap(a, b);
    b = b % a; swap(s, u);
    s = s - u * q; swap(t, v);
    t = t - v * q;
  }
  return {u, v};
}

int fexp(int x, int p) {
  int r = 1;
  while (p > 0) {
    if (p & 1) p -= 1, r = (r * x);
    else p >>= 1, x = (x * x);
  }
  return r;
}

int fexp(int x, int p, int mod) {
  int r = 1;
  while (p > 0) {
    if (p & 1) p -= 1, r = (r * x) % mod;
    else p >>= 1, x = (x * x) % mod;
  }
  return r;
}

// Prime power factors (p, a) where p^a || n.
set<pair<int, int>> factor(int n) {
  set<pair<int, int>> S;
  for (int i = 2, a = 0; i <= n; i++) {
    while (n % i == 0) n /= i, a++;
    if (a > 0) S.insert({i, a}), a = 0;
  }
  return S;
}

// Set of all divisors of n, except 1.
set<int> divisors(int n) {
  set<int> S = {n};
  for (int i = 2; i*i <= n; i++) {
    if (n % i == 0) {
      S.insert(i);
      S.insert(n/i);
    }
  }
  return S;
}

// # of {x such that 0 < x < n and gcd(n, x) = 1}.
int phi(int n) {
  int ans = 1;
  set<pair<int, int>> S = factor(n);
  for (auto& s : S) ans *= (fexp(s.fst, s.snd) - fexp(s.fst, s.snd-1));
  return ans;
}

// primitive root of n.
// Assuming (Z/nZ)* is cyclic, this function returns an element g which generates the group.
int primroot(int n) {
  int p = phi(n);
  set<int> D = divisors(p);
  for (int g = 2; g < n; g++) {
    if (gcd(g, n) == 1) {
      for (auto& d : D) {
        if (fexp(g, d, n) == 1 and d != p) break;
        if (fexp(g, d, n) == 1 and d == p) return g;
      }
    }
  }
}

// Chinese Remainder Theorem
// m_i are pair wise co-prime
// return solution to system x == a_i mod m_i
int crt(vector<int>& A, vector<int>& M, int n) {
  int m = accumulate(all(M), 1, multiplies<int>());
  int x = 0;
  for (int i = 0; i < n; i++) {
    int k = m/M[i];
    int d = egcd(M[i], k % M[i]).snd;
    x += A[i] * k * d % m;
    x += m;
  }
  return x % m;
}

// For small m this functions are accurate and fast.
// For large m use the boost implementation.
ld psi(ld r, ld m) {
  ld ans = eg + log(2.0*m) + pi/(2.0*tan(r*pi/m));
  for (ld n = 1; n <= (m-1)/2; n += 1.0)
    ans -= 2*cos(2*pi*n*r/m) * log(sin(pi*n/m));
  
  return -ans;
}

// terms defined above determines when we truncate the infinite sum.
// Note: sum small to large for higher precision
ld T(ld r, ld q) {
  ld n = terms, s = 0.0;
  while (n >= -eps) {
    s += (log(n + r/q)/(n + r/q) - log(n + 1.0)/(n + 1.0));
    n -= 1.0;
  }
  return s;
}

// vector of row-major order indices in multi-dimensional array.
vector<int> indices(int* n, int i, int m) {
  vector<int> I(m);
  for (int j = m - 1, k = i; j >= 0; j--) {
    I[j] = k % n[j];
    k -= I[j];
    k /= n[j];
  }
  return I;
}
  
// ---------------------------------------------------------------------------

// input precomp arrays (coefficients in formula)
ld p[maxq], t[maxq];

// fft inputs and outputs
fftwl_complex a[maxq], b[maxq], at[maxq], bt[maxq];

// dimension sizes of the multi-dimensional array.
int n[maxd];

/*  Compute the Euler-Kronecker constant of the cyclotomic field Q(\mu_q).
 *  \mu_q denotes the q'th roots of unity.
 *  The total running time using the boost::gamma function is O(terms * q * log q).
 *  where terms is the number of terms in the sum T.
 */

ld gammaq(int q) {
  if (q <= 1 or q >= maxq) return -1;

  complex<ld> gq = eg;
  set<int> D = divisors(q);
  for (auto& d : D) {

    for (int i = 0; i < q; i++) a[i][0] = a[i][1] = b[i][0] = b[i][1] = at[i][0] = at[i][1] = bt[i][0] = bt[i][1] = 0;

    // if 8 divides d, then a slightly different approach is neccesary.
    if (d % 8 != 0) {

      // Setup fft for each divisor.
      // For each unique prime factor p_i^a_i the fft will have an additional dimension.
      set<pair<int, int>> S = factor(d);
      vector<int> F, P;
      for (auto& s : S) F.push_back(fexp(s.fst, s.snd)), P.push_back(s.fst);

      // Compute the size of each dimension, and find a primitive root for each p_i^a_i.
      int m = S.size(); 
      vector<int> g(m);
      for (int i = 0; i < m; i++) {
        n[i] = phi(F[i]);
        g[i] = primroot(F[i]);
      }

      // fftw backwards is used since our exponents are positive.
      fftwl_plan pa = fftwl_plan_dft(m, n, a, at, FFTW_BACKWARD, FFTW_ESTIMATE);
      fftwl_plan pb = fftwl_plan_dft(m, n, b, bt, FFTW_BACKWARD, FFTW_ESTIMATE);

      // Prepare psi and T values.
      // Use boost for faster run time.
      for (int i = 1; i <= d-1; i++) {
        if (gcd(i, d) == 1) {
          // p[i] = boost::math::digamma((ld) i / d);
          p[i] = psi(i, d);
          t[i] = T(i, d);
        }
      }

      // Our sum has phi(d) terms total.
      // Here k is the product of the n_i, in other words, the multi-dim array has k cells.
      int k = phi(d);
      for (int i = 0; i < k; i++) {
        // Retrieve indices for cell i. (\vec k)
        vector<int> I = indices(n, i, m);

        /*  We must re-arrange the coefficients p[i] and t[i] we computed above.
         *  This is the transition from one sum over r to m sums over n_i, described in the report.
         *  Create the system of equations: c[j] = g[j] ^ I[j] mod p_j^a_j.
         *  a[I] = T(Solution to above system mod d, by CRT / d)
         *  b[I] = psi(Solution to above system mod d, by CRT / d)
         */

        vector<int> c(m);
        for (int j = 0; j < m; j++) c[j] = fexp(g[j], I[j], F[j]);
        a[i][0] = t[crt(c, F, m)];
        b[i][0] = p[crt(c, F, m)];
      }

      // Execute the dft on a and b.
      fftwl_execute(pa);
      fftwl_execute(pb);

      // Sum over appropriate entries of tranformed arrays.
      for (int i = 1; i < k; i++) {
        vector<int> I = indices(n, i, m);

        // Chi is primitive iff it is the product of primitive characters.
        bool prim = true;
        for (int j = 0; j < m; j++)
          if (I[j] % P[j] == 0)
            prim = false;

        if (prim) {
          ld logd = log(d);
          complex<ld> x = {at[i][0], at[i][1]};
          complex<ld> y = {bt[i][0], bt[i][1]};
          gq += x/y;
          gq -= logd;
        }
      }
      
      fftwl_destroy_plan(pa);
      fftwl_destroy_plan(pb);
    
    // ---------------------------------------------------------------------------
    // 8 divides d.
    } else {
      
      // The power of 2 will have 2 dimensions associated to it.
      // n[0] = 2, and n[1] = 2^(a - 2)
      set<pair<int, int>> S = factor(d);
      vector<int> F, P;
      for (auto& s : S) {
        if (s.fst == 2) n[0] = 2, n[1] = fexp(2, s.snd-2);
        F.push_back(fexp(s.fst, s.snd)), P.push_back(s.fst);
      }

      // The rest of the powers are handled as before.
      // There is no primitive root of 2^a when a > 2.
      int m = S.size() + 1;
      vector<int> g(m-2);
      for (int i = 0; i < F.size()-1; i++) {
        n[2+i] = phi(F[i+1]);
        g[i] = primroot(F[i+1]);
      }

      fftwl_plan pa = fftwl_plan_dft(m, n, a, at, FFTW_BACKWARD, FFTW_ESTIMATE);
      fftwl_plan pb = fftwl_plan_dft(m, n, b, bt, FFTW_BACKWARD, FFTW_ESTIMATE);

      for (int i = 1; i <= d-1; i++) {
        if (gcd(i, d) == 1) {
          // p[i] = boost::math::digamma((ld) i / d);
          p[i] = psi(i, d);
          t[i] = T(i, d);
        }
      }

      int k = phi(k);
      for (int i = 0; i < k; i++) {
        vector<int> I = indices(n, i, m);
        
        /*  The group of units mod 2^a (a > 3) is generated by 2 elements, we use (-1) and 5.
         *  We still have m equations, the equation corresponding to 2 is:
         *  (-1)^I[0] 5^I[1] mod 2^a.
         *  The rest are computed as before.
         */

        vector<int> c(m-1);
        c[0] = (fexp((-1 + F[0]) % F[0], I[0], F[0]) * fexp(5, I[1], F[0])) % F[0];
        for (int j = 1; j < m-1; j++) c[j] = fexp(g[j-1], I[j+1], F[j]);
        a[i][0] = t[crt(c, F, m-1)];
        b[i][0] = p[crt(c, F, m-1)];
      }

      fftwl_execute(pa);
      fftwl_execute(pb);

      for (int i = 1; i < k; i++) {
        vector<int> I = indices(n, i, m);

        // Primitivity of chi mod 2^a is checked seperately.
        bool prim = true;
        if (I[1] % 2 == 0)
          prim = false;

        for (int j = 2; j < m; j++)
          if (I[j] % P[j-1] == 0)
            prim = false;

        if (prim) {
          ld logd = log(d);
          complex<ld> x = {at[i][0], at[i][1]};
          complex<ld> y = {bt[i][0], bt[i][1]};
          gq += x/y;
          gq -= logd;
        }
      }
      
      fftwl_destroy_plan(pa);
      fftwl_destroy_plan(pb);
    }
  }
  
  return real(gq);
}

int main() {
  int q;
  cin >> q;
  cout << q << endl;
  cout << gammaq(q) << endl;
}