algorithm - セット {1,2,3,..,n} の互いに素な部分集合の数を計算する方法

Question

私はこのタスクを解決しています (問題 I) . ステートメントは次のとおりです。

{1, 2, 3, ..., n}セットの互いに素な部分集合はいくつありますか? 整数の集合は、要素の 2 つすべてが互いに素である場合、互いに素であると呼ばれます。最大公約数が 1 の場合、2 つの整数は互いに素です。

入力

入力の最初の行には、2 つの整数nとm( 1 <= n <= 3000, 1 <= m <= 10^9 + 9)が含まれます。

出力

{1, 2, 3, ..., n}moduloの余素部分集合の数を出力しmます。

例

入力: 4 7 出力: 5

、、、、、 、 、、、の12の互いに素な部分集合があり{1,2,3,4}ます。{}{1}{2}{3}{4}{1,2}{1,3}{1,4}{2,3}{3,4}{1,2,3}{1,3,4}

---

素数を使えば解けると思います。（各素数を使用したかどうかを追跡しています）..しかし、よくわかりません。

このタスクを解決するためのヒントを得ることができますか?

• このシーケンスは、http: //oeis.org/A084422にあります。

Answer 1

これは、シーケンスの最初の 200 要素を 1秒未満で通過する答えであり、正しい答え 200 → 374855124868136960 が得られます。最適化 (編集 1 を参照) を使用すると、最初の 500 エントリを 90 秒未満で計算できます。これは迅速です - - @David Eisenstatの答えは、開発できればより良いものになる可能性があります。私自身の独自の回答を含め、これまでに与えられたアルゴリズムに対して別のアプローチを取っていると思うので、個別に投稿します。

最適化した後、私は本当にグラフの問題をコーディングしていることに気付き、ソリューションをグラフの実装として書き直しました (編集 2 を参照)。グラフの実装は、いくつかの最適化を可能にし、はるかに洗練されており、1 桁以上高速であり、スケーリングも優れていますf(600)。27 秒と比較して、1.5 秒で計算されます。

ここでの主なアイデアは、再帰関係を使用することです。任意のセットについて、基準を満たすサブセットの数は次の合計です。

• 1 つの要素が削除されたサブセットの数。と
• その要素が確実に含まれているサブセットの数。

2 番目のケースでは、要素が確実に含まれている場合、それと互いに素でない他の要素はすべて削除する必要があります。

効率の問題:

• 最大の要素を削除して、その要素が他のすべての要素とすでに互いに素である可能性を最大化することにしました。この場合、再帰呼び出しを 2 回ではなく 1 回だけ行う必要があります。
• キャッシング/メモ化が役立ちます。

以下のコード。

#include <cassert>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include <iostream>
#include <ctime>

const int PRIMES[] = // http://rlrr.drum-corps.net/misc/primes1.shtml
    { 2, 3, 5, ...
      ..., 2969, 2971, 2999 };
const int NPRIMES = sizeof(PRIMES) / sizeof(int);

typedef std::set<int> intset;
typedef std::vector<intset> intsetvec;

const int MAXCALC = 200; // answer at http://oeis.org/A084422/b084422.txt
intsetvec primeFactors(MAXCALC +1);

typedef std::vector<int> intvec;

// Caching / memoization
typedef std::map<intvec, double> intvec2dbl;
intvec2dbl set2NumCoPrimeSets;

double NumCoPrimeSets(const intvec& set)
{
    if (set.empty())
        return 1;

    // Caching / memoization
    const intvec2dbl::const_iterator i = set2NumCoPrimeSets.find(set);
    if (i != set2NumCoPrimeSets.end())
        return i->second;

    // Result is the number of coprime sets in:
    //      setA, the set that definitely has the first element of the input present
    // +    setB, the set the doesn't have the first element of the input present

    // Because setA definitely has the first element, we remove elements it isn't coprime with
    // We also remove the first element: as this is definitely present it doesn't make any
    // difference to the number of sets
    intvec setA(set);
    const int firstNum = *setA.begin();
    const intset& factors = primeFactors[firstNum];
    for(int factor : factors) {
        setA.erase(std::remove_if(setA.begin(), setA.end(),
            [factor] (int i) { return i % factor == 0; } ), setA.end());
    }

    // If the first element was already coprime with the rest, then we have setA = setB
    // and we can do a single call (m=2). Otherwise we have two recursive calls.
    double m = 1;
    double c = 0;
    assert(set.size() - setA.size() > 0);
    if (set.size() - setA.size() > 1) {
        intvec setB(set);
        setB.erase(setB.begin());
        c = NumCoPrimeSets(setB);
    }
    else {
        // first elt coprime with rest
        m = 2;
    }
    const double numCoPrimeSets = m * NumCoPrimeSets(setA) + c;

    // Caching / memoization
    set2NumCoPrimeSets.insert(intvec2dbl::value_type(set, numCoPrimeSets));
    return numCoPrimeSets;
}


int main(int argc, char* argv[])
{
    // Calculate prime numbers that factor into each number upto MAXCALC
    primeFactors[1].insert(1); // convenient
    for(int i=2; i<=MAXCALC; ++i) {
        for(int j=0; j<NPRIMES; ++j) {
            if (i % PRIMES[j] == 0) {
                primeFactors[i].insert(PRIMES[j]);
            }
        }
    }

    const clock_t start = clock();

    for(int n=1; n<=MAXCALC; ++n) {
        intvec v;
        for(int i=n; i>0; --i) { // reverse order to reduce recursion
            v.push_back(i);
        }
        const clock_t now = clock();
        const clock_t ms = now - start;
        const double numCoPrimeSubsets = NumCoPrimeSets(v);
        std::cout << n << ", " << std::fixed << numCoPrimeSubsets << ", " << ms << "\n";
    }

    return 0;
}

