Otimização de Custos em Permutações, Fatoração de Polinômios Ciclótomicos e Cobertura Mínima em Árvores

Problema 1: Otimização de Custos em Permutações

Considere um jogo onde o custo de uma permutação é determinado pelo seu número de inversões. Uma operação de reordenação está disponível, que incorre em um custo fixo c e altera o custo da permutação para o custo médio ponderado X de todas as permutações possíveis. O objetivo é minimizar o custo total. Primeiramente, é crucial observar que a resposta para uma permutação depende unicamente do seu número de inversões.

Se nenhuma operação de reordenação for realizada, podemos calcular o custo para cada número possível de inversões, danotado por fi. Se a operação de reordenação for utilizada, o custo passa a ser a média ponderada X de todos os fi, mais o custo c. Logicamente, a reordenação será aplicada às permutações que originalmente teriam um fi maior do que X + c. Ao iterar sobre possíveis valores para X + c, transformamos todos os fi que excedem este valor para X + c. A partir disso, um novo valor médio ponderado X' pode ser calculado. Resolvendo X' = X, encontramos o valor ótimo de X e, em seguida, verificamos se X + c está dentro do intervalo considerado. Além disso, é necessário calcular o número de permutações de comprimento n que possuem k inversões, o que pode ser eficientemente determinado via programação dinâmica.

Exemplo de Implementação para Permutações

#include <iostream>
#include <algorithm>
#include <vector>
#include <numeric> // Para std::gcd
#include <string.h> // Para memset

// Usando __int128 para lidar com valores grandes em cálculos de fração.
// Esta é uma extensão GNU GCC e pode não ser universalmente suportada.
typedef __int128_t int128_t;

// Função GCD para int128_t
int128_t gcd(int128_t a, int128_t b) {
    while (b) {
        a %= b;
        std::swap(a, b);
    }
    return a;
}

// Estrutura para representar frações
struct Fracao {
    int128_t numerador, denominador;

    Fracao(int128_t num = 0, int128_t den = 1) : numerador(num), denominador(den) {
        simplificar();
    }

    void simplificar() {
        if (denominador == 0) return; // Erro: divisão por zero
        int128_t divisor_comum = gcd(numerador, denominador);
        numerador /= divisor_comum;
        denominador /= divisor_comum;
        if (denominador < 0) { // Garantir denominador positivo
            numerador = -numerador;
            denominador = -denominador;
        }
    }

    // Operador de adição
    Fracao operator+(const Fracao& outra) const {
        int128_t novo_denominador = (denominador / gcd(denominador, outra.denominador)) * outra.denominador;
        int128_t novo_numerador = numerador * (novo_denominador / denominador) + outra.numerador * (novo_denominador / outra.denominador);
        return Fracao(novo_numerador, novo_denominador);
    }

    // Adição com inteiro
    Fracao operator+(int128_t val_inteiro) const {
        return Fracao(numerador + val_inteiro * denominador, denominador);
    }

    // Multiplicação
    Fracao operator*(const Fracao& outra) const {
        return Fracao(numerador * outra.numerador, denominador * outra.denominador);
    }

    // Multiplicação com inteiro
    Fracao operator*(int128_t val_inteiro) const {
        return Fracao(numerador * val_inteiro, denominador);
    }

    // Função para comparar Fração com inteiro (min)
    static Fracao min_val(int128_t val_inteiro, const Fracao& fracao_val) {
        Fracao int_como_fracao(val_inteiro * fracao_val.denominador, fracao_val.denominador);
        return (int_como_fracao.numerador < fracao_val.numerador) ? int_como_fracao : fracao_val;
    }
};

// Sobrecarga do operador << para imprimir int128_t e Fracao
std::ostream& operator<<(std::ostream& out, int128_t x) {
    if (x == 0) return out << "0";
    std::string s = "";
    bool negativo = false;
    if (x < 0) {
        negativo = true;
        x = -x;
    }
    while (x > 0) {
        s += (x % 10) + '0';
        x /= 10;
    }
    if (negativo) s += '-';
    std::reverse(s.begin(), s.end());
    return out << s;
}

std::ostream& operator<<(std::ostream& out, const Fracao& f) {
    out << f.numerador;
    if (f.denominador != 1) {
        out << "/" << f.denominador;
    }
    return out;
}


int128_t dp_permutacao_inversoes[305]; // dp[k] = custo para k inversões
int128_t contagens_inversoes[20][305]; // contagens_inversoes[n][k] = número de permutações de tamanho n com k inversões
int128_t sequencia_atual[305]; // A[i] na entrada

// Estrutura para armazenar pares (custo, num_inversoes)
struct CustoInversoes {
    int128_t custo;
    int num_inversoes;

    bool operator<(const CustoInversoes& outra) const {
        return custo < outra.custo;
    }
};

