Segunda parte de la pregunta de divisibilidad de la suma factorial

Question

Qué primos $p$ dividir la suma de los factoriales $1! + 2! + 3! + 4! + 5! + \cdots + (p-1)!$ ? Esto está relacionado con mi pregunta anterior .

Answer 1

Obsérvese que heurísticamente, si estos números estuvieran equidistribuidos mod $p$ esperaríamos que cada uno fuera cero (es decir, esperaríamos $p$ para dividir el resultado) con probabilidad $\approx 1/p$ suponiendo que todos los valores son independientes (lo que parece una suposición razonable), deberíamos esperar tener $\sum_{p\lt n}\frac{1}{p}\approx M+\ln \ln n$ de ellos menos de $n$ , donde $M$ es el Constante de Meissel-Mertens - esta suma es aproximadamente 2,887 para $n=10^6$ Así que, aunque heurísticamente "esperemos" otra en ese rango, no es sorprendente que no haya ninguna, y a pesar de la aparente falta de más soluciones, debido a la divergencia de la serie, ¡el número "esperado" de soluciones sigue siendo infinito!

Añadido: He montado mi propia implementación de software en C++, confirmando los resultados (que sólo $p=3$ y $p=11$ trabajo) para $p\lt 5\times 10^5$ pero también observando la distribución de los resultados $\mod p$ para ver si surgía algún patrón. Calculé el valor de $\dfrac{\sum_{1}^{p-1}i!\bmod p}{p}$ (la división por $p$ que sirve para "normalizar" los valores en el rango $[0, 1)$ ) para todos los $p$ en ese rango y luego los agrupó en bins; estos son los gráficos con 100 bins y 256 bins. (Estos recuentos de bins deben ser lo suficientemente pequeños en relación con los valores primos en consideración que no debería haber ningún efecto de aliasing sustancial en el binning de los datos). Los resultados son un poco desiguales, pero no parece haber ninguna desviación hacia valores concretos (por ejemplo $\frac{p-1}2$ ) que aparecerá como un pico en los contenedores. Si tengo más tiempo, puedo echar un vistazo al cuadrado de la suma y, sobre todo, a la inversa de la suma $\pmod p$ pero como esto último requeriría que escribiera un algoritmo GCD rápido para encontrar la inversa, tendrá que esperar para más adelante.

(Nota: la caída brusca en el último elemento es un artefacto de mi adición de un valor 0 al final de los datos de la bandeja para que mi hoja de cálculo no normalizara automáticamente el rango; no es un valor real).

Answer 2

Tengo algunos pequeños comentarios analíticos y una búsqueda numérica relativamente agresiva.

En primer lugar, he comprobado todos los impar $p$ (no necesariamente primo) hasta $10^6$ y encontré sólo los mismos cinco que Yoni arriba, es decir, 3, 9, 11, 33 y 99. Esto tomó ~1 hora usando 1 núcleo de mi MacBook Pro, con una versión temprana no optimizada de mi código.

En segundo lugar, por supuesto que sí: \begin{equation} p\nmid\sum_{n=1}^{p-1}n!\quad \forall p\in\mathbb{N}, p \textrm{ even} \end{equation} porque todos los términos de la suma son pares excepto el primero, por lo que la suma es siempre impar.

En tercer lugar, la búsqueda numérica se acelera ligeramente utilizando la identidad: \begin{equation} \sum_{n=1}^{p-1} n! = 1+2[1+3[1+4[...[1+(p-1)]]]] \end{equation} He escrito un software para comprobar los primos $p$ para $p|\sum_{n=1}^{p-1} n!$ utilizando la identidad anterior, y configurarlo para utilizar múltiples hilos de forma vergonzosamente paralela. Usando este software he mostrado: \begin{equation} \textrm{3 and 11 are the only primes }p<10^7\textrm{ that divide }\sum_{n=1}^{p-1} n! \end{equation} Si quiere compilarlo usted mismo, tendrá que enlazarlo con la biblioteca primesieve.

#include <stdio.h>
#include <stdlib.h>
#include <primesieve/soe/PrimeSieve.h>

// To compile: g++ -O3 factorialSum.cpp -lprimesieve -o factorialSum
//         or: g++ -O3 -static factorialSum.cpp -lprimesieve -o factorialSum

// We perform arithmetic at u1 precision; p is restricted to u2 precision
typedef uint_fast64_t u1;
typedef uint_fast32_t u2;