時間特性は、私の最初の回答よりもずっと良く見えます。しかし、それでも 5 秒で 3000 には達しません。

---

編集 1

このメソッドに対して実行できる興味深い最適化がいくつかあります。全体として、これにより、より大きなn.

• セット内のすでに互いに素であるすべての数値は、1 つの前処理ステップで削除できます。m数値が削除された場合、元のセットには、削減されたものよりも 2 ^m因子倍の組み合わせがあります (それぞれの互いに素であるため、または他の要素とは独立してセット外に)。
• 最も重要なことは、セット内の任意の場所にある削除する要素を選択できることです。最も接続された要素を削除するのが最も効果的であることがわかりました。
• 以前に使用された再帰的な関係を一般化して、削除されたすべての要素が同じ素因数を持つ複数の要素を削除できます。たとえば、 set{2, 3, 15, 19, 45}の場合、数値 15 と 45 は同じ素因数 3 と 5 を持ちます。一度に2 つの数値が削除されるため、サブセットの数{2, 3, 15, 19, 45}= 15 または 45 のいずれかの組み合わせの数の2 倍( 15 または 45 のいずれかが存在する場合、3 が存在しない必要があるため、セット{2, 19}) + 15 および 45 が存在しない場合のサブセットの数 (セット の場合{2, 3, 19})
• number 型に使用shortすると、パフォーマンスが約 10% 向上しました。
• 最後に、セットを標準化することでより良いキャッシュ ヒットを得ることを期待して、セットを同等の素因数を持つセットに変換することも可能です。たとえば、{ 3, 9, 15}は と同等 (同形)2, 4, 6です。これは最も急進的なアイデアでしたが、パフォーマンスへの影響はおそらく最も少なかったでしょう。

具体的な例を示した方が理解しやすいかもしれません。私はセット {1..12} を選択しました。これは、どのように機能するかを理解するのに十分な大きさですが、理解できるように十分に小さいものです。

NumCoPrimeSets({ 1 2 3 4 5 6 7 8 9 10 11 12 })
Removed 3 coprimes, giving set { 2 3 4 5 6 8 9 10 12 } multiplication factor now 8
Removing the most connected number 12 with 8 connections
To get setA, remove all numbers which have *any* of the prime factors { 2 3 }
setA = { 5 }
To get setB, remove 2 numbers which have *exactly* the prime factors { 2 3 }
setB = { 2 3 4 5 8 9 10 }
**** Recursing on 2 * NumCoPrimeSets(setA) + NumCoPrimeSets(setB)

NumCoPrimeSets({ 5 })
Base case return the multiplier, which is 2
NumCoPrimeSets({ 2 3 4 5 8 9 10 })
Removing the most connected number 10 with 4 connections
To get setA, remove all numbers which have *any* of the prime factors { 2 5 }
setA = { 3 9 }
To get setB, remove 1 numbers which have *exactly* the prime factors { 2 5 }
setB = { 2 3 4 5 8 9 }
**** Recursing on 1 * NumCoPrimeSets(setA) + NumCoPrimeSets(setB)

NumCoPrimeSets({ 3 9 })
Transformed 2 primes, giving new set { 2 4 }
Removing the most connected number 4 with 1 connections
To get setA, remove all numbers which have *any* of the prime factors { 2 }
setA = { }
To get setB, remove 2 numbers which have *exactly* the prime factors { 2 }
setB = { }
**** Recursing on 2 * NumCoPrimeSets(setA) + NumCoPrimeSets(setB)

NumCoPrimeSets({ })
Base case return the multiplier, which is 1
NumCoPrimeSets({ })
Base case return the multiplier, which is 1
**** Returned from recursing on 2 * NumCoPrimeSets({ }) + NumCoPrimeSets({ })
Caching for{ 2 4 }: 3 = 2 * 1 + 1
Returning for{ 3 9 }: 3 = 1 * 3

NumCoPrimeSets({ 2 3 4 5 8 9 })
Removed 1 coprimes, giving set { 2 3 4 8 9 } multiplication factor now 2
Removing the most connected number 8 with 2 connections
To get setA, remove all numbers which have *any* of the prime factors { 2 }
setA = { 3 9 }
To get setB, remove 3 numbers which have *exactly* the prime factors { 2 }
setB = { 3 9 }
**** Recursing on 3 * NumCoPrimeSets(setA) + NumCoPrimeSets(setB)

NumCoPrimeSets({ 3 9 })
Transformed 2 primes, giving new set { 2 4 }
Cache hit, returning 3 = 1 * 3
NumCoPrimeSets({ 3 9 })
Transformed 2 primes, giving new set { 2 4 }
Cache hit, returning 3 = 1 * 3
**** Returned from recursing on 3 * NumCoPrimeSets({ 3 9 }) + NumCoPrimeSets({ 3 9 })
Caching for{ 2 3 4 8 9 }: 12 = 3 * 3 + 3
Returning for{ 2 3 4 5 8 9 }: 24 = 2 * 12

**** Returned from recursing on 1 * NumCoPrimeSets({ 3 9 }) + NumCoPrimeSets({ 2 3 4 5 8 9 })
Caching for{ 2 3 4 5 8 9 10 }: 27 = 1 * 3 + 24
Returning for{ 2 3 4 5 8 9 10 }: 27 = 1 * 27

**** Returned from recursing on 2 * NumCoPrimeSets({ 5 }) + NumCoPrimeSets({ 2 3 4 5 8 9 10 })
Caching for{ 2 3 4 5 6 8 9 10 12 }: 31 = 2 * 2 + 27
Returning for{ 1 2 3 4 5 6 7 8 9 10 11 12 }: 248 = 8 * 31