int main() {
    std::ios_base::sync_with_stdio(false);
    std::cin.tie(NULL);

    // Pré-calcular o número de permutações com k inversões
    // contagens_inversoes[n][k] = número de permutações de tamanho n com k inversões
    contagens_inversoes[0][0] = 1; // 0 elementos, 0 inversões, 1 permutação (a vazia)
    for (int n_elem = 0; n_elem < 16; ++n_elem) {
        for (int k_inv = 0; k_inv <= n_elem * (n_elem - 1) / 2; ++k_inv) {
            if (contagens_inversoes[n_elem][k_inv] == 0) continue;
            for (int novo_elem_pos = 0; novo_elem_pos <= n_elem; ++novo_elem_pos) {
                // Adicionar o (n_elem+1)-ésimo elemento pode criar 0 a n_elem novas inversões
                contagens_inversoes[n_elem + 1][k_inv + novo_elem_pos] += contagens_inversoes[n_elem][k_inv];
            }
        }
    }

    int num_testes;
    std::cin >> num_testes;
    while (num_testes--) {
        long long n_ll, a_ll, b_ll, c_ll, d_ll;
        std::cin >> n_ll >> a_ll >> b_ll >> c_ll >> d_ll;
        int128_t n = n_ll;
        int128_t custo_a = a_ll;
        int128_t custo_b = b_ll;
        int128_t custo_c = c_ll;
        int128_t num_consultas = d_ll;

        int128_t max_inversoes = n * (n - 1) / 2;

        // Calcular dp_permutacao_inversoes[k] (custo mínimo para k inversões sem reordenação)
        // O custo da i-ésima inversão é a
        // Se i > N - i, podemos inverter a ordem e ter N - i inversões, com custo b
        // Isso parece ser uma variação de custo para alcançar um certo número de inversões
        // ou uma estratégia para reduzir o número de inversões por um custo fixo.
        memset(dp_permutacao_inversoes, 0, sizeof(dp_permutacao_inversoes));
        for (int i = 1; i <= max_inversoes; ++i) {
            if (i <= max_inversoes - i) {
                dp_permutacao_inversoes[i] = dp_permutacao_inversoes[i - 1] + custo_a;
            } else {
                dp_permutacao_inversoes[i] = Fracao::min_val(dp_permutacao_inversoes[i - 1] + custo_a, dp_permutacao_inversoes[max_inversoes - i] + custo_b).numerador; // A min_val compara int com Fracao
            }
        }
        
        std::vector<CustoInversoes> custos_por_inversoes(max_inversoes + 1);
        for (int i = 0; i <= max_inversoes; ++i) {
            custos_por_inversoes[i] = {dp_permutacao_inversoes[i], i};
        }
        std::sort(custos_por_inversoes.begin(), custos_por_inversoes.end());

        Fracao valor_X_otimo(std::numeric_limits<int128_t>::max(), 1); // Inicializa com um valor muito grande
        
        int128_t fatorial_n = 1;
        for (int i = 1; i <= n; ++i) fatorial_n *= i;

        int128_t soma_ponderada_custos = 0; // Somatório de f[n][p[i].second] * p[i].first
        int128_t soma_contagens = 0;      // Somatório de f[n][p[i].second]

        // Encontrar o X ótimo iterando sobre as possibilidades de cortes
        // sp: soma das contagens de permutações
        // sd: soma das contagens * custo (f_i)
        // Isso encontra o ponto de corte ideal para aplicar a reordenação
        // onde X + c é menor do que o custo original f_i.
        for (int i = 0; i <= max_inversoes; ++i) {
            soma_contagens += contagens_inversoes[n][custos_por_inversoes[i].num_inversoes];
            soma_ponderada_custos += contagens_inversoes[n][custos_por_inversoes[i].num_inversoes] * custos_por_inversoes[i].custo;

            // Condição para X+c estar no intervalo [custos_por_inversoes[i].custo, custos_por_inversoes[i+1].custo]
            // A ideia é que X + c é o valor limite. Se f_k > X + c, ele se torna X + c.
            // A equação para X é: (soma_dos_f_i_originais + (Fatorial_N - soma_contagens_ate_corte) * (X+c)) / Fatorial_N = X
            // Ou rearranjando: X = (soma_ponderada_custos + custo_c * (Fatorial_N - soma_contagens)) / soma_contagens_total
            // Aqui, o corte é p[i].first, e estamos testando X + c = p[i].first
            // (soma_ponderada_custos + (fatorial_n - soma_contagens) * custos_por_inversoes[i].custo + custo_c * soma_contagens) / fatorial_n = X
            // Na verdade, a lógica é: X é a média dos custos (onde alguns podem ser "tampados" por (X+c)).
            // X = (soma_dos_custos_originais_abaixo_do_corte + soma_dos_f_i_acima_do_corte_com_limite_X_mais_c) / Fatorial_N
            // X = (sd_abaixo_corte + (Fatorial_N - sp_abaixo_corte) * (X + c)) / Fatorial_N
            // Fatorial_N * X = sd_abaixo_corte + (Fatorial_N - sp_abaixo_corte) * X + (Fatorial_N - sp_abaixo_corte) * c
            // X * (Fatorial_N - (Fatorial_N - sp_abaixo_corte)) = sd_abaixo_corte + (Fatorial_N - sp_abaixo_corte) * c
            // X * sp_abaixo_corte = sd_abaixo_corte + (Fatorial_N - sp_abaixo_corte) * c
            // X = (sd_abaixo_corte + (Fatorial_N - sp_abaixo_corte) * c) / sp_abaixo_corte
            // O código original é: X = (Fraction) { sd + c * (fac - sp), sp } + c;
            // Parece que o X na verdade é o valor que se compara a F_i, ou seja, o ponto de corte.
            // Para encontrar o ponto de corte Y = X + C, precisamos que Y seja tal que os valores menores que Y não mudam, e os maiores viram Y.
            // X = (sum_{f_k <= Y} f_k * count_k + sum_{f_k > Y} Y * count_k) / Fatorial_N
            // Y é o candidato a X + c. O loop itera sobre p[i].first (que é um f_k) como candidato a Y.
            // Seja Y = p[i].first.
            // soma_ponderada_custos é sum_{k=0 to i} f_k * count_k
            // soma_contagens é sum_{k=0 to i} count_k
            // O resto das permutações (Fatorial_N - soma_contagens) terão seu custo limitado a Y (ou p[i+1].first, depende da interpretação).
            // O original parece estar calculando o X efetivo quando o corte está entre p[i].first e p[i+1].first
            // Se o corte for K, X = (sum(fk*countk for fk < K) + sum(K*countk for fk >= K)) / Fatorial_N
            // X_avg = (soma_ponderada_custos + (Fatorial_N - soma_contagens) * K) / Fatorial_N
            // Onde K é o valor de corte. A ideia é que K = X_avg + custo_c
            // K - custo_c = (soma_ponderada_custos + (Fatorial_N - soma_contagens) * K) / Fatorial_N
            // (K - custo_c) * Fatorial_N = soma_ponderada_custos + (Fatorial_N - soma_contagens) * K
            // K * Fatorial_N - custo_c * Fatorial_N = soma_ponderada_custos + K * Fatorial_N - K * soma_contagens
            // K * soma_contagens = soma_ponderada_custos + custo_c * Fatorial_N
            // K = (soma_ponderada_custos + custo_c * Fatorial_N) / soma_contagens
            //
            // No código: sp = soma_contagens_ate_p[i], sd = soma_ponderada_custos_ate_p[i].
            // A linha: X = (Fraction) { sd + c * (fac - sp), sp } + c;
            // Isto parece estar calculando o X no caso onde os custos > p[i].first são limitados a p[i].first
            // E o X final é X_avg + c.
            // Reinterpretando a linha original: `X = (Fraction) { sd + c * (fac - sp), sp } + c;`
            // Isso parece ser uma solução para `X_avg = (sum_original + (total_perm - sp) * (X_avg + c)) / total_perm`
            // `X_avg * total_perm = sum_original + (total_perm - sp) * X_avg + (total_perm - sp) * c`
            // `X_avg * (total_perm - (total_perm - sp)) = sum_original + (total_perm - sp) * c`
            // `X_avg * sp = sum_original + (total_perm - sp) * c`
            // `X_avg = (sum_original + (total_perm - sp) * c) / sp`
            // Se `sum_original` = `sd` (soma ponderada dos f_i onde f_i <= p[i].first),
            // e `sp` = `sum_{f_i <= p[i].first} count_i`.
            // Então, o `X_avg` calculado é quando o ponto de corte é *p[i].first*.
            // A variável `X` final que o código busca é o *custo final da operação de reordenação*, ou seja, `X_avg + c`.
            // Então `X_final = ((sd + (fatorial_n - soma_contagens) * custo_c) / soma_contagens) + custo_c`.
            // Ou seja, o termo `+ c` no final da linha do código é o custo fixo.
            // O `(sd + c * (fac - sp), sp)` calcula o `X_avg`.
            // A condição `max((__int128)c, p[i].first) * sp <= sd + c * fac && sd + c * fac <= p[i+1].first * sp`
            // Verifica se o `X_avg + c` calculado (que seria o ponto de corte `Y`) realmente cai entre `p[i].first` e `p[i+1].first`.
            // Simplificando a condição: `Y_candidato = (soma_ponderada_custos + custo_c * (fatorial_n - soma_contagens)) / soma_contagens + custo_c`
            // A verificação é se `custos_por_inversoes[i].custo <= Y_candidato <= custos_por_inversoes[i+1].custo`.
            // Multiplicando por soma_contagens:
            // `custos_por_inversoes[i].custo * soma_contagens <= soma_ponderada_custos + custo_c * (fatorial_n - soma_contagens) + custo_c * soma_contagens`
            // `custos_por_inversoes[i].custo * soma_contagens <= soma_ponderada_custos + custo_c * fatorial_n`
            // E similarmente para o limite superior.
            if (i == max_inversoes) { // Caso especial para o último elemento
                Fracao X_avg_calc(soma_ponderada_custos + custo_c * (fatorial_n - soma_contagens), soma_contagens);
                valor_X_otimo = X_avg_calc + custo_c;
                break;
            }

            int128_t lim_inf_Y = custos_por_inversoes[i].custo;
            int128_t lim_sup_Y = custos_por_inversoes[i+1].custo;

            // X_final_candidato = X_avg + custo_c, onde X_avg é calculado com ponto de corte = lim_inf_Y
            // Ou seja, para todas as permutações com custo original > lim_inf_Y, o custo se torna lim_inf_Y.
            // X_avg = (soma_ponderada_custos + (fatorial_n - soma_contagens) * lim_inf_Y) / fatorial_n
            // X_final_candidato = X_avg + custo_c
            // O código usa um jeito diferente, onde `X` é a média *final* de custo
            // X * soma_contagens = soma_ponderada_custos + custo_c * (fatorial_n - soma_contagens)
            // Este X é o X_avg. O custo total é X_avg + custo_c.
            // Se Y_corte = X_avg + c, então:
            // Y_corte = (soma_ponderada_custos + (fatorial_n - soma_contagens) * c) / soma_contagens + c
            // Vamos chamar o Y_corte calculado de Y_calc.
            // Y_calc = (soma_ponderada_custos + (fatorial_n - soma_contagens) * custo_c) / soma_contagens + custo_c;
            // Precisamos que lim_inf_Y <= Y_calc <= lim_sup_Y.
            // (lim_inf_Y - custo_c) * soma_contagens <= soma_ponderada_custos + (fatorial_n - soma_contagens) * custo_c
            // lim_inf_Y * soma_contagens - custo_c * soma_contagens <= soma_ponderada_custos + fatorial_n * custo_c - soma_contagens * custo_c
            // lim_inf_Y * soma_contagens <= soma_ponderada_custos + fatorial_n * custo_c (Condição 1)

            // lim_sup_Y * soma_contagens - custo_c * soma_contagens >= soma_ponderada_custos + (fatorial_n - soma_contagens) * custo_c
            // lim_sup_Y * soma_contagens >= soma_ponderada_custos + fatorial_n * custo_c (Condição 2)

            // Portanto, a condição é:
            // (soma_ponderada_custos + fatorial_n * custo_c) deve estar entre (lim_inf_Y * soma_contagens) e (lim_sup_Y * soma_contagens)
            // (inclusive no limite inferior, exclusivo no superior para o loop i < N).

            if (lim_inf_Y * soma_contagens <= soma_ponderada_custos + fatorial_n * custo_c && 
                soma_ponderada_custos + fatorial_n * custo_c < lim_sup_Y * soma_contagens) {
                
                // Se a condição é verdadeira, encontramos o intervalo onde X_avg + c se encaixa
                // O valor X_avg é (soma_ponderada_custos + custo_c * (fatorial_n - soma_contagens)) / soma_contagens
                // E o custo final da reordenação é X_avg + custo_c
                valor_X_otimo = Fracao(soma_ponderada_custos + custo_c * (fatorial_n - soma_contagens), soma_contagens) + custo_c;
                break; // Encontramos o X ótimo
            }
        }
        
        while (num_consultas--) {
            int128_t contador_inversoes_atual = 0;
            for (int i = 1; i <= n; ++i) {
                std::cin >> sequencia_atual[i];
                for (int j = 1; j < i; ++j) {
                    if (sequencia_atual[j] > sequencia_atual[i]) {
                        contador_inversoes_atual++;
                    }
                }
            }
            std::cout << Fracao::min_val(dp_permutacao_inversoes[contador_inversoes_atual], valor_X_otimo) << "\n";
        }
    }
    return 0;
}

