resueltos - divisiones sucesivas para descomponer cada numero en factores primos

Chudnovsky división binaria y factorización

Question

En este artículo , se proporciona una formulación recursiva rápida de la fórmula de Chudnovsky pi utilizando división binaria. En Python:

C = 640320
C3_OVER_24 = C**3 // 24
def bs(a, b):
    if b - a == 1:
        if a == 0:
            Pab = Qab = 1
        else:
            Pab = (6*a-5)*(2*a-1)*(6*a-1)
            Qab = a*a*a*C3_OVER_24
        Tab = Pab * (13591409 + 545140134*a) # a(a) * p(a)
        if a & 1:
            Tab = -Tab
    else:
        m = (a + b) // 2
        Pam, Qam, Tam = bs(a, m)
        Pmb, Qmb, Tmb = bs(m, b)

        Pab = Pam * Pmb
        Qab = Qam * Qmb
        Tab = Qmb * Tam + Pam * Tmb
    return Pab, Qab, Tab

N = int(digits/DIGITS_PER_TERM + 1)
# Calclate P(0,N) and Q(0,N)
P, Q, T = bs(0, N)
one = 10**digits
sqrtC = sqrt(10005*one, one)
return (Q*426880*sqrtC) // T

Este método ya es muy rápido, pero se menciona que una implementación en el sitio web de la biblioteca GMP, gmp-chudnovsky.c , también determina el numerador y el denominador en la división binaria. Como el código está optimizado y es difícil de entender para mí, ¿cuál es la idea general de cómo se hace esto? No puedo decir si las fracciones se simplifican, los números se mantienen en forma factorizada en lugar de multiplicarse por completo, o ambos.

Aquí hay una muestra de código de la división binaria:

/* binary splitting */
void
bs(unsigned long a, unsigned long b, unsigned gflag, long int level)
{
  unsigned long i, mid;
  int ccc;

  if (b-a==1) {
    /*
      g(b-1,b) = (6b-5)(2b-1)(6b-1)
      p(b-1,b) = b^3 * C^3 / 24
      q(b-1,b) = (-1)^b*g(b-1,b)*(A+Bb).
    */
    mpz_set_ui(p1, b);
    mpz_mul_ui(p1, p1, b);
    mpz_mul_ui(p1, p1, b);
    mpz_mul_ui(p1, p1, (C/24)*(C/24));
    mpz_mul_ui(p1, p1, C*24);

    mpz_set_ui(g1, 2*b-1);
    mpz_mul_ui(g1, g1, 6*b-1);
    mpz_mul_ui(g1, g1, 6*b-5);

    mpz_set_ui(q1, b);
    mpz_mul_ui(q1, q1, B);
    mpz_add_ui(q1, q1, A);
    mpz_mul   (q1, q1, g1);
    if (b%2)
      mpz_neg(q1, q1);

    i=b;
    while ((i&1)==0) i>>=1;
    fac_set_bp(fp1, i, 3);  /*  b^3 */
    fac_mul_bp(fp1, 3*5*23*29, 3);
    fp1[0].pow[0]--;

    fac_set_bp(fg1, 2*b-1, 1);  /* 2b-1 */
    fac_mul_bp(fg1, 6*b-1, 1);  /* 6b-1 */
    fac_mul_bp(fg1, 6*b-5, 1);  /* 6b-5 */

    if (b>(int)(progress)) {
      printf("."); fflush(stdout);
      progress += percent*2;
    }

  } else {
    /*
      p(a,b) = p(a,m) * p(m,b)
      g(a,b) = g(a,m) * g(m,b)
      q(a,b) = q(a,m) * p(m,b) + q(m,b) * g(a,m)
    */
    mid = a+(b-a)*0.5224;     /* tuning parameter */
    bs(a, mid, 1, level+1);

    top++;
    bs(mid, b, gflag, level+1);
    top--;

    if (level == 0)
      puts ("");

    ccc = level == 0;

    if (ccc) CHECK_MEMUSAGE;

    if (level>=4) {           /* tuning parameter */
#if 0
      long t = cputime();
#endif
      fac_remove_gcd(p2, fp2, g1, fg1);
#if 0
      gcd_time += cputime()-t;
#endif
    }

    if (ccc) CHECK_MEMUSAGE;
    mpz_mul(p1, p1, p2);

    if (ccc) CHECK_MEMUSAGE;
    mpz_mul(q1, q1, p2);

    if (ccc) CHECK_MEMUSAGE;
    mpz_mul(q2, q2, g1);

    if (ccc) CHECK_MEMUSAGE;
    mpz_add(q1, q1, q2);

    if (ccc) CHECK_MEMUSAGE;
    fac_mul(fp1, fp2);

    if (gflag) {
      mpz_mul(g1, g1, g2);
      fac_mul(fg1, fg2);
    }
  }

  if (out&2) {
    printf("p(%ld,%ld)=",a,b); fac_show(fp1);
    if (gflag)
      printf("g(%ld,%ld)=",a,b); fac_show(fg1);
  }
}