以下のコード

#include <cassert>
#include <vector>
#include <set>
#include <map>
#include <unordered_map>
#include <queue>
#include <algorithm>
#include <fstream>
#include <iostream>
#include <ctime>

typedef short numtype;

const numtype PRIMES[] = // http://rlrr.drum-corps.net/misc/primes1.shtml
    ...
const numtype NPRIMES = sizeof(PRIMES) / sizeof(numtype);

typedef std::set<numtype> numset;
typedef std::vector<numset> numsetvec;

const numtype MAXCALC = 200; // answer at http://oeis.org/A084422/b084422.txt
numsetvec primeFactors(MAXCALC +1);

typedef std::vector<numtype> numvec;

// Caching / memoization
typedef std::map<numvec, double> numvec2dbl;
numvec2dbl set2NumCoPrimeSets;

double NumCoPrimeSets(const numvec& initialSet)
{
    // Preprocessing step: remove numbers which are already coprime
    typedef std::unordered_map<numtype, numvec> num2numvec;
    num2numvec prime2Elts;
    for(numtype num : initialSet) {
        const numset& factors = primeFactors[num];
        for(numtype factor : factors) {
             prime2Elts[factor].push_back(num);
        }
    }

    numset eltsToRemove(initialSet.begin(), initialSet.end());
    typedef std::vector<std::pair<numtype,int>> numintvec;
    numvec primesRemaining;
    for(const num2numvec::value_type& primeElts : prime2Elts) {
        if (primeElts.second.size() > 1) {
            for (numtype num : primeElts.second) {
                eltsToRemove.erase(num);
            }
            primesRemaining.push_back(primeElts.first);
        }
    }

    double mult = pow(2.0, eltsToRemove.size());
    if (eltsToRemove.size() == initialSet.size())
        return mult;

    // Do the removal by creating a new set
    numvec set;
    for(numtype num : initialSet) {
        if (eltsToRemove.find(num) == eltsToRemove.end()) {
            set.push_back(num);
        }
    }

    // Transform to use a smaller set of primes before checking the cache
    // (beta code but it seems to work, mostly!)
    std::sort(primesRemaining.begin(), primesRemaining.end());
    numvec::const_iterator p = primesRemaining.begin();
    for(int j=0; p!= primesRemaining.end() && j<NPRIMES; ++p, ++j) {
        const numtype primeRemaining = *p;
        if (primeRemaining != PRIMES[j]) {
            for(numtype& num : set) {
                while (num % primeRemaining == 0) {
                    num = num / primeRemaining * PRIMES[j];
                }
            }
        }
    }

    // Caching / memoization
    const numvec2dbl::const_iterator i = set2NumCoPrimeSets.find(set);
    if (i != set2NumCoPrimeSets.end())
        return mult * i->second;

    // Remove the most connected number
    typedef std::unordered_map<numtype, int> num2int;
    num2int num2ConnectionCount;
    for(numvec::const_iterator srcIt=set.begin(); srcIt!=set.end(); ++srcIt) {
        const numtype src = *srcIt;
        const numset& srcFactors = primeFactors[src];
        for(numvec::const_iterator tgtIt=srcIt +1; tgtIt!=set.end(); ++tgtIt) {
            for(numtype factor : srcFactors) {
                const numtype tgt = *tgtIt;
                if (tgt % factor == 0) {
                    num2ConnectionCount[src]++;
                    num2ConnectionCount[tgt]++;
                }
            }
        }
    }
    num2int::const_iterator connCountIt = num2ConnectionCount.begin();

    numtype numToErase = connCountIt->first;
    int maxConnCount = connCountIt->second;
    for (; connCountIt!=num2ConnectionCount.end(); ++connCountIt) {
        if (connCountIt->second > maxConnCount || connCountIt->second == maxConnCount && connCountIt->first > numToErase) {
            numToErase = connCountIt->first;
            maxConnCount = connCountIt->second;
        }
    }

    // Result is the number of coprime sets in:
    //      setA, the set that definitely has a chosen element of the input present
    // +    setB, the set the doesn't have the chosen element(s) of the input present

    // Because setA definitely has a chosen element, we remove elements it isn't coprime with
    // We also remove the chosen element(s): as they are definitely present it doesn't make any
    // difference to the number of sets
    numvec setA(set);
    const numset& factors = primeFactors[numToErase];
    for(numtype factor : factors) {
        setA.erase(std::remove_if(setA.begin(), setA.end(),
            [factor] (numtype i) { return i % factor == 0; } ), setA.end());
    }

    // setB: remove all elements which have the same prime factors
    numvec setB(set);
    setB.erase(std::remove_if(setB.begin(), setB.end(),
        [&factors] (numtype i) { return primeFactors[i] == factors; }), setB.end());
    const size_t numEltsWithSamePrimeFactors = (set.size() - setB.size());

    const double numCoPrimeSets = 
        numEltsWithSamePrimeFactors * NumCoPrimeSets(setA) + NumCoPrimeSets(setB);

    // Caching / memoization
    set2NumCoPrimeSets.insert(numvec2dbl::value_type(set, numCoPrimeSets));
    return mult * numCoPrimeSets;
}

int main(int argc, char* argv[])
{
    // Calculate prime numbers that factor into each number upto MAXCALC
    for(numtype i=2; i<=MAXCALC; ++i) {
        for(numtype j=0; j<NPRIMES; ++j) {
            if (i % PRIMES[j] == 0) {
                primeFactors[i].insert(PRIMES[j]);
            }
        }
    }

    const clock_t start = clock();

    std::ofstream fout("out.txt");

    for(numtype n=0; n<=MAXCALC; ++n) {
        numvec v;
        for(numtype i=1; i<=n; ++i) {
            v.push_back(i);
        }
        const clock_t now = clock();
        const clock_t ms = now - start;
        const double numCoPrimeSubsets = NumCoPrimeSets(v);
        fout << n << ", " << std::fixed << numCoPrimeSubsets << ", " << ms << "\n";
        std::cout << n << ", " << std::fixed << numCoPrimeSubsets << ", " << ms << "\n";
    }

    return 0;
}

