#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <time.h>
#include <vector>
#include <algorithm>

typedef std::vector<uint32_t> uint32_v;


static double timediff(const timespec *tp0, const timespec *tp1)
{
    return double (tp1->tv_sec - tp0->tv_sec) +
        1e-9 * (tp1->tv_nsec - tp0->tv_nsec);
}

static inline uint32_t mul(uint32_t a, uint32_t b, uint32_t p)
{
    return uint64_t(a) * b % p;
}

static inline uint32_t pow(uint32_t a, uint32_t n, uint32_t p)
{
    uint32_t res = 1, b = a;
    while (n) {
        if (n & 1)
            res = mul(res, b, p);
        b = mul(b, b, p);
        n >>= 1;
    }
    return res;
}

struct Factor
{
    uint32_t l;
    uint32_t ps[9];
    uint32_t ks[9];
};


void factor(uint32_t a, Factor &f)
{
    uint32_t l = 0;

    assert (a);
    if (a % 2 == 0) {
        uint32_t k = 0;
        while (a % 2 == 0) {
            a /= 2;
            k++;
        }
        f.ps[l] = 2;
        f.ks[l] = k;
        l++;
    }

    for (uint32_t p = 3; ; p += 2) {
        uint32_t q = a / p, r = a % p;
        if (q < p)
            break;
        if (r == 0) {
            uint32_t k = 0;
            while (a % p == 0) {
                a /= p;
                k++;
            }
            f.ps[l] = p;
            f.ks[l] = k;
            l++;
        }
    }

    if (a > 1) {
        f.ps[l] = a;
        f.ks[l] = 1;
        l++;
    }

    f.l = l;
}

uint32_t find_deg(uint32_t a, uint32_t p, Factor fq)
{
    uint32_t q = p - 1;
    assert (pow(a, q, p) == 1);
    for (uint32_t l = 0; l < fq.l; l++) {
        uint32_t pl = fq.ps[l], kl = fq.ks[l];
        while (kl) {
            uint32_t qq = q / pl;
            if (pow(a, q / pl, p) != 1)
                break;
            q /= pl;
            kl--;
        }
    }
    //printf("a=%u, p=%u, deg=%u\n", a, p, q);
    assert (q);
    return q;
}

void solve(uint32_t a, uint32_t p, uint32_t l, uint32_t r, uint32_v &v)
{
    assert (a < p);
    assert (l <= r);
    uint32_t q = p - 1;
    Factor fp, fq;
    factor(p, fp);
    assert (fp.l == 1 && fp.ks[0] == 1);
    factor(q, fq);
    
    v.clear();
    if (a <= 1) {
        if (l <= a && a <= r)
            v.push_back(a);
        return;
    }
    
    uint32_t m = find_deg(a, p, fq);
    for (uint32_t x = l; x <= r; x++) {
        if (pow(x, m, p) == 1)
            v.push_back(x);
    }
    printf("size=%u\n", (uint32_t) v.size());
    std::sort(v.begin(), v.end());
}

void test(uint32_t a, uint32_t p, uint32_t l, uint32_t r)
{
    uint32_v v;
    v.reserve(100);
    
    printf("a=%u, p=%u, l=%u, r=%u\n", a, p, l, r);
    timespec tp0, tp1;
    clock_gettime(CLOCK_MONOTONIC, &tp0);
    solve(a, p, l, r, v);
    clock_gettime(CLOCK_MONOTONIC, &tp1);
    printf("%.9fs\n", timediff(&tp0, &tp1));    
}


int main()
{
    test(29, 999999937, 204411352, 204412600); // 124999993 / 1248
    test(7, 999999937, 450141522, 450142959);  // 111111105 / 1437
    test(32, 999999751, 239615046, 239616639); //  99999976 / 1593
    test(41, 999999937, 395909863, 395911989); //  76923073 / 2126
    test(6, 999999751, 217939683, 217942142);  //  66666651 / 2459
    test(6, 999999883, 724460982, 724463816);  //  55555550 / 2834
    test(35, 999999937, 399778083, 399780936); //  55555553 / 2853
    test(2, 999999937, 432617324, 432621146);  //  41666665 / 3822
    //test(13, 815730721, 0, 815730721);  //  10 / 815730721
    return 0;
}