traceformula/dirichlet.h at master · rfurman/traceformula · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
// Return a list of pairs (g,e) where (Z/NZ)* is generated by the order e elements g
void gen_and_pow_of_znstar(int N, vector<int>& gens, vector<int>& exps) {
    gens.clear();
    exps.clear();
    pari_sp av = avma;
    GEN zns = znstar(stoi(N));
    int phi = itos(gel(zns,1));
    int ngen = lg(gel(zns,2));
    vector<pair<int,int> > ret;
    for(int i=1; i<ngen; i++) {
	    gens.push_back(itos(gel(gel(gel(zns,3),i),2)));
            exps.push_back(itos(gel(gel(zns,2),i)));
    }
    avma = av;
}

struct enum_znstar {
    vector<int> exps;
    vector<int> gens;
    vector<int> cur_exp;
    int cur_n;
    int N;

    enum_znstar(int _N) : N(_N) {
        gen_and_pow_of_znstar(N, gens, exps);
        cur_exp = vector<int>(exps.size());
        cur_n=1;
    }
    bool next() {
        cur_exp[0]++;
        cur_n = (cur_n * gens[0]) % N;
        for(int i=0; cur_exp[i]>=exps[i]; i++) {
            cur_exp[i]-=exps[i];
            if(i==exps.size()-1) return false;
            cur_exp[i+1]++;
            cur_n = (cur_n * gens[i+1]) % N;
        }
        return true;
    }
    int get() { return cur_n; }
    operator int() { return cur_n; }
    int operator[](int i) {
        return cur_exp[i];
    }
    friend ostream& operator<<(ostream& os, const enum_znstar& v) {
        os << "(" << v.cur_n << " = ";
        for(int i=0; i<v.cur_exp.size(); i++) {
            os << (i?"*":"") << v.gens[i] << "^" << v.cur_exp[i];
        }
        os << " mod " << v.N << ")";
        return os;
    }
};


CCMatrix matrix_of_characters_by_parity(int N, int k) {
    int phiN = phi(N);
    enum_znstar e1(N), e2(N);
    vector<vector<int> > ret(phiN, vector<int>(N));
    int chi=0;
    int nchi=0;
    do {
        do {
            int& cur = ret[chi][e2];
            for(int i=0; i<e1.gens.size(); i++) {
                cur += ((e1.cur_exp[i] * e2.cur_exp[i])%e1.exps[i])*(phiN/e1.exps[i]);
                cur %= phiN;
            }
            cur %= phiN;
        } while(e2.next());
        if(ret[chi][N-1]==(k%2?phiN/2:0)) nchi++;
        chi++;
    } while(e1.next());

    CCMatrix ret2 = CCMatrix::Zero(nchi, N);
    RR pi = mpfr::const_pi();
    for(int i=0, ii=0; i<phiN; i++) if(ret[i][N-1]==(k%2?phiN/2:0)) {
        for(int j=0, jj=0; j<N; j++)
            if(__gcd(N,j)==1)
                ret2.coeffRef(ii,jj++)=polar((RR)1, -2*pi*(1.*ret[i][j]/phiN));
            else
                ret2.coeffRef(ii,jj++)=0;
        ii++;
    }
    return ret2;
}