n=6005分程度で処理可能です。時間は依然として指数関数的に見えますが、50 から 60 程度ごとに 2 倍にnなります。1つだけ計算したグラフをn以下に示します。

---

編集 2

このソリューションは、グラフに関してより自然に実装されます。さらに 2 つの最適化が行われました。

• 最も重要なことは、グラフ G を 2 つの集合 A と B に分割して、A と B の間に接続がない場合、coprimes(G) = coprimes(A) * coprimes(B) です。

• 次に、一連の素因数のすべての数値を 1 つのノードにまとめることができるため、ノードの値はその素因数の数値のカウントになります。

以下のコードでは、Graphクラスは元の隣接行列とノード値をFilteredGraph保持し、クラスは残りのノードの現在のリストを として保持するbitsetため、ノードが削除されると、新しい隣接行列をビット マスキングによって計算できます (比較的再帰で渡すデータはほとんどありません)。

#include "Primes.h"

#include <cassert>
#include <bitset>
#include <vector>
#include <set>
#include <map>
#include <unordered_map>
#include <algorithm>
#include <iostream>
#include <ctime>

// Graph declaration

const int MAXGROUPS = 1462; // empirically determined

class Graph
{
    typedef std::bitset<MAXGROUPS> bitset;
    typedef std::vector<bitset> adjmatrix;
    typedef std::vector<int> intvec;
public:
    Graph(int numNodes)
        : m_nodeValues(numNodes), m_adjMatrix(numNodes) {}
    void SetNodeValue(int i, int v) { m_nodeValues[i] = v; }
    void SetConnection(int i, int j)
    {
        m_adjMatrix[i][j] = true;
        m_adjMatrix[j][i] = true;
    }
    int size() const { return m_nodeValues.size(); }
private:
    adjmatrix m_adjMatrix;
    intvec m_nodeValues;
    friend class FilteredGraph;
};

class FilteredGraph
{
    typedef Graph::bitset bitset;
public:
    FilteredGraph(const Graph* unfiltered);

    int FirstNode() const;
    int RemoveNode(int node);
    void RemoveNodesConnectedTo(int node);

    double RemoveDisconnectedNodes();

    bool AttemptPartition(FilteredGraph* FilteredGraph);

    size_t Hash() const { return std::hash<bitset>()(m_includedNodes); }
    bool operator==(const FilteredGraph& x) const
    { return x.m_includedNodes == m_includedNodes && x.m_unfiltered == m_unfiltered; }

private:
    bitset RawAdjRow(int i) const {
        return m_unfiltered->m_adjMatrix[i];
    }
    bitset AdjRow(int i) const {
        return RawAdjRow(i) & m_includedNodes;
    }
    int NodeValue(int i) const {
        return m_unfiltered->m_nodeValues[i];
    }

    const Graph* m_unfiltered;
    bitset m_includedNodes;
};

// Cache
namespace std {
    template<>
    class hash<FilteredGraph> {
    public:
        size_t operator()(const FilteredGraph & x) const { return x.Hash(); }
    };
}
typedef std::unordered_map<FilteredGraph, double> graph2double;
graph2double cache;

// MAIN FUNCTION

double NumCoPrimesSubSets(const FilteredGraph& graph)
{
    graph2double::const_iterator cacheIt = cache.find(graph);
    if (cacheIt != cache.end())
        return cacheIt->second;

    double rc = 1;

    FilteredGraph A(graph);
    FilteredGraph B(graph);

    if (A.AttemptPartition(&B)) {
        rc = NumCoPrimesSubSets(A);
        A = B;
    }

    const int nodeToRemove = A.FirstNode();
    if (nodeToRemove < 0) // empty graph
        return 1;

    // Graph B is the graph with a node removed
    B.RemoveNode(nodeToRemove);

    // Graph A is the graph with the node present -- and hence connected nodes removed
    A.RemoveNodesConnectedTo(nodeToRemove);
    // The number of numbers in the node is the number of times it can be reused
    const double removedNodeValue = A.RemoveNode(nodeToRemove);

    const double A_disconnectedNodesMult = A.RemoveDisconnectedNodes();
    const double B_disconnectedNodesMult = B.RemoveDisconnectedNodes();

    const double A_coprimes = NumCoPrimesSubSets(A);
    const double B_coprimes = NumCoPrimesSubSets(B);

    rc *= removedNodeValue * A_disconnectedNodesMult * A_coprimes
                           + B_disconnectedNodesMult * B_coprimes;

    cache.insert(graph2double::value_type(graph, rc));
    return rc;
}

// Program entry point
int Sequence2Graph(Graph** ppGraph, int n);

int main(int argc, char* argv[])
{
    const clock_t start = clock();

    int n=800; // runs in approx 6s on my machine

    Graph* pGraph = nullptr;

    const int coPrimesRemoved = Sequence2Graph(&pGraph, n);
    const double coPrimesMultiplier = pow(2,coPrimesRemoved);
    const FilteredGraph filteredGraph(pGraph);
    const double numCoPrimeSubsets = coPrimesMultiplier * NumCoPrimesSubSets(filteredGraph);
    delete pGraph;
    cache.clear(); // as it stands the cache can't cope with other Graph objects, e.g. for other n

    const clock_t now = clock();
    const clock_t ms = now - start;
    std::cout << n << ", " << std::fixed << numCoPrimeSubsets << ", " << ms << "\n";

    return 0;
}

