#1181460 (C++23) No.3663 LCM Decomposition

• 問題: No.3663 LCM Decomposition
• ユーザー: noya2
• 提出日時: 2026-08-30 15:59:39
• 言語: C++23 (gcc 15.3.0 + boost 1.92.0)
• 結果: AC
• 実行時間: 235 ms / 1,000 ms
• コード長: 17,351 bytes
• コンパイル時間: 3,024 ms
• コンパイル使用メモリ: 360,804 KB
• 実行使用メモリ: 6,272 KB
• 最終ジャッジ日時: 2026-08-30 15:59:46
• 合計ジャッジ時間: 5,697 ms
• ジャッジサーバーID （参考情報）: judge3_0 / judge1_0

ファイルパターン	結果
sample	AC * 1
other	AC * 14

ソースコード

#line 2 "/Users/noya2/Desktop/Noya2_library/template/template.hpp"
using namespace std;

#include<bits/stdc++.h>
#line 1 "/Users/noya2/Desktop/Noya2_library/template/inout_old.hpp"
namespace noya2 {

template <typename T, typename U>
ostream &operator<<(ostream &os, const pair<T, U> &p){
    os << p.first << " " << p.second;
    return os;
}
template <typename T, typename U>
istream &operator>>(istream &is, pair<T, U> &p){
    is >> p.first >> p.second;
    return is;
}

template <typename T>
ostream &operator<<(ostream &os, const vector<T> &v){
    int s = (int)v.size();
    for (int i = 0; i < s; i++) os << (i ? " " : "") << v[i];
    return os;
}
template <typename T>
istream &operator>>(istream &is, vector<T> &v){
    for (auto &x : v) is >> x;
    return is;
}

void in() {}
template <typename T, class... U>
void in(T &t, U &...u){
    cin >> t;
    in(u...);
}

void out() { cout << "\n"; }
template <typename T, class... U, char sep = ' '>
void out(const T &t, const U &...u){
    cout << t;
    if (sizeof...(u)) cout << sep;
    out(u...);
}

template<typename T>
void out(const vector<vector<T>> &vv){
    int s = (int)vv.size();
    for (int i = 0; i < s; i++) out(vv[i]);
}

struct IoSetup {
    IoSetup(){
        cin.tie(nullptr);
        ios::sync_with_stdio(false);
        cout << fixed << setprecision(15);
        cerr << fixed << setprecision(7);
    }
} iosetup_noya2;

} // namespace noya2
#line 1 "/Users/noya2/Desktop/Noya2_library/template/const.hpp"
namespace noya2{

const int iinf = 1'000'000'007;
const long long linf = 2'000'000'000'000'000'000LL;
const long long mod998 =  998244353;
const long long mod107 = 1000000007;
const long double pi = 3.14159265358979323;
const vector<int> dx = {0,1,0,-1,1,1,-1,-1};
const vector<int> dy = {1,0,-1,0,1,-1,-1,1};
const string ALP = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
const string alp = "abcdefghijklmnopqrstuvwxyz";
const string NUM = "0123456789";

void yes(){ cout << "Yes\n"; }
void no(){ cout << "No\n"; }
void YES(){ cout << "YES\n"; }
void NO(){ cout << "NO\n"; }
void yn(bool t){ t ? yes() : no(); }
void YN(bool t){ t ? YES() : NO(); }

} // namespace noya2
#line 2 "/Users/noya2/Desktop/Noya2_library/template/utils.hpp"

#line 6 "/Users/noya2/Desktop/Noya2_library/template/utils.hpp"

