8
\$\begingroup\$

I wrote a program in Python 3 to check all tuples of the form (s,t) where s and t are positive integers satisfying three conditions:

  1. The inequalities s <= t^2 and t <= s^2 hold, and the value s+t divides st(s+1)(t+1);
  2. The value of s+1 is prime;
  3. The inequalities both s * [math.ceil( math.ceil( (t^2)/(s+1) ) * ((s+1)/t) )] > s(s+t) and s > t hold.

I'm interested in how often a tuple (s,t) satisfies each of these conditions in turn; that is, how often a tuple meets condition (1), versus meeting both (1) and (2), versus meeting (1) - (3).

import math

def IsPrime(n):
# This returns True if a number is prime, and False if the number is composite.
    if n % 2 == 0 and n > 2: 
        return False
    return all(n % i for i in range(3, int(math.sqrt(n)) + 1, 2)) 

def IsGQ(s,t):
# This checks the divisibility condition.
    return (s*t*(s+1)*(t+1)) % (s+t) == 0 and s <= t**2 and t <= s**2

def IsNotTransitive(s,t):
    n = math.ceil((t**2)/(s+1))
    k = math.ceil(n*((s+1)/t))
    return (s*k) > t*(t+s) and s > t

rng = 1000 # The upper limit that `t` will iterate up to
quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) and (2) 
prime_quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) - (3) 
intransitive_quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) - (4)

# The next 5 lines place all prime numbers up to rng^2 into a list. 
# Ideally, this cuts down on how many times s+1 has to be checked to be a prime.
primes = [] 
for i in range(1, rng**2):
    if IsPrime(i):
        primes.append(i)

for t in range(4, rng + 1): # Due to project details, I don't care about t<4. 
    if t % 50 == 0: 
        print("We have reached t = " + str(t) + ".")
    for s in range(2, t**2 + 1): # To satisfy condition (1), I don't want s to iterate above t^2
        if IsGQ(s,t): # Condition (1) and (2)?
            quads += 1
            if s+1 in primes:            # Condition (3)?
                prime_quads += 1
                if IsNotTransitive(s,t): # Condition (4)?
                    intransitive_quads += 1

print(str(quads))
print(str(prime_quads))
print(str(intransitive_quads))  

Currently, this runs to completion in over 10 minutes for rng = 1000. Ideally, running this for rng = 10000 in a reasonable amount of time is my current goal; computing these numbers for the value rng = 100000 is my most optimistic goal.

Is there any way to reorder my iteration or change my method of checking for primeness that could be more efficient? It seems that the time for s+1 in primes to complete grows quadratically with rng. Since, after a certain point, s+1 will be less than every value in primes that hasn't been checked yet, would it be better to write my own "in" function along the lines of

def CheckPrime(s, primes):
    for n in primes:
        if s+1 == n:
            return True
        elif s+1 < n:
            return False

or would this still be slower than the current code?

Note: IsPrime comes from this Stack Overflow answer.

\$\endgroup\$
  • \$\begingroup\$ you said it takes 10 minutes to run the program . but your program consists of two parts. calculating the primes list and testing the numbers. How long takes each part? \$\endgroup\$ – miracle173 Apr 11 '18 at 5:39
2
\$\begingroup\$

You are correct that the largest time sink in this function is your IsPrime function. Your CheckPrime function should in theory be faster, however, there is a better way.

First, a relatively small optimization: Instead of using the IsPrime from the linked SO question, use the Sieve of Eratosthenes to generate a list of of True / False values depending on the index. Here's a simple implementation:

def sieve(limit):
  primes = [True] * limit
  primes[0] = primes[1] = False

  for i in range(limit):
    if primes[i]:
      for n in range(i ** 2, limit, i):
        primes[n] = False

  return primes

If you cache the result of this function, checking if a number is prime becomes as simple as cache[n], O(1) instead of roughly O(n).

This change certainly helps, with rng = 1000, on my system it took 4 minutes, 3 seconds to run. *

The next best change I see is reducing the range that you iterate over s. If condition 1 does not pass, there's no need to even touch that value for s. With this in mind, I changed s in range(2, t**2 + 1) to s in range(int(t ** 0.5)) and ran the code again. While this will certainly help with larger ranges, it didn't help too much with rng = 1000, resulting in a time of (again) 4 minutes and 3 seconds. I suspect that this is because at the upper bounds of the limit, we only get to skip ~100 iterations. It will make a difference if you run this for rng = 100000.

