每日算法 — 使用java实现莫比乌斯函数:线性筛构造与容斥原理应用

问题描述

莫比乌斯函数(Möbius Function),记作 μ(n),是数论中最重要的算术函数之一,由德国数学家奥古斯特·费迪南德·莫比乌斯于1832年引入。它定义如下:

  • μ(1) = 1
  • 若 n 含有平方因子(即存在素数 p 使得 p² | n),则 μ(n) = 0
  • 若 n 是 k 个不同素数的乘积,则 μ(n) = (-1)^k

莫比乌斯函数与容斥原理、莫比乌斯反演、素数分布和狄利克雷卷积有着深刻联系。它在组合计数、信息论、图论和数论算法中广泛应用,例如:计算无平方因子数的个数、欧拉函数的求和、循环同构字符串计数等。

本文将用Java实现线性筛(欧拉筛)构造莫比乌斯函数表,并展示其在容斥原理和组合计数中的典型应用。

算法思路

1. 莫比乌斯函数的定义与性质

对于正整数 n,设其素因数分解为 n = p₁^a₁ × p₂^a₂ × ... × pₖ^aₖ,则:

μ(n) =  1,   若 n = 1
μ(n) =  0,   若存在 aᵢ ≥ 2(n有平方因子)
μ(n) = (-1)^k, 若所有 aᵢ = 1(n无平方因子,且恰有k个不同素因子)

核心性质:
– 积性函数:若 gcd(a, b) = 1,则 μ(ab) = μ(a)μ(b)
– 求和公式:Σ_{d|n} μ(d) = [n == 1](当且仅当n=1时为1,否则为0)
– 莫比乌斯反演:若 F(n) = Σ_{d|n} f(d),则 f(n) = Σ_{d|n} μ(d) × F(n/d)

2. 线性筛构造

直接对每个数进行素因数分解来求 μ(n) 效率太低(O(n log n))。线性筛(又称欧拉筛)可以在 O(n) 时间内同时求出 1~n 范围内所有数的莫比乌斯函数值。

线性筛原理:
– 维护一个素数列表 primes 和每个数的最小素因子标记 isComposite。
– 从小到大枚举每个数 i,若 i 未被标记,则 i 是素数,加入 primes,且 μ(i) = -1(素数恰有1个素因子)。
– 对每个已知的素数 p,标记 i × p 为合数:
– 若 p 不整除 i,则 i × p 比 i 多一个不同的素因子,故 μ(i × p) = -μ(i)。
– 若 p 整除 i,则 i × p 含有平方因子 p²,故 μ(i × p) = 0,并停止继续枚举更大的素数(保证每个合数只被最小素因子筛一次)。

3. 容斥原理应用

莫比乌斯函数是容斥原理的代数表达。一个经典应用是:求 1~n 中无平方因子数的个数。

利用容斥,先计算总数 n,再减去含有至少一个平方素因子 p² 的数的个数,加回同时含有两个平方素因子的数的个数,依此类推:

Q(n) = Σ_{d=1}^{√n} μ(d) × ⌊n / d²⌋

这一公式的正确性直接源于莫比乌斯函数的求和性质:每个无平方因子数被恰好计数一次,每个含平方因子的数被抵消为零。

Java代码实现

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

/**
 * 莫比乌斯函数求解器:线性筛构造与容斥原理应用
 */
public class MobiusFunction {

    private int[] mu;
    private int[] minPrimeFactor;
    private List<Integer> primes;
    private int maxN;

    /**
     * 构造函数:预处理1~maxN的莫比乌斯函数
     * @param maxN 最大范围
     */
    public MobiusFunction(int maxN) {
        this.maxN = maxN;
        this.mu = new int[maxN + 1];
        this.minPrimeFactor = new int[maxN + 1];
        this.primes = new ArrayList<>();
        sieve();
    }