namespace noya2{

unsigned long long inner_binary_gcd(unsigned long long a, unsigned long long b){
    if (a == 0 || b == 0) return a + b;
    int n = __builtin_ctzll(a); a >>= n;
    int m = __builtin_ctzll(b); b >>= m;
    while (a != b) {
        int mm = __builtin_ctzll(a - b);
        bool f = a > b;
        unsigned long long c = f ? a : b;
        b = f ? b : a;
        a = (c - b) >> mm;
    }
    return a << std::min(n, m);
}

template<typename T> T gcd_fast(T a, T b){ return static_cast<T>(inner_binary_gcd(std::abs(a),std::abs(b))); }

long long sqrt_fast(long long n) {
    if (n <= 0) return 0;
    long long x = sqrt(n);
    while ((x + 1) * (x + 1) <= n) x++;
    while (x * x > n) x--;
    return x;
}

template<typename T> T floor_div(const T n, const T d) {
    assert(d != 0);
    return n / d - static_cast<T>((n ^ d) < 0 && n % d != 0);
}

template<typename T> T ceil_div(const T n, const T d) {
    assert(d != 0);
    return n / d + static_cast<T>((n ^ d) >= 0 && n % d != 0);
}

template<typename T> void uniq(std::vector<T> &v){
    std::sort(v.begin(),v.end());
    v.erase(unique(v.begin(),v.end()),v.end());
}

template <typename T, typename U> inline bool chmin(T &x, U y) { return (y < x) ? (x = y, true) : false; }

template <typename T, typename U> inline bool chmax(T &x, U y) { return (x < y) ? (x = y, true) : false; }

template<typename T> inline bool range(T l, T x, T r){ return l <= x && x < r; }

} // namespace noya2
#line 8 "/Users/noya2/Desktop/Noya2_library/template/template.hpp"

#define rep(i,n) for (int i = 0; i < (int)(n); i++)
#define repp(i,m,n) for (int i = (m); i < (int)(n); i++)
#define reb(i,n) for (int i = (int)(n-1); i >= 0; i--)
#define all(v) (v).begin(),(v).end()

using ll = long long;
using ld = long double;
using uint = unsigned int;
using ull = unsigned long long;
using pii = pair<int,int>;
using pll = pair<ll,ll>;
using pil = pair<int,ll>;
using pli = pair<ll,int>;

namespace noya2{

/*　~ (. _________ . /)　*/

}

using namespace noya2;


#line 2 "c.cpp"

#line 2 "/Users/noya2/Desktop/Noya2_library/math/factorize.hpp"

#line 6 "/Users/noya2/Desktop/Noya2_library/math/factorize.hpp"
#include <initializer_list>
#line 10 "/Users/noya2/Desktop/Noya2_library/math/factorize.hpp"

namespace fast_factorize {

/*
    See : https://judge.yosupo.jp/submission/189742
*/

// ---- gcd ----

uint64_t gcd_stein_impl( uint64_t x, uint64_t y ) {
    if( x == y ) { return x; }
    const uint64_t a = y - x;
    const uint64_t b = x - y;
    const int n = __builtin_ctzll( b );
    const uint64_t s = x < y ? a : b;
    const uint64_t t = x < y ? x : y;
    return gcd_stein_impl( s >> n, t );
}

uint64_t gcd_stein( uint64_t x, uint64_t y ) {
    if( x == 0 ) { return y; }
    if( y == 0 ) { return x; }
    const int n = __builtin_ctzll( x );
    const int m = __builtin_ctzll( y );
    return gcd_stein_impl( x >> n, y >> m ) << ( n < m ? n : m );
}

// ---- is_prime ----

uint64_t mod_pow( uint64_t x, uint64_t y, uint64_t mod ) {
    uint64_t ret = 1;
    uint64_t acc = x;
    for( ; y; y >>= 1 ) {
        if( y & 1 ) {
            ret = __uint128_t(ret) * acc % mod;
        }
        acc = __uint128_t(acc) * acc % mod;
    }
    return ret;
}

bool miller_rabin( uint64_t n, const std::initializer_list<uint64_t>& as ) {
    return std::all_of( as.begin(), as.end(), [n]( uint64_t a ) {
        if( n <= a ) { return true; }

        int e = __builtin_ctzll( n - 1 );
        uint64_t z = mod_pow( a, ( n - 1 ) >> e, n );
        if( z == 1 || z == n - 1 ) { return true; }

        while( --e ) {
            z = __uint128_t(z) * z % n;
            if( z == 1 ) { return false; }
            if( z == n - 1 ) { return true; }
        }

        return false;
    });
}

bool is_prime( uint64_t n ) {
    if( n == 2 ) { return true; }
    if( n % 2 == 0 ) { return false; }
    if( n < 4759123141 ) { return miller_rabin( n, { 2, 7, 61 } ); }
    return miller_rabin( n, { 2, 325, 9375, 28178, 450775, 9780504, 1795265022 } );
}

// ---- Montgomery ----

class Montgomery {
    uint64_t mod;
    uint64_t R;
public:
    Montgomery( uint64_t n ) : mod(n), R(n) {
       for( size_t i = 0; i < 5; ++i ) {
          R *= 2 - mod * R;
       }
    }

