/* ============================================================================
УНИВЕРСАЛЬНЫЙ ГЕНЕРАТОР CRT-УСЛОВИЙ v61 (gp2c-strict, исправлен Этап 2 для T>=48)
============================================================================ */
is_power_of_2(n) = {
if(n <= 0, return(0));
my(m = n);
while(m % 2 == 0, m = m \ 2);
return(m == 1);
}
analyze_chains_v61(res, k_max, verbose) = {
my(rows, row, L, T, x0, M, conds, tau_R_vec, total_trials, full_chains);
my(pos_hits, match_dist, matches, i, k, n, tau_val, emp_density, avg_indiv);
my(naive_theory, corr_est, r_total, r_p, pct, r_m, update_freq, expected_prob, p_i);
rows = matsize(res)[1];
if(rows == 0 || res[1,4] == 0,
if(verbose >= 1, print("Пустой пул."));
return(0)
);
for(row = 1, rows,
if(res[row,5] <= 0, next);
L = res[row,1]; T = res[row,2]; x0 = res[row,3]; M = res[row,4];
conds = vector(L, j, res[row, 5+j]);
tau_R_vec = vector(L, j, T \ numdiv(conds[j]));
total_trials = 0; full_chains = 0;
pos_hits = vector(L); match_dist = vector(L+1);
update_freq = max(1, (k_max + 1) \ 100);
for(k = 0, k_max,
n = x0 + k * M;
total_trials = total_trials + 1;
matches = 0;
for(i = 0, L-1,
tau_val = numdiv(n + i);
if(tau_val == T,
matches = matches + 1;
pos_hits[i+1] = pos_hits[i+1] + 1
)
);
match_dist[matches+1] = match_dist[matches+1] + 1;
if(matches == L, full_chains = full_chains + 1);
if(verbose >= 1 && k % update_freq == 0,
printf("\r[Анализ системы %d: %3d%%]", row, 100 * k \ (k_max + 1));
);
); if(verbose >= 1, print(""));
if(verbose >= 1,
print("\n========================================");
printf("СИСТЕМА #%d | x0=%d | M=%d | logM=%.2f\n", row, x0, M, log(M));
printf("tau_R = %s\n", Str(tau_R_vec));
r_total = real(total_trials);
emp_density = if(total_trials > 0, real(full_chains) / r_total, 0.0);
printf("Проверено: %d | Полных: %d | Плотность: %.3e\n", total_trials, full_chains, emp_density);
print("Частоты по позициям (τ=T):");
for(i = 1, L,
if(verbose >= 1,
r_p = real(pos_hits[i]);
pct = r_p / r_total * 100.0;
printf(" pos %d (τ_R=%d): %5.1f%%\n", i-1, tau_R_vec[i], pct)
)
);
expected_prob = 1.0;
print("Расчёт ожидаемой вероятности полной цепочки:");
for(i = 1, L,
if(r_total > 0, p_i = real(pos_hits[i]) / r_total, p_i = 0.0);
expected_prob = expected_prob * p_i;
);
printf(" Ожидаемая (prod pos): %.4e\n", expected_prob);
printf(" Эмпирическая (full/trials): %.4e\n", emp_density);
if(expected_prob > 1e-50 && emp_density > 0,
printf(" Коэффициент корреляции CRT: %.2f\n", emp_density / expected_prob);
, print(" Коэффициент корреляции CRT: N/A");
);
if(full_chains == 0, print("⚠ Нет цепочек."))
)
);
return(1)
}
generate_crt_pool_v61(L, T, max_logLcm, max_mod, K, seed, strict_composite, max_tau_R, verbose) = {
my(pool_conds, pool_lcm, pool_scores, primes, num_primes);
my(i, j, k, p_idx, r, conflicts, best_cnt, best_off_cnt, best_offs, best_off);
my(i_pos, target_pos, pos_count, p_big, e, needed, new_val, new_lcm, new_tau, new_tau_R, found, p_used, skip, worst_bad);
my(x, M_crt, crt_ok, g, mi, ai, score_sum, tau_R_i, cur_tau, result, tmp, c, tau_R, all_ok, w, logM_factor, is_ok, p_pow, a);
my(unique_cnt, is_dup, tmp_M, tmp_x0, base_i, step_i, trap_tau, reject_reason, budget_phase2, max_a, max_a_2, max_a_3);
my(best_e, score, best_score, v2_global, v2_mi);
if(!strict_composite, strict_composite = 1);
if(!max_tau_R, max_tau_R = 8);
if(max_tau_R > T \ 2, max_tau_R = T \ 2); while(max_tau_R > 1 && !is_power_of_2(max_tau_R), max_tau_R = max_tau_R - 1);
if(max_tau_R < 2, max_tau_R = 2);
/* Паритетная безопасность: 2^5 только для T<48 */
max_a_2 = if(T < 48, 5, 3);
max_a_3 = 3;
budget_phase2 = max_logLcm * 0.85;
num_primes = primepi(max(3000, L+500));
primes = primes(num_primes);
setrand(seed);
pool_conds = matrix(K, L, i, j, 1);
pool_lcm = vector(K, j, 1);
pool_scores = vector(K, j, 0.0);
best_offs = vector(60);
if(verbose >= 1, printf("v61: L=%d, T=%d, max_logLcm=%.1f, K=%d, strict=%d, max_tau_R=%d, 2^a<= %d\n", L, T, max_logLcm, K, strict_composite, max_tau_R, max_a_2));
/* Этап 2: Малые простые (ИСПРАВЛЕНО: проверка только на делимость T) */
if(verbose >= 1, print("Этап 2: Малые простые"));
for(p_idx = 1, num_primes,
if(primes[p_idx] > min(L, 60), break);
p = primes[p_idx];
for(k = 1, K,
best_off_cnt = 0; best_cnt = L+1;
for(r = 0, p-1,
conflicts = 0;
forstep(i_pos = r, L-1, p,
cur_tau = numdiv(pool_conds[k,i_pos+1]*p);
/* Убрана проверка is_power_of_2, разрешены промежуточные tau_R */
if(T % cur_tau != 0, conflicts = conflicts + 1);
);
if(conflicts < best_cnt, best_cnt = conflicts);
if(conflicts < best_cnt, best_off_cnt = 0);
if(conflicts == best_cnt, best_off_cnt = best_off_cnt + 1);
if(conflicts == best_cnt, best_offs[best_off_cnt] = r);
);
if(best_off_cnt == 0, next);
best_off = best_offs[random(best_off_cnt)+1];
forstep(i_pos = best_off, L-1, p,
cur_tau = numdiv(pool_conds[k,i_pos+1]*p);
tau_R = T \ cur_tau;
is_ok = (T % cur_tau == 0);
if(strict_composite && tau_R == 2, is_ok = 0);
if(is_ok,
pool_conds[k,i_pos+1] = pool_conds[k,i_pos+1] * p;
pool_lcm[k] = lcm(pool_lcm[k], p);
);
); pos_count = 0;
forstep(i_pos = best_off, L-1, p, pos_count = pos_count + 1);
if(pos_count == 1, target_pos = best_off + 1,
pos_count > 1, target_pos = best_off + p * (pos_count \ 2) + 1,
target_pos = 0);
if(target_pos > 0,
if(p == 2, max_a = max_a_2,
if(p == 3, max_a = max_a_3,
max_a = 2));
p_pow = 1;
for(a = 1, max_a,
p_pow = p_pow * p;
new_val = pool_conds[k,target_pos] * p_pow;
new_tau = numdiv(new_val);
new_tau_R = T \ new_tau;
if((max_mod == 0 || new_val <= max_mod) && (T % new_tau == 0),
if(!strict_composite || new_tau_R != 2,
new_lcm = lcm(pool_lcm[k], p_pow);
if(log(new_lcm) <= budget_phase2,
pool_conds[k,target_pos] = new_val;
pool_lcm[k] = new_lcm;
);
);
);
);
);
);
);
/* Этап 3: Добивка (жестко требует степени двойки) */
if(verbose >= 1, print("Этап 3: Добивка"));
for(k = 1, K,
setrand(seed + k * 7919);
while(log(pool_lcm[k]) < max_logLcm,
target_pos = 0; worst_bad = 0;
for(i = 1, L,
cur_tau = numdiv(pool_conds[k,i]);
if(T % cur_tau == 0,
tau_R = T\cur_tau;
if(!is_power_of_2(tau_R) || tau_R > max_tau_R,
if(tau_R > worst_bad, worst_bad = tau_R; target_pos = i);
);
);
);
if(target_pos == 0, break);
needed = T \ numdiv(pool_conds[k, target_pos]);
found = 0;
skip = random(2);
p_big = nextprime(L+1);
while(skip > 0, p_big = nextprime(p_big+1); skip = skip - 1);
while(p_big < 50000 && found == 0,
p_used = 0;
for(j = 1, L, if(pool_conds[k,j] % p_big == 0, p_used = 1));
if(p_used == 0,
best_e = 0; best_score = -1e9;
for(e = 1, 3,
if(needed % (e+1) == 0,
new_val = pool_conds[k,target_pos] * p_big^e;
new_lcm = lcm(pool_lcm[k], p_big^e);
new_tau = numdiv(new_val);
new_tau_R = T \ new_tau;
if((max_mod == 0 || new_val <= max_mod) && log(new_lcm) <= max_logLcm,
score = 0.0;
if(is_power_of_2(new_tau_R) && new_tau_R <= max_tau_R, score = score + 5000.0);
if(new_tau_R == 4, score = score + 2000.0);
if(new_tau_R == 8, score = score + 800.0);
if(new_tau_R == 16, score = score + 300.0);
score = score - real(log(new_lcm)) * 10.0;
if(score > best_score, best_score = score; best_e = e);
);
);
);
if(best_e > 0,
e = best_e;
new_val = pool_conds[k,target_pos] * p_big^e;
new_lcm = lcm(pool_lcm[k], p_big^e);
new_tau = numdiv(new_val);
new_tau_R = T \ new_tau;
pool_conds[k,target_pos] = new_val;
pool_lcm[k] = new_lcm;
found = 1;
);
);
if(found == 0, p_big = nextprime(p_big+1));
);
if(found == 0, pool_scores[k] = -1.0; break);
);
);
/* Этап 4: CRT + Жесткий фильтр ловушек + Бюджет-скор */
if(verbose >= 1, print("Этап 4: CRT + стабильность + фильтр"));
result = matrix(K, 5+L);
for(k = 1, K,
if(pool_scores[k] < 0, next);
reject_reason = "";
x = Mod(0,1); M_crt = 1; crt_ok = 1; all_ok = 1; score_sum = 0.0;
for(i = 1, L,
mi = pool_conds[k,i];
if(mi == 1, next); ai = Mod(-i+1, mi); g = gcd(M_crt, mi);
if(g > 1, if((lift(x)-lift(ai))%g != 0, crt_ok = 0; reject_reason = "CRT incompatible"));
if(crt_ok == 0, break);
x = chinese(x, ai); M_crt = lcm(M_crt, mi));
if(crt_ok == 0 || M_crt == 0,
pool_scores[k] = -1.0; next;
);
/* Жесткая проверка паритетных ловушек */
if(T >= 36,
v2_global = valuation(M_crt, 2);
for(i = 1, L,
mi = pool_conds[k,i];
if(mi % 2 == 0,
v2_mi = valuation(mi, 2);
if(v2_mi >= 5 && v2_global == v2_mi,
reject_reason = "parity trap (v2_Mi == v2_LCM)";
all_ok = 0; break;
);
);
);
);
if(all_ok == 0, pool_scores[k] = -1.0; next);
/* Финальный скоринг (только степени двойки) */
for(i = 1, L,
cur_tau = numdiv(pool_conds[k,i]);
if(T % cur_tau != 0, all_ok = 0; break);
if(cur_tau >= T, all_ok = 0; break);
tau_R_i = T\cur_tau;
if(!is_power_of_2(tau_R_i) || tau_R_i > max_tau_R, all_ok = 0; break);
w = 0.0;
if(tau_R_i == 4, w = 1.0);
if(tau_R_i == 8, w = 0.6);
if(tau_R_i == 16, w = 0.25);
score_sum = score_sum + w;
);
if(all_ok == 0, pool_scores[k] = -1.0; next);
if(M_crt > 1,
pool_scores[k] = real(score_sum) / real(log(M_crt));
if(log(M_crt) > 45.0, pool_scores[k] = pool_scores[k] * (1.0 - 0.02*(log(M_crt)-45.0)));
, pool_scores[k] = 0.0);
result[k,1] = L; result[k,2] = T; result[k,3] = lift(x); result[k,4] = M_crt; result[k,5] = pool_scores[k];
for(i = 1, L, result[k,5+i] = pool_conds[k,i]));
/* Этап 5: Сортировка и дедуп */
if(verbose >= 1, print("Этап 5: Сортировка и дедуп")); for(i = 1, K-1, for(j = i+1, K,
if(result[j,5] > result[i,5],
tmp = vector(5+L);
for(c = 1, 5+L, tmp[c] = result[i,c]);
for(c = 1, 5+L, result[i,c] = result[j,c]);
for(c = 1, 5+L, result[j,c] = tmp[c]))));
if(verbose >= 1, print("Топ-5 (tau_R):"));
for(k = 1, 5,
if(result[k,5] > 0,
printf(" #%d: M=%-12d logM=%.2f tau_R=%s\n", k, result[k,4], log(result[k,4]), Str(vector(L, i, T\numdiv(result[k,5+i]))))));
unique_cnt = 0;
for(i = 1, K,
if(result[i,5] > 0,
is_dup = 0; tmp_M = result[i,4]; tmp_x0 = result[i,3];
for(j = 1, unique_cnt,
if(result[j,4] == tmp_M && result[j,3] == tmp_x0, is_dup = 1; break));
if(is_dup == 0,
unique_cnt = unique_cnt + 1;
for(c = 1, 5+L, result[unique_cnt, c] = result[i, c])));
);
for(i = unique_cnt+1, K, for(c = 1, 5+L, result[i, c] = 0));
if(verbose >= 1, printf("Уникальных систем: %d из %d\n", unique_cnt, K));
return(result)
}
search_chains_fast_v61(res, target, max_k, verbose) = {
my(rows, found, k, n, i, ok, L, T, x0, M, checked, batch, row);
rows = matsize(res)[1];
found = 0; checked = 0;
if(max_k > 0, batch = max_k \ 20, batch = 1000);
if(batch < 1000, batch = 1000);
for(row = 1, rows,
if(res[row,5] <= 0, next);
L = res[row,1]; T = res[row,2]; x0 = res[row,3]; M = res[row,4];
if(verbose >= 1, printf("Система %d (score=%.4f): ", row, res[row,5]));
for(k = 0, max_k,
checked = checked + 1;
n = x0 + k*M; ok = 1;
for(i = 0, L-1,
if(numdiv(n+i) != T, ok = 0);
if(ok == 0, break);
);
if(ok,
found = found + 1;
if(verbose >= 1, printf("\n [✓ #%d] k=%7d n=%d", found, k, n));
if(found >= target, return(found));
);
if(verbose >= 1 && checked % batch == 0, printf("."));
);
if(verbose >= 1, print("")); );
if(verbose >= 1, printf("Завершено. Найдено: %d из %d\n", found, target));
return(found)
}
Я её ещё не запускал, кстати. Та с которой хорошие паттерны выпали для D(48,6), была 59 версия, а это 63-я, исправленная.
Тот который получше (где частоты получше) это был вот такой:
