Suma de la función totiente $\sum_{k=1}^n \varphi(k)$

Question

Exploré algunos sitios de desafíos de ciencias de la computación/teoría de los números por diversión, y presentaron el siguiente problema, exactamente como sigue:

Dejemos que $$P(n) = \sum_{k=1}^n \varphi(k)$$

Encuentre $P(10^{16})$

He buscado durante bastante tiempo sobre esto y he probado diferentes enfoques:

1. Utilizando la fórmula de $$\varphi(n)= n \cdot \prod_{i=1}^k \frac{p_i-1}{p_i}$$ Intenté calcular cada $\varphi(n)$ en el rango, pero esto se vuelve muy ineficiente para grandes $n$ . Podría llegar hasta $10^7$ con este enfoque. Más allá de esto, se vuelve demasiado lento.

2. Intenté una diferente, más directa. Wikipedia y Wolfram Alpha sugieren fórmulas similares para calcular directamente $P(n)$ : $$P(n) = \sum_{k=1}^n \varphi(k)= \frac12 \cdot \biggl (1+ \sum_{k=1}^n\mu (k)\cdot \lfloor {\frac nk} \rfloor^2\biggl)$$ Esta fórmula parecía mucho más prometedora. La probé y conseguí llegar mucho más lejos que $10^7$ pero aún está lejos del objetivo. Calculando previamente un tamiz de la función Moebius, pude llegar a un poco menos de $10^9$ . Mi memoria era insuficiente y no podía calcular más valores en un colador. E incluso si pudiera, todavía toma mucho tiempo y está muy lejos de $10^{16}$ .

Aquí está parte del código que he utilizado para mi segundo enfoque escrito en Java:

public static BigInteger PhiSummatoryFunction (long limit)
{
    BigInteger sum = BigInteger.ZERO;
    int [] m = MoebiusSieve(limit);
    for (int i=1;i<m.length;i++)
        sum=sum.add(BigInteger.valueOf((long) (m[i]*Math.floor(limit/i)*Math.floor(limit/i))));
    return sum.add(BigInteger.ONE).divide(BigInteger.ONE.add(BigInteger.ONE));
}

Donde MoebiusSieve es una función que calcula los valores de la función Moebius hasta un determinado límite en un tamiz, utilizando un método similar al de Eratóstenes.

---

3. Después de entender e implementar el método recursivo sugerido en un enlace proporcionado en los comentarios: $$P(n) = \frac {n(n+1)}{2} - \sum_{i=2}^\sqrt n P(\lfloor \frac ni \rfloor) - \sum_{j=1}^\sqrt n P(j) \cdot (\lfloor \frac nj \rfloor - \lfloor \frac n{j+1} \rfloor)$$

Puedo calcular valores hasta $P(10^{11})$ y con la máxima asignación de memoria, precalculando tantos $\varphi(n)$ como sea posible y, en consecuencia, todos los $P(n)$ que puedo para la memoización, puedo calcular $P(10^{12})$ en poco más de 20 minutos. Una gran mejora, pero todavía un poco lejos de $P(10^{16})$ . No pasa nada si el cálculo tarda un poco más, pero me temo que $P(10^{16})$ tardaría un tiempo exponencialmente mayor, a juzgar por el "salto" en el tiempo de cálculo entre $P(10^{11})$ y $P(10^{12})$ . Mi memoria me permite "guardar" hasta $350,000,000 \space (n)$ valores O hasta $700,000,000 \space (k)$ valores. ¿Quizás haya una forma de realizar la suma utilizando los valores (k) en lugar de los (n)?

Todos mis cálculos sugieren y demuestran que mi recursión es la que más tiempo consume. Esto es obvio, pero estoy seguro de que consume más tiempo del que debería, como señala qwr . Así que a continuación pongo el código de recursión, con algo de documentación. Me parece que esta es la forma correcta de hacer este cálculo, pero mi implementación no es óptima.

    public static BigInteger phiR (long limit, long [] s) // limit is 10^t, s is the sieve of precomputed values of `P(n)`. Can store maximum 350,000,000 values
{                                                                                                                                               
    if (limit<s.length)                                 
        return BigInteger.valueOf(s[(int) limit]);
    BigInteger sum = BigInteger.valueOf(limit).multiply(BigInteger.valueOf(limit).add(BigInteger.ONE)).divide(BigInteger.valueOf(2)); // this corresponds to the n'th triangular number
    BigInteger midsum1=BigInteger.ZERO;
    BigInteger midsum2=BigInteger.ZERO;
    for (long m=2;m*m<=limit;m++) // computing the first sum, first for changing floor(limit/m) values
        midsum1=midsum1.add(phiR((long) Math.floor(limit/m),s));
    for (long d=1;d*d<=limit;d++) // computing the second sum
        if ((double)d!=Math.floor(limit/d))
            midsum2=midsum2.add(BigInteger.valueOf((long) (Math.floor(limit/d)-Math.floor(limit/(d+1)))).multiply(phiR(d,s)));
    sum=sum.subtract(midsum1).subtract(midsum2);
    return sum;
}

Me sugirieron que usara dictámenes por qwr además de la matriz, para los valores grandes de $n$ pero no sé nada al respecto. Se puede hacer otra mejora para que el plazo sea de un día más o menos?

---

Answer 1

La solución recursiva que mostró toma $O(n^{2/3})$ tiempo (ver mi respuesta ). Hasta donde yo sé, esto es lo más rápido que se puede conseguir para obtener respuestas exactas utilizando enfoques recursivos. Mi implementación en PyPy calcula $P(10^{12})$ en 22 segundos (extrapolando $P(10^{16})$ tardará entre 1 y 2 días con suficiente memoria), así que tendrás que trabajar en la eficiencia de tu implementación.

Mi algoritmo es casi una traducción directa de mi respuesta enlazada. Actualmente tengo un límite fijado en $n^{2/3}$ . Para todos los $n$ menos de este límite, las sumas de los totientes/totales se calculan previamente utilizando un Enfoque del tamiz de Eratóstenes . Para los valores de $n$ mayor que el corte, calcular recursivamente y almacenar en un diccionario. Para $P(10^{16})$ es probable que se quede sin memoria al intentar almacenar todos los valores de los totales y que tenga que cargar y almacenar desde el disco.

Puede buscar las referencias en OEIS A002088 para obtener recursos que a menudo me resultan útiles.

En la OEIS, sólo los términos hasta $P(10^{18})$ son conocidos ( A064018 ). Los valores más altos fueron calculados por Hiroaki Yamanouchi (una investigación superficial indica que puede ser el legendario Min_25 en el Proyecto Euler, SPOJ, CodeChef, etc.) en OEIS; puede intentar ponerse en contacto con él para ver cómo lo hizo.

Answer 2

Asintóticamente, es posible que pueda hacerlo mejor utilizando el enfoque descrito en Lagarias y Odlyzko, Informática $\pi(x)$ Un método analítico (1987) J. Algorithms vol.8 pp.173-191. En la sección 2 se discuten las condiciones que bastan para adaptar el enfoque a funciones distintas de $\pi(x)$ y $\Phi(x)$ parece que los conoce. Sin embargo, hay que señalar que aunque asintóticamente su algoritmo es $O(x^{1/2 + \varepsilon})$ La constante oculta puede ser bastante más alta que el algoritmo que está utilizando actualmente.