// Graph implementation

FilteredGraph::FilteredGraph(const Graph* unfiltered)
    : m_unfiltered(unfiltered)
{
    for(int i=0; i<m_unfiltered->size(); ++i) {
        m_includedNodes.set(i);
    }
}

int FilteredGraph::FirstNode() const
{
    int firstNode=0;
    for(; firstNode<m_unfiltered->size() && !m_includedNodes.test(firstNode); ++firstNode) {
    }
    if (firstNode == m_unfiltered->size())
        return -1;
    return firstNode;
}

int FilteredGraph::RemoveNode(int node)
{
    m_includedNodes.set(node, false);
    return NodeValue(node);
}

void FilteredGraph::RemoveNodesConnectedTo(const int node)
{
    const bitset notConnected = ~RawAdjRow(node);
    m_includedNodes &= notConnected;
}

double FilteredGraph::RemoveDisconnectedNodes()
{
    double mult = 1.0;
    for(int i=0; i<m_unfiltered->size(); ++i) {
        if (m_includedNodes.test(i)) {
            const int conn = AdjRow(i).count();
            if (conn == 0) {
                m_includedNodes.set(i, false);;
                mult *= (NodeValue(i) +1);
            }
        }
    }
    return mult;
}

bool FilteredGraph::AttemptPartition(FilteredGraph* pOther)
{
    typedef std::vector<int> intvec;
    intvec includedNodesCache;
    includedNodesCache.reserve(m_unfiltered->size());
    for(int i=0; i<m_unfiltered->size(); ++i) {
        if (m_includedNodes.test(i)) {
            includedNodesCache.push_back(i);
        }
    }

    if (includedNodesCache.empty())
        return false;

    const int startNode= includedNodesCache[0];

    bitset currFoundNodes;
    currFoundNodes.set(startNode);
    bitset foundNodes;
    do {
        foundNodes |= currFoundNodes;
        bitset newFoundNodes;
        for(int i : includedNodesCache) {
            if (currFoundNodes.test(i)) {
                newFoundNodes |= AdjRow(i);
            }
        }
        newFoundNodes &= ~ foundNodes;
        currFoundNodes = newFoundNodes;
    } while(currFoundNodes.count() > 0);

    const size_t foundCount = foundNodes.count();
    const size_t thisCount = m_includedNodes.count();

    const bool isConnected = foundCount == thisCount;

    if (!isConnected) {
        if (foundCount < thisCount) {
            pOther->m_includedNodes = foundNodes;
            m_includedNodes &= ~foundNodes;
        }
        else {
            pOther->m_includedNodes = m_includedNodes;
            pOther->m_includedNodes &= ~foundNodes;
            m_includedNodes = foundNodes;
        }
    }

    return !isConnected;
}


// Initialization code to convert sequence from 1 to n into graph
typedef short numtype;
typedef std::set<numtype> numset;

bool setIntersect(const numset& setA, const numset& setB)
{
    for(int a : setA) {
        if (setB.find(a) != setB.end())
            return true;
    }
    return false;
}

int Sequence2Graph(Graph** ppGraph, int n)
{
    typedef std::map<numset, int> numset2int;
    numset2int factors2count;

    int coPrimesRemoved = n>0; // for {1}
    // Calculate all sets of prime factors, and how many numbers belong to each set
    for(numtype i=2; i<=n; ++i) {
        if ((i > n/2) && (std::find(PRIMES, PRIMES+NPRIMES, i) !=PRIMES+NPRIMES)) {
            ++coPrimesRemoved;
        }
        else {
            numset factors;
            for(numtype j=0; j<NPRIMES && PRIMES[j]<n; ++j) {
                if (i % PRIMES[j] == 0) {
                    factors.insert(PRIMES[j]);
                }
            }
            factors2count[factors]++;
        }
    }

    // Create graph
    Graph*& pGraph = *ppGraph;
    pGraph = new Graph(factors2count.size());
    int srcNodeNum = 0;
    for(numset2int::const_iterator i = factors2count.begin(); i!=factors2count.end(); ++i) {
        pGraph->SetNodeValue(srcNodeNum, i->second);
        numset2int::const_iterator j = i;
        int tgtNodeNum = srcNodeNum+1;
        for(++j; j!=factors2count.end(); ++j) {
            if (setIntersect(i->first, j->first)) {
                pGraph->SetConnection(srcNodeNum, tgtNodeNum);
            }
            ++tgtNodeNum;
        }
        ++srcNodeNum;
    }

    return coPrimesRemoved;
}

coprimes( n) を計算するためのグラフを以下に赤で示します (古いアプローチを黒で示しています)。

観測された (指数関数的な) 増加率に基づくn=3000と、プログラムが爆発しないと仮定すると、予測は 30 時間になります。これは、特により多くの最適化により、計算上実行可能に見え始めていますが、必要な 5 秒にはほど遠いです! 必要な解決策は間違いなく短くて甘いものですが、これは楽しいものでした...

Answer 2

これは、指定されたシーケンスを 5秒未満で取得するアプローチですn=62(最適化を使用すると、5 秒で実行されますが、この問題に対する 2 回目の試みの方がうまくいくn=75ことに注意してください)。問題のモジュロ部分は、関数が大きくなるにつれて数値エラーを回避することと関係があると想定しているため、今のところ無視しています。

このアプローチは、素数ごとにサブセット内の最大 1 つの数値を選択できるという事実に基づいています。

