Чтоб видеть сходимость и вклад каждой константы, последовательно считаю интегралы со всё большим количеством констант.
Погодите, вы ведь при этом заменяете интеграл суммы на сумму интегралов? Извините, если что...
-- добавлено через 56 минут --Новая версия программы.
Добавил опцию --bias
Можно устанавливать наклон второго прохода этой опцией.
При --bias 0 будет фактически равномерный проход с удвоенной выборкой.
Изменения по выводимой инфе - убрал непонятное.
Дисперсии первого прогона считаются по другому (должно быть немного точнее).
Ссылка
https://disk.yandex.ru/d/BLYFw4hJuTc2PQТекст справки
Код:
hl1_mc5_1 — Sequential Importance Sampling для оценки коэффициентов A[j].
ОПИСАНИЕ:
Коэффициенты A[j] входят в формулу Харди–Литтлвуда для количества
простых кортежей заданного диаметра d в интервале [1, x]: оно
асимптотически равно Σⱼ A[j] · ∫₂ˣ dt / (ln t)^{|H₀|+j}, где |H₀| —
длина базового паттерна, а j — число дополнительных простых позиций
(«загрязнений») внутри того же диаметра. Каждый коэффициент A[j] равен
сумме локальных констант по всем допустимым конфигурациям с j
загрязнениями и определяет совокупный вклад кортежей длины |H₀|+j
в эту сумму.
Программа использует Sequential Importance Sampling с двумя прогонами:
- Прогон 1 (смещение=0): Равномерное сэмплирование, точное для малых j
- Прогон 2 (смещение=bias): Importance Sampling, лучшее для больших j
Результаты объединяются через обратные дисперсии.
ИСПОЛЬЗОВАНИЕ:
./hl1_mc5 [опции] MAXN h1 h2 ... hk
АРГУМЕНТЫ:
MAXN Максимальный индекс j для вывода (A[0]..A[MAXN]).
h1 h2 ... hk Элементы множества H₀ (обязательные позиции кортежа).
Диаметр d = max(h_i). H₀ должно содержать 0 и d.
ОПЦИИ:
--size N Количество сэмплов на проход (K, M, G, T). По умолчанию: 5M
--bias X Параметр наклона второго прогона (по умолчанию: 2.0).
Увеличение смещает выборку в сторону больших s,
уменьшение делает распределение весов равномернее.
--full Выводить числа целиком (%.0f)
--plain Вывод без квадратных скобок (для парсинга)
--nostdev Не печатать массив относительных погрешностей (СКО)
--help Показать эту справку
ПРИМЕРЫ:
./hl1_mc5 --size 10M 99 0 210 420 # d=420, 10M сэмплов на проход, bias=2.0
./hl1_mc5 --size 1M --bias 1.0 30 0 210 420 630 840 1050 1260 # мягкий наклон
-- добавлено через 37 минут --Можно вызвать прямо из pari
Функция:
Код:
hl1(maxn, H0, size="1M", bias=2.0) =
{
my(cmd, lines, h0_str, i);
h0_str = Str(H0[1]);
for(i = 2, #H0, h0_str = Str(h0_str, " ", H0[i]));
cmd = Str("./hl1_mc5 --size ", size, " --bias ", bias,
" ", maxn, " ", h0_str);
lines = externstr(cmd);
[eval(lines[1]), eval(lines[2])]
}
Запуск:
Код:
? [A, ERR] = hl1(90, [0, 210, 420], "1M", 2.0);
Кортеж: диаметр 420, обязательных 3, возможных загрязнений 208, запрещенных позиций 0, простых 45
Счёт: выборка на проход 1M, смещение 2.0, потоков доступно 7, расчёт до j=90
Прогон 1. Смещение 0.0
[100%] Прошло: 00:00:13 Осталось: 00:00:00 потоков используется: 1
Прогон 2. Смещение 2.0
[100%] Прошло: 00:00:14 Осталось: 00:00:00 потоков используется: 1
Результат: 59 достоверных коэффициентов (err<100%)
? #A
%58 = 59
? A[1]
%59 = 1.5034947902255200000000000000000000000 E80
?
Сообщения и прогресс "пробросятся" в терминал pari (чтобы видеть как дела и когда закончится), а результаты попадут в переменные.
