divrep/coulfact.c at master · hvds/divrep · 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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#include <stdlib.h>
#include "coulfact.h"
#include "gmp_main.h"   /* prime_count */

/* for sorting */
int _mpz_comparator(const void *va, const void *vb) {
    return mpz_cmp(*(mpz_t *)va, *(mpz_t *)vb);
}

void init_fact(t_fact *f) {
    f->count = 0;
    f->size = 16;
    f->ppow = malloc(f->size * sizeof(t_ppow));
}
void free_fact(t_fact *f) {
    free(f->ppow);
}
void add_fact(t_fact *f, t_ppow pp) {
    uint count = f->count++;
    if (f->count > f->size) {
        uint size = f->size * 2;
        f->ppow = realloc(f->ppow, size * sizeof(t_ppow));
        f->size = size;
    }
    f->ppow[count] = pp;
}
void reverse_fact(t_fact *f) {
    t_ppow pp;
    uint c = f->count;
    for (uint i = 0; i + i + 1 < c; ++i) {
        uint j = c - i - 1;
        pp = f->ppow[i];
        f->ppow[i] = f->ppow[j];
        f->ppow[j] = pp;
    }
}

void init_zfact(t_zfact *f) {
    f->count = 0;
    f->size = 16;
    f->ppow = malloc(f->size * sizeof(t_zpow));
    for (int i = 0; i < f->size; ++i)
        mpz_init(f->ppow[i].p);
}
void free_zfact(t_zfact *f) {
    for (int i = 0; i < f->size; ++i)
        mpz_clear(f->ppow[i].p);
    free(f->ppow);
}
void add_zfact(t_zfact *f, t_zpow pp) {
    uint count = f->count++;
    if (f->count > f->size) {
        uint size = f->size * 2;
        f->ppow = realloc(f->ppow, size * sizeof(t_zpow));
        for (int i = f->count; i < size; ++i)
            mpz_init(f->ppow[i].p);
        f->size = size;
    }
    mpz_set(f->ppow[count].p, pp.p);
    f->ppow[count].e = pp.e;
}

uint try_simple_fact(uint n, uint d, t_fact *f) {
    uint e = 0;
    while ((n % d) == 0) {
        n /= d;
        ++e;
    }
    if (e) {
        t_ppow pp;
        pp.p = d;
        pp.e = e;
        add_fact(f, pp);
    }
    return n;
}

void simple_fact(uint n, t_fact *f) {
    uint d = 3;
    if (n > 1)
        n = try_simple_fact(n, 2, f);
    while (n > 1) {
        n = try_simple_fact(n, d, f);
        d += 2;
    }
    return;
}

uint simple_tau(t_fact *f) {
    uint t = 1;
    for (uint i = 0; i < f->count; ++i)
        t *= f->ppow[i].e + 1;
    return t;
}

uint simple_valuation(ulong n, ulong p) {
    uint v = 0;
    while ((n % p) == 0) {
        ++v;
        n /= p;
    }
    return v;
}

uint simple_prime_count(ulong n) {
    mpz_t zn, zc;
    mpz_init_set_ui(zn, n);
    mpz_init(zc);
    prime_count(zc, zn);
    uint c = mpz_get_ui(zc);
    mpz_clear(zn);
    mpz_clear(zc);
    return c;
}

uint tiny_gcd(uint a, uint b) {
    if (a > b)
        return tiny_gcd(b, a);
    if (a == 0)
        return b;
    return tiny_gcd(b % a, a);
}

ulong simple_gcd(ulong a, ulong b) {
    if (a > b)
        return simple_gcd(b, a);
    if (a == 0)
        return b;
    return simple_gcd(b % a, a);
}