/* Returns:
* k = omega(n) if omega(n) <= m,
* 0 if omega(n) > m,
* 0 if s == 0 and n is not squarefree.
*
* Arguments:
* n: positive integer (t_INT)
* m: maximum number of distinct prime factors allowed (t_INT)
* s: squarefree flag (t_INT): s == 0 → require squarefree; s != 0 → ignore exponents
*/
GEN
omega_upto(GEN n, GEN m, GEN s)
{
pari_sp av = avma;
long M, nb = 0;
long squarefree = (typ(s) == t_INT && signe(s) == 0); /* s == 0 */
ulong p, B;
if (typ(n) != t_INT || signe(n) <= 0) pari_err_TYPE("omega_upto", n);
if (typ(m) != t_INT || signe(m) <= 0) pari_err_TYPE("omega_upto", m);
if (cmpiu(n, 1) <= 0) return gen_0;
M = itos(m);
if (M <= 0) return gen_0;
n = absi(n); /* работает с |n| и возвращает новый объект */
/* --- Этап 1: обработка малых простых (2, 3, 5, 7) --- */
/* 2 */
if (!mod2(n)) /* если чётное */
{
long v = vali(n);
if (squarefree && v > 1) { set_avma(av); return gen_0; }
nb++;
if (nb > M) { set_avma(av); return gen_0; }
n = shifti(n, -v); /* n /= 2^v */
if (is_pm1(n)) goto done;
}
/* 3, 5, 7 */
{
const ulong small_primes[] = {3, 5, 7};
long j;
for (j = 0; j < 3; j++)
{
p = small_primes[j];
if (dvdii(n, utoipos(p)))
{
long v = Z_lvalrem(n, p, &n);
if (squarefree && v > 1) { set_avma(av); return gen_0; }
nb++;
if (nb > M) { set_avma(av); return gen_0; }
if (is_pm1(n)) goto done;
}
}
}
/* Граница trial division как в factor() */
B = tridiv_bound(n);
/* Trial division по глобальной primetab */
{
long i, l = lg(primetab);
for (i = 1; i < l; i++)
{
ulong p_val = itou(gel(primetab, i));
if (p_val > B) break;
if (p_val <= 7) continue; /* уже обработаны */
if (dvdii(n, utoipos(p_val)))
{
long v = Z_lvalrem(n, p_val, &n);
if (squarefree && v > 1) {
set_avma(av);
return gen_0;
}
nb++;
if (nb > M) {
set_avma(av);
return gen_0;
}
if (is_pm1(n)) goto done;
}
}
}
/* --- Этап 2: остаток --- */
if (is_pm1(n)) goto done;
/* Если остаток простой — учтём его */
if (ifac_isprime(n))
{
nb++;
if (nb > M) { set_avma(av); return gen_0; }
goto done;
}
/* Если уже найдено M делителей, а остаток >1 → будет >M */
if (nb >= M)
{
set_avma(av);
return gen_0;
}
/* --- Этап 3: ifac_start только для составного остатка --- */
{
GEN part = ifac_start(n, squarefree ? 1 : 0);
GEN here;
for (;;)
{
here = ifac_main(&part);
if (!here) break; /* факторизация завершена */
if (here == gen_0)
{
set_avma(av);
return gen_0; /* квадрат найден в Moebius-режиме */
}
nb++;
if (nb > M)
{
set_avma(av);
return gen_0;
}
ifac_delete(here);
}
}
done:
set_avma(av);
return (nb >= 1) ? stoi(nb) : gen_0;
}