Next, I remembered reading that the cost of a function call in Python is fairly high [source] - 150ns per call. That might not seem like much, but when you are calling a function roughly 300 hundred million times (IsGQ function, rng = 1000), it can be significant. With this in mind, I inlined the IsGQ function, and dropped the range checks since they are now handled by the range function, and moved the s > t check for condition 3 (4?) into the if statement. This helped a bit, though less than I expected, resulting in a run time of 3 minutes and 54 seconds.

There's probably more optimization that can be done here, but I am no expert in Python, and nothing else is jumping out to me.

Here's the code I used for the final run:

import math

rng = 1000

def sieve(limit):
  primes = [True] * limit
  primes[0] = primes[1] = False

  for i in range(limit):
    if primes[i]:
      for n in range(i ** 2, limit, i):
        primes[n] = False

  return primes

print('Populating prime cache')
primeCache = sieve(rng ** 2 + 2)
print('Prime cache populated')

def IsNotTransitive(s,t):
  n = math.ceil((t**2)/(s+1))
  k = math.ceil(n*((s+1)/t))
  return (s*k) > t*(t+s)

quads = 0
prime_quads = 0
intransitive_quads = 0

for t in range(4, rng + 1):
  if t % 50 == 0: 
    print("We have reached t = " + str(t))
  for s in range(int(t ** 0.5), t**2 + 1):
    if s * t * (s + 1) * (t + 1) % (s + t) == 0:
      quads += 1
      if primeCache[s + 1]:
        prime_quads += 1
        if s > t and IsNotTransitive(s,t):
          intransitive_quads += 1

print(str(quads))
print(str(prime_quads))
print(str(intransitive_quads))  

* All times given are CPU times, not real time. Real time was generally only a few seconds higher. 4m7s for the first run.

\$\endgroup\$
2
\$\begingroup\$

You are running into several performance problems.

First: if you want to generate a list of primes, you should use a sieve which is much more efficient.

Second: you are checking s+1 in primes, which is inefficient for a list. Instead you should make primes be a set, so that your lookup time is constant. Third: Instead of for s in range(2, t**2 + 1):, you should use for s in range(floor(t**.5), t**2 + 1):, as this prevents wasting time checking values where s^2<t.

Between these two changes you should see a pretty large speed boost.

\$\endgroup\$
  • \$\begingroup\$ I also noticed the range check - I was surprised that changing it actually had no noticeable effect! I suspect this is because 1000 ** 0.5 is only 100, and the inner loop is run so many times that 100 is nothing in the grand scheme of things. \$\endgroup\$ – Gerrit0 Apr 10 '18 at 21:42
  • \$\begingroup\$ Makes sense, it still is cleaner, and lets you remove the check from IsGq. \$\endgroup\$ – Oscar Smith Apr 10 '18 at 21:54
  • \$\begingroup\$ @Gerrit0 1000**0.5=31622... So the effect is even less \$\endgroup\$ – miracle173 Apr 11 '18 at 5:57
  • \$\begingroup\$ if one times the Code one sees that None of your "improvement" has an Impact on the overall performance of the code \$\endgroup\$ – miracle173 Apr 11 '18 at 8:29
  • 1
    \$\begingroup\$ @OscarSmith The difference is irrelevant, I show some measurements in my post \$\endgroup\$ – miracle173 Apr 16 '18 at 17:06
2
\$\begingroup\$

If you want to improve Performance you should know where you spend your time. You can either use a profiler or sin this case it is sufficient if you insert some statements that retrieve the actual time, like timer.timer() The difference of the function return values between this function is the elapsed time in seconds.

percentages
  t   time   %pr   %lp      cnt %quad     quad
 50   0.05    20    80    42864   4.1     1766
100   0.32   4.1    96   338239   1.7     5769
150   1.26   2.7    97  1136114  0.94    10675
200   2.84   2.4    98  2686489  0.63    16994
250   6.02   2.3    98  5239364  0.46    24122
300  10.43   2.1    98  9044739  0.36    32359

So for t=300 the elapsed time was 10.43 seconds, 2.1% of this time where used to calculate the primes, 98% of the time your program spent in the Loop. cntis the number, how often the function IsGQ was called and quadsshows how often IsGQ was true and quad was incremented. %quad Shows you, how many of These cnt numbers where actually pass the IsGQ test.

the %pr/%lp numbers Show you, that optimizing the first part of your algorithn, were the prime numbers are found, does not Change anything. Even if the creation of the primt takes 0 seconds, your algorithm will not get faster.

