Computing n-th power residue symbols

[geoff]

The current 1.2.x versions of sr5sieve checks whether terms of a sequence are quadratic residues before putting the sequence into the sieve. For version 1.3.x I am trying to extend this to check whether they are cubic, quartic or quintic residues too. For this to work I need to find a fast way to test, given a sequence k*b^n+c and that r divides p-1, whether -c*k*b^d is an r-th power residue wrt p, for 0 <= d < r.

When r=2 this can be done with a single division then a lookup into a precomputed table of Legendre symbols to exclude those terms with Legendre symbol (-ckb^d/p) = -1.

For r in general it can be done by computing x_i=(1/b^i)^((p-1)/r) (mod p) for 0 <= i < r once for each p, then computing y_k=k^((p-1)/r) (mod p) once for each k that is a quadratic residue wrt p, then excluding those terms k*b^n+c with n=i (mod r) which don't satisfy x_i=y_k.

I have a working prototype that implements this for r in {3,5}, but it actually results in a small slowdown (it takes longer to check whether sequences are cubic residues than the time saved by excluding 2/3 of them from the sieve). I am now looking for a faster algorithm to calculate the r-th power symbol of a wrt p than computing a^((p-1)/r) (mod p) with powmod(), but a faster powmod() would also help a little and might be enough to get a small speedup. (If the r-th power symbol could be computed as fast as the Legendre symbol then something like 35% speedup could be achieved, but even doubling the speed of powmod would probably only get a 10% speedup at best).

It is probably not possible to construct lookup tables for the cubic and quintic residue symbols, but the references I have seen when searching for cubic residues all seem to work in a ring of Eisenstein integers, and that sounds like I will need to do a bit of study to figure out how they relate to plain integers modulo a prime that the sieve works with.

If anyone can help with ideas or references for making fast versions of these functions it would be a big help:

uint64_t powmod(uint64_t a, uint64_t n, uint64_t p)
{
  /* Return a^n (mod p). Assume 0 < a < p and n=(p-1)/r, where p is a
     large prime and r is a small integer.
     If it helps then assume that the function will be called about 100 times
     with the same values of n and p. */
}

int is_cubic_residue(uint64_t a, uint64_t p)
{
  /* Return 1 if x^3=a (mod p) has a solution, 0 otherwise.
     Assume 0 < a < p where p is a large prime, 3|(p-1), and (a/p)=1.
     If it helps then it is OK to return a false 1 with low probability, but
     never a false 0. */
}

int is_quartic_residue(uint64_t a, uint64_t p) /* As above, replace 3 with 4 */
int is_quintic_residue(uint64_t a, uint64_t p) /* As above, replace 3 with 5 */

[R. Gerbicz]

I'm also writing a sieve that will sieving k*b^n+c numbers. See your problems:

Originally Posted by geoff

For r in general it can be done by computing x_i=(1/b^i)^((p-1)/r) (mod p) for 0 <= i < r once for each p, then computing y_k=k^((p-1)/r) (mod p) once for each k that is a quadratic residue wrt p, then excluding those terms k*b^n+c with n=i (mod r) which don't satisfy x_i=y_k.

Supposing that you have converted your sequences k*b^n-1.
Suppose also that k*b^n==1 mod p, where p is prime.
There are two cases for each r|(p-1) prime value: if b^((p-1)/r)!=1 mod p then you can find n mod r, suppose that n=x*r+y, where 0<=y<r
k*b^n==1 mod p then / ^((p-1)/r) you get:
k^((p-1)/r)*b^(y*(p-1)/r)==1 mod p so
k^((p-1)/r)*(b^((p-1)/r))^y==1 mod p.
This equation is solvable in all cases, and there is only one y solution (mod r),
because b^((p-1)/r)!=1 mod p means that the order of b modulo p is divisible by r.
It means that by baby step giant step method you can find the required y value for each remaining k value by O(sqrt(remaining*r)) mulmods, where remaining is the number of remaining k value. ( not counting the computation time of the required b^((p-1)/r) and k^((p-1)/r) mod p values!). If r is small then I'm using a table to see: if there is no n value for some k that n mod r=y then I eliminate that k value.

If b^((p-1)/r)==1 mod p then the equation: k^((p-1)/r)==1 mod p, so if you compute this and this isn't true then you can eliminate that k value.

You can extend this for primepowers r=q^h, where q is prime and h>1.

Originally Posted by geoff

I have a working prototype that implements this for r in {3,5}, but it actually results in a small slowdown (it takes longer to check whether sequences are cubic residues than the time saved by excluding 2/3 of them from the sieve). I am now looking for a faster algorithm to calculate the r-th power symbol of a wrt p than computing a^((p-1)/r) (mod p) with powmod(), but a faster powmod() would also help a little and might be enough to get a small speedup. (If the r-th power symbol could be computed as fast as the Legendre symbol then something like 35% speedup could be achieved, but even doubling the speed of powmod would probably only get a 10% speedup at best).

There is a much faster way for this if you are using more prime(power) divisors of p-1. See if 3 and 5 divides p-1 then first compute x1=a^((p-1)/15) mod p then x2=x1^5 mod p and x3=x1^3 mod p, the x2 give the a^((p-1)/3) mod p and x3 give a^((p-1)/5) mod p. You can easily generalize this for more terms by *binary tree* method.

At the end you know some informations: n mod q^i mod p for each k value, you can solve these by Chinese remainder theorem: you will know n mod prod for each k value, where prod is the product of all q^i.

[geoff]

Thanks for this explaination, I still have a way to go before I can implement it but I get the idea. I think it will supercede my approach, unless there is a fundamentally faster way to decide whether k is an r-th power residue than computing k^((p-1)/r).