// This evaluates $\sum_{n=1}^{p-1} n!$ based on the definition.
// We don't actually use this in our program below.
u2 factorialSum_modp_method1(u2 p){
  u1 factorial_modp=1;
  u1 factorialSum_modp=0;
  for(u1 n=1; n<=p-1; n++){
    factorial_modp*=n;
    factorial_modp%=p;
    factorialSum_modp+=factorial_modp;
    factorialSum_modp%=p;
  }
  return factorialSum_modp;
}

// This evaluates $\sum_{n=1}^{p-1} n! = 1+2[1+3[1+4[...[1+(p-1)]]]]$ via the RHS of the expression.
// This is somewhat faster than 'method1' above.
u2 factorialSum_modp_method2(u2 p){
  u1 factorialSum_modp=1;
  for(u1 n=p-1; n>1; n--)
    factorialSum_modp=(factorialSum_modp*n+1)%p;
  return factorialSum_modp;
}

// This prints out p iff p divides $\sum_{n=1}^{p-1} n!$.
// We pass this as an arument to PrimeSieve.generatePrimes().
void check_factorialSum_modp(u2 p){
  if(factorialSum_modp_method2(p)==0)
    printf("%u\n", (unsigned)p);
}

void error_incorrectUsage(char *args[]){
  printf("Usage: %s thread numThreads blockSize pMin pMax\n", args[0]);
  printf("       thread is a positive integer 0, ..., numThreads-1\n");
  printf("       numThreads is a natural number, 1, 2, ...\n");
  printf("       blockSize is the interleave size between threads\n");
  printf("       pMin and pMax are integers, the upper and lower limits on candidate primes p\n\n");
  printf("e.g.:  %s 0 1 10000 0 100000\n", args[0]);
  printf("       will check [0, 100000] for primes p that divide $\\sum_{n=1}^{p-1} n!$ \n");
  printf("       (expression written LaTeX).\n\n");
  printf("In order to use multiple threads, run multiple instances with different values of 'thread'.\n");
  exit(1);
}

int main(int nargs, char *args[]){

  // Parse command line arguments and check for proper usage
  if( !(nargs==6) ) error_incorrectUsage(args);
  char *error1, *error2, *error3, *error4, *error5;
  u2 thread     = strtoul(args[1], &error1, 10);
  u2 numThreads = strtoul(args[2], &error2, 10);
  u2 blockSize  = strtoul(args[3], &error3, 10);
  u2 pMin       = strtoul(args[4], &error4, 10);
  u2 pMax       = strtoul(args[5], &error5, 10);
  if(*error1 || *error2 || *error3 || *error4 || *error5) error_incorrectUsage(args);

  // Sieve intervals [pLower, pUpper] that belong to this thread
  PrimeSieve ps;
  for(u2 pLower=pMin+thread*blockSize; pLower<=pMax-blockSize; pLower+=numThreads*blockSize){
    u2 pUpper=pLower+blockSize;
    ps.generatePrimes(pLower, pUpper, check_factorialSum_modp);
    printf("interval [%u,%u] completed on thread %u\n", pLower, pUpper, thread);
  }

}

Answer 3

Lo siguiente es más un comentario que una respuesta utilizable, pero es demasiado largo para la caja de comentarios

Otra observación que tiene una simetría sorprendente, que parece ocurrir sólo si el módulo es un primo (y esto es, por tanto, una observación en sí misma).

Definimos la matriz triangular inferior de los números eulerianos:

$ \qquad \qquad \displaystyle E= \Tiny \left[ \begin{array} {rrrrr} 1 & . & . & . & . & . \\ 1 & 0 & . & . & . & . \\ 1 & 1 & 0 & . & . & . \\ 1 & 4 & 1 & 0 & . & . \\ 1 & 11 & 11 & 1 & 0 & . \\ 1 & 26 & 66 & 26 & 1 & 0 \end{array} \right] $

Observamos, que las sumas de las filas se acumulan a los factoriales por lo que $E \cdot V(1) = F $ donde $V(1)$ es la columna-vector de unos, da la columna-vector de factoriales $F$ . Si ahora sumamos las filas hasta el índice del primo $p$ en cuestión y tomar este módulo $p$ entonces tenemos otra expresión de la suma. Si ahora cambiamos el orden de cálculo: calculamos primero las sumas de columnas (módulo p) y luego obtenemos un vector fila con entradas simétricas, aparentemente si p es primo. Por ejemplo, para $p=13$ (y tenga en cuenta que en Pari/GP el índice comienza en 1 y no en 0)

$ \qquad \qquad V(1)^T_{p+1} ~ \cdot E_{p+1,p+1} \pmod p = \\\qquad \qquad \left[ \begin{array} {rrrrrrrrrrrrr} 1 & 2 & 3 & 11 & 6 & 12 & 5 & 12 & 6 & 11 & 3 & 2 & 1 & 0 \end{array} \right]$