Problema 2: Fatoração de Polinômios Ciclótomicos

Dada a identidade fundamental para polinômios ciclótomicos, \(x^n - 1 = \prod_{d | n} \Phi_d(x)\), onde \(\Phi_n(x)\) representa o n-ésimo polinômio ciclotômico. Aplicando o logaritmo natural em ambos os lados e utilizando a Inversão de Möbius, podemos derivar uma fórmula explícita para \(\Phi_n(x)\): \(\Phi_n(x) = \prod_{d | n}(x^d - 1)^{\mu(\frac{n}{d})}\). Isso significa que cada \(\Phi_n(x)\) pode ser expresso como um produto e quociente de binômios da forma \((x^d - 1)\), onde os expoentes são determinados pelos valores da função de Möbius \(\mu\).

A multiplicação de polinômios pode ser vista como uma variação do problema da mochila 0/1, enquanto a divisão de polinômios, especialmente por termos como \((x-1)\), pode ser modelada como uma "mochila infinita" usando a série geométrica \(\frac{1}{1-x} = \sum_{i=0}^\infty x^i\). No entanto, o grau de \(\Phi_n(x)\) é \(\varphi(n)\) (função totiente de Euler), o que torna a abordagem direta com polinômios de grau arbitrário computacionalmente intensiva para alguns n. Uma alternativa mais eficiente, e a que foi implementada, envolve a construção de \(\Phi_n(x)\) recursivamente a partir de seus divisores, explorando propriedades conhecidas.

Exemplo de Implementação para Fatoração

#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric> // Para std::iota, std::gcd

typedef long long int_ll;

const int MAX_N = 2000005;
std::vector<int_ll> divisores_de_n[MAX_N]; // divisores_de_n[j] armazena divisores 'i' de 'j' tal que 'n % i == 0'
std::vector<int_ll> polinomios_ciclotomicos[MAX_N]; // Armazena coeficientes de Phi_k(x)
int_ll mobius[MAX_N]; // Valores da função de Möbius
int_ll euler_phi[MAX_N]; // Valores da função totiente de Euler
int_ll primos[MAX_N]; // Lista de números primos
int_ll maior_divisor_primo[MAX_N]; // Maior divisor primo de i (ou o próprio i se for primo)
int_ll contador_primos;
bool eh_composto[MAX_N]; // Marca números compostos