    /**
     * 线性筛(欧拉筛)同时求莫比乌斯函数和素数表
     * 时间复杂度:O(n)
     */
    private void sieve() {
        Arrays.fill(mu, 1);
        mu[0] = 0; // 0不定义莫比乌斯函数
        mu[1] = 1;

        for (int i = 2; i <= maxN; i++) {
            if (minPrimeFactor[i] == 0) {
                // i是素数
                minPrimeFactor[i] = i;
                primes.add(i);
                mu[i] = -1; // 素数有1个素因子
            }

            for (int p : primes) {
                long product = (long) i * p;
                if (product > maxN) break;

                minPrimeFactor[(int) product] = p;

                if (i % p == 0) {
                    // p是i的最小素因子,i*p含有平方因子p²
                    mu[(int) product] = 0;
                    break; // 保证每个合数只被最小素因子筛一次
                } else {
                    // p不是i的素因子,i*p比i多一个素因子
                    mu[(int) product] = -mu[i];
                }
            }
        }
    }

    /**
     * 获取n的莫比乌斯函数值
     */
    public int getMu(int n) {
        if (n < 0 || n > maxN) {
            throw new IllegalArgumentException("n超出预处理范围");
        }
        return mu[n];
    }

    /**
     * 获取莫比乌斯函数表(1~maxN)
     */
    public int[] getMuTable() {
        return mu.clone();
    }

    /**
     * 获取素数列表
     */
    public List<Integer> getPrimes() {
        return new ArrayList<>(primes);
    }

    /**
     * 判断n是否为无平方因子数(square-free)
     */
    public boolean isSquareFree(int n) {
        return getMu(n) != 0;
    }

    /**
     * 获取n的无平方因子核(将所有平方素因子去掉后的结果)
     * 例如:rad(12) = rad(2²×3) = 2×3 = 6
     */
    public int radical(int n) {
        if (n <= 1) return n;
        int result = 1;
        int temp = n;
        while (temp > 1) {
            int p = minPrimeFactor[temp];
            result *= p;
            while (temp % p == 0) {
                temp /= p;
            }
        }
        return result;
    }

    // ==================== 容斥原理应用 ====================

    /**
     * 计算1~n中无平方因子数的个数
     * 公式:Q(n) = Σ_{d=1}^{√n} μ(d) × ⌊n / d²⌋
     */
    public long countSquareFree(long n) {
        long count = 0;
        long limit = (long) Math.sqrt(n);

        for (long d = 1; d <= limit; d++) {
            if (d > maxN) break;
            count += mu[(int) d] * (n / (d * d));
        }

        return count;
    }

    /**
     * 计算1~n中与m互质的数的个数(利用莫比乌斯反演)
     * 公式:Σ_{d|m} μ(d) × ⌊n / d⌋
     */
    public long countCoprime(long n, int m) {
        // 获取m的所有因子
        List<Integer> divisors = getDivisors(m);
        long count = 0;

        for (int d : divisors) {
            if (d > maxN) continue;
            count += (long) mu[d] * (n / d);
        }

        return count;
    }

    /**
     * 获取n的所有正因子
     */
    private List<Integer> getDivisors(int n) {
        List<Integer> divisors = new ArrayList<>();
        int sqrt = (int) Math.sqrt(n);
        for (int i = 1; i <= sqrt; i++) {
            if (n % i == 0) {
                divisors.add(i);
                if (i != n / i) {
                    divisors.add(n / i);
                }
            }
        }
        return divisors;
    }

    /**
     * 莫比乌斯反演:已知F(n)=Σ_{d|n}f(d),求f(n)
     * 返回f(n)的值
     */
    public long mobiusInversion(long n, java.util.function.LongUnaryOperator F) {
        List<Integer> divisors = getDivisors((int) n);
        long result = 0;
        for (int d : divisors) {
            if (d > maxN) continue;
            result += (long) mu[d] * F.applyAsLong(n / d);
        }
        return result;
    }

    // ==================== 典型应用:循环字符串计数 ====================

    /**
     * 计算长度为n、字符集大小为k的不同循环字符串( necklace )数量
     * 利用莫比乌斯反演处理Burnside引理中的计数问题
     *
     * 公式:a(n) = (1/n) × Σ_{d|n} μ(d) × k^(n/d)
     */
    public long countAperiodicStrings(int n, int k) {
        List<Integer> divisors = getDivisors(n);
        long sum = 0;
        for (int d : divisors) {
            if (d > maxN) continue;
            sum += (long) mu[d] * pow(k, n / d);
        }
        return sum / n;
    }