    uint64_t fma( uint64_t a, uint64_t b, uint64_t c ) const {
        const __uint128_t d = __uint128_t(a) * b;
        const uint64_t    e = c + mod + ( d >> 64 );
        const uint64_t    f = uint64_t(d) * R;
        const uint64_t    g = ( __uint128_t(f) * mod ) >> 64;
        return e - g;
    }

    uint64_t mul( uint64_t a, uint64_t b ) const {
        return fma( a, b, 0 );
    }
};

// ---- Pollard's rho algorithm ----

uint64_t pollard_rho( uint64_t n ) {
    if( n % 2 == 0 ) { return 2; }
    const Montgomery m( n );

    constexpr uint64_t C1 = 1;
    constexpr uint64_t C2 = 2;
    constexpr uint64_t M = 512;

    uint64_t Z1 = 1;
    uint64_t Z2 = 2;
retry:
    uint64_t z1 = Z1;
    uint64_t z2 = Z2;
    for( size_t k = M; ; k *= 2 ) {
        const uint64_t x1 = z1 + n;
        const uint64_t x2 = z2 + n;
        for( size_t j = 0; j < k; j += M ) {
            const uint64_t y1 = z1;
            const uint64_t y2 = z2;

            uint64_t q1 = 1;
            uint64_t q2 = 2;
            z1 = m.fma( z1, z1, C1 );
            z2 = m.fma( z2, z2, C2 );
            for( size_t i = 0; i < M; ++i ) {
                const uint64_t t1 = x1 - z1;
                const uint64_t t2 = x2 - z2;
                z1 = m.fma( z1, z1, C1 );
                z2 = m.fma( z2, z2, C2 );
                q1 = m.mul( q1, t1 );
                q2 = m.mul( q2, t2 );
            }
            q1 = m.mul( q1, x1 - z1 );
            q2 = m.mul( q2, x2 - z2 );

            const uint64_t q3 = m.mul( q1, q2 );
            const uint64_t g3 = gcd_stein( n, q3 );
            if( g3 == 1 ) { continue; }
            if( g3 != n ) { return g3; }

            const uint64_t g1 = gcd_stein( n, q1 );
            const uint64_t g2 = gcd_stein( n, q2 );

            const uint64_t C = g1 != 1 ? C1 : C2;
            const uint64_t x = g1 != 1 ? x1 : x2;
            uint64_t       z = g1 != 1 ? y1 : y2;
            uint64_t       g = g1 != 1 ? g1 : g2;

            if( g == n ) {
                do {
                    z = m.fma( z, z, C );
                    g = gcd_stein( n, x - z );
                } while( g == 1 );
            }
            if( g != n ) {
                return g;
            }

            Z1 += 2;
            Z2 += 2;
            goto retry;
        }
    }
}

void factorize_impl( uint64_t n, std::vector<uint64_t>& ret ) {
    if( n <= 1 ) { return; }
    if( is_prime( n ) ) { ret.push_back( n ); return; }

    const uint64_t p = pollard_rho( n );

    factorize_impl( p, ret );
    factorize_impl( n / p, ret );
}

std::vector<uint64_t> factorize( uint64_t n ) {
    std::vector<uint64_t> ret;
    factorize_impl( n, ret );
    std::sort( ret.begin(), ret.end() );
    return ret;
}

} // namespace fast_factorize

namespace noya2 {

std::vector<std::pair<long long, int>> factorize(long long n){
    std::vector<std::pair<long long, int>> ans;
    auto ps = fast_factorize::factorize(n);
    int sz = ps.size();
    for (int l = 0, r = 0; l < sz; l = r){
        while (r < sz && ps[l] == ps[r]) r++;
        ans.emplace_back(ps[l], r-l);
    }
    return ans;
}

std::vector<long long> divisors(long long n){
    auto ps = fast_factorize::factorize(n);
    int sz = ps.size();
    std::vector<long long> ans = {1};
    for (int l = 0, r = 0; l < sz; l = r){
        while (r < sz && ps[l] == ps[r]) r++;
        int e = r - l;
        int len = ans.size();
        ans.reserve(len*(e+1));
        long long mul = ps[l];
        while (true){
            for (int i = 0; i < len; i++){
                ans.emplace_back(ans[i]*mul);
            }
            if (--e == 0) break;
            mul *= ps[l];
        }
    }
    return ans;
}

std::vector<long long> divisors(const std::vector<std::pair<long long, int>> &pes){
    std::vector<long long> ans = {1};
    for (auto [p, e] : pes){
        int len = ans.size();
        ans.reserve(len*(e+1));
        long long mul = p;
        while (true){
            for (int i = 0; i < len; i++){
                ans.emplace_back(ans[i]*mul);
            }
            if (--e == 0) break;
            mul *= p;
        }
    }
    return ans;
}

bool is_prime(long long n){
    if (n <= 1) return false;
    return fast_factorize::is_prime(n);
}

} // namespace noya2
#line 2 "/Users/noya2/Desktop/Noya2_library/math/bitwise_convolution.hpp"