// Crivo para calcular Möbius e Euler Phi
void Crivo(int n_max) {
    euler_phi[1] = mobius[1] = 1;
    for (int_ll i = 2; i <= n_max; ++i) {
        if (!eh_composto[i]) {
            mobius[i] = -1;
            euler_phi[i] = i - 1;
            primos[++contador_primos] = i;
        }
        for (int_ll j = 1; j <= contador_primos && primos[j] * i <= n_max; ++j) {
            eh_composto[primos[j] * i] = true;
            if (i % primos[j] == 0) {
                mobius[primos[j] * i] = 0;
                euler_phi[primos[j] * i] = euler_phi[i] * primos[j];
                break;
            }
            mobius[primos[j] * i] = -mobius[i];
            euler_phi[primos[j] * i] = euler_phi[i] * euler_phi[primos[j]];
        }
    }
    // Para Phi_n(x) onde n não é livre de quadrados, usamos Phi_n(x) = Phi_{n/p}(x^p)
    // Precisamos do maior divisor 'p' tal que p|n e p é primo.
    // Ou, se n é da forma p^k * m, onde p não divide m, então Phi_n(x) = Phi_{p*m}(x^{p^(k-1)})
    // A implementação usa uma propriedade diferente para !mu[x] casos.
    // Para simplificar, vou adaptar a lógica do mx original para maior_divisor_primo
    // O original usa mx[j] = max(mx[j], i) se !mu[i]. Isto encontra o maior divisor 'i' de 'j' tal que mu[i] != 0
    // Isso é útil para Phi_{p^k}(x) = Phi_p(x^{p^{k-1}}).
    // O mx no código original parece ser o maior divisor 'd' de 'x' tal que mu[d] != 0.
    // Isso não é trivial e é a parte mais obscura do código original.
    // Vou usar maior_divisor_primo para o maior fator primo de i, pois mu[i] != 0 se i é livre de quadrados.
    // Se mu[x] == 0, significa que x não é livre de quadrados, i.e., x tem um fator primo ao quadrado.
    // Se x = p*m onde p|m, então Phi_x(y) = Phi_{x/p}(y^p). O maior divisor quadrado de x é crucial.
    // A linha `mx[j] = max(mx[j], i);` dentro do `if (!mu[i]) continue;` é o que confunde.
    // Parece que para cada `j`, `mx[j]` é o maior `i` tal que `i|j` e `mu[i]` é não-zero.
    // Isso é o maior divisor livre de quadrados de `j`.
    std::iota(maior_divisor_primo, maior_divisor_primo + n_max + 1, 0); // Inicializa com o próprio número
    for (int_ll i = 1; i <= n_max; ++i) {
        if (!mobius[i]) continue; // Se i não é livre de quadrados, não é o tipo de divisor que estamos procurando.
        for (int_ll j = i; j <= n_max; j += i) {
            // Se j tem um divisor 'i' que é livre de quadrados, 'i' é um candidato para mx[j].
            // O mx[j] deve ser o maior tal divisor.
            // Para simplificar a lógica e manter a equivalência com a original:
            // mx[j] = o maior divisor de j que é livre de quadrados.
            // Isso é o que a linha original parece fazer:
            // for (int i = 1; i <= n; i++) { if (!mu[i]) continue; for (int j = i; j <= n; j += i) mx[j] = max(mx[j], i); }
            maior_divisor_primo[j] = std::max(maior_divisor_primo[j], i);
        }
    }
}

// Compara dois vetores de coeficientes (polinômios)
bool comparar_polinomios(const std::vector<int_ll>& p1, const std::vector<int_ll>& p2) {
    if (p1.size() != p2.size()) {
        return p1.size() < p2.size();
    }
    for (int i = p1.size() - 1; i >= 0; --i) {
        if (p1[i] != p2[i]) {
            return p1[i] < p2[i];
        }
    }
    return false; // São iguais
}

// Imprime um polinômio formatado
void imprimir_polinomio(std::ostream& os, std::vector<int_ll>& p) {
    os << "(";
    int grau_max = p.size() - 1;

    // Ajustar o sinal do primeiro coeficiente (de maior grau) se negativo
    // para garantir que o termo líder seja positivo e o resto se ajuste.
    bool inverter_sinais = false;
    if (grau_max >= 0 && p[grau_max] < 0) {
        inverter_sinais = true;
        for (int i = 0; i <= grau_max; ++i) {
            p[i] = -p[i];
        }
    }

    bool primeiro_termo_impresso = false;
    for (int i = grau_max; i >= 0; --i) {
        if (p[i] == 0) continue;

        if (primeiro_termo_impresso) {
            os << (p[i] > 0 ? "+" : "");
        }
        
        if (p[i] == 1) {
            if (i == 0) os << "1"; // Coeficiente 1, grau 0
            else os << "x"; // Coeficiente 1, grau > 0
        } else if (p[i] == -1) {
            if (i == 0) os << "-1";
            else os << "-x";
        } else {
            os << p[i];
            if (i > 0) os << "x";
        }
        
        if (i > 1) os << "^" << i;
        
        primeiro_termo_impresso = true;
    }
    // Se o polinômio é zero, imprime "(0)"
    if (!primeiro_termo_impresso) os << "0";
    os << ")";
}