y luego sumamos el vector resultante modulo(p) y restamos 1 porque tenemos también el 0! en la suma.
No veo en realidad lo que esto podría darnos -y no puedo profundizar en esto- pero también puede ser significativo (y posiblemente en el camino hacia una forma -más- cerrada) que las columnas del triángulo euleriano sean expresables como secuencias de términos principales de series geométricas y sus derivadas y así podrían tener una forma cerrada para sus segmentos principales.
De todos modos -quizás podamos hacer que la fila resultante y su suma modular sean más simples mediante algún patrón.

---

Aquí están los vectores resultantes para los primeros primos: $ \Tiny \begin{matrix} 2: & 1 & 1 & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 3: & 1 & 2 & 1 & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 5: & 1 & 2 & 3 & 2 & 1 & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 7: & 1 & 2 & 3 & 1 & 3 & 2 & 1 & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 11: & 1 & 2 & 3 & 10 & 8 & 8 & 8 & 10 & 3 & 2 & 1 & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 13: & 1 & 2 & 3 & 11 & 6 & 12 & 5 & 12 & 6 & 11 & 3 & 2 & 1 & . & . & . & . & . & . & . & . & . & . & . \\ 17: & 1 & 2 & 3 & 13 & . & 13 & 15 & 9 & 3 & 9 & 15 & 13 & . & 13 & 3 & 2 & 1 & . & . & . & . & . & . & . \\ 19: & 1 & 2 & 3 & 14 & 15 & 1 & 7 & 8 & 15 & 10 & 15 & 8 & 7 & 1 & 15 & 14 & 3 & 2 & 1 & . & . & . & . & . \\ 23: & 1 & 2 & 3 & 16 & 9 & . & 15 & 12 & 18 & 12 & 6 & 17 & 6 & 12 & 18 & 12 & 15 & . & 9 & 16 & 3 & 2 & 1 & . \\ 29: & 1 & 2 & 3 & 19 & 24 & 9 & 18 & . & 28 & 23 & 7 & 24 & 22 & 14 & 6 & 14 & 22 & 24 & 7 & 23 & 28 & . & 18 & 9 \end{matrix} $