#line 6 "/Users/noya2/Desktop/Noya2_library/math/bitwise_convolution.hpp"

namespace noya2 {

// f mapsto F
// F(S) = sum[T subset S] f(T)
//      = f(S) o ...
template<typename T, class Op = std::plus<T>>
void bitwise_zeta_transform(std::vector<T> &a){
    uint32_t sz = a.size();
    assert(std::has_single_bit(sz));
    Op op{};
    for (uint32_t p2 = 1; p2 < sz; p2 <<= 1){
        for (uint32_t u = 0; u < sz; u += (p2 << 1)){
            for (uint32_t d = 0; d < p2; d++){
                a[d ^ u ^ p2] = op(a[d ^ u ^ p2], a[d ^ u]);
            }
        }
    }
}


// inverse transform of zeta : F mapsto f
template<typename T, class Opinv = std::minus<T>>
void bitwise_mobius_transform(std::vector<T> &a){
    uint32_t sz = a.size();
    assert(std::has_single_bit(sz));
    Opinv opinv{};
    for (uint32_t p2 = 1; p2 < sz; p2 <<= 1){
        for (uint32_t u = 0; u < sz; u += (p2 << 1)){
            for (uint32_t d = 0; d < p2; d++){
                a[d ^ u ^ p2] = opinv(a[d ^ u ^ p2], a[d ^ u]);
            }
        }
    }
}

template<typename T>
std::vector<T> bitwise_or_convolution(std::vector<T> a, std::vector<T> b){
    assert(a.size() == b.size());
    uint32_t sz = a.size();
    bitwise_zeta_transform<T, std::plus<T>>(a);
    bitwise_zeta_transform<T, std::plus<T>>(b);
    for (uint32_t i = 0; i < sz; i++){
        a[i] *= b[i];
    }
    bitwise_mobius_transform<T, std::minus<T>>(a);
    return a;
}

template<typename T>
std::vector<T> bitwise_and_convolution(std::vector<T> a, std::vector<T> b){
    assert(a.size() == b.size());
    uint32_t sz = a.size();
    std::reverse(a.begin(), a.end());
    std::reverse(b.begin(), b.end());
    bitwise_zeta_transform<T, std::plus<T>>(a);
    bitwise_zeta_transform<T, std::plus<T>>(b);
    for (uint32_t i = 0; i < sz; i++){
        a[i] *= b[i];
    }
    bitwise_mobius_transform<T, std::minus<T>>(a);
    std::reverse(a.begin(), a.end());
    return a;
}

template<typename T>
void hadamard_transform(std::vector<T> &a){
    uint32_t sz = a.size();
    assert(std::has_single_bit(sz));
    for (uint32_t p2 = 1; p2 < sz; p2 <<= 1){
        for (uint32_t u = 0; u < sz; u += (p2 << 1)){
            for (uint32_t d = 0; d < p2; d++){
                T s = a[d ^ u], t = a[d ^ u ^ p2];
                a[d ^ u] = s + t;
                a[d ^ u ^ p2] = s - t;
            }
        }
    }
}

template<typename T>
std::vector<T> bitwise_xor_convolution(std::vector<T> a, std::vector<T> b){
    assert(a.size() == b.size());
    uint32_t sz = a.size();
    hadamard_transform(a);
    hadamard_transform(b);
    T iv = T(1) / T(sz);
    for (uint32_t i = 0; i < sz; i++){
        a[i] *= b[i] * iv;
    }
    hadamard_transform(a);
    return a;
}

template<typename T, int LOG>
std::vector<T> subset_convolution(const std::vector<T> &a, const std::vector<T> &b){
    assert(a.size() == b.size());
    using fps = std::array<T, LOG + 1>;
    constexpr int s = LOG;
    uint32_t sz = a.size();
    auto bitwise_or_zeta_transform = [&]<typename is_zeta>(std::vector<fps> &c, is_zeta z){
        for (uint32_t p2 = 1; p2 < sz; p2 <<= 1){
            for (uint32_t u = 0; u < sz; u += (p2 << 1)){
                for (uint32_t d = 0; d < p2; d++){
                    for (int i = 0; i <= s; i++){
                        if constexpr (is_zeta::value){
                            c[d ^ u ^ p2][i] += c[d ^ u][i];
                        }
                        else {
                            c[d ^ u ^ p2][i] -= c[d ^ u][i];
                        }
                    }
                }
            }
        }
    };
    auto lift = [&](const std::vector<T> &c){
        std::vector<fps> C(sz, fps{});
        for (uint32_t i = 0; i < sz; i++){
            C[i][std::popcount(i)] = c[i];
        }
        return C;
    };
    auto A = lift(a), B = lift(b);
    bitwise_or_zeta_transform(A, std::true_type{});
    bitwise_or_zeta_transform(B, std::true_type{});
    for (uint32_t i = 0; i < sz; i++){
        fps c = {};
        for (int x = 0; x <= s; x++){
            for (int y = 0; y <= s - x; y++){
                c[x + y] += A[i][x] * B[i][y];
            }
        }
        std::swap(A[i], c);
    }
    bitwise_or_zeta_transform(A, std::false_type{});
    std::vector<T> ab(sz);
    for (uint32_t i = 0; i < sz; i++){
        ab[i] = A[i][std::popcount(i)];
    }
    return ab;
}

} // namespace noya2
#line 5 "c.cpp"