int main() {
    std::ios_base::sync_with_stdio(false);
    std::cin.tie(NULL);

    Crivo(MAX_N - 1);

    int num_testes;
    std::cin >> num_testes;
    while (num_testes--) {
        std::vector<std::vector<int_ll>> resultados_fatoracao;
        int_ll n_val;
        std::cin >> n_val;

        // Limpar e preencher divisores_de_n para a entrada atual n_val
        for (int_ll i = 1; i * i <= n_val; ++i) {
            if (n_val % i == 0) {
                divisores_de_n[n_val].push_back(i);
                if (i * i != n_val) {
                    divisores_de_n[n_val].push_back(n_val / i);
                }
            }
        }
        std::sort(divisores_de_n[n_val].begin(), divisores_de_n[n_val].end());

        // Para cada divisor 'x' de n_val, calcular Phi_x(t)
        for (int_ll x : divisores_de_n[n_val]) {
            polinomios_ciclotomicos[x].assign(euler_phi[x] + 1, 0); // Redimensionar e zerar
            polinomios_ciclotomicos[x][0] = 1; // Termo constante é 1 para começar
            
            // Este bloco lida com casos onde x não é livre de quadrados (mu[x] == 0)
            // Usando a propriedade Phi_{p*m}(x) = Phi_m(x^p) se p|m
            // No código original, se !mu[x] (x não livre de quadrados), ele usa mx[x]
            // que é o maior divisor de x que é livre de quadrados.
            // Aparentemente, a propriedade usada é: se x = p * m, onde p é o menor primo tal que p^2 | x,
            // então Phi_x(t) = Phi_{x/p}(t^p).
            // O mx[x] no original é o maior divisor livre de quadrados de x.
            // Para !mu[x], a linha é: `for (int j = 0; j <= phi[mx[x]]; j++) poly[x][j * x / mx[x]] = poly[mx[x]][j];`
            // Isso significa que Phi_x(t) é basicamente Phi_{mx[x]}(t^(x/mx[x])).
            // Exemplo: Phi_4(x) = Phi_2(x^2) = x^2 + 1. Aqui mx[4]=2, x/mx[x]=2.
            // Phi_8(x) = Phi_4(x^2) = x^4+1. Aqui mx[8]=2, x/mx[x]=4.
            // O que se encaixa é que se x = p^k * m, onde m é livre de quadrados,
            // então Phi_x(t) = Phi_m(t^{p^k}). Aqui m é o mx[x].
            if (!mobius[x]) {
                int_ll fator_potencia = x / maior_divisor_primo[x];
                // Precisamos garantir que maior_divisor_primo[x] já tenha sido calculado (no sentido de Phi_{maior_divisor_primo[x]})
                // Isso é garantido pela ordem de processamento (x é divisor, e maior_divisor_primo[x] é divisor de x e <= x).
                for (int_ll j = 0; j <= euler_phi[maior_divisor_primo[x]]; ++j) {
                    if (j * fator_potencia <= euler_phi[x]) {
                        polinomios_ciclotomicos[x][j * fator_potencia] = polinomios_ciclotomicos[maior_divisor_primo[x]][j];
                    }
                }
            } else { // x é livre de quadrados (mu[x] != 0)
                // Os laços seguintes implementam a divisão de polinômios para obter Phi_x(t)
                // Baseado na relação x^n - 1 = product_{d|n} Phi_d(x)
                // e Phi_n(x) = (x^n - 1) / product_{d|n, d!=n} Phi_d(x).
                // Isso é feito como uma operação de "mochila" com sinais.
                // Na verdade, Phi_n(x) = product_{d|n} (x^d - 1)^mu(n/d)
                // O código parece implementar de outra forma: Phi_n(x) = x^n - 1 / prod_{d|n, d!=n} Phi_d(x)
                // Para Phi_x, ele constrói o polinômio de grau euler_phi[x].
                // É uma forma de calcular Phi_x(t) = (x^x - 1) / Prod_{d|x, d<x} Phi_d(x)
                // O código faz isso de forma iterativa:
                // Inicialmente, poly[x][0] = 1.
                // Depois, para cada divisor 'v' de 'x':
                // Se mu[v] == mu[x], significa que 'v' e 'x' têm a mesma paridade de fatores primos (livre de quadrados).
                // Isso corresponde a termos no denominador que são (x^v - 1).
                // O original subtrai, simulando (1 - x^v)^-1 = 1 + x^v + x^(2v) + ...
                // Mas parece ser algo mais complexo.
                // A lógica é provavelmente Phi_x(t) = Prod_{d|x} (1 - t^d)^{mu(x/d)}.
                // Os laços no código parecem estar expandindo os termos (1 - x^v) com base nos divisores.
                // O termo mu[v] == mu[x] ou mu[v] != mu[x] controla se é um termo (1 - x^v) ou (1 - x^v)^(-1).
                // Se mu[x] = 1, e mu[v] = 1, então mu[x/v] pode ser -1.
                // Se mu[x] = -1, e mu[v] = -1, então mu[x/v] pode ser 1.
                // Essa é uma forma de fazer convolução polinomial com base na inversão de Möbius.
                // É um DP baseado na relação recursiva de Phi_n(x).
                // Por exemplo, para Phi_6(x) = x^2 - x + 1.
                // Divisores de 6: 1, 2, 3, 6.
                // mu[1]=1, mu[2]=-1, mu[3]=-1, mu[6]=1.
                // poly[6][0]=1.
                // v=1: mu[1]=1, mu[6]=1 (iguais). Subtrai poly[6][j-1]. poly[6] = (1, -1, 0, ...)
                // v=2: mu[2]=-1, mu[6]=1 (diferentes). Soma poly[6][j-2]. poly[6] = (1, -1, 1, -1, ...)
                // v=3: mu[3]=-1, mu[6]=1 (diferentes). Soma poly[6][j-3].
                // É uma construção de Phi_x(t) a partir de (1-t^v) ou (1+t^v) termos.
                // O original usa:
                // for (auto v : fac[x]) { if (mu[v] == mu[x]) { for (int j = phi[x]; j >= v; j--) poly[x][j] -= poly[x][j - v]; } }
                // for (auto v : fac[x]) { if (mu[v] != mu[x]) { for (int j = v; j <= phi[x]; j++) poly[x][j] += poly[x][j - v]; } }
                // Esta é uma implementação da fórmula: Phi_n(x) = prod_{d|n, mu(n/d)=1} (x^d - 1) / prod_{d|n, mu(n/d)=-1} (x^d - 1)
                // Que é equivalente a Phi_n(x) * prod_{d|n, mu(n/d)=-1} (x^d-1) = prod_{d|n, mu(n/d)=1} (x^d-1).
                // O código está fazendo um DP para resolver esta equação.
                // A forma como mu[v] é comparado a mu[x] implica que ele está usando a relação:
                // Phi_n(x) = \prod_{d|n} (x^{n/d}-1)^{\mu(d)} (esta é a formula de Möbius)
                // ou Phi_n(x) = (x^n-1) / \prod_{d|n, d<n} Phi_d(x)
                // O loop `j >= v` para subtração e `j <= v` para adição é o que faz o "polinômio de mochila".
                // `poly[x][j] -= poly[x][j-v]` significa `poly[x]` é dividido por `(1 - x^v)`.
                // `poly[x][j] += poly[x][j-v]` significa `poly[x]` é multiplicado por `(1 + x^v)` ou dividido por `(1 - x^v)` (se os sinais estiverem invertidos).
                
                // Primeiro loop para termos com mu[v] == mu[x]
                for (int_ll v : divisores_de_n[x]) {
                    if (v == 0) continue; // v deve ser positivo
                    if (mobius[v] == mobius[x]) {
                        for (int_ll j = euler_phi[x]; j >= v; --j) {
                            polinomios_ciclotomicos[x][j] -= polinomios_ciclotomicos[x][j - v];
                        }
                    }
                }
                // Segundo loop para termos com mu[v] != mu[x]
                for (int_ll v : divisores_de_n[x]) {
                    if (v == 0) continue;
                    if (mobius[v] != mobius[x]) {
                        for (int_ll j = v; j <= euler_phi[x]; ++j) {
                            polinomios_ciclotomicos[x][j] += polinomios_ciclotomicos[x][j - v];
                        }
                    }
                }
            }
            resultados_fatoracao.push_back(polinomios_ciclotomicos[x]);
        }
        std::sort(resultados_fatoracao.begin(), resultados_fatoracao.end(), comparar_polinomios);

        for (auto& poly_vec : resultados_fatoracao) {
            imprimir_polinomio(std::cout, poly_vec);
        }
        std::cout << "\n";

        // Limpar estruturas para o próximo caso de teste
        for (int_ll x : divisores_de_n[n_val]) {
            divisores_de_n[x].clear(); // Limpa os divisores pré-calculados para x
            polinomios_ciclotomicos[x].clear(); // Limpa o polinômio para x
        }
        divisores_de_n[n_val].clear(); // Limpa para n_val
        polinomios_ciclotomicos[n_val].clear();
    }
    return 0;
}

Problema 3: Visualização em Árvores

Este problema explora conceitos de geometria em árvores, especificamente "círculos" e diâmetros. Dada uma coleção de pontos S em uma árvore, definimos seu diâmetro como o caminho mais longo entre quaisquer dois pontos em S. O ponto médio do diâmetro, mid(S), é o centro de seu círculo de cobertura mínimo. Um "círculo" em uma árvore, c(u, r), é o conjunto de todos os pontos x tais que dist(x, u) \(\le\) r. O círculo de cobertura mínima para um conjunto S, c(S), é então c(mid(S), len(S)/2).

A fusão de círculos de cobertura mínima segue regras específicas: se c(S) está contido em c(T), então c(S \(\cup\) T) = c(T). Similarmente, se c(T) está contido em c(S), então c(S \(\cup\) T) = c(S). Caso contrário (os círculos se intersectam ou são disjuntos), o diâmetro do conjunto combinado len(S \(\cup\) T) é dado por dist(mid(S), mid(T)) + len(S)/2 + len(T)/2.

A solução emprega uma abordagem de "dividir para conquistar" (divide and conquer) sobre os intervalos de pontos [l, r]. No passo de fusão, para um ponto médio mid, precisamos considerar a contribuição de pares de subproblemas [l, mid] e (mid, r]. À medida que r aumenta para o intervalo direito, o círculo de cobertura c((mid, r]) também se expande. Isso significa que a relação entre c([l, mid]) e c((mid, r]) muda monotonicamente: de c([l, mid]) contendo c((mid, r]), para eles se intersectando/sendo disjuntos, e finalmente c((mid, r]) contendo c([l, mid]). Essas fronteiras podem ser mantidas com dois ponteiros, p e q.