• 最初の素数 2 から始めます。2 を含めない場合、この素数の組み合わせは 1 つになります。2 を含めると、2 で割り切れる数と同じ数の組み合わせがあります。
• 次に、2 番目の素数 3 に進み、それを含めるかどうかを決定します。含めない場合、この素数には 1 つの組み合わせがあります。2 を含めると、3 で割り切れる数と同じ数の組み合わせがあります。
• ... 等々。

例を挙げる{1,2,3,4}と、セット内の数値を次のように素数にマッピングします。この段階で説明を簡単にするため、素数として 1 を含めました。

1 → {1}
2 → {2,4}
3 → {3}

「素数」1 には 2 つの組み合わせ (含まないか 1)、素数 2 には 3 つの組み合わせ (含まないか 2 か 4)、3 には 2 つの組み合わせ (含まないか 1) があります。 3)。したがって、サブセットの数は です2 * 3 * 2 = 12。

同様{1,2,3,4,5}に

1 → {1}
2 → {2,4}
3 → {3}
5 → {5}

与える2 * 3 * 2 * 2= 24。

しかし、 の場合{1,2,3,4,5,6}、物事はそれほど単純ではありません。我々は持っています

1 → {1}
2 → {2,4,6}
3 → {3}
5 → {5}

しかし、素数 2 に数 6 を選択した場合、素数 3 に数を選択することはできません (脚注として、私の最初のアプローチでは、これに戻るかもしれませんが、3 の選択肢が6 を選んだときに半分に減ったので、素数 2 の組み合わせの数として 4 ではなく 3.5 を使用し2 * 3.5 * 2 * 2 = 28、正しい答えを出しました。n=17ただし、このアプローチを を超えて機能させることはできませんでした。)

私がこれに対処した方法は、各レベルで素因数の各セットの処理を分割することです。素因数{2,4}が{2}あり{6}ますが、素因数があり{2,3}ます。1 (素数ではない) の誤ったエントリを省略すると、次のようになります。

2 → {{2}→{2,4}, {2,3}→6}
3 → {{3}→{3}}
5 → {{5}→{5}}

組み合わせの数を計算するための 3 つのパスがあります。素数 2 を選択しないパスが 1 つと、実行するパスが 2 つ ( through{2}→{2,4}と through ) あり{2,3}→6ます。

• 最初のパスには1 * 2 * 2 = 4、3 を選択するかしないか、5 を選択するかどうかを選択できるため、組み合わせがあります。
• 同様に、2 番目には 2 * 2 * 2 = 8、2 または 4 のいずれかを選択できることに注意して組み合わせがあります。
• 素数31 * 1 * 2 = 2の選択肢は 1 つしかなく、それを使用しないためです。

全体として、これにより4 + 8 + 2 = 14組み合わせが得られます (最適化のため、1 番目と 2 番目のパスをまとめて を取得することもできます3 * 2 * 2 = 12)。1 を選択するか選択しないかの選択肢もあるので、組み合わせの総数は です2 * 14 = 28。

パスを再帰的に実行する C++ コードを以下に示します。(Visual Studio 2012 で記述された C++11 ですが、プラットフォーム固有のものは何も含めていないため、他の gcc でも動作するはずです)。

#include <cassert>
#include <vector>
#include <set>
#include <algorithm>
#include <iterator>
#include <iostream>
#include <ctime>

const int PRIMES[] = // http://rlrr.drum-corps.net/misc/primes1.shtml
    { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
        53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
        103, 107, 109, 113, 127, 131, 137, 139, 149,
        151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199 };
const int NPRIMES = sizeof(PRIMES) / sizeof(int);

typedef std::vector<int> intvec;
typedef std::set<int> intset;
typedef std::vector<std::set<int>> intsetvec;

struct FactorSetNumbers
{
    intset factorSet;
    intvec numbers; // we only need to store numbers.size(), but nice to see the vec itself
    FactorSetNumbers() {}
    FactorSetNumbers(const intset& factorSet_, int n)
        : factorSet(factorSet_)
    {
        numbers.push_back(n);
    }
};
typedef std::vector<FactorSetNumbers> factorset2numbers;
typedef std::vector<factorset2numbers> factorset2numbersArray;

double NumCoPrimeSubsets(
    const factorset2numbersArray& factorSet2Numbers4FirstPrime,
    int primeIndex, const intset& excludedPrimes)
{
    const factorset2numbers& factorSet2Numbers = factorSet2Numbers4FirstPrime[primeIndex];
    if (factorSet2Numbers.empty())
        return 1;

    // Firstly, we may choose not to use this prime number at all
    double numCoPrimeSubSets = NumCoPrimeSubsets(factorSet2Numbers4FirstPrime,
        primeIndex + 1, excludedPrimes);
    // Optimization: if we're not excluding anything, then we can collapse
    // the above call and the first call in the loop below together
    factorset2numbers::const_iterator i = factorSet2Numbers.begin();
    if (excludedPrimes.empty()) {
        const FactorSetNumbers& factorSetNumbers = *i;
        assert(factorSetNumbers.factorSet.size() == 1);
        numCoPrimeSubSets *= (1 + factorSetNumbers.numbers.size());
        ++i;
    }
    // We are using this prime number. The number of subsets for this prime number is the sum of
    // the number of subsets for each set of integers whose factors don't include an excluded factor
    for(; i!=factorSet2Numbers.end(); ++i) {
        const FactorSetNumbers& factorSetNumbers = *i;
        intset intersect;
        std::set_intersection(excludedPrimes.begin(),excludedPrimes.end(),
            factorSetNumbers.factorSet.begin(),factorSetNumbers.factorSet.end(),
            std::inserter(intersect,intersect.begin()));
        if (intersect.empty()) {
            intset unionExcludedPrimes;
            std::set_union(excludedPrimes.begin(),excludedPrimes.end(),
                factorSetNumbers.factorSet.begin(),factorSetNumbers.factorSet.end(),
                std::inserter(unionExcludedPrimes,unionExcludedPrimes.begin()));
            // Optimization: don't exclude on current first prime,
            // because can't possibly occur later on
            unionExcludedPrimes.erase(unionExcludedPrimes.begin());
            numCoPrimeSubSets += factorSetNumbers.numbers.size() *
                NumCoPrimeSubsets(factorSet2Numbers4FirstPrime,
                    primeIndex + 1, unionExcludedPrimes);
        }
    }
    return numCoPrimeSubSets;
}

