Numerical Error

Question

I wrote a quick sieve to test if a number is prime or not. I have two questions:

1) I tested a 200 digit prime, and it incorrectly says it wasn't prime. I believe this is down to floating point error (or something like that). How can I make this more accurate?

2) Is there a better way of writing this? I used decimal to deal with larger numbers. Is this the best way?

import math
from decimal import *
def isprime(n):
    i = 2
    a = 1
    if n == 1:
        return 0
    if n == 2 or n == 3:
        return 1
    while i < n**0.5 + 1:
        if Decimal(math.fmod(n,i)) == 0:
            a = 0
            i = n
        if Decimal(math.fmod(n,i)) != 0:
            i += 1
            a = 1
    return a

This is a duplicate question that has already been answered. [Python and primes](http://stackoverflow.com/questions/18833759/python-prime-number-checker) — Dresden, Apr 06 '16 at 00:30
python have integer of arbitrary precision, you don't need to use floating point number to make trial division, but trial division is to slow to big numbers, better check other more specialized algorithm like [Miller-Rabin test](https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants_of_the_test) — Copperfield, Apr 06 '16 at 00:48
This is **not** a duplicate of [Python and primes](http://stackoverflow.com/questions/18833759/python-prime-number-checker); there is virtually no overlap since the two main thrusts here are issues of numeric precision and primality testing for middlish to large numbers whereas [Python and primes](http://stackoverflow.com/questions/18833759/python-prime-number-checker) - and oodles of similar topics - deal only with piddling small numbers that fit in a machine word. — DarthGizka, Apr 06 '16 at 05:18
First: what you did is trial division, not a sieve; many people fail to understand the difference. Second: Floating-point inaccuracy is definitely the cause of your problem, as your algorithm is correctly stated. Third: A 200-digit number is far too large to be tested for primality in this way. Fourth: You will need a Miller-Rabin or Baillie-Wagstaff primality tester; let me know if you want code. — user448810, Apr 06 '16 at 18:43
Python 3 has infinite precision integer math, you should use it instead of `Decimal`. — Antti Haapala -- Слава Україні, Apr 07 '16 at 06:37

score 2 · Answer 1 · answered Apr 06 '16 at 16:58

The standard double floating point format can only represent integers exactly up to 2^53 (9007199254740992, 16 digits); from that point on there are gaps between representable whole numbers, and these gaps get bigger as the numbers get bigger.

64-bit builds of python use 64-bit integers natively, on other platforms you could use numpy's int64. That doesn't get you anywhere near 200 digits but it gets you far enough away from the 32-bit integer ranges for naive code to get painfully slow.

For example, trial division needs to consider only up to pi(sqrt(2^32)) = 6542 potential small prime divisors when working with integers up to 2^32, or sqrt(2^32)/2 = 32767 odd candidate divisors plus the number 2, or somewhere between these two extremes when using a wheel with an order higher than 2. Near 2^64, the number of small prime divisors to be tested is pi(sqrt(2^64)) = 203,280,221 already...

Deterministic primality testing beyond 2^32 is the domain of algorithms like Miller-Rabin or Baillie-PSW. Miller-Rabin is deterministic up to certain thresholds - up to somewhere around 2^64 - when used with certain carefully selected sets of bases; see The best known SPRP bases sets. Baillie-PSW is also known to be deterministic up to at least 2^64.

This means that beyond 2^64 you need to use a big integer type of some sort, and you have to make do with probabilistic algorithms for primality testing (giving you 'industrial-grade' primes rather than proven ones). Or plan on spending insanely huge amounts of time on actually proving the primality of a 200-digit number... A 200-digit number has a 100-digit square root (some 300 bits), so even computing all potential small prime divisors for fueling a naive trial division test is not feasible anymore.

Numerical Error

1 Answers1