A soma das distâncias \(\sum dist(x, y)\) em uma árvore é um problema comum. Esta implementação o resolve eficientemente usando Decomposição Heavy-Light (HLD) e uma Árvore de Índices Binários (BIT) para consultas de soma de alcance e atualizações de ponto. A complexidade desta abordagem é \(\mathcal{O}(N \log^3 N)\), que é prática devido à pequena constante. Antes de aplicar HLD, a estrutura da árvore é modificada: cada aresta \((u, v)\) é substituída por dois nós e duas arestas, i.e., \((u, \text{novo\_no})\) e \((\text{novo\_no}, v)\). Isso efetivamente dobra todas as distâncias de caminho e garante que todos os "midpoints" do diâmetro caiam em nós existentes, simplificando os cálculos.

Exemplo de Implementação para Visualização

#include <iostream>
#include <vector>
#include <algorithm>

// Macro para obter o bit menos significativo
#define lowbit(x) ((x) & (-(x)))

const int MAX_NOS = 200005; // Ajustado para 2*N devido aos nós intermediários

int num_nos_original; // N original, o total de nós na árvore expandida é 2*N - 1
std::vector<int> adj[MAX_NOS * 2]; // Lista de adjacência para a árvore expandida
int profundidade[MAX_NOS * 2]; // profundidade[x] = profundidade do nó x
int ordem_dfs[MAX_NOS * 2]; // ordem_dfs[x] = tempo de entrada na DFS para o nó x
int ordem_dfs_inversa[MAX_NOS * 2]; // ordem_dfs_inversa[t] = nó com tempo de entrada t
int pai[MAX_NOS * 2]; // pai[x] = pai do nó x na DFS
int cabeca_caminho[MAX_NOS * 2]; // cabeca_caminho[x] = nó líder do caminho pesado que x pertence
int tamanho_subarvore[MAX_NOS * 2]; // tamanho_subarvore[x] = tamanho da subárvore de x
int filho_pesado[MAX_NOS * 2]; // filho_pesado[x] = filho pesado de x
int contador_dfs; // Contador de tempo para DFS

// Adiciona uma aresta bidirecional
void adicionar_aresta(int u, int v) {
    adj[u].push_back(v);
    adj[v].push_back(u);
}

// DFS1 para calcular tamanhos de subárvore, pais e filhos pesados
void dfs1_hld(int u, int p, int d) {
    tamanho_subarvore[u] = 1;
    pai[u] = p;
    profundidade[u] = d;
    filho_pesado[u] = 0; // Resetar filho pesado

    for (int v : adj[u]) {
        if (v == p) continue;
        dfs1_hld(v, u, d + 1);
        tamanho_subarvore[u] += tamanho_subarvore[v];
        if (tamanho_subarvore[v] > tamanho_subarvore[filho_pesado[u]]) {
            filho_pesado[u] = v;
        }
    }
}

// DFS2 para atribuir cabeças de caminho e ordens DFS
void dfs2_hld(int u, int h) {
    ordem_dfs[u] = ++contador_dfs;
    ordem_dfs_inversa[contador_dfs] = u;
    cabeca_caminho[u] = h;

    if (filho_pesado[u]) {
        dfs2_hld(filho_pesado[u], h);
    }
    for (int v : adj[u]) {
        if (v == pai[u] || v == filho_pesado[u]) continue;
        dfs2_hld(v, v); // Novos caminhos pesados começam com v
    }
}

// Calcula o LCA (Lowest Common Ancestor) de dois nós
int calcular_lca(int u, int v) {
    while (cabeca_caminho[u] != cabeca_caminho[v]) {
        if (profundidade[cabeca_caminho[u]] < profundidade[cabeca_caminho[v]]) {
            v = pai[cabeca_caminho[v]];
        } else {
            u = pai[cabeca_caminho[u]];
        }
    }
    return (profundidade[u] < profundidade[v]) ? u : v;
}

// Calcula a distância entre dois nós
int calcular_distancia(int u, int v) {
    return profundidade[u] + profundidade[v] - 2 * profundidade[calcular_lca(u, v)];
}

// Encontra o k-ésimo ancestral de um nó
int k_esimo_ancestral(int u, int k) {
    while (k > 0) {
        if (profundidade[u] - profundidade[cabeca_caminho[u]] < k) {
            k -= (profundidade[u] - profundidade[cabeca_caminho[u]] + 1);
            u = pai[cabeca_caminho[u]];
        } else {
            return ordem_dfs_inversa[ordem_dfs[u] - k];
        }
    }
    return u; // Retorna u se k=0 ou se for um ancestral acima da raiz (não deveria acontecer)
}

// Encontra o nó no caminho de u para v que está a distância k de u
int encontrar_no_no_caminho(int u, int v, int k) {
    int ancestral_comum = calcular_lca(u, v);
    if (profundidade[u] - profundidade[ancestral_comum] >= k) {
        return k_esimo_ancestral(u, k);
    } else {
        return k_esimo_ancestral(v, calcular_distancia(u, v) - k);
    }
}

// Estrutura para armazenar informações do diâmetro
struct InfoDiametro {
    int p1, p2, comprimento; // p1 e p2 são os endpoints do diâmetro, comprimento é a distância
    void adicionar_ponto(int novo_ponto) {
        int d1 = calcular_distancia(novo_ponto, p1);
        int d2 = calcular_distancia(novo_ponto, p2);
        // Garantir que d1 <= d2 para simplificar a lógica
        if (d1 > d2) {
            std::swap(d1, d2);
            std::swap(p1, p2);
        }
        if (d2 > comprimento) {
            comprimento = d2;
            p1 = novo_ponto; // O novo_ponto se torna um dos endpoints do novo diâmetro
        }
    }
};

// Árvore de Índices Binários (BIT) para soma de alcance e atualização de ponto
// Usamos duas BITs para implementar atualização de alcance e consulta de alcance em O(logN)
struct BIT {
    long long arvore_indices[MAX_NOS * 2]; // O tamanho deve ser grande o suficiente para o dfn
    int tamanho_max;

    BIT(int max_val) : tamanho_max(max_val) {
        std::fill(arvore_indices, arvore_indices + tamanho_max + 1, 0);
    }

    void atualizar(int indice, long long valor) {
        for (; indice <= tamanho_max; indice += lowbit(indice)) {
            arvore_indices[indice] += valor;
        }
    }

    long long consultar(int indice) {
        long long soma = 0;
        for (; indice > 0; indice -= lowbit(indice)) {
            soma += arvore_indices[indice];
        }
        return soma;
    }
};

BIT bit_soma_valores(MAX_NOS * 2); // Para soma(v)
BIT bit_soma_indices_x_valores(MAX_NOS * 2); // Para soma(i*v)

// Atualiza um alcance [l, r] com um valor 'val'
void atualizar_alcance(int l, int r, long long val) {
    bit_soma_valores.atualizar(l, val);
    bit_soma_valores.atualizar(r + 1, -val);
    bit_soma_indices_x_valores.atualizar(l, - (long long)l * val);
    bit_soma_indices_x_valores.atualizar(r + 1, (long long)(r + 1) * val);
}

// Consulta a soma de um alcance [l, r]
long long consultar_alcance_soma(int l, int r) {
    // query(r) = sum(val_i) * (r+1) + sum(i*val_i)
    long long soma_r = bit_soma_valores.consultar(r) * (r + 1) + bit_soma_indices_x_valores.consultar(r);
    long long soma_l_menos_1 = bit_soma_valores.consultar(l - 1) * l + bit_soma_indices_x_valores.consultar(l - 1);
    return soma_r - soma_l_menos_1;
}