int main(int argc, char* argv[])
{
    const int MAXCALC = 80;

    intsetvec primeFactors(MAXCALC +1);
    // Calculate prime numbers that factor into each number upto MAXCALC
    for(int i=2; i<=MAXCALC; ++i) {
        for(int j=0; j<NPRIMES; ++j) {
            if (i % PRIMES[j] == 0) {
                primeFactors[i].insert(PRIMES[j]);
            }
        }
    }

    const clock_t start = clock();

    factorset2numbersArray factorSet2Numbers4FirstPrime(NPRIMES);
    for(int n=2; n<=MAXCALC; ++n) {
        {
            // For each prime, store all the numbers whose first prime factor is that prime
            // E.g. for the prime 2, for n<=20, we store
            // {2},       { 2,  4,  8, 16 }
            // {2, 3},    { 6, 12, 18 }
            // {2, 5},    { 5, 10, 20 }
            // {2, 7},    { 14 }
            const int firstPrime = *primeFactors[n].begin();
            const int firstPrimeIndex = std::find(PRIMES, PRIMES + NPRIMES, firstPrime) - PRIMES;
            factorset2numbers& factorSet2Numbers = factorSet2Numbers4FirstPrime[firstPrimeIndex];
            const factorset2numbers::iterator findFactorSet = std::find_if(factorSet2Numbers.begin(), factorSet2Numbers.end(),
                [&](const FactorSetNumbers& x) { return x.factorSet == primeFactors[n]; });
            if (findFactorSet == factorSet2Numbers.end()) {
                factorSet2Numbers.push_back(FactorSetNumbers(primeFactors[n], n));
            }
            else {
                findFactorSet->numbers.push_back(n);
            }

            // The number of coprime subsets is the number of coprime subsets for the first prime number,
            // starting with an empty exclusion list
            const double numCoPrimeSubSetsForNEquals1 = 2;
            const double numCoPrimeSubsets = numCoPrimeSubSetsForNEquals1 *
                NumCoPrimeSubsets(factorSet2Numbers4FirstPrime,
                0, // primeIndex
                intset()); // excludedPrimes
            const clock_t now = clock();
            const clock_t ms = now - start;
            std::cout << n << ", " << std::fixed << numCoPrimeSubsets << ", " << ms << "\n";
        }
    }

    return 0;
}

タイミング: 0.1 秒未満で 40 までのシーケンス、0.5 秒で 50 までのシーケンス、2.5 秒で 60 まで、20 秒で 70 まで、157 秒で 80 までのシーケンスを計算します。

これは確かに正しい数値を出力しているように見えますが、予想どおり、コストがかかりすぎます。特に、少なくとも指数関数的な時間 (そしておそらく組み合わせ時間) がかかります。

明らかに、このアプローチは必要に応じて拡張できません。ただし、他の人にアイデアを与える (またはこのアプローチを失敗として除外する) 何かがここにある可能性があります。次の 2 つの可能性があるようです。

1. このアプローチは、次のいくつかの組み合わせにより機能します。
   • 私が見つけていない巧妙な数学的最適化がいくつかあり、計算を完全に排除します。
   • bitsetではなく使用など、いくつかの効率の最適化がありますset。
   • キャッシング。これは、再帰呼び出し構造をツリー構造に変更し、段階的に更新できる可能性があるという点で、最も有望に思えます(親子関係は乗算を示し、兄弟関係は加算を示します)。
     
2. このアプローチを機能させることはできません。
   • これとはほとんど関係のないアプローチがいくつかあります。
   • 私が使用した最初のアプローチが機能する可能性があります。これははるかに「閉じた形式」のソリューションであり、upto では非常に効率的に機能n=17し、at 以上では失敗しn=18、少数で終了しました。私は長い時間をかけてパターンを書き上げ、なぜ突然 が失敗したのかを理解しようとしましたがn=18、それを見ることができませんでした。これに戻るか、誰かが興味を持っている場合は別の回答として含めます。

編集: いくつかのトリックを使用していくつかの最適化を行い、可能な場合は既存の計算をやり直さないようにしました。コードは約 10 倍高速です。良さそうに聞こえますが、それは絶え間ない改善にすぎません。#subsets(n+1)本当に必要なのは、この問題に対する洞察#subsets(n)です。

Answer 3

これは私がそれを行う方法です：

1. mod mまでの数の素因数を求めるn
2. セットのキューを作成しq、それに空のセットを追加し、カウンターを 1 に設定します。
3. キューが空でない間、キューXから要素をポップします
4. kからmax(X)までのすべての数値についてn、集合の因数が数値の因数と交差するかどうかを確認します。そうでない場合は、キューに追加し、X U kカウンターを 1 増やします。そうでない場合は、次の k に進みます。
5. 返却カウンター

2 つの重要な点を指摘する必要があります。

• までの数の因数分解は必要ありませんが、nそれらの素因数だけが必要です。つまり、12 の場合は 2 と 3 だけが必要です。このようにして、2 つの数が互いに素であるかどうかをチェックすることは、2 つの集合の交点が空であるかどうかをチェックすることになります。
• 作成する新しいセットごとの「因数のセット」を追跡できます。このようにして、セット内の他のすべての数に対してすべての新しい数をテストする必要はありませんが、素因数セットを次の 1 つに対して交差させるだけです。セット全体。