I compared your prime creation algorithm to the sieving algorithm proposed by @OscarSmith or @Gerrit0 and got the following results

prime creation
  t      t*t    org    imp
 50     2500   0.01   0.00
100    10000   0.01   0.00
150    22500   0.03   0.01
200    40000   0.07   0.01
250    62500   0.14   0.02
300    90000   0.22   0.02

So your algorithm is rather bad compared to the sieve algorithm, but nevertheless ist impact on the program execution time is negligible.

I compared your algorithm (org) to the improved method to check primes proposed by @Gerrit0 and a method, where the primes are stored as Python set, and to a Loop where condition 2 and 3 are not tested at all, and finally to an empty loop where no condition is tested at all.

loop time 
  t      cnt    org    imp    set  nochk  empty
 50    42864   0.04   0.03   0.03   0.03  0.013
100   338239   0.31   0.22   0.21   0.20  0.028
150  1136114   1.23   0.83   0.72   0.71  0.098
200  2686489   2.77   1.62   1.63   1.66   0.22
250  5239364   5.88   3.50   3.57   3.41   0.44
300  9044739  10.21   6.27   6.08   6.10   0.75

You see that your method to find a prime in your least takes almost double time of @Gerrit0's method. Using set does not improve the Speed and there is not other method that can improve your program Performance because even if you skip the prime test your program will not run faster as if you use @Gerrit0's or my method. I think there is one Problem with @Gerrit0's method that it needs considerable more Memory than my set. For large t that can make a difference. If you 'set' to 'empty', you will loose most of your time in the IsGQ which most of the time Returns Falseas %quad shows.

Based on this data one can try to model the time needed for the algorithm and I got

$$ \log(\text{time})=3,107 \log(t) -6,691$$

'log' is the logarithm of base 10. I used Excel to calculate a linear regression model. From this I estimated the running time of the loop part of your program as

loop  
     t     est. time
    50          0,04
   100          0,33
   150          1,18
   200          2,87
   250          5,75
   300         10,13
   400         24,76
   500         49,52
  1000        426,71
  2000       3676,76
 10000     546057,82
100000  698785591,19

Even if you can half the running time by using a better method to check if a number is a prime you will need more than 3 days to run your program for t=1000 and about 11 years to run it for t=100000.

The running time of your program is proportional to t^3. This seems to be true because how often is IsGQ called? It is approximatly

$$\sum_{s=1}^ts^2=\frac{(2t+1)(t+1)t}{6}\approx\frac{t^3}{3}$$

(see A000330)

So it seems you cant make your algorithm considerable faster with the medhods proposed in the posts of @OscarSmith or @Gerrit0.

The run time of the empty loop compared to the runtime of the nochk look Shows that you can avoid IsGQ-testing numbers.


Maybe the following idea can improve performance: in IsGQ you check if s+t divides st(s+1)(t+1). If you divide st(s+1)(t+1) by s+t the remainder is

$$s^2(s+1)(s-1)$$ or, because of the symmetry of s and t

$$t^2(t+1)(t-1)$$