// Adiciona um valor `v` a todos os nós no caminho de `u` até a raiz (ou até a cabeça do caminho 1)
void adicionar_no_caminho(int u, long long v) {
    while (cabeca_caminho[u] != 1) { // Enquanto não chegarmos ao caminho que contém a raiz
        atualizar_alcance(ordem_dfs[cabeca_caminho[u]], ordem_dfs[u], v);
        u = pai[cabeca_caminho[u]];
    }
    // Lidar com o caminho que contém a raiz (excluir o nó 1, se necessário, ou incluir tudo)
    // Se u == 1, apenas o nó 1 está sendo atualizado.
    // O original usa dfn[1]+1, dfn[x] se x != 1, isso parece pular a raiz 1 em algumas operações.
    // Vamos usar ordem_dfs[1] para a raiz
    if (u != 0) { // u != 0 porque pai[1] é 0, profundidade[0] = -1.
        atualizar_alcance(ordem_dfs[1], ordem_dfs[u], v);
    }
}

// Consulta a soma dos valores no caminho da raiz até `u`
long long consultar_caminho_soma(int u) {
    long long soma_total = 0;
    while (cabeca_caminho[u] != 1) {
        soma_total += consultar_alcance_soma(ordem_dfs[cabeca_caminho[u]], ordem_dfs[u]);
        u = pai[cabeca_caminho[u]];
    }
    if (u != 0) {
        soma_total += consultar_alcance_soma(ordem_dfs[1], ordem_dfs[u]);
    }
    return soma_total;
}

long long soma_diametros_globais; // Resposta final

int comprimentos_diametros[MAX_NOS]; // len[i]
int midpoints_diametros[MAX_NOS]; // md[i]

// Função principal de dividir e conquistar
void resolver_subproblema(int l, int r) {
    if (l == r) return;
    if (r == l + 1) {
        soma_diametros_globais += calcular_distancia(l, r);
        return;
    }

    int mid = (l + r) / 2;
    
    // Cálculo para o lado direito (mid+1 a r)
    long long soma_comprimentos_direita = 0;
    InfoDiametro info_direita = {mid + 1, mid + 1, 0};
    for (int i = mid + 1; i <= r; ++i) {
        info_direita.adicionar_ponto(i);
        comprimentos_diametros[i] = info_direita.comprimento;
        midpoints_diametros[i] = encontrar_no_no_caminho(info_direita.p1, info_direita.p2, info_direita.comprimento / 2);
        soma_comprimentos_direita += comprimentos_diametros[i];
    }

    // Cálculo para o lado esquerdo (l a mid) e fusão
    InfoDiametro info_esquerda = {mid, mid, 0};
    int ptr_q = mid + 1; // Ponteiro para o início da região "conter" (c(esquerda) contém c(direita))
    int ptr_p = mid + 1; // Ponteiro para o início da região "intersectar/disjunto" (c(direita) contém c(esquerda))

    for (int i = mid; i >= l; --i) {
        info_esquerda.adicionar_ponto(i);
        int diametro_esq = info_esquerda.comprimento;
        int midpoint_esq = encontrar_no_no_caminho(info_esquerda.p1, info_esquerda.p2, info_esquerda.comprimento / 2);

        // Movimentar ptr_q: Expandir a região onde c(esquerda) contém c(direita)
        // Condição: dist(midpoint_esq, midpoint_dir) + diametro_esq/2 <= diametro_dir/2
        // Isso significa que o círculo direito é pequeno o suficiente para ser contido no esquerdo
        while (ptr_q <= r && !(calcular_distancia(midpoint_esq, midpoints_diametros[ptr_q]) + diametro_esq / 2 <= comprimentos_diametros[ptr_q] / 2)) {
            // Este ponto já não está na região de "contido" (c(esquerda) contém c(direita))
            // mas está entrando na região de "intersectar/disjunto" (ou já estava na região de c(direita) contém c(esquerda))
            // Remove da soma de comprimentos_direita e adiciona para cálculos da região do meio
            soma_comprimentos_direita -= comprimentos_diametros[ptr_q];
            // Para a região intermediária, a contribuição é dist(mid(S), mid(T)) + r(S) + r(T)
            // r(T) = comprimentos_diametros[ptr_q] / 2
            adicionar_no_caminho(midpoints_diametros[ptr_q], 1); // Adiciona 1 ao count do nó md[ptr_q]
            ptr_q++;
        }

        // Movimentar ptr_p: Expandir a região onde c(direita) contém c(esquerda)
        // Condição: dist(midpoint_esq, midpoint_dir) + diametro_dir/2 <= diametro_esq/2
        // Isso significa que o círculo esquerdo é pequeno o suficiente para ser contido no direito
        while (ptr_p < ptr_q && calcular_distancia(midpoint_esq, midpoints_diametros[ptr_p]) + comprimentos_diametros[ptr_p] / 2 <= diametro_esq / 2) {
            // Este ponto já não está na região de "intersectar/disjunto"
            // mas está entrando na região de c(direita) contém c(esquerda)
            // Remove das contas da região do meio
            adicionar_no_caminho(midpoints_diametros[ptr_p], -1); // Remove 1 do count do nó md[ptr_p]
            ptr_p++;
        }
        
        // Contribuição da região onde c(direita) contém c(esquerda) (l <= i <= mid, mid+1 <= j <= ptr_p-1)
        // Para cada j em [mid+1, ptr_p-1], c(midpoints_diametros[j]) contém c(info_esquerda)
        // A soma dos diâmetros é diametro_esq. Então (ptr_p - (mid + 1)) * diametro_esq.
        soma_diametros_globais += (long long)(ptr_p - (mid + 1)) * diametro_esq;

        // Contribuição da região de interseção/disjunção (l <= i <= mid, ptr_p <= j <= ptr_q-1)
        // A soma dos diâmetros é sum_{j=ptr_p to ptr_q-1} (dist(midpoint_esq, midpoints_diametros[j]) + diametro_esq/2 + comprimentos_diametros[j]/2)
        // = (ptr_q - ptr_p) * diametro_esq/2 + sum_{j=ptr_p to ptr_q-1} (dist(midpoint_esq, midpoints_diametros[j]) + comprimentos_diametros[j]/2)
        // dist(x,y) = dep[x] + dep[y] - 2*dep[LCA(x,y)]
        // dist(midpoint_esq, md[j]) = prof[midpoint_esq] + prof[md[j]] - 2*prof[LCA(midpoint_esq, md[j])]
        // sum (prof[midpoint_esq] + prof[md[j]]) = (ptr_q - ptr_p)*prof[midpoint_esq] + sum prof[md[j]]
        // sum (2*prof[LCA(midpoint_esq, md[j])])
        // `consultar_caminho_soma(midpoint_esq)` retorna a soma de valores ao longo do caminho.
        // `query(x)` no código original é `consultar_caminho_soma(midpoint_esq)`.
        // A lógica de `sl` e `Query(xmid)` no código original é um pouco densa.
        // sl = sum (len[p]/2 + dep[md[p]])
        // A contribuição é: (ptr_q - ptr_p) * (diametro_esq / 2) + sl + (long long)profundidade[midpoint_esq] * (ptr_q - ptr_p) - 2 * consultar_caminho_soma(midpoint_esq);
        // O termo `sl` parece ser o `sum (len[p]/2 + dep[md[p]])` para `p` de `ptr_p` a `ptr_q-1`.
        // Mas a `adicionar_no_caminho` e `consultar_caminho_soma` está trabalhando com o `dfn` dos nós.
        // A consulta para `dist(x, y)` usando HLD e BIT pode ser reescrita para ser mais clara.
        // Para cada nó `md[j]` no intervalo [ptr_p, ptr_q-1], queremos calcular `dist(midpoint_esq, md[j])`.
        // `dist(A, B) = dep[A] + dep[B] - 2 * dep[LCA(A, B)]`
        // `sum dist(midpoint_esq, md[j]) = (ptr_q - ptr_p) * dep[midpoint_esq] + sum dep[md[j]] - 2 * sum dep[LCA(midpoint_esq, md[j])]`
        // A `bit_soma_valores` e `bit_soma_indices_x_valores` estão sendo usadas para calcular `sum dep[md[j]]` e `sum dep[LCA]` indiretamente.
        // O `bit_soma_valores.query(r)` retorna a soma de todos os `v` adicionados no intervalo `[1, r]`
        // O `bit_soma_indices_x_valores.query(r)` retorna a soma de `i*v` no intervalo `[1, r]`
        // Este é um método padrão de BIT para soma de alcance. A `adicionar_no_caminho` e `consultar_caminho_soma` são consultas de caminho.
        // A consulta de caminho para `sum_{j=ptr_p to ptr_q-1} dist(midpoint_esq, midpoints_diametros[j])` é complexa.
        // Ela é calculada como `(profundidade[midpoint_esq] * count + sum_{j=ptr_p to ptr_q-1} profundidade[md[j]] - 2 * sum_{j=ptr_p to ptr_q-1} profundidade[LCA(midpoint_esq, md[j])])`
        // `sum profundidade[md[j]]` pode ser obtido usando uma BIT auxiliar nos nós `md[j]` marcados.
        // `sum profundidade[LCA(midpoint_esq, md[j])]` é a parte mais difícil e é onde a BIT para caminho HLD entra.
        // O `sl` no código original parece ser a soma das profundidades. `sl += len[q] / 2 + dep[md[q]]`
        // `sl` mantém a soma de `comprimentos_diametros[j]/2 + profundidade[midpoints_diametros[j]]` para `j` em `[ptr_p, ptr_q-1]`.
        // A expressão `(long long)(ptr_q - ptr_p) * (diametro_esq / 2) + sl + (long long)profundidade[midpoint_esq] * (ptr_q - ptr_p) - 2 * consultar_caminho_soma(midpoint_esq);`
        // Soma_diametros_globais += (qtd_nos_meio) * raio_esq + sum(raio_dir) + sum(profundidade[md_dir]) + (qtd_nos_meio)*profundidade[mid_esq] - 2*sum(profundidade[LCA(mid_esq, md_dir)])
        // O `sl` está somando os `len[q]/2 + dep[md[q]]` para os nós que estão em `[ptr_p, ptr_q-1]`.
        // A parte `sum_{j=ptr_p to ptr_q-1} (dist(midpoint_esq, midpoints_diametros[j]) + comprimentos_diametros[j]/2)`
        // é `(profundidade[midpoint_esq] * (ptr_q - ptr_p) + sum_{j=ptr_p to ptr_q-1} profundidade[midpoints_diametros[j]] - 2 * sum_{j=ptr_p to ptr_q-1} profundidade[LCA(midpoint_esq, midpoints_diametros[j])]) + sum_{j=ptr_p to ptr_q-1} comprimentos_diametros[j]/2`
        // No código, `sl` tem `sum_{j=ptr_p to ptr_q-1} (profundidade[midpoints_diametros[j]] + comprimentos_diametros[j]/2)`
        // `Query(xmid)` retorna `sum_{j=ptr_p to ptr_q-1} profundidade[LCA(midpoint_esq, midpoints_diametros[j])]`.
        // Então, a expressão é `(ptr_q - ptr_p) * diametro_esq / 2 + (ptr_q - ptr_p) * profundidade[midpoint_esq] + sl - 2 * consultar_caminho_soma(midpoint_esq)`
        long long current_sl_sum = 0; // Para sumarizar profundidade[md[j]] + comprimentos_diametros[j]/2
        for(int j = ptr_p; j < ptr_q; ++j) {
            current_sl_sum += (long long)profundidade[midpoints_diametros[j]] + (comprimentos_diametros[j] / 2);
        }
        
        soma_diametros_globais += (long long)(ptr_q - ptr_p) * (diametro_esq / 2);
        soma_diametros_globais += (long long)(ptr_q - ptr_p) * profundidade[midpoint_esq];
        soma_diametros_globais += current_sl_sum;
        soma_diametros_globais -= 2 * consultar_caminho_soma(midpoint_esq);


        // Contribuição da região onde c(esquerda) contém c(direita) (ptr_q a r)
        // Para cada j em [ptr_q, r], c(info_esquerda) contém c(midpoints_diametros[j])
        // A soma dos diâmetros é diametro_esq. Então (r - ptr_q + 1) * diametro_esq.
        soma_diametros_globais += soma_comprimentos_direita;
    }
    // Limpar as BITs para os nós que foram adicionados e removidos da região do meio
    for (int j = ptr_p; j < ptr_q; ++j) {
        adicionar_no_caminho(midpoints_diametros[j], -1); // Desfaz as operações para o próximo passo
    }

    resolver_subproblema(l, mid);
    resolver_subproblema(mid + 1, r);
}