void solve(){
    ll n; in(n);
    auto pes = factorize(n);
    auto ds = divisors(pes);
    ll l, r; in(l,r);
    int okcnt = ranges::count_if(ds, [&](ll x){
        return l <= x && x <= r;
    });
    if (okcnt < 3){
        out(-1);
        return ;
    }
    int sz = pes.size();
    vector<ll> pws(sz);
    rep(i,sz){
        auto [p, e] = pes[i];
        ll prd = 1;
        while (e--){
            prd *= p;
        }
        pws[i] = prd;
    }
    vector<ll> inner(1<<sz);
    rep(s,1<<sz){
        ll prd = 1;
        rep(i,sz){
            if (s >> i & 1){
                prd *= pws[i];
            }
        }
        for (ll d : ds){
            if (l <= d && d <= r && d % prd == 0){
                inner[s] = d;
            }
        }
    }
    vector<bool> dp1(1<<sz), dp2(1<<sz);
    {
        vector<uint> a(1<<sz);
        rep(s,1<<sz){
            dp1[s] = (inner[s] != 0);
            a[s] = uint(inner[s] != 0);
        }
        auto aa = bitwise_or_convolution(a, a);
        rep(s,1<<sz){
            dp2[s] = (aa[s] != 0);
        }
    }
    auto [st, u] = [&]{
        int tot = (1<<sz)-1;
        rep(s,1<<sz){
            if (dp2[s] && dp1[tot^s]){
                return pair(s,tot^s);
            }
        }
        return pair(-1,-1);
    }();
    if (u == -1){
        out(-1);
        return ;
    }
    auto [s, t] = [&]{
        rep(sub,1<<sz){
            if (dp1[sub] && dp1[st^sub]){
                return pair(sub,st^sub);
            }
        }
        abort();
    }();
    vector<ll> ans = {inner[s],inner[t],inner[u]};
    uniq(ans);
    set<ll> sts(all(ans));
    while (ans.size() < 3u){
        for (ll d : ds){
            if (l <= d && d <= r){
                if (!sts.contains(d)){
                    ans.emplace_back(d);
                    sts.insert(d);
                    break;
                }
            }
        }
    }
    sort(all(ans));
    assert(ans.size() == 3u);
    // out(n);
    // out(ds);
    // out(inner);
    // out(s,t,u);
    out(ans);
}

int main(){
    int t = 1; in(t);
    while (t--) { solve(); }
}