Skip to content

Instantly share code, notes, and snippets.

@JSmol

JSmol/euler.cc Secret

Created Jan 13, 2021
Embed
What would you like to do?
Euler-Kronecker Constants
/* 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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment