1
$\begingroup$

I'm trying to speedup my primality test where every $n \in \mathbb{Z}$ is prime if all lattice points on $x+y=n$ are visible from the origin.

lattice

It's really inefficient at the moment because it requires $\approx\frac {n}2$ $gcd$ computations.

For example, to determine primality of $n = 249999991$, it takes $\approx 12$ seconds on my computer!

To speed things up, I tried reducing the $\frac {n}2$ steps by running the C program below to see if there's a way to calculate a maximum bound on the number of steps needed to determine primality.

The max steps needed for each $10,000,000$ range follows:

Range [0,10000000] Max step 1568
Range [10000000,20000000] Max step 2231
Range [20000000,30000000] Max step 2738
Range [30000000,40000000] Max step 3161
Range [40000000,50000000] Max step 3534
Range [50000000,60000000] Max step 3870
Range [60000000,70000000] Max step 4181
Range [70000000,80000000] Max step 4470
Range [80000000,90000000] Max step 4739
Range [90000000,100000000] Max step 4986
Range [100000000,110000000] Max step 5243
Range [110000000,120000000] Max step 5474
Range [120000000,130000000] Max step 5699
Range [130000000,140000000] Max step 5915
Range [140000000,150000000] Max step 6120
Range [150000000,160000000] Max step 6323
Range [160000000,170000000] Max step 6518
Range [170000000,180000000] Max step 6705
Range [180000000,190000000] Max step 6890
Range [190000000,200000000] Max step 7053
Range [200000000,210000000] Max step 7244
Range [210000000,220000000] Max step 7415
Range [220000000,230000000] Max step 7580
Range [230000000,240000000] Max step 7736
Range [240000000,250000000] Max step 7904
Range [250000000,260000000] Max step 8055
Range [260000000,270000000] Max step 8213
Range [270000000,280000000] Max step 8364
Range [280000000,290000000] Max step 8514
Range [290000000,300000000] Max step 8658
Range [300000000,310000000] Max step 8799
Range [310000000,320000000] Max step 8940
Range [320000000,330000000] Max step 9074
Range [330000000,340000000] Max step 9219
Range [340000000,350000000] Max step 9350
Range [350000000,360000000] Max step 9486
Range [360000000,370000000] Max step 9615
Range [370000000,380000000] Max step 9744
Range [380000000,390000000] Max step 9869
Range [390000000,400000000] Max step 9998
Range [400000000,410000000] Max step 10116
Range [410000000,420000000] Max step 10241
Range [420000000,430000000] Max step 10365
Range [430000000,440000000] Max step 10481
Range [440000000,450000000] Max step 10605
Range [450000000,460000000] Max step 10716
Range [460000000,470000000] Max step 10836
Range [470000000,480000000] Max step 10946
Range [480000000,490000000] Max step 11066
Range [490000000,500000000] Max step 11174
Range [500000000,510000000] Max step 11286
Range [510000000,520000000] Max step 11393
Range [520000000,530000000] Max step 11510
Range [530000000,540000000] Max step 11613
Range [540000000,550000000] Max step 11723
Range [550000000,560000000] Max step 11831
Range [560000000,570000000] Max step 11936
Range [570000000,580000000] Max step 12041
Range [580000000,590000000] Max step 12140
Range [590000000,600000000] Max step 12240
Range [600000000,610000000] Max step 12348
Range [610000000,620000000] Max step 12444
Range [620000000,630000000] Max step 12548
Range [630000000,640000000] Max step 12630
Range [640000000,650000000] Max step 12735
Range [650000000,660000000] Max step 12839
Range [660000000,670000000] Max step 12936
Range [670000000,680000000] Max step 13026
Range [680000000,690000000] Max step 13133
Range [690000000,700000000] Max step 13224
Range [700000000,710000000] Max step 13320
Range [710000000,720000000] Max step 13410
...

Notice the range $[240000000,250000000]$ has max steps $7904$.
I modified the program to impose that $7904$ bound.

Running the primality test again for $n = 249999991$, it now takes $\approx \frac {1}{1000}$ second. (compared to $12$ seconds)

Question: Is there a way to calculate a maximum bound on the number of steps required to determine primality?

For example, I could apply the "bound method" (if such a thing exists) to say it will require $123456$ steps to determine primality of

$5311379928167670986895882065524686273295931177$
$27031923199444138200403559860852242739162502265$
$22928566888932948624650101534657933765270723940$
$9519978766587351943831270835393219031728127$

Thanks.

#include <stdio.h>
#include <stdint.h>
#include <math.h>


uint64_t gcd(uint64_t a, uint64_t b)
{
    return (b != 0) ? gcd(b, a % b) : a;
}


int isPrimeNumber(uint64_t num)
{
    uint64_t i;

    for (i = 2; i <= sqrt(num); i++)
    {
        if (num % i == 0)
        {
            // not prime
            return 0;
        }
    }

    // prime
    return 1;
}


int main(int argc, char* argv)
{
    uint64_t x, y, n;
    uint64_t i, mPrime;
    uint64_t step, maxStep = 0;
    uint64_t begin, end;
    uint64_t interval = 10000000;

    for (i = 0; i < 200; i++)
    {
        begin = i * interval;
        end = begin + interval;

        printf("Range [%llu,%llu] ", begin, end);

        for (n= begin; n < end; n++)
        {
            mPrime = n % 10;

            if (mPrime == 1 || mPrime == 3 || mPrime == 7 || mPrime == 9)
            {
                if (isPrimeNumber(n) != 1)
                {
                    // Start near line x=y
                    x = (n / 2) + 2;
                    y = n - x;

                    step = 0;

                    do {
                        step++;

                        if (gcd(x, y) != 1)
                        {
                            // Lattice point not visible
                            break;
                        }

                        x++;
                        y--;
                    } while (x < n - 1);

                    if (step > maxStep) { maxStep = step; }
                }
            }
        }

        printf("Max step %llu\n", maxStep);
    }

    return 0;
}

EDIT:

Many thanks to Ravi's answer!

He was able to determine the test just needs to step thru this interval $[n/2, n/2 + \sqrt{n}/2]$ to determine primality.

wow

I've updated the C program with his suggested interval.

// Primality Test
// Every n is prime if all lattice points on x+y=n are visible from the origin.

#include <stdio.h>
#include <stdint.h>
#include <math.h>


uint64_t gcd(uint64_t a, uint64_t b)
{
    return (b != 0) ? gcd(b, a % b) : a;
}


int isPrimeNumber(uint64_t n)
{
    if (n == 1) return 0;
    if (n == 2 || n == 3) return 1;
    if (n % 2 == 0) return 0;

    // Start near line x=y.
    uint64_t x = (n / 2) + 2;
    uint64_t y = n - x;
    uint64_t count = sqrt(n) / 2;

    for (uint64_t i = 0; i < count; ++i) {
        // Check lattice point visibility...
        if (gcd(x, y) != 1) return 0;
        x++; y--;
    }

    return 1;
}


int main(int argc, char* argv)
{
    uint64_t n = 1000000007;

    if (isPrimeNumber(n) == 1)
    {
        printf("%llu prime.", n);
    }
    else
    {
        printf("%llu not prime.", n);
    }

    return 0;
}
$\endgroup$
3
  • $\begingroup$ It's sufficient to check $x \leq \sqrt{n}$. It's also sufficient to check whether $x$ itself is a divisor of $n$, instead of checking that they have a common factor. You can't do substantially better than that with this strategy, because of integers $n = pq$ where $p$ and $q$ are primes of approximately the same size. $\endgroup$ Jan 9, 2021 at 20:01
  • 1
    $\begingroup$ Actually I'm confused: you're asking for optimizations to your primality test, but the code you give already includes a much faster primality test (the isPrimeNumber routine)--which is in fact exactly the algorithm that would result from my previous comment. So what are you looking for? $\endgroup$ Jan 9, 2021 at 20:45
  • $\begingroup$ I'm looking for a pattern that reduces the number of steps needed to determine primality. I omit all the primes using that function you noted, to see what non-primes are taking for steps. Then, whenever a prime is encountered in that range, I know the max steps needed to apply. The steps seem to increase by $\approx 200$ per $10,000,000$, so can I extrapolate that going forward to say, $1$ trillion? $\endgroup$
    – vengy
    Jan 9, 2021 at 20:55

1 Answer 1

2
$\begingroup$

Claim: if $n$ is composite, it shares a common factor with some $x$ in the interval $[n/2, n/2 + \sqrt{n}/2]$.

Proof: Let $p$ be the smallest prime divisor of $n$. Then $p \leq \sqrt{n}$, since otherwise $n$ would be prime. Any multiple of $p$ shares a common factor (namely $p$) with $n$. The closest multiple of $p$ to $n/2$ is either $n/2$ itself (if $n$ is even) or $\frac{n \pm p}{2}$ (if $n$ is odd). Since $0 < p \leq \sqrt{n}$, both cases produce an integer in $[n/2, n/2 + \sqrt{n}/2]$ that shares the factor $p$ with $n$.

So it's enough to go from $x = n/2$ up to $n/2 + \sqrt{n}/2$. With this improvement, your algorithm will still be somewhat slower than your isPrimeNumber routine--just because calculating a gcd is slower than calculating divisibility. This will be feasible for testing numbers with maybe 15-20 digits, but you would need completely different methods to test something much larger, like the 183-digit number you mention.

$\endgroup$
3
  • $\begingroup$ Nice! So basically, the test just needs to step thru this interval $[n/2, n/2 + \sqrt{n}/2]$ in the worst case to determine primality. For large $n$, how does it compare to $gcd()$ computations in $O(log$ $n)$ versus $O(\sqrt{n})$? Thanks! $\endgroup$
    – vengy
    Jan 9, 2021 at 22:07
  • $\begingroup$ This algorithm would use $O(\sqrt n \log n)$ divisions--that is, $O(\sqrt n)$ gcd computations, each of which involves $O(\log n)$ divisions. For comparison, isPrimeNumber requires $O(\sqrt n)$ divisions. (State-of-the-art deterministic primality tests typically run in $O((\log n)^c)$ time for various small constants $c$; this makes it feasible to test the primality of arbitrary numbers with many thousands of digits, and numbers of special forms with millions of digits.) $\endgroup$ Jan 10, 2021 at 1:01
  • $\begingroup$ Nice explanation. Thanks Ravi! $\endgroup$
    – vengy
    Jan 10, 2021 at 1:34

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .