Функция Эйлера

[Geen]

Verbery в сообщении #1553976 писал(а):

Вот конечный рабочий вариант:

Есть у меня ощущение, что этот вариант категорически неверен (несмотря на тесты)....

[mihaild]

Я причин, почему это должно работать неправильно, не вижу. Функция Эйлера может посчитаться некорректно (например от $1913^3$ получится отрицательным), но вроде бы ни на каких таких аргументах мы её значение вычислять не будем. Еще в некоторых местах лишние модули, но там они ничего не делают.

[Geen]

mihaild в сообщении #1553981 писал(а):

Еще в некоторых местах лишние модули, но там они ничего не делают.

Что бы было где проверять результат (например, в WolframAlpha), надо брать $n$ и $k$ меньше 10000. И брать по модулю 10007. И сразу станет понятно - лишние модули или неправильные.

-- 06.05.2022, 13:22 --

И, на самом деле, вообще непонятно почему третий цикл лучше взятия целой части от деления (или целочисленного деления)...

[mihaild]

Geen в сообщении #1553984 писал(а):

Что бы было где проверять результат (например, в WolframAlpha), надо брать $n$ и $k$ меньше 10000.

Неужели WolframAlpha не умеет в длинную арифметику??
И тут для корректности реализации важно, что модуль больше $2n$ (без этого даже с вычислением простых будет непонятно что).

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

[Geen]

mihaild в сообщении #1553987 писал(а):

Неужели WolframAlpha не умеет в длинную арифметику??

Я думаю он просто отказывается факторизовать числа длиной в десятки тысяч десятичных цифр
По крайней мере, в бесплатой версии

-- 06.05.2022, 14:09 --

mihaild в сообщении #1553987 писал(а):

И тут для корректности реализации важно, что модуль больше $2n$

На мой взгляд, это противоестественное ограничение...

[B@R5uk]

Кто-нибудь может сказать столько простых сомножителей (с учётом кратности) в числе

$$500\;000\;!$$

?

Хочется прикинуть оценку сложности разложения биномиального коэффициента на простые сомножители.

[mihaild]

B@R5uk, в смысле сумма степеней в разложении на простые множители? 1786531. Но нужно учитывать, что простые множители нетипично маленькие для числа такого размера, так что алгоритм, который будет пробовать искать разложение, не проверив предварительно почти все нужные маленькие (до 500000) простые числа вряд ли имеет шансы на успех, а вот просто перебирающий все простые подряд - должен работать нормально.

[Dmitriy40]

B@R5uk
Очевидно $\pi(500000)=41538$. Причём ровно $\pi(500000)-\pi(500000/2)=41538-22044=19494$ из них ровно в первой степени, $\pi(500000/2)-\pi(500000/3)=22044-15225=6819$ из них ровно в квадрате, и так далее вплоть до $\pi(\sqrt{500000})$, после которых степени простых будут увеличиваться сложнее (до $\pi(\sqrt[3]{500000})$, а потом до $\pi(\sqrt[4]{500000})$ и так далее вплоть до корня 19-й степени, больше не дающего простых).
Кстати программа подсчёта суммы всех степеней:

? n=0;forprime(p=2,500000,for(k=1,oo,if(p^k>500000,break);n+=500000\p^k));n
1786531

[B@R5uk]

mihaild в сообщении #1554077 писал(а):

B@R5uk, в смысле сумма степеней в разложении на простые множители? 1786531.

Да. Всего в три с половиной раза больше самого числа. Я думал ситуация хуже будет.

Я просто тут размышляю: нельзя ли обойтись без разложения вообще? Вся загвоздка с функцией Эйлера в том, что нам нужно представление аргумента в виде произведения небольших чисел, а не в виде частного двух гигантских: $A_n^k$ и $k!$. Ну так если памяти не жалко, то одно на другое можно поделить и без разложения (результат останется в виде произведения "маленьких" чисел). А потом найти в частном все простые сомножители и воспользоваться прекрасной формулой:

$$\varphi(x)=\frac{x}{p_1\cdot\ldots\cdot p_i}\times(p_1-1)\times\ldots\times(p_i-1)$$

Вот только трудоёмкость, я боюсь, хуже окажется, чем через разложение.

[Geen]

$n!$ разложить на простые сомножители несложно (алгоритм был дан выше), а вот $n!+1$, к примеру, практически нереально.

-- 07.05.2022, 20:31 --

B@R5uk в сообщении #1554080 писал(а):

воспользоваться прекрасной формулой

Эту формулу тоже обсуждали выше - деление по модулю.... несколько затратно

[B@R5uk]

Geen в сообщении #1554083 писал(а):

деление по модулю.... несколько затратно

Именно по этому я предлагаю обойтись без деления по модулю.

[mihaild]

B@R5uk в сообщении #1554080 писал(а):

Я думал ситуация хуже будет

Простое число $p$ входит в факториал со степенью $\sum\limits_{i=1}^\infty \lfloor \frac{n}{p^i}\rfloor$, что примерно равно $\frac{n}{p}$ при $p > 2$ и примерно равно $n$ при $p = 2$. Используя примерную оценку "есть $\frac{n}{\log n}$ простых чисел меньших $n$, и $k$-е простое число примерно равно $k \cdot \log k$", получаем, что сумма всех степеней оценивается сверху примерно как $n + \sum\limits_{k = 2}^{n / \log n} \frac{n}{k \log k} \sim n \cdot (\log \log n + 1)$, что даже неплохо согласуется с получившимся результатом.

(отношение суммы степеней в $n!$ к $n$, посчитанное точно и приближением выше)

[B@R5uk]

Вот такая порнография получается, если не использовать стандартное разложение на простые множители:

Используется синтаксис Java

import java.util.*;

public class Test_Binomial_3 {
    
    private static final int LIMIT = 500_000;
//    private static final int LIMIT = 200;
    private static final int [] mulsPower = new int [LIMIT + 1];
    private static final int [] divsPower = new int [LIMIT / 2 + 1];
    
    public static void main (String [] args) {
        int k, l, m, setSize, subSize, val, divPower;
        double logResult = 0.0;
        Random rand = new Random ();
        
        val = (int) (2 * Math .sqrt (LIMIT));
        setSize = LIMIT - rand .nextInt (val);
        subSize = rand .nextInt (setSize / 2 - val) + val + 1;
        //setSize = 50;
        //subSize = 22;
        //setSize = 43;
        //subSize = 20;
        //setSize = 49;
        //subSize = 16;
        System .out .println ("n = " + setSize + ", k = " + subSize + "\n");
        
        Arrays .fill (mulsPower, 0);
        for (k = setSize - subSize + 1; setSize >= k; ++k) {
            mulsPower [k] = 1;
        }
        Arrays .fill (divsPower, 0);
        for (k = 2; subSize >= k; ++k) {
            divsPower [k] = 1;
        }
        outerLoop:
        for (k = subSize; 2 <= k; --k) {
            //display (setSize, subSize);
            //System .out .println ("div = " + k);
            divPower = divsPower [k];
            if (0 == divPower) {
                continue;
            }
            val = Integer .min (divPower, mulsPower [k]);
            mulsPower [k] -= val;
            divPower      -= val;
            if (0 == divPower) {
                divsPower [k] = 0;
                continue;
            }
            for (l = 2; k >= l * l; ++l) {
                if (0 != k % l) {
                    continue;
                }
                m = k / l;
                val = Integer .min (divPower, mulsPower [m]);
                mulsPower [m] -= val;
                divsPower [l] += val;
                divPower      -= val;
                if (0 == divPower) {
                    divsPower [k] = 0;
                    continue outerLoop;
                }
                if (m == l) {
                    continue;
                }
                val = Integer .min (divPower, mulsPower [l]);
                mulsPower [l] -= val;
                divsPower [m] += val;
                divPower      -= val;
                if (0 == divPower) {
                    divsPower [k] = 0;
                    continue outerLoop;
                }
            }
            for (l = setSize / k, m = l * k; 2 <= l; --l, m -= k) {
                val = Integer .min (divPower, mulsPower [m]);
                mulsPower [m] -= val;
                mulsPower [l] += val;
                divPower      -= val;
                if (0 == divPower) {
                    divsPower [k] = 0;
                    continue outerLoop;
                }
            }
            val = Integer .min (divPower, mulsPower [k]);
            mulsPower [k] -= val;
            divPower      -= val;
            if (0 == divPower) {
                divsPower [k] = 0;
                continue;
            }
            System .out .print ("GCD for " + k);
            for (l = 4; setSize > l; ++l) {
                if (0 == mulsPower [l]) {
                    continue;
                }
                val = gcd (k, l);
                if (1 == val) {
                    continue;
                }
                val = Integer .min (val, k / val);
                System .out .println (" with " + l + " by " + val);
                divsPower [k] = divPower;
                for (l = k / val, m = k; 2 <= l; --l, m -= val) {
                    divPower = divsPower [m];
                    divsPower [m] = 0;
                    divsPower [l]   += divPower;
                    divsPower [val] += divPower;
                }
                continue outerLoop;
            }
            System .out .println (" - ERROR !");
        }
        //display (setSize, subSize);
        System .out .println ("\nDone.\n");
        
        for (k = 2, l = 0; setSize > k; ++k) {
            val = mulsPower [k];
            l += val;
            logResult += val * Math .log10 (k);
        }
        val = (int) Math .floor (logResult);
        System .out .println ("Result:                " + Math .pow (10, logResult - val) + " * 10 ^ " + val);
        System .out .println ("Number of multipliers: " + l);
    }
    
    @SuppressWarnings ("unused")
    private static void display (int setSize, int subSize) {
        int k;
        for (k = 2; subSize >= k; ++k) {
            System .out .println (String .format ("%6d%6d%6d", k, mulsPower [k], divsPower [k]));
        }
        for (; setSize >= k; ++k) {
            System .out .println (String .format ("%6d%6d",    k, mulsPower [k]));
        }
    }
    
    private static int gcd (int a, int b) {
        while (0 != b) {
            int c = b;
            b = a % b;
            a = c;
        }
        return a;
    }
}
 

Числитель $A_n^k$, представленный в виде множителей последовательных чисел в степенях, хранимых в переменной mulsPower, сокращается со знаменателем $k!$, степени которого хранятся в divsPower. Формально разложение на простые множители отсутствует, но оно неявно происходит, когда очередной обрабатываемых множитель делителя не сокращается ни с одним множителем делимого. Это означает, что этот делящий множитель составной, и, с помощью поиска наибольшего общего делителя (функция gcd), этот делящий множитель раскладывается на сомножители (а за одно и все остальные, меньшее его и имеющие тот же самый делитель).

Программа не находит функцию Эйлера, только представление её аргумента в виде "маленьких" множителей. Для завершения, нужно добавить таблицу простых чисел, поиск простых сомножителей и собственно расчёт по формуле

$$\varphi(x)=\frac{x}{p_1\cdot\ldots\cdot p_i}\times(p_1-1)\times\ldots\times(p_i-1)$$

Результат работы программы сейчас приблизительно выглядит так:

Используется синтаксис Text

n = 498691, k = 241303

GCD for 128690 with 257390 by 10
GCD for 128686 with 257390 by 2
GCD for 64347 with 257391 by 3
GCD for 20699 with 14 by 7
GCD for 13211 with 22 by 11
GCD for 11219 with 26 by 13
GCD for 7633 with 34 by 17
GCD for 6289 with 38 by 19
GCD for 5267 with 46 by 23
GCD for 4321 with 58 by 29
GCD for 4061 with 62 by 31

Done.

Result:                4.233477958222776 * 10 ^ 149999
Number of multipliers: 28221

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

[Geen]

B@R5uk в сообщении #1554262 писал(а):

Для завершения

Зачем всё это? Игнорируя некотрые "помарки", программа была приведена выше.

B@R5uk в сообщении #1554262 писал(а):

и собственно расчёт по формуле

Я так и не понял, как Вы собираетесь считать по этой формуле без деления?

-- 09.05.2022, 18:43 --

И да, время работы программы должно быть порядка 10 миллисекунд...

[B@R5uk]

Geen, как, как? Ну включите голову. Первый множитель в формуле (который "с делением") — это значок, обозначающий величину, а не то, как она считается. И я не говорит, что деление не нужно, я говорил, что не нужно деление по модулю.