Instead of iterating thru all pairs (t,s) (there are O(t^3)) you iterate over all possible modules m=t+s (there are only O(t^2) and calculate all possible s (there shouldn't be too many, it depends on the primefactorization of m). Maybe you get an algorithm with a running time proportional to t^2 instead of t^3. But this is a rather complex Task and there are still a lot of numbers where IsGQ results in Falsebecause t is to large.


Finally the implementation of IsNotTransitive is strange. Why do you use float calculation. Now you have to consider if floating point errors have an impact on your logic. ceil can be implemented in integer arithmetic:

ceil(a/b) = -((-a)//b)

This seems to be even a little bit faster.

\$\endgroup\$
1
\$\begingroup\$

First things first: there's an out-by-one bug in your code, which made me suspect an out-by-one bug in mine until I tracked it down:

for i in range(1, rng**2):
    if IsPrime(i):
        primes.append(i)

for t in range(4, rng + 1):
    for s in range(2, t**2 + 1): # To satisfy condition (1), I don't want s to iterate above t^2
        ...
            if s+1 in primes:            # Condition (3)?

When \$s = t^2\$, \$s + t = t(t+1)\$ so the first condition always succeeds. But when \$t\$ is rng that means we're testing rng**2+1 in primes, and that can give false negatives.


Is there any way to reorder my iteration or change my method of checking for primeness that could be more efficient?

Wrong question. It's not a matter of re-ordering the iteration: it's a matter of iterating over something much smaller.

You're looking for \$s\$ such that \$(s+t) \mid st(s+1)(t+1)\$. But \$s(s+1) = (s+t)(s+1-t) + (t-1)t\$, so that's equivalent to looking for \$s\$ such that \$(s+t) \mid (t-1)t^2(t+1)\$. It's pretty easy to factor \$(t-1)t^2(t+1)\$, since it's already partly factored. To make it really efficient you can use a sieve (which you already need for the primality test anyway).

I adapted some code I wrote previously to enumerate factors given a prime factorisation, and was able to run with rng=1000 in about 1.4 seconds. Online demo.

The bottleneck (about 55% for rng=10000) was the sieve generation, and that would be a serious memory problem for rng=100000. Since we use it for factorisation up to t+1, and only for primality testing beyond that, I changed the sieve to only run up to rng+2 and IsPrime to case split. This gave a 70% speedup (presumably beating the obvious limit of 55% by improving the cache locality), down to just over 30 seconds for rng=10000. Online demo.

Then there are small speedups available by pushing the range tests on s into the generator and not pushing onto the heap values of s+t which are too large. I did think that there would be a speedup by observing that when you increment t you can reuse the factorisations of t and t+1, but that actually doesn't seem to work. My fastest version takes 20 seconds to go up to rng=10000:

import math
import heapq

rng = 10000 # The upper limit that `t` will iterate up to
quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) and (2) 
prime_quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) - (3) 
intransitive_quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) - (4)

# Pre-calculate sieve
rr = rng + 2
max_prime = [0] * rr
for p in range(2, rr):
    if max_prime[p] == 0:
        max_prime[p] = p
        for q in range(p * p, rr, p):
            max_prime[q] = p

small_primes = [p for p in range(2, rr) if max_prime[p] == p]

def IsPrime(n):
    if n < rng:
        return max_prime[n] == n
    for p in small_primes:
        if n % p == 0:
            return False
        if p * p > n:
            break
    return True

def CandidateS(t):
    # Factor (t-1)t^2(t+1)
    primeFactors = dict()
    for q in [t-1, t, t, t+1]:
        r = q
        while r > 1:
            p = max_prime[r]
            primeFactors[p] = primeFactors.get(p, 0) + 1
            r //= p

    primes = [(1, 1)] + sorted(primeFactors.items())
    q = [(1, 0, 1)]
    lower = t + int(math.ceil(math.sqrt(t)))
    cutoff = t * (t + 1)
    while len(q) > 0:
        # d is the divisor
        # i is the index of its largest "prime" in primes
        # a is the exponent of that "prime"
        (d, i, a) = heapq.heappop(q)
        if d >= lower:
            yield d - t
        if a < primes[i][1]:
            next = d * primes[i][0]
            if next <= cutoff:
                heapq.heappush(q, (next, i, a + 1))
        if i + 1 < len(primes):
            next = d * primes[i + 1][0]
            if next <= cutoff:
                heapq.heappush(q, (next, i + 1, 1))
            # The condition i > 0 is to avoid duplicates arising because
            # d == d // primes[0][0]
            if i > 0 and a == 1:
                next = d // primes[i][0] * primes[i + 1][0]
                if next <= cutoff:
                    heapq.heappush(q, (next, i + 1, 1))

def IsNotTransitive(s,t):
    n = math.ceil((t**2)/(s+1))
    k = math.ceil(n*((s+1)/t))
    return (s*k) > t*(t+s) and s > t

for t in range(4, rng + 1): # Due to project details, I don't care about t<4. 
    if t % 500 == 0: 
        print("We have reached t = " + str(t) + ".")
    for s in CandidateS(t):
        quads += 1
        if IsPrime(s+1):            # Condition (3)?
            prime_quads += 1
            if IsNotTransitive(s,t): # Condition (4)?
                intransitive_quads += 1

print(str(quads))
print(str(prime_quads))
print(str(intransitive_quads))

Online demo

\$\endgroup\$
  • \$\begingroup\$ +1 Yes, calculating the divisors of (t-1)t^2(t+1) is the right approach. I think even rng=1000000 will be possible (maybe one has to improve the factorization method). \$\endgroup\$ – miracle173 Apr 23 '18 at 5:40
  • \$\begingroup\$ @miracle173, the next place I'd look for a major speedup if pushing that far would be in the primality test. I'd switch to using Baillie-PSW rather than trial division by primes. \$\endgroup\$ – Peter Taylor Apr 23 '18 at 20:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.