Answer 4

これが、先ほど述べた別のアプローチです。
以前使用していたものよりもはるかに高速です。coprime_subsets(117)オンラインインタプリタ(イデオン)を使えば5秒以内に計算できます。

このコードは、空のセットから開始し、すべての数値をすべての可能なサブセットに挿入して、可能なセットを構築します。

primes_to_3000 = set([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999])
# primes up to sqrt(3000), used for factoring numbers
primes = set([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53])

factors = [set() for _ in xrange(3001)]
for p in primes:
    for n in xrange(p, 3001, p):
        factors[n].add(p)


def coprime_subsets(highest):
    count = 1
    used = {frozenset(): 1}
    for n in xrange(1, highest+1):
        if n in primes_to_3000:
            # insert the primes into all sets
            count <<= 1
            if n < 54:
                used.update({k.union({n}): v for k, v in used.iteritems()})
            else:
                for k in used:
                    used[k] *= 2
        else:
            for k in used:
                # only insert into subsets that don't share any prime factors
                if not factors[n].intersection(k):
                    count += used[k]
                    used[k.union(factors[n])] += used[k]
    return count

---

これが私のアイデアとPythonでの実装です。遅いようですが、(オンライン インタープリターを使用して) テストした方法と同じかどうかはわかりません...
「実際の」コンピューターで実行すると違いが生じる可能性がありますが、できます。現時点でそれをテストします。

このアプローチには 2 つの部分があります。

• 素因数のリストを事前に生成する
• 可能なサブセットの数を決定するためのキャッシュされた再帰関数を作成します。
  • 各数値について、その要因をチェックして、サブセットに追加できるかどうかを確認します
  • 追加できる場合は、再帰ケースのカウントを取得し、合計に追加します

その後、モジュロを取るだけだと思います...

これが私のpython実装です（改良版）：

# primes up to 1500
primes = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499

factors = [set() for _ in xrange(3001)]
for p in primes:
    for n in xrange(p, 3001, p):
        factors[n].add(p)


def coprime_subsets(highest, current=1, factors_used=frozenset(), cache={}):
    """
    Determine the number of possible coprime subsets of numbers,
    using numbers starting at index current.
    factor_product is used for determining if a number can be added
    to the current subset.
    """
    if (current, factors_used) in cache:
        return cache[current, factors_used]
    count = 1
    for n in xrange(current, highest+1):
        if factors_used.intersection(factors[n]):
            continue
        count += coprime_subsets(highest, n+1, factors_used.union(factors[n]))
    cache[current, factors_used] = count
    return count

他にもアイデアがあり、時間があれば実装してみたいと思います。別のアプローチがかなり高速になる可能性があると思います。

Answer 5

提案された回答と質問の前文は、尋ねられた質問とは異なる質問に対処しているようです。

質問は次のとおりです。

Output the number of coprime subsets of {1, 2, 3, ..., n} modulo m.

提案された回答は、別のものに対処しようとします。

Output the number of coprime subsets of {1, 2, 3, ..., n}.

これらの質問は同等ではありません。1 つ目は有限環 Z_m を扱い、2 つ目は整数環 Z を扱い、算術演算がまったく異なります。

さらに、質問の前文にある「最大公約数が 1 の場合、2 つの整数は互いに素である」という定義は Z_m には適用されません。除数。

同じ反論が質問の例にも当てはまります。3 と 4 は 7 を法とする互いに素ではありません。なぜなら、どちらも 7 を法とする 2 で割り切れるからです: 4=(2*2)%7 と 3=(2*5)%7.

実際、Z_m 算術は非常に奇妙で、少なくとも素数 m の場合、O(1) 時間で答えを出すことができます: 任意の n と素数 m について、m を法とする互いに素なペアは存在しません。理由は次のとおりです。Z_m の非ゼロ要素は次数 m-1 の巡回群を形成します。これは、Z_m の非ゼロ要素 a および b に対して、Z_m の c に対して a=bc と書けることを意味します。これは、素数 m に対して Z_m に互いに素なペアがないことを証明しています。

質問の例から: {2, 3} mod 7 と {3, 4} mod 7: 2=(3*3)%7 と 3=(4*6)%7 を見てみましょう。したがって、{2,3} は Z_7 で互いに素ではなく (どちらも 3 で割り切れます)、{3,4} は Z_7 で互いに素ではありません (どちらも 4 で割り切れます)。

次に、非素数 m の場合を考えてみましょう。ma を素数ベキ m=p_1^i_1*...*p_k^i_k の積として書きます。a と b に共通の素因数がある場合、それらは明らかに互いに素ではありません。それらの少なくとも 1 つ、たとえば b が素数 p_1、...、p_k のいずれも割り切らない場合、素数 m の場合とほぼ同じ理由で、a と b は共通の因数を持ちます。b は乗法となります。 Z_m の単位であるため、a は Z_m の b で割り切れます。

したがって、m が複合であり、a と b が m の異なる素因数で割り切れる場合を考慮する必要があります。たとえば、a が p で割り切れ、b が q で割り切れるとします。この場合、それらは実際に余素になる可能性があります。たとえば、6 を法とする 2 と 3 は互いに素です。

したがって、余素ペアの問題は、次の手順に要約されます。

1. n より小さい m の素因数を見つけます。存在しない場合、互いに素なペアはありません。

2. n より小さい m の因数のままであるこれらの素因数のベキの積を列挙します。

3. それらの中の Z-comprime ペアの数を計算します。