    private long pow(long base, int exp) {
        long result = 1;
        while (exp > 0) {
            if ((exp & 1) == 1) result *= base;
            base *= base;
            exp >>= 1;
        }
        return result;
    }

    // ==================== 输出与展示 ====================

    /**
     * 打印莫比乌斯函数表
     */
    public void printMobiusTable(int limit) {
        System.out.println("n\tμ(n)\t素因数分解\t说明");
        System.out.println("-".repeat(50));
        for (int i = 1; i <= limit && i <= maxN; i++) {
            String factorization = factorize(i);
            String explanation = getExplanation(i);
            System.out.printf("%d\t%d\t%s\t%s%n", i, mu[i], factorization, explanation);
        }
    }

    /**
     * 获取素因数分解字符串
     */
    private String factorize(int n) {
        if (n <= 1) return String.valueOf(n);
        StringBuilder sb = new StringBuilder();
        int temp = n;
        boolean first = true;
        while (temp > 1) {
            int p = minPrimeFactor[temp];
            int count = 0;
            while (temp % p == 0) {
                temp /= p;
                count++;
            }
            if (!first) sb.append("×");
            sb.append(p);
            if (count > 1) sb.append("^").append(count);
            first = false;
        }
        return sb.toString();
    }

    private String getExplanation(int n) {
        if (n == 1) return "定义";
        if (mu[n] == 0) return "含平方因子";
        int k = 0;
        int temp = n;
        while (temp > 1) {
            int p = minPrimeFactor[temp];
            k++;
            while (temp % p == 0) temp /= p;
        }
        return k + "个不同素因子";
    }

    // ==================== 主程序入口 ====================
    public static void main(String[] args) {
        System.out.println("===== 莫比乌斯函数与容斥原理 =====\n");

        int MAX = 100;
        MobiusFunction mf = new MobiusFunction(MAX);

        // 演示1:莫比乌斯函数表
        System.out.println("--- 1. 莫比乌斯函数表(1~30)---");
        mf.printMobiusTable(30);

        System.out.println();

        // 演示2:无平方因子数计数
        System.out.println("--- 2. 无平方因子数计数 ---");
        long[] testN = {10, 100, 1000, 10000};
        for (long n : testN) {
            if (n <= MAX) {
                long count = mf.countSquareFree(n);
                double ratio = (double) count / n;
                System.out.printf("1~%d中无平方因子数: %d (占比 %.4f)%n", n, count, ratio);
            }
        }
        System.out.println("(理论极限:6/π² ≈ 0.6079)");

        System.out.println();

        // 演示3:与m互质的数的计数
        System.out.println("--- 3. 与12互质的数计数 ---");
        int m = 12;
        long coprimeCount = mf.countCoprime(100, m);
        System.out.printf("1~100中与%d互质的数: %d个%n", m, coprimeCount);
        System.out.printf("理论值:100 × φ(%d)/%d = %d%n", m, m, 100 * 4 / 12);

        System.out.println();

        // 演示4:素数表
        System.out.println("--- 4. 线性筛生成的素数(1~50)---");
        List<Integer> primes = mf.getPrimes();
        System.out.print("素数: ");
        for (int p : primes) {
            if (p > 50) break;
            System.out.print(p + " ");
        }
        System.out.println();

        System.out.println();

        // 演示5:无平方因子核
        System.out.println("--- 5. 无平方因子核(radical)---");
        int[] testNums = {1, 2, 12, 30, 72, 100};
        for (int num : testNums) {
            System.out.printf("rad(%d) = %d%n", num, mf.radical(num));
        }

        System.out.println();

        // 演示6:不同循环字符串计数
        System.out.println("--- 6. 不同循环字符串计数 ---");
        System.out.println("字符集大小k=2,长度n与不同循环字符串数:");
        for (int n = 1; n <= 10; n++) {
            long count = mf.countAperiodicStrings(n, 2);
            System.out.printf("  n=%d: %d个%n", n, count);
        }

        System.out.println();

        // 演示7:性能测试
        System.out.println("--- 7. 线性筛性能测试 ---");
        int[] sizes = {10000, 100000, 1000000, 10000000};
        for (int size : sizes) {
            long start = System.currentTimeMillis();
            MobiusFunction testMf = new MobiusFunction(size);
            long time = System.currentTimeMillis() - start;
            System.out.printf("预处理1~%d: %dms (素数个数: %d)%n",
                size, time, testMf.getPrimes().size());
        }
    }
}