[Si descomponemos las columnas de la matriz euleriana en sus componentes-series nos involucramos con los segmentos principales de las series geométricas y sus derivadas - y todo ello modulo p. Como al final necesitamos la suma de las columnas-sumas he expresado la suma parcial final en términos de sus componentes-series (truncadas), aquí para el ejemplo primo p=5 . Obtenemos como resultado final la suma de todos los elementos de la siguiente matriz $$\small \begin{bmatrix} 1\cdot 5^0 & 0\cdot 4^0 & 0\cdot 3^0 & 0\cdot 2^0 & 0\cdot 1^0 \\ 1\cdot 5^1 & -1\cdot 4^1 & 0\cdot 3^1 & 0\cdot 2^1 & 0\cdot 1^1 \\ 1\cdot 5^2 & -2\cdot 4^2 & 1\cdot 3^2 & 0\cdot 2^2 & 0\cdot 1^2 \\ 1\cdot 5^3 & -3\cdot 4^3 & 3\cdot 3^3 & -1\cdot 2^3 & 0\cdot 1^3 \\ 1\cdot 5^4 & -4\cdot 4^4 & 6\cdot 3^4 & -4\cdot 2^4 & 1\cdot 1^4\\ 1\cdot 5^5 & -5\cdot 4^5 & 10\cdot 3^5 & -10\cdot 2^5 & 5\cdot 1^5 \end{bmatrix} \pmod 5 $$ donde -porque p es un primo- tenemos ceros sistemáticos y congruencias expresables en números negativos $\pmod p$ y consigue esto:

$$\small \left[ \begin{array} {rrrrrrrrrrr} 1 & . & .& . & . \\ . & 1\cdot 1^1 & . & . & . \\ .& -2\cdot 1^2 & 1\cdot 2^2 & . & . \\ . & 3\cdot 1^3 & -3\cdot 2^3 & 1\cdot 3^3 & . \\ . & -4\cdot 1^4 & 6\cdot 2^4 & -4\cdot 3^4 & 1\cdot 4^4\\ . & . & . & . & . \end{array} \right] \pmod 5 $$ que se puede generalizar a otros módulos primos p de forma obvia.
La idea con todo esto es tener alguna ventaja para una mejor/más intuitiva expresión de la suma quizás por formas cerradas de la serie geométrica truncada (módulo p). Todavía podemos reconocer las sumas de filas como factoriales (mientras que el módulo p), pero el cambio de orden de cálculo y las formas cerradas de las sumas de columnas, que esperamos sean ventajosas, podrían ser útiles.

Podemos hacer aún más desde la modularidad. La siguiente representación se puede hacer análogamente incluso para primos p más altos, pero utilizo sólo el ejemplo anterior. La matriz anterior se puede reescribir en congruencia completa como $$\small \left[ \begin{array} {rrrrrrrrrrr} 1 & . & .& . & . \\ . & 1\cdot 1^1 & . & . & . \\ .& 3\cdot 1^2 & 1\cdot 2^2 & . & . \\ . & 3\cdot 1^3 & 2\cdot 2^3 & 1\cdot 3^3 & . \\ . & 1\cdot 1^4 & 1\cdot 2^4 & 1\cdot 3^4 & 1\cdot 4^4\\ . & . & . & . & . \end{array} \right] \pmod 5 $$ y esos binomios simples dan entonces como columnas $$\small \left[ \begin{array} {rrrrrrrrrrr} 1 & 1^1\cdot 2^3 & 2^2\cdot 3^2 & 3^3\cdot 4^1 & 1\\ \end{array} \right] \pmod 5 $$ y si junto cada dos términos leyendo de izquierda a derecha entonces esto es
$$\small \begin{array} {rrrrrrrrrrr} 2 \cdot 1 & + 2 \cdot 2^3 & + 2^2\cdot 2^2 \\ \end{array} \equiv 4\pmod 5 $$ donde por la definición en el OP tuvimos que decrementar en 1 para borrar el $0!$ plazo.

Si escribimos el primo impar $p=2q+1$ entonces la fórmula general (incluyendo el $0!$ ) se convierte ahora: $$ S_p \equiv \sum_{k=0}^{p-1} k! \equiv 2 \sum _{k=0}^{q-1} k^k (k+1)^{p-1-k} + (q \cdot -q)^q \pmod p $$

que podría ser finalmente accesible con un poco más de conocimiento sobre las raíces primitivas módulo p.

Answer 4

La respuesta de Gottfried y las respuestas de Steven ponen de manifiesto que hay que elegir en una búsqueda numérica directa: o bien se puede calcular $(p-1)!$ $(\textrm{mod }p)$ por separado para cada $p$ o puede almacenar $(p-1)!$ y luego calcular $p! = p(p-1)!$ para probar el siguiente $p$ . Este mismo compromiso aparece en otros problemas similares, como la búsqueda de Números marrones .

Para comprobar cada $p<N$ El primer enfoque requiere $O(N^2)$ multiplicaciones de números de orden $p$ El segundo enfoque requiere $N$ multiplica, pero cada multiplicación es muy costosa, estando entre un número de orden $p$ y una de orden $p!$ . Otra consideración es que la segunda solución, almacenar $(p-1)!$ sin reducir, requiere una memoria lineal en $p$ . Realmente no quiero ver lo que pasa cuando $(p-1)!$ ya no cabe en la caché. Además, no es obvio cómo paralelizar el segundo algoritmo.

Mi impresión es que si quiere comprobar todo $p<N$ y no sólo los primos $p$ es más rápido para almacenar $(p-1)!$ para los pequeños $N \lesssim 10^6$ pero para los más grandes $N$ es más rápido de calcular $(p-1)!$ $(\textrm{mod }p)$ por separado para cada $p$ . Si quiere comprobar sólo los primos $p$ entonces la "primera aproximación" obtiene un $O(\log(n))$ y el "segundo enfoque" apenas se beneficia, lo que significa que el primer enfoque es más rápido: véase el " Actualización: " debajo del código para una explicación.

Aquí hay un código C que utiliza el GMP para implementar este "segundo enfoque"; yo implementé el "primer enfoque" en mi primera respuesta anterior.

#include <gmp.h>
#include <time.h>

#define NMAX 100000

int main(){

    // Print begin time & date
    time_t beginTime;
    struct tm * timeinfo;
    time ( &beginTime );
    timeinfo = localtime ( &beginTime );
    gmp_printf("# Begin time/date is: %s", asctime (timeinfo) ); 

    // Stores n! and sum(n!)
    mpz_t cur_factorial;
    mpz_t sum_factorial;
    mpz_t sum_factorial_modn;
    mpz_init_set_ui(cur_factorial, 1);
    mpz_init_set_ui(sum_factorial, 1);
    mpz_init(sum_factorial_modn);

    for(unsigned long int n=2; n<=NMAX; n++){

      // Check to see if n divides sum(m=1..n-1) m!
      if( mpz_divisible_ui_p(sum_factorial, n) )
    gmp_printf("%u\n", n);

      // Update n!
      mpz_mul_ui(cur_factorial, cur_factorial, n);

      // Update sum(n!)
      mpz_add(sum_factorial, sum_factorial, cur_factorial);

    }

    // Print end time & date; total run time
    time_t endTime;
    time ( &endTime );
    timeinfo = localtime ( &endTime );
    gmp_printf("# End time/date is:   %s", asctime (timeinfo) ); 
    gmp_printf("# Total runtime of %.0f seconds.\n", difftime(endTime, beginTime) ); 

    // Release and return
    mpz_clear(cur_factorial);
    mpz_clear(sum_factorial);
}

---

Actualización:

Este código que implementa el "segundo enfoque" es estrictamente más lento que el de mi respuesta anterior. El código anterior tarda ~9s para un límite superior de $10^5$ , que sólo tarda ~5s usando el código de mi respuesta anterior. Una búsqueda hasta $10^6$ tarda unos 19 minutos con este código, frente a unos 8 minutos con el otro. En parte, esto se debe a que, con el otro código, se obtiene un enorme aumento de velocidad si sólo se prueban los primos $p$ . Este código está configurado para comprobar todos los $n$ incluso las que no tienen sentido. Puede que no parezca una comparación justa, pero modificando el código anterior para que sólo compruebe impar $n$ no proporciona ninguna mejora en la velocidad; la cuestión es que de todos modos hay que calcular $(n-1)!$ por cada $n$ con este método, así que ya has pagado el coste; la comprobación de la divisibilidad no lleva ningún tiempo en comparación con las multiplicaciones.

Answer 5

Además del post de @Douglas, aquí hay una implementación muy corta (y creo que: eficiente) en Pari/GP - sin embargo, yo mismo no lo ejecuté a N alto: [Actualización] He añadido un método de cálculo mejorado al final[/actualización]

MaxN=101   \\ set some upper limit for the computation
S=[1,1]    \\ contains current factorial and current sum
for(k=3,MaxN, S *= [k-1,k-1;0,1]; if(isprime(k),print([k,S[1] % k,S[2]%k])))

Resultado :

p | (p-1)! %p | S %p
--+-----------+-----
[3, 2, 0]
[5, 4, 3]
[7, 6, 5]
[11, 10, 0]
[13, 12, 9]
[17, 16, 12]
[19, 18, 8]
[23, 22, 20]
[29, 28, 16]
[31, 30, 1]
[37, 36, 4]
[41, 40, 3]
...    

Hasta el primo 1000 MaxN=7919

MaxN = prime(1000)
...
...
[7919, 7918, 2882]

esto necesitó 281 mseg.

Una versión que utiliza el cómputo del módulo es la siguiente

gettime()
{for(j=1,1000,p=prime(j); 
   S=[1,1];
   for(k=3,p,
        S*=[k-1,k-1;0,1];
         S=S % p
       );
   print([p,S[1] % p,S[2]%p])
   )}
gettime()

%1124 = 9219   \\ this is msecs

---

[Lo siguiente sólo es relevante, si queremos una rutina, que produzca la lista sobre todos los primos hasta algún límite.
Tiene entonces alguna ventaja hacer un compromiso entre el primer método en el que necesitamos todos los factoriales y sumas una sola vez pero que necesita muchos bits para la representación de los números y el modulo-computación completo, que opera sobre números pequeños pero debe hacer todo el recorrido de cálculo de los factoriales (residuales) para cada primo de forma que tenemos un doble bucle.
Si calculamos y almacenamos los factoriales y las sumas parciales necesarias en trozos de, digamos, $t=10$ términos, de manera que tengamos una lista de valores precalculados $$ \small \begin{array} {r|l|l|} & f & S \\ \hline \\ c_1=& 10! & \sum_{k=1}^{10} k! \\ c_2=& {20!\over 10!} & {\sum_{k=11}^{20} k! \over 10!} \\ c_3=& {30!\over 20!} & {\sum_{k=21}^{30} k! \over 20!} \\ \ldots \end{array} $$ entonces podemos componer nuestros resultados para cada primo mediante la multiplicación modular de esas constantes y un resto corto. Afinando el parámetro $t$ para mayor permite entonces encontrar un compromiso entre el almacenamiento/tamaño de los números almacenados y el número de pasos en el bucle interno, que contiene el cálculo modular por primo.
Utilizando algunos datos empíricos he hecho una estimación con una curva cuadrática en Excel, que dice, para la lista de primos hasta 1000000 necesito 90 min con tamaño de bloque $t=20$ 30 minutos con $t=100$ y 15 minutos con $t=500$ .
(Si se desea, puedo adjuntar aquí el código Pari/GP)