int main() {
    std::ios_base::sync_with_stdio(false);
    std::cin.tie(NULL);
    std::cout.tie(NULL);

    std::cin >> num_nos_original;

    // A árvore tem N nós. Para cada aresta (u,v), um nó (N+i) é adicionado
    // A nova árvore terá N + (N-1) = 2N-1 nós.
    // As distâncias originais são 1, mas após a transformação, cada aresta original se torna 2.
    // Ex: u -- v (dist 1) torna-se u -- (N+i) -- v (dist 1+1=2)
    // Isso garante que os midpoints de qualquer caminho cairão em um nó existente.
    // profundidade[0] = -1 é uma convenção para pai da raiz.
    profundidade[0] = -1;
    for (int i = 1; i < num_nos_original; ++i) {
        int u, v;
        std::cin >> u >> v;
        // Adiciona um nó intermediário para cada aresta
        adicionar_aresta(u, num_nos_original + i);
        adicionar_aresta(v, num_nos_original + i);
    }
    
    // Total de nós na árvore expandida é num_nos_original + (num_nos_original - 1)
    // O nó raiz é 1.
    // num_nos_total_expandida = num_nos_original + (num_nos_original - 1); // máximo possível é 2*N-1
    int total_nos_expandida = 2 * num_nos_original - 1; // Para o BIT tamanho maximo

    // Inicializa a BIT com o tamanho máximo do dfn (que pode ir até total_nos_expandida)
    bit_soma_valores = BIT(total_nos_expandida);
    bit_soma_indices_x_valores = BIT(total_nos_expandida);

    dfs1_hld(1, 0, 0); // Raiz 1, pai 0, profundidade 0
    dfs2_hld(1, 1);    // Raiz 1, cabeça de caminho 1

    resolver_subproblema(1, num_nos_original); // Resolvendo para os nós originais de 1 a N
    
    // Cada diâmetro é contado uma vez. A transformação de arestas duplica as distâncias.
    // Então, a soma final dos comprimentos dos diâmetros é dividida por 2.
    std::cout << soma_diametros_globais / 2 << "\n";

    return 0;
}

Tags: programação dinâmica permutações inversões Frações Polinômios Ciclótomicos

Publicado em 8-19 22:52