Dmitriy40Смотрю исходный код forprime
Там крайне заморочено всё. Разрабы очень старались, чтобы было безопасно (пока возможно, используется тип
ulong, потом переключается в
GEN) , чтобы таблица по которой итерируется (или её куски) помещалась в L2 кэш то есть под
forprime() решетом готовится своя таблица простых по которым он потом идёт, и все такое сильно запутанное.
А вот в движке факторизации, когда используется тот же
factor(n,D) с ограничением по множителям, таких заморочек нет - там они просто ходят по готовой primetab, но это непубличная таблица и как я понял, сделано так специально.
Тоже самое по
ispseudoprime(). В движке факторизации есть функция
ifac_isprime() которая используется для проверки полученных делителей и частного на простоту, но к этому моменту уже ясно что маленькие делители отсутсвуют. А "пользовательская"
ispseudoprime() запутанней, т.к. в ней может ещё быть флаг, а в тестируемом числе могут быть маленькие делители и т.п..
Думаю, что было бы вполне безопасно использовать
primetab и
ifac_isprime(), но моих знаний не хватает пока как сделать это в пользовательской функции, т.е. довести до ума результат трансляции
gp2c не влезая (не перекомпилируя) в ядро.
Ранее сделанная
omega_upto() как раз ходит по
primetab и использует
ifac_isprime(0)...
YadryaraЯ кстати выжал из ИИ ещё одну функцию,
omega_range() которая считает количество различных простых множителей не с двойки до лимита, а в заданном диапазоне, но пока не уверен что она вам нужна
/* Returns:
* k = number of distinct prime factors of n in [l, u], if k <= f
* and n is squarefree,
* 0 otherwise.
*
* Arguments:
* n: positive integer (t_INT)
* f: maximum allowed number of distinct primes in [l, u] (t_INT)
* l: lower bound for primes (t_INT, >= 2)
* u: upper bound for primes (t_INT, >= l)
*/
GEN
omega_range(GEN n, GEN f, GEN l, GEN u)
{
pari_sp av = avma;
long F, nb = 0;
ulong L, U;
if (typ(n) != t_INT || signe(n) <= 0) pari_err_TYPE("omega_range", n);
if (typ(f) != t_INT || signe(f) <= 0) pari_err_TYPE("omega_range", f);
if (typ(l) != t_INT || typ(u) != t_INT) pari_err_TYPE("omega_range", l);
if (cmpii(l, u) > 0 || cmpiu(l, 1) <= 0) pari_err_TYPE("omega_range", l);
if (cmpiu(n, 1) <= 0) return gen_0;
F = itos(f);
if (F <= 0) return gen_0;
L = itou(l);
U = itou(u);
if (L > U || L < 2) return gen_0;
n = absi(n);
/* Обработка 2, если в диапазоне */
if (L <= 2 && 2 <= U)
{
if (!mod2(n))
{
long v = vali(n);
if (v > 1) { set_avma(av); return gen_0; }
nb++;
if (nb > F) { set_avma(av); return gen_0; }
n = shifti(n, -v);
if (is_pm1(n)) goto done;
}
}
/* Trial division по primetab */
{
long i, lim = lg(primetab);
for (i = 1; i < lim; i++)
{
ulong p = itou(gel(primetab, i));
if (p > U) break;
if (p < L) continue;
if (dvdii(n, utoipos(p)))
{
long v = Z_lvalrem(n, p, &n);
if (v > 1) { set_avma(av); return gen_0; }
nb++;
if (nb > F) { set_avma(av); return gen_0; }
if (is_pm1(n)) goto done;
}
}
}
/* Проверка остатка: если он простой и в [l, u] — учитываем */
if (!is_pm1(n))
{
if (cmpiu(n, L) >= 0 && cmpiu(n, U) <= 0)
{
if (ifac_isprime(n))
{
nb++;
if (nb > F) { set_avma(av); return gen_0; }
goto done;
}
}
}
done:
set_avma(av);
return (nb >= 1) ? stoi(nb) : gen_0;
}
04.12.2025, 04:38 
.
