1

I'm trying to maximize the Euler Totient function on Python given it can use large arbitrary numbers. The problem is that the program gets killed after some time so it doesn't reach the desired ratio. I have thought of increasing the starting number into a larger number, but I don't think it's prudent to do so. I'm trying to get a number when divided by the totient gets higher than 10. Essentially I'm trying to find a sparsely totient number that fits this criteria.

Here's my phi function:

def phi(n):
    amount = 0

    for k in range(1, n + 1):
        if fractions.gcd(n, k) == 1:
            amount += 1

    return amount
Jonathan Leffler
  • 730,956
  • 141
  • 904
  • 1,278
Marorin
  • 167
  • 3
  • 12
  • my first thought is to use the multiplicative property of "if `gcd(m, n) = 1` then `φ(mn) == φ(m) φ(n)`" and try to do a logarithmic decomposition through regression – Aaron Feb 28 '17 at 15:29
  • Don't just use a `for` loop over all the numbers up to `n`, this can get rather slow for large `n`. Instead, 1) get all the prime numbers up to `n` (or `sqrt(n)`? not sure); 2) use those to get the prime factors of `n`, 3) use those to get the totient. More work, but should be faster. – tobias_k Feb 28 '17 at 15:42
  • Possible duplicate of [Computing Eulers Totient Function](http://stackoverflow.com/questions/18114138/computing-eulers-totient-function) – Juan T Feb 28 '17 at 17:27
  • Its not a duplicate, there's a difference in just finding Totient and finding a ratio with said toitient over large numbers. Finding a number that gets the ratio farther is a bit difficult and you aren't guaranteed to find a ratio like this simply. – Marorin Feb 28 '17 at 17:31

3 Answers3

1

The most likely candidates for high ratios of N/phi(N) are products of prime numbers. If you're just looking for one number with a ratio > 10, then you can generate primes and only check the product of primes up to the point where you get the desired ratio

def totientRatio(maxN,ratio=10):
    primes    = []
    primeProd = 1
    isPrime   = [1]*(maxN+1)
    p         = 2
    while p*p<=maxN:
        if isPrime[p]:
            isPrime[p*p::p] = [0]*len(range(p*p,maxN+1,p))
            primes.append(p)
            primeProd *= p
            tot = primeProd
            for f in primes:
                tot -= tot//f
            if primeProd/tot >= ratio:
                return primeProd,primeProd/tot,len(primes)
        p += 1 + (p&1)

output:

totientRatio(10**6)
16516447045902521732188973253623425320896207954043566485360902980990824644545340710198976591011245999110, 
10.00371973209101, 
55

This gives you the smallest number with that ratio. Multiples of that number will have the same ratio.

n = 16516447045902521732188973253623425320896207954043566485360902980990824644545340710198976591011245999110

n*2/totient(n*2) = 10.00371973209101

n*11*13/totient(n*11*13) = 10.00371973209101

No number will have a higher ratio until you reach the next product of primes (i.e. that number multiplied by the next prime).

n*263/totient(n*263) = 10.041901868473037

Removing a prime from the product affects the ratio by a proportion of (1-1/P).
For example if m = n/109, then m/phi(m) = n/phi(n) * (1-1/109)

(n//109) / totient(n//109)    = 9.91194248684247

10.00371973209101 * (1-1/109) = 9.91194248684247

This should allow you to navigate the ratios efficiently and find the numbers that meed your need.

For example, to get a number with a ratio that is >= 10 but closer to 10, you can go to the next prime product(s) and remove one or more of the smaller primes to reduce the ratio. This can be done using combinations (from itertools) and will allow you to find very specific ratios:

m = n*263/241

m/totient(m) = 10.000234225865265

m = n*(263...839) / (7 * 61 * 109 * 137)  # 839 is 146th prime 

m/totient(m) = 10.000000079805726
Alain T.
  • 40,517
  • 4
  • 31
  • 51
0

I have a partial solution for you, but the results don't look good.. (this solution may not give you an answer with modern computer hardware (amount of ram is limiting currently)) I took an answer from this pcg challenge and modified it to spit out ratios of n/phi(n) up to a particular n

import numba as nb
import numpy as np
import time

n = int(2**31)

@nb.njit("i4[:](i4[:])", locals=dict(
    n=nb.int32, i=nb.int32, j=nb.int32, q=nb.int32, f=nb.int32))

def summarum(phi):
    #calculate phi(i) for i: 1 - n
    #taken from <a>https://codegolf.stackexchange.com/a/26753/42652</a>
    phi[1] = 1
    i = 2
    while i < n:
        if phi[i] == 0:
            phi[i] = i - 1
            j = 2
            while j * i < n:
                if phi[j] != 0:
                    q = j
                    f = i - 1
                    while q % i == 0:
                        f *= i
                        q //= i
                    phi[i * j] = f * phi[q]
                j += 1
        i += 1
    #divide each by n to get ratio n/phi(n)
    i = 1
    while i < n: #jit compiled while loop is faster than: for i in range(): blah blah blah
        phi[i] = i//phi[i]
        i += 1
    return phi

if __name__ == "__main__":
    s1 = time.time()
    a = summarum(np.zeros(n, np.int32))
    locations = np.where(a >= 10)
    print(len(locations))

I only have enough ram on my work comp. to test about 0 < n < 10^8 and the largest ratio was about 6. You may or may not have any luck going up to larger n, although 10^8 already took several seconds (not sure what the overhead was... spyder's been acting strange lately)

Community
  • 1
  • 1
Aaron
  • 10,133
  • 1
  • 24
  • 40
0

p55# is a sparsely totient number satisfying the desired condition.

Furthermore, all subsequent primorial numbers are as well, because pn# / phi(pn#) is a strictly increasing sequence:

p1# / phi(p1#) is 2, which is positive. For n > 1, pn# / phi(pn#) is equal to pn-1#pn / phi(pn-1#pn), which, since pn and pn-1# are coprime, is equal to (pn-1# / phi(pn-1#)) * (pn/phi(pn)). We know pn > phi(pn) > 0 for all n, so pn/phi(pn) > 1. So we have that the sequence pn# / phi(pn#) is strictly increasing.

I do not believe these to be the only sparsely totient numbers satisfying your request, but I don't have an efficient way of generating the others coming to mind. Generating primorials, by comparison, amounts to generating the first n primes and multiplying the list together (whether by using functools.reduce(), math.prod() in 3.8+, or ye old for loop).

As for the general question of writing a phi(n) function, I would probably first find the prime factors of n, then use Euler's product formula for phi(n). As an aside, make sure to NOT use floating-point division. Even finding the prime factors of n by trial division should outperform computing gcd n times, but when working with large n, replacing this with an efficient prime factorization algorithm will pay dividends. Unless you want a good cross to die on, don't write your own. There's one in sympy that I'm aware of, and given the ubiquity of the problem, probably plenty of others around. Time as needed.

Speaking of timing, if this is still relevant enough to you (or a future reader) to want to time... definitely throw the previous answer in the mix as well.

Walker
  • 23
  • 4
Walker
  • 1
  • 1