/* Written by Robert Hutchins  Latest version 7/15/2021 */

#include <stdio.h>
#include <ctype.h>
#include <string.h>

#define CONGMAX 100000

int                             congruarr[CONGMAX];
int                             sieve[CONGMAX];
int                             primRootArr[CONGMAX];
int                             quadResidueArr[CONGMAX];

int                             i;
int                             j;
int                             inx;
int                             priminx;
int                             sawEven;
int                             sawZeroPrimRoots;
int                             primRootCnt;
int                             quadResCnt;
int                             totient;

int gcd(int x, int y)
{
    while (x != y)
    {
        if (x > y)
            x = x - y;
        else
            y = y - x;
    }
    return x;
}

int isCoprime(int x, int y)
{
    return (gcd(x, y) == 1);
}

int primroots(int num)
{
    int                         preroot;
    int                         primroot;

    memset(congruarr, 0, sizeof congruarr);
    memset(primRootArr, 0, sizeof primRootArr);
    inx = 0;
    for (i = 1; i < num; i++)
    {
        if (isCoprime(i, num))
        {
            congruarr[inx] = i;
            inx++;
        }
    }
    totient = inx;
    primRootCnt = 0;
    for (i = 1; i < totient; i++)
    {
        preroot = congruarr[i];
        priminx = 0;
        while (1)
        {
            priminx++;
            if (preroot == 1)
                break;
            preroot = (preroot * congruarr[i]) % num;
        }
        if (priminx == totient)
        {
            primRootArr[primRootCnt] = congruarr[i];quadResCnt;
            primRootCnt++;
            //printf(" %d ", congruarr[i]);
        }
    }

    sawEven = 0;
    sawZeroPrimRoots = 0;
    if (primRootCnt == 0)  // Composite but no primroot
    {
        //sawEven = 0;  // 
        sawEven = 1;
        sawZeroPrimRoots = 1;
    }
    else if (primRootCnt > 0)
    {
        for (i = 0; i < primRootCnt; i++)
        {
            primroot = primRootArr[i];
            if ((primroot % 2) == 0)
            {
                sawEven = 1;
                break;
            }
        }
    }

    //if (sawEven) // print out composites with even prim roots
    if (!sawEven)  // to print out composites with only odd prim roots
    //if (sawZeroPrimRoots)  to print out only zero prim roots
#if 0
    {
        const int    MAXPRINT = 17;
        int maxprint;

        printf("%d ", num);
        printf("totient %d  ", totient);
        /*for (i = 0; i < totient; i++)
        {
            printf(" %d ", congruarr[i]);
        }
        printf(")  "); */

        //maxprint = (primRootCnt > MAXPRINT) ? MAXPRINT : primRootCnt;
        maxprint =  primRootCnt;
        printf("primRootCnt %d (", primRootCnt);
        //for (i = 0; i < primRootCnt; i++)
        for (i = 0; i < maxprint; i++)
        {
            printf(" %d ", primRootArr[i]);
        }
        printf(")  ");
        printf("\n");
    }
#endif
    // Only print number, no other info
    {
        printf("%d,\n", num);
    }
}

int quadraticResidue(int num)
{
    int                         preQuad;
    int                         quadRes;
    int                         j;
    int                         alreadySawThisQuad;

    memset(congruarr, 0, sizeof congruarr);
    memset(quadResidueArr, 0, sizeof quadResidueArr);
    inx = 0;
    for (i = 1; i < num; i++)
    {
        if (isCoprime(i, num))
        {
            congruarr[inx] = i;
            inx++;
        }
    }
    totient = inx;
    quadResCnt = 0;
    alreadySawThisQuad = 0;
    /* This loop only uses the coprimes to calculate quadratic residues */
    /*for (i = 0; i < totient; i++)
    {
        preQuad = congruarr[i];
        preQuad = ( congruarr[i] * congruarr[i]) % num; */
    /* This loop uses all numbers < n */
    for (i = 1; i < num; i++)
    {
        preQuad = i;
        preQuad = ( i * i) % num;
        if (preQuad == 0)
            continue; // Don't use zero because only non-zero allowed
        for (j = 0; j < quadResCnt; j++)
        {
            if (preQuad == quadResidueArr[j])
            {
                alreadySawThisQuad = 1;
            }
        }
        if (!alreadySawThisQuad)
        {
            quadResidueArr[quadResCnt] = preQuad;
            quadResCnt++;
        }
    }
    printf("%d ", num);
    printf("totient %d  ", totient);
    for (i = 0; i < totient; i++)
    {
        printf(" %d ", congruarr[i]);
    }
    printf(")  ");

    printf("quad residue count %d  ", quadResCnt);
    for (j = 0; j < quadResCnt; j++)
    {
        printf(" %d ", quadResidueArr[j]);
    }
    printf("\n");
}

void genSieve()
{
    int                         j;
    int                         totj;
    
    memset(sieve, 0, CONGMAX);

    for (i = 2; i < CONGMAX; i++)
    {
        if (sieve[i] == 0)
        {
            //printf("%d is prime\n", i);
            totj = i;
            j = i;
            while ((totj += j) < CONGMAX)
                sieve[totj] = 1;
        }
    }
    /*for (i = 2; i < CONGMAX; i++)
        printf("i = %d %d\n", i, sieve[i]); */
    printf("sieve is initialized\n");
}

int genPrimeComposite(int * num, int prime)
{
    int n;

    n = *num;
    if (prime)
    {
        n++;
        if (sieve[n] != 0)
        {
            while (sieve[n] != 0)
            {
                n++;
            }
        }
        *num = n;
        return *num;
    }
    else  // composite
    {
        n++;
        if (sieve[n] == 0)
        {
            while (sieve[n] == 0)
            {
                n++;
            }
        }
        *num = n;
        return *num;
    }
}

int genInc(int * num)
{
    return (*num)++;
}

int main(
  int     argc,
  char *  argv[])
{
    int                         n;

    genSieve();
    n = 2;  // first primroot is 3 if prime, 4 if composite
    while (n < CONGMAX)
    {
        /*genInc(&n);
        quadraticResidue(n); */
        genPrimeComposite(&n, 0);
        primroots(n);
    }
}