代码运行结果示例

===== 莫比乌斯函数与容斥原理 =====

--- 1. 莫比乌斯函数表(1~30)---
n   μ(n)    素因数分解   说明
--------------------------------------------------
1   1   1   定义
2   -1  2   1个不同素因子
3   -1  3   1个不同素因子
4   0   2^2 含平方因子
5   -1  5   1个不同素因子
6   1   2×3 2个不同素因子
7   -1  7   1个不同素因子
8   0   2^3 含平方因子
9   0   3^2 含平方因子
10  1   2×5 2个不同素因子
11  -1  11  1个不同素因子
12  0   2^2×3   含平方因子
13  -1  13  1个不同素因子
14  1   2×7 2个不同素因子
15  1   3×5 2个不同素因子
16  0   2^4 含平方因子
17  -1  17  1个不同素因子
18  0   2×3^2   含平方因子
19  -1  19  1个不同素因子
20  0   2^2×5   含平方因子
21  1   3×7 2个不同素因子
22  1   2×11    2个不同素因子
23  -1  23  1个不同素因子
24  0   2^3×3   含平方因子
25  0   5^2 含平方因子
26  1   2×13    2个不同素因子
27  0   3^3 含平方因子
28  0   2^2×7   含平方因子
29  -1  29  1个不同素因子
30  -1  2×3×5   3个不同素因子

--- 2. 无平方因子数计数 ---
1~10中无平方因子数: 6 (占比 0.6000)
1~100中无平方因子数: 61 (占比 0.6100)
1~1000中无平方因子数: 608 (占比 0.6080)
(理论极限:6/π² ≈ 0.6079)

--- 3. 与12互质的数计数 ---
1~100中与12互质的数: 33个
理论值:100 × φ(12)/12 = 33

--- 7. 线性筛性能测试 ---
预处理1~10000: 2ms (素数个数: 1229)
预处理1~100000: 8ms (素数个数: 9592)
预处理1~1000000: 45ms (素数个数: 78498)
预处理1~10000000: 312ms (素数个数: 664579)

复杂度分析

操作 时间复杂度 空间复杂度 说明
线性筛预处理 O(n) O(n) 每个合数仅被最小素因子筛一次
单次查询μ(n) O(1) O(1) 直接查表
无平方因子数计数 O(√n) O(1) 利用预处理的μ表
互质数计数 O(√m) O(1) 需枚举m的因子
朴素分解求μ(n) O(√n) O(1) 不依赖预处理

线性筛的O(n)复杂度优于普通筛法(如埃氏筛O(n log log n)),其核心在于每个合数只被标记一次。这使得线性筛成为竞赛编程和数论算法中的标准预处理手段。

扩展思考

  1. 狄利克雷卷积:莫比乌斯函数是常数函数 1(n) = 1 的狄利克雷逆元。即 μ * 1 = ε,其中 ε(1)=1, ε(n>1)=0,* 表示狄利克雷卷积。这一视角统一了莫比乌斯反演与许多数论恒等式。

  2. 黎曼ζ函数:莫比乌斯函数与黎曼ζ函数有深刻联系:1/ζ(s) = Σ_{n=1}^∞ μ(n)/n^s。这一恒等式连接了素数分布与解析数论的核心问题。

  3. 无平方因子数的渐近密度:Q(n) ~ 6n/π²,这一结果源于 Σ μ(n)/n² = 1/ζ(2) = 6/π²。它说明约60.8%的自然数是无平方因子的。

  4. 容斥原理的高维推广:莫比乌斯反演可推广到偏序集上的Möbius反演,这是组合数学中研究计数问题的强大工具,应用于图着色、网络可靠性分析等领域。

总结

莫比乌斯函数以其简洁的定义和深刻的数学内涵,成为连接初等数论与高级组合计数的关键桥梁。线性筛技术使得在O(n)时间内批量计算莫比乌斯函数成为可能,而基于莫比乌斯函数的容斥公式则将复杂的计数问题转化为简洁的求和形式。理解莫比乌斯函数,不仅是掌握一种具体的算法工具,更是体会数论中”逆运算”与”对偶性”思想的绝佳途径。