sqrtmod

#include "ulong_extras.h"
#include "profiler.h"

mp_limb_t n_sqrtmod_brute(mp_limb_t a, mp_limb_t p)
{
    mp_limb_t t, t2;

    if (a <= 1)
    {
        return a;
    }

/*
    if (n_jacobi_unsigned(a, p) == -1)
        return 0;
*/

    t = t2 = 1;

    while (t <= (p - 1) / 2)
    {
        if (t2 == a)
            return t;

        t2 = n_addmod(t2, 2*t + 1, p);
        t++;
    }

    return 0;
}

int main()
{
    mp_limb_t a, p, t1, t2;

    for (p = 2; p < 10000; p = n_nextprime(p, 0))
    {
        for (a = 0; a < p; a++)
        {
            t1 = n_sqrtmod(a, p);
            t2 = n_sqrtmod_brute(a, p);

            if (t1 > p / 2)
                t1 = p - t1;

            if (t1 != t2)
                abort();
        }
    }

    for (p = 2; p < 10000; p = n_nextprime(p, 0))
    {
        printf("p = %lu\n", p);

        TIMEIT_START
        for (a = 0; a < p; a++)
            t1 = n_sqrtmod(a, p);
        TIMEIT_STOP

        TIMEIT_START
        for (a = 0; a < p; a++)
            t2 = n_sqrtmod_brute(a, p);
        TIMEIT_STOP
    }
}