Answer 1

Los siguientes comentarios son clave:

/*
    g(b-1,b) = (6b-5)(2b-1)(6b-1)
    p(b-1,b) = b^3 * C^3 / 24
    q(b-1,b) = (-1)^b*g(b-1,b)*(A+Bb).
*/

/*
    p(a,b) = p(a,m) * p(m,b)
    g(a,b) = g(a,m) * g(m,b)
    q(a,b) = q(a,m) * p(m,b) + q(m,b) * g(a,m)
*/

/*
          p*(C/D)*sqrt(C)
    pi = -----------------
             (q+A*p)
*/

Observar:

1. Si escalamos p y q por un factor de k en la división final, no afecta el resultado.
2. Al calcular la operación de fusión que corresponde al segundo conjunto de comentarios, si tuviéramos que escalar cada uno de p(a,m) , g(a,m) q(a,m) por k entonces ese factor simplemente sería llevado a p(a,b) , g(a,b) , q(a,b) ; de manera similar, si tuviéramos que escalar cada uno de p(m,b) , g(m,b) q(m,b) entonces ese factor se llevaría a cabo.

3. La línea

   fac_remove_gcd(p2, fp2, g1, fg1);

   es efectivamente

   k = gcd(p2, g1);
   p2 /= k; // p(m,b) /= k
   g1 /= k; // g(a,m) /= k

   Esto tiene el efecto neto de downscaling p(a,b) , g(a,b) q(a,b) por ese gcd . En las dos observaciones anteriores, esta reducción de escala se lleva a cabo limpiamente hasta llegar al resultado final.

Posdata

He intentado tres formas de implementar el factoring en Python.

• Rastrear las factorizaciones usando listas inmutables ralentiza las cosas horrendamente porque hay demasiado trabajo ocupado manteniendo las listas.
• Usar gcd de gmpy no da una aceleración
• El cálculo previo de una lista de números primos hasta 6 * N (es decir, que pueden dividir g ) y la prueba de cada uno de esos números primos reduce la velocidad de los factores en un factor de 2 a 3.

Concluyo que obtener una aceleración con este enfoque requiere el uso de estado mutable para rastrear las factorizaciones, por lo que es un gran éxito para la facilidad de mantenimiento.

Answer 2

No miré el código completo, pero lo eché un vistazo rápido para comprender mejor el extracto que me proporcionó en su pregunta.

Con el fin de responder a algunos puntos de su pregunta, primero mire este fragmento de código:

typedef struct {
  unsigned long max_facs;
  unsigned long num_facs;
  unsigned long *fac;
  unsigned long *pow;
} fac_t[1];

Define un nuevo tipo como una estructura (si no conoces a C en absoluto, digamos que es como un objeto Pyhton rudimentario que incorpora variables pero no tiene ningún método). Este tipo permite manejar números enteros como dos valores enteros y dos matrices (por ejemplo: dos listas):

• factor más grande
• numero de factores
• lista de factores
• lista de poderes (correspondientes) para estos factores

Al mismo tiempo, el código mantiene los mismos números que los enteros grandes del tipo libgmp (esto es lo que significan mpz_t p y mpz_t g en los argumentos de la función).

Ahora, ¿qué pasa con la función que está mostrando? Se llama fac_remove_gcd ; el fac inicial tiene que ver con el nombre del tipo descrito anteriormente; las dos palabras siguientes son fáciles de entender: divida dos números enteros de tipo fac por su gcd .

El código C itera sobre las dos listas de factores en ambas listas; es fácil sincronizar ambas listas ya que los factores están ordenados (sección del código alrededor de las sentencias else if y else ); siempre que se detecten dos factores comunes (enunciado inicial if ), dividir es una cuestión de simple sustracción: la menor potencia se sustrae en ambas listas de potencias para este factor (por ejemplo, con a = 2 * 5 ^ 3 * 17 yb = 3 * 5 ^ 5 * 19, el valor 3 se restará en las listas de potencias para ayb en la posición correspondiente al factor 5 que lleva a a = 2 * 5 ^ 0 * 17 yb = 3 * 5 ^ 2 * 19 )

Durante esta operación, se crea un número (del mismo tipo de fac ) y se llama fmul ; obviamente es el gcd de ambos números.

Después de este paso, el gcd llamado fmul y el ser del tipo fac se convierte en un entero grande GMP con la función (también en el código del programa) llamada bs_mul . Esto permite calcular el GCD como un número entero grande para sincronizar los nuevos valores de los números divididos en ambas formas: enteros grandes y tipo de fac especial. Una vez que el GCD se computa como un número entero grande, es fácil dividir ambos números iniciales por el GCD.

Por lo tanto, las funciones actúan sobre dos versiones de cada número inicial.

Espero que pueda ayudar.