Подтвердить что ты не робот

Чудновское двоичное расщепление и факторинг

В этой статье дается краткая рекурсивная формулировка формулы Чудновского pi с использованием бинарного расщепления. В 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

Этот метод уже очень быстрый, но он упомянул, что реализация на веб-сайте библиотеки GMP, gmp-chudnovsky.c, также влияет на числитель и знаменатель при двоичном расщеплении. Поскольку код оптимизирован и трудно понять, какова общая идея, как это делается? Я не могу сказать, упрощены ли дроби, числа хранятся в факторизованной форме, а не полностью умножаются, или и то, и другое.

Вот пример кода двоичного разбиения:

/* 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);
  }
}
4b9b3361

Ответ 1

Я не смотрел полный код, но я быстро взглянул на него, чтобы лучше понять выдержку, которую вы предоставили в своем вопросе.

Чтобы ответить на некоторые вопросы из вашего вопроса, сначала посмотрите на этот фрагмент кода:

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

Он определяет новый тип как структуру (если вы вообще не знаете C, допустим, что он похож на рудиментарный объект ввода объектов Pyhton, но без методов). Этот тип позволяет обрабатывать целые числа как два целых значения и два массива (скажем: два списка):

  • наибольший фактор
  • количество факторов
  • список факторов
  • список (соответствующих) полномочий для этих факторов

В то же время код сохраняет те же числа, что и большие целые числа из типа libgmp (это то, что подразумевается под mpz_t p и mpz_t g в аргументах функции).

Теперь, как насчет функции, которую вы показываете. Он называется fac_remove_gcd; начальная fac связана с именем ранее описанного типа; два следующих слова легко понять: разделите два целых числа типа fac на их gcd.

Код C выполняет итерацию по двум спискам факторов в обоих списках; легко синхронизировать оба списка, поскольку факторы упорядочены (раздел кода вокруг операторов else if и else); когда обнаруживаются два общих фактора (начальный оператор if), деление является вопросом простой подстановки: наименьшая мощность вычитается в обоих списках степеней для этого фактора (например, при a = 2 * 5 ^ 3 * 17 и b = 3 * 5 ^ 5 * 19, значение 3 будет вычитано в списках степеней для a и b в позиции, соответствующей фактору 5, приводящему к a = 2 * 5 ^ 0 * 17 и b = 3 * 5 ^ 2 * 19).

Во время этой операции создается номер (того же типа fac) и называется fmul; это, очевидно, gcd обоих чисел.

После этого шага gcd, называемый fmul, и тип fac преобразуется в большое целое число GMP с функцией (также в коде программы), называемой bs_mul. Это позволяет вычислить GCD как большое целое, чтобы синхронизировать новые значения разделенных чисел в обеих формах: большие целые числа и специальный тип fac. Когда GCD вычисляется как большое целое число, легко делить оба начальных числа на GCD.

Таким образом, функции действуют на две версии каждого начального числа.

Надеюсь, это поможет.

Ответ 2

Ключевыми являются следующие комментарии:

/*
    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)
*/

Заметим:

  • Если мы масштабируем как p, так и q в k в конечном делении, это не влияет на результат.
  • При вычислении операции слияния, которая соответствует второму набору комментариев, если бы мы масштабировали каждый из p(a,m), g(a,m) и q(a,m) на k, то этот коэффициент просто переносится на p(a,b), g(a,b), q(a,b); аналогично, если бы мы масштабировали каждый из p(m,b), g(m,b) и q(m,b), то этот коэффициент несли бы.
  • Линия

    fac_remove_gcd(p2, fp2, g1, fg1);
    

    эффективно

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

    Это имеет чистый эффект уменьшения масштаба p(a,b), g(a,b) и q(a,b) тем, что gcd. По предыдущим двум наблюдениям это масштабирование полностью переносится до конечного результата.

Postscript

Я пробовал три способа реализации факторинга в Python.

  • Отслеживание факторизаций с использованием неизменяемых списков замедляет ситуацию ужасно, потому что слишком много занятой работы, поддерживающей списки.
  • Использование gmpy gcd не дает ускорения
  • Предварительный подсчет списка простых чисел до 6 * N (т.е. который может делить g), и тестирование каждого из этих простых чисел замедляет работу в 2 - 3 раза.

Я пришел к выводу, что для ускорения такого подхода требуется использование изменяемого состояния для отслеживания факторизации, поэтому это большой успех для удобства обслуживания.