[Shitpost] How to Teach Kids Basic OI Number Theory Using the Methods of a Mathematics Textbook

minstdfx · 2026-10-28 01:18:33 · 算法·理论

我看隔壁数学系搞的不错，术语极大丰富，各种人话基本消灭，逻辑严谨、体系顺序也受重视，如果加上讲的是 oi 内容，数学课本就是我们理想中的 oi 讲稿。

:::info[英文版/English]

Local Structure of the Integers and Finite Summation of Arithmetic Functions

Conventions

This text rests on ordinary set theory and the constructions of the integers and rational numbers. Finite sets, maps, equivalence relations, finite sequences, induction, and the theorem of finite recursion are taken as established. Real logarithms, integrals, and elementary asymptotic estimates are used only in complexity proofs. Write

When l>r, set [l,r]=\varnothing. An empty sum is 0, and an empty product is 1. The notation [P] denotes the indicator of the proposition P: it is 1 when the proposition holds and 0 otherwise. Unless otherwise specified, divisors are positive and primes are positive primes; moduli are always positive integers, and remainders are always the least nonnegative representatives. The notation p^e\parallel n means p^e\mid n and p^{e+1}\nmid n. Integer square roots and other root expressions used as bounds are interpreted with the explicitly indicated rounding. An unsubscripted \log is the natural logarithm; changing to any fixed base in an asymptotic operation count affects only the constant.

All rings are commutative and have an identity; the particular uses of groups, units, and residue classes are defined in the text. The degenerate case of modulus 1 is allowed. A general coefficient ring is denoted by R, and integers are mapped into R by the canonical unital homomorphism. Division in a ring uses the inverses of units. Ordinary multiplication, pointwise multiplication, and Dirichlet convolution are distinguished throughout; the last is consistently denoted by *. Equalities of formal series mean coefficientwise equality unless convergence conditions are supplied separately.

Computational conventions. Time bounds count elementary arithmetic operations, and storage bounds count integers or ring elements. If integers, addresses, moduli, and intermediate results can be represented exactly in a constant number of machine words, the corresponding additions, subtractions, multiplications, divisions, comparisons, rounding operations, and random accesses are counted as elementary operations; the bit complexity of arbitrary-precision arithmetic is accounted for separately. The costs of operations on ring elements and their storage are specified in each construction. Algorithmic complexity includes binary exponentiation, inversion, modulus factorization, local function evaluation, and preprocessing.

Hash-table lookups are charged at their expected cost; error probabilities and orbit models for randomized algorithms are described in the corresponding chapters.

Within each chapter, definitions, lemmas, propositions, theorems, corollaries, constructions, and examples share one sequence numbered “chapter.item”; equations are numbered separately. Chapters follow the order of mathematical dependence: divisibility and residue-class structures come first, followed by convolution of arithmetic functions, and finally finite representations, sieves, and lattice-point sums. Algorithm names are retained at the corresponding constructions for comparison with the original notes; a page correspondence, editorial corrections, and supplementary sources appear at the end.

Chapter I — Divisibility and Ideals of the Integers

§1 Divisibility and Division with Remainder

Definition 1.1. For a,b\in\mathbf Z, one says that a divides b, and writes a\mid b, if there exists q\in\mathbf Z such that b=aq. For n>0, put

By definition, every integer divides 0, whereas 0\mid b if and only if b=0. For n>0, the pairing of factors d\mapsto n/d is an involution on D(n).

Proposition 1.2. For integers a,b,c, divisibility is transitive; if a\mid b and a\mid c, then a\mid ub+vc for all u,v\in\mathbf Z; if a\mid b and b\mid a, then |a|=|b|. For t\ne0,

Proof. The first two assertions follow by substitution in the definition of divisibility. If a,b are both nonzero, the equalities b=aq and a=br give qr=1, so q,r\in\{1,-1\}; if either is zero, both are zero. The reverse implication in (1.2) uses cancellation for multiplication in the integers. \square

Theorem 1.3 (Division with remainder). Let m>0 and a,d\in\mathbf Z. There is a unique pair of integers (q,r) satisfying

Proof. Taking q=\lfloor(a-d)/m\rfloor and r=a-qm proves existence. If (q',r') is another such pair, then m(q-q')=r'-r and |r'-r|<m; hence both sides are zero. \square

When d=0, the resulting r is denoted by a\bmod m. Any m consecutive integers have distinct remainders, which therefore run through 0,\ldots,m-1: the difference of two of these integers has absolute value less than m, so they can be congruent only if they are equal.

§2 Ideals of the Integers and Greatest Common Divisors

Definition 1.4. An ideal of the ring of integers is a nonempty subset I\subseteq\mathbf Z with the following properties: if x,y\in I, then x-y\in I; if x\in I and t\in\mathbf Z, then tx\in I. For a finite family of integers a_1,\ldots,a_r, put

This is the smallest ideal containing all the a_i. The ideal generated by the empty family is, by convention, \{0\}.

Theorem 1.5. Every ideal of the integers has a unique expression d\mathbf Z, where d\in\mathbf N.

Proof. If I=\{0\}, take d=0. Otherwise, I contains a positive integer; let d be the least such integer. For any a\in I, division with remainder gives a=qd+r with 0\le r<d. Since r\in I, minimality forces r=0. Thus I\subseteq d\mathbf Z; the reverse inclusion follows from the ideal property. If d\mathbf Z=e\mathbf Z, the generators divide one another; since d,e\ge0, one has d=e. \square

Definition 1.6. The greatest common divisor of a finite family of integers is the nonnegative generator of the ideal they generate; it is denoted by \gcd(a_1,\ldots,a_r). Two integers are relatively prime if their greatest common divisor is 1. In particular, \gcd(0,0)=0.

Theorem 1.7 (Bézout). Let d=\gcd(a_1,\ldots,a_r). Then d\mid a_i for every i, and there exist integers u_i such that

Every common divisor divides d. Consequently, unless the family consists entirely of zeros, d is the greatest positive common divisor in the usual order; the equation \sum_i a_i x_i=c has an integer solution if and only if d\mid c.

Proof. The equality (a_1,\ldots,a_r)=d\mathbf Z gives the first two assertions. A common divisor divides every integer linear combination, and hence divides d. The last assertion is equivalent to c\in d\mathbf Z, including the case d=0. \square

Proposition 1.8. For integers a,b,q,

If \gcd(a,b)=1 and b\mid ac, then b\mid c. If also a\mid c, then ab\mid c.

Proof. The first equality follows from (a,b)=(b,a-qb), and the second from multiplying the generated ideal by t. Choose ua+vb=1 and multiply by c to obtain b\mid c. If c=at, then b\mid at implies b\mid t, whence ab\mid c. \square

Construction 1.9 (Extended Euclidean algorithm). Apply division with remainder successively to |a|,|b| until the remainder is zero. The last nonzero remainder is \gcd(a,b); if both inputs are initially zero, return 0 directly. If the recursive call has obtained

return x=y' and y=x'-qy', so that ax+by=d. Finally, adjust the coefficients according to the signs of the original inputs.

Proof. Equation (1.6) preserves the greatest common divisor, and the strictly decreasing nonnegative remainders ensure termination. Back-substitution gives the stated coefficients. Let u>v be consecutive positive remainders. If v\le u/2, a halving has already occurred; otherwise, the next remainder is u-v<u/2. Thus the size is at least halved every two steps, and the number of divisions is O(1+\log(\min(|a|,|b|)+1)). \square

Remark. This bound counts integer divisions; it does not include coefficient sizes or the bit costs of integer multiplication and division. Successively taking greatest common divisors of several integers gives a result independent of the grouping, since the generated ideal is the same.

§3 Primes and Unique Factorization

Definition 1.10. A prime is a positive integer with exactly two positive divisors. An integer greater than 1 that is not prime is called composite. For a prime p and a positive integer n, put

Lemma 1.11. If a prime p divides ab, then p\mid a or p\mid b.

Proof. If p\nmid a, then \gcd(p,a) is a positive divisor of p different from p, and hence equals 1. Apply Proposition 1.8. \square

Theorem 1.12 (Unique factorization). Every n>0 has a unique expression

in which only finitely many exponents are nonzero; n=1 corresponds to the empty product.

Proof. Proceed by induction on n. If n>1 is not prime, write n=uv with 1<u,v<n and factor u,v separately to obtain existence. Given two factorizations, any prime factor on the left divides some prime factor on the right by Lemma 1.11, so these factors are equal. Cancel this factor and apply induction to obtain uniqueness. \square

Definition 1.13. The least common multiple of a nonempty finite family of positive integers is its least positive common multiple; it is denoted by \operatorname{lcm}(a_1,\ldots,a_r).

Proposition 1.14. For positive integers a,b,

The corresponding minimum and maximum formulas hold for every nonempty finite family of positive integers.

Proof. The divisibility criterion follows directly from unique factorization. Each exponent of a common divisor is at most the corresponding minimum, and these minima do give a common divisor; the argument for common multiples is analogous. Finally, use \min(s,t)+\max(s,t)=s+t. \square

Corollary 1.15. If positive integers a_1,\ldots,a_r are pairwise relatively prime and their product is the kth power of a positive integer, where k>0, then each a_i is the kth power of a positive integer.

Proof. Any prime occurs in at most one a_i. Its exponent in the full product is divisible by k, so its exponent in that a_i is also divisible by k. \square

§4 Congruences and Linear Equations

Definition 1.16. For m>0, integers a,b are congruent modulo m if m\mid a-b; one writes a\equiv b\pmod m. The equivalence class of an integer a is denoted by [a]_m. The set of these classes is denoted by \mathbf Z/m\mathbf Z, with addition and multiplication defined by

Proposition 1.17. Congruence is an equivalence relation. The operations in (1.11) are independent of the representatives and make \mathbf Z/m\mathbf Z a commutative ring with identity. If a positive integer d satisfies d\mid a,b,m, then

For any integer c, let g=\gcd(c,m). Then

Proof. The equivalence relation and the ring operations follow from the compatibility of divisibility with addition and multiplication. Equation (1.12) is cancellation in the integer equality a-b=tm. Write c=gc' and m=gm', so that \gcd(c',m')=1; Proposition 1.8 gives m'\mid c'(a-b) if and only if m'\mid a-b. \square

Remark. Cancelling a nonunit factor in a congruence generally requires changing the modulus as well. The ring modulo 1 has one element, [0]_1=[1]_1; all integers are congruent modulo 1.

Theorem 1.18 (Linear congruences). Let m>0 and d=\gcd(a,m). The equation

has a solution if and only if d\mid b. If ua+vm=d and d\mid b, then all integer solutions are

Consequently, there are exactly d solution classes modulo m.

Proof. The original equation is equivalent to the existence of an integer solution of ax+my=b. Theorem 1.7 therefore gives the solvability criterion and a particular solution. The difference of two solutions satisfies a(x-x_0)\equiv0\pmod m; apply (1.13). A complete residue system of length m contains exactly d representatives spaced by m/d. \square

§5 Combining Congruence Conditions

Theorem 1.19 (Chinese remainder theorem). Let m_1,\ldots,m_r be pairwise relatively prime positive integers, and let M=\prod_i m_i. The map

is a ring isomorphism. In particular, for any prescribed a_i, the system x\equiv a_i\pmod{m_i} has exactly one solution class modulo M.

Proof. Put M_i=M/m_i and use Bézout's theorem to choose u_i such that u_iM_i\equiv1\pmod{m_i}. Then

satisfies every equation: in the ith component, the ith term gives a_i and all other terms vanish. If two integers satisfy the same conditions, each m_i divides their difference; pairwise relative primality implies that M divides their difference. The map plainly preserves the ring operations and is therefore an isomorphism. A component with modulus 1 imposes no condition and may be omitted. \square

Theorem 1.20 (Combining two arbitrary moduli). Let M,N>0 and d=\gcd(M,N). The system

has a solution if and only if d\mid B-A. When it is solvable, choose uM+vN=d and put

so that all solutions are given by x\equiv A+Mt_0\pmod L.

Proof. Write x=A+Mt. The problem becomes Mt\equiv B-A\pmod N. Theorem 1.18 gives the solvability criterion and the unique class of t modulo N/d. Changing t by one period changes x by L, and L=\operatorname{lcm}(M,N). \square

Construction 1.21 (ExCRT). Start with (A,M)=(0,1), and successively combine the current class x\equiv A\pmod M with the next congruence condition by Theorem 1.20. After each step, normalize A to [0,M-1]. If a step is incompatible, the system has no solution; otherwise, the final class is the set of all integer solutions.

Proof. The induction invariant is that the current residue class is exactly the solution set of the equations already processed. Theorem 1.20 preserves this invariant. The empty system has all integers as its solution set, in agreement with the initial state. \square

Remark. The combined modulus is the least common multiple of the moduli already processed. Divide by the greatest common divisor before multiplying, using an integer representation that accommodates the resulting product. If all intermediate moduli are at most L, then r combinations use O(r\log(L+1)) integer divisions; the bit costs of arbitrary-precision integers are counted separately.

§6 Two Existence Problems on Residue Classes

Example 1.22 (Xiaokai's Doubt). Let a,b>1 and \gcd(a,b)=1. The greatest nonnegative integer that cannot be written as ax+by with x,y\in\mathbf N is

Proof. The remainders of 0,b,\ldots,(a-1)b modulo a are distinct, so every class contains exactly one of these terms jb. Within that class, the representable nonnegative integers are exactly jb+ta with t\in\mathbf N: if n=ax+by, write y=j+sa with 0\le j<a; then n=jb+a(x+sb)\ge jb. The converse is immediate. Thus the integers in that class below jb are not representable, and the greatest candidate is jb-a. Taking the maximum over all classes gives (a-1)b-a. Since a,b>1 are relatively prime, this number is nonnegative and hence belongs to the required set. \square

Remark. If a=1 or b=1, every nonnegative integer is representable, so there is no greatest nonrepresentable nonnegative integer. Equation (1.20) then gives -1, which can only be used as a separately prescribed sentinel value. If a,b are not relatively prime, infinitely many nonnegative integers are not representable.

Example 1.23 (CF338D, GCD Table). Given positive integers a_1,\ldots,a_k and upper bounds n,Y, determine whether there exist

Put L=\operatorname{lcm}(a_1,\ldots,a_k). First check L\le n, and then use ExCRT to solve

If this system has no solution, neither does the original problem. Otherwise, all its solutions are y\equiv y_0\pmod L. Take the least positive representative y_*=y_0 if y_0>0, or y_*=L if the normalized representative is y_0=0. The original problem is solvable if and only if

Proof. Every original solution satisfies L\mid x and (1.22). Since a_i divides both L and y+i-1, and L\mid x,

Thus, whenever an original solution exists, x may be replaced by L. Moreover, adding or subtracting L from y preserves all these greatest common divisors. Hence the least positive representative suffices to test both the bound and the equalities. The converse is verified by (x,y)=(L,y_*). \square

Remark. For a contiguous segment of length k within an n\times m table, the boundary condition is y+k-1\le m, so Y=m-k+1; if Y<1, there is no solution. The condition y\le m in the abbreviated slide statement corresponds to Y=m. Computation of L may stop as soon as L>n is established. Test the resulting candidate by the greatest-common-divisor equalities in (1.23).

Chapter II — Finite Operations in Residue Class Rings

§1 Units, Inverses, and Powers

Definition 2.1. A unit of a commutative ring with identity R is an element having a multiplicative inverse. The set of units is denoted by R^\times; it is an abelian group under ring multiplication. For m>0, put

For an element a of any monoid, set a^0=1, including the case of the zero ring and the case a=0 in a ring.

Proposition 2.2. The class [a]_m is a unit if and only if \gcd(a,m)=1. Its inverse, when it exists, is unique; if the extended Euclidean algorithm gives ua+vm=1, the inverse is [u]_m.

Proof. Invertibility is equivalent to solvability of ax\equiv1\pmod m; apply Theorem 1.18. If ax=ay=1, then x=x(ay)=(xa)y=y. \square

Remark. The group U_1 has one element, and \varphi(1)=1. For a prime p, every nonzero residue class is a unit, so \mathbf F_p=\mathbf Z/p\mathbf Z is a field. For a general composite modulus, a nonzero element need not be invertible.

Construction 2.3 (Binary exponentiation). For a monoid element a and e\in\mathbf N, maintain (z,b,t), initially (1,a,e). While t>0, put z\leftarrow zb if t is odd, and then put b\leftarrow b^2 and t\leftarrow\lfloor t/2\rfloor. Return z.

Proof. The invariant is zb^t=a^e. Write t=2q+\epsilon, where \epsilon\in\{0,1\}. After one update, the product is zb^\epsilon(b^2)^q, equal to its previous value. At termination, t=0, so z=a^e. The number of multiplications is O(1+\log(e+1)). \square

Remark. The algorithm requires only associativity and therefore also applies to matrices. Modular multiplication uses exact integer products followed by reduction, or an equivalent implementation without overflow. Exponent 0 gives the identity element, whose normalized representative modulo 1 is still 0.

Construction 2.4 (Batch inversion). Let a_1,\ldots,a_r be units of R. Put P_0=1 and P_i=P_{i-1}a_i. Compute P_r^{-1} once, and then, for i=r down to 1, compute

Proof. Multiplying the two formulas by a_i and P_{i-1}, respectively, gives the identity. The procedure uses one inversion and O(r) multiplications. \square

Proposition 2.5. For a prime p and 1<i<p, put q=\lfloor p/i\rfloor and r=p\bmod i. Then

Thus the inverses of 1,\ldots,p-1 may be computed recursively from 1^{-1}=1.

Proof. Since p=qi+r, one has qi\equiv-r\pmod p. Here 1\le r<i<p because p is prime. Multiplication by i^{-1}r^{-1} gives the result. \square

Remark. This recurrence uses invertibility of nonzero residues modulo a prime. For a composite modulus, even if i is a unit, the remainder m\bmod i may not be a unit; for example, when m=8,i=3, the remainder is 2.

§2 Unit Groups and Reduction of Exponents

Theorem 2.6 (Euler). If \gcd(a,m)=1, then

In particular, if p is prime and p\nmid a, then a^{p-1}\equiv1\pmod p and a^{-1}\equiv a^{p-2}\pmod p.

Proof. Multiplication by [a]_m permutes the finite set U_m. Taking the product of all its elements gives

The common product is a unit and may be cancelled. For a prime modulus, \varphi(p)=p-1; then multiply by a^{-1}. \square

Proposition 2.7. If m=\prod_{p\mid m}p^{\alpha_p}, then

The empty product gives \varphi(1)=1.

Proof. The units modulo p^\alpha are precisely the residue classes not divisible by p, of which there are p^\alpha-p^{\alpha-1}. The Chinese remainder theorem identifies units modulo m with families of units in the prime-power components, so their cardinalities multiply. \square

Theorem 2.8 (Reduction of exponents for nonunits). Let m>0, a\in\mathbf Z, and e\in\mathbf N. If \gcd(a,m)=1, the exponent may be replaced by e\bmod\varphi(m). For arbitrary a, if e\ge\varphi(m), then

If e<\varphi(m), the original exponent is retained.

Proof. The assertion is automatic for m=1. Consider each full prime-power component p^\alpha\mid m. If p\nmid a, then \varphi(p^\alpha)\mid\varphi(m), and Euler's theorem permits changes of the exponent by multiples of \varphi(m). If p\mid a, then a^t is zero modulo p^\alpha for t\ge\alpha. Moreover,

where the last inequality follows from the inductive inequality 2^{\alpha-1}\ge\alpha. Hence both sides of (2.7) vanish in this component. Combine the componentwise conclusions by the Chinese remainder theorem. The unit case follows directly from Theorem 2.6. \square

Remark. In the noncoprime case, choose the exponent according to the threshold in Theorem 2.8. For example, 2^2\equiv0\pmod4, whereas 2^{2\bmod\varphi(4)}=2^0\equiv1\pmod4. In a recursive computation of a power tower, retain both the exponent modulo \varphi(m) and whether it has reached \varphi(m); treat a zero exponent according to a^0=1.

§3 Local Computation of Binomial Coefficients

Definition 2.9. For n\in\mathbf N and k\in\mathbf Z, let \binom nk be the number of k-element subsets of an n-element set, with value 0 when k<0 or k>n. Thus, in the polynomial ring over the integers,

For 0\le k\le n, ordering the elements of the subset and of its complement separately gives the integer identity n!=\binom nk k!(n-k)!.

Lemma 2.10. In a commutative ring of prime characteristic p,

Proof. For 0<j<p, the integer identity j!\binom pj=p(p-1)\cdots(p-j+1) implies p\mid j!\binom pj. Since p\nmid j!, Lemma 1.11 gives p\mid\binom pj. Thus all intermediate terms in the binomial expansion vanish. \square

Theorem 2.11 (Lucas). Let p be prime, and let n=\sum_i n_ip^i and k=\sum_i k_ip^i be the base-p expansions, with 0\le n_i,k_i<p. Then

Proof. Repeated application of Lemma 2.10 in \mathbf F_p[X] gives

The exponent selected from each factor on the right has the form j_ip^i, where 0\le j_i\le n_i<p. These exponents add without carrying. Hence the only way to obtain X^k is to take j_i=k_i in every position; if k_i>n_i in any position, the coefficient is zero. Comparison of coefficients proves the assertion. \square

Remark. The factorials and inverse factorials of 0,\ldots,p-1 can be precomputed in O(p) time and space, after which each evaluation uses O(1+\log_p(n+1)) digit operations. This method applies when a table of length p can be stored. For n\ge p, the factorial n! is zero modulo p and has no multiplicative inverse.

Theorem 2.12 (Legendre). For a prime p and n\in\mathbf N, let s_p(n) denote the sum of the base-p digits of n. Then

Here 0!=1, so both sides are zero when n=0.

Proof. Each 1\le t\le n contributes \nu_p(t) to the total exponent, namely one for every j\ge1 satisfying p^j\mid t. Interchanging these two finite summations gives the first equality. If n=\sum_i n_ip^i, then

which is the second equality. \square

Proposition 2.13. Let t_1+\cdots+t_r=n, where t_i\in\mathbf N. The multinomial coefficient

is odd if and only if the binary addition of the t_i involves no carries: their sets of binary 1 positions are pairwise disjoint and have union equal to the set of 1 positions of n.

Proof. By Theorem 2.12, its 2-adic valuation is \sum_i s_2(t_i)-s_2(n). Combining two copies of 2^j into one copy of 2^{j+1} decreases the total digit sum by exactly 1. After all carries are performed, the digit sum is s_2(n). Thus the difference is zero if and only if no carry is needed. \square

§4 Prime-Power Moduli and Reduced Factorials

Definition 2.14. Fix a prime p and define positive integers

Both are 1 when n=0, and both are relatively prime to p.

Proposition 2.15. The following integer identity holds:

Proof. Partition the factors of n! into those not divisible by p and the factors p,2p,\ldots,\lfloor n/p\rfloor p. Removing one factor p from each factor in the latter part leaves \lfloor n/p\rfloor!. Removing all its remaining factors p gives exactly U_p(\lfloor n/p\rfloor). \square

Construction 2.16. Let q=p^\alpha, where \alpha\ge1. Precompute

If n=vq+r with 0\le r<q, then

Substituting this formula into (2.17) gives a recursive computation of U_p(n)\bmod q.

Proof. In every complete interval of length q, the factors not divisible by p correspond modulo q to those in the first block. The final partial block corresponds to 1,\ldots,r. The recursive argument is divided by p at each step and therefore eventually reaches 0. \square

Theorem 2.17 (Binomial coefficients modulo a prime power). Let 0\le k\le n and q=p^\alpha. Put

If e\ge\alpha, then \binom nk\equiv0\pmod q; otherwise,

Proof. Since the binomial coefficient is an integer, e\ge0. Separate the powers of p in the three factorials and cancel the common factors p in the integers to obtain

The two U_p factors on the left are relatively prime to q and are therefore invertible modulo q, giving (2.21). If e\ge\alpha, the right-hand side is zero. \square

Construction 2.18 (Binomial coefficients for an arbitrary modulus). For M>1, first obtain the factorization M=\prod_i p_i^{\alpha_i}, compute each prime-power component by Theorem 2.17, and combine them by the Chinese remainder theorem. For M=1, return the normalized remainder 0 directly. This construction is commonly called ex-Lucas.

Remark. The cost of prime factorization is counted separately. For a component q, the table (2.18) requires O(q) time and space. A direct implementation using binary exponentiation at each recursive level computes one corresponding reduced factorial in O((1+\log_p(n+1))(1+\log(n+1))) modular multiplications. If a table of length q exceeds available storage, another method for computing G_p may be used. To evaluate the binomial coefficient, first separate the prime powers in the integers, then invert the units in the denominator.

§5 Parity Counting and Binary Recurrences

Example 2.19 (gym102586E). Given distinct nonnegative integers A_1,\ldots,A_k, count modulo 2 the ordered sequences of length n whose entries are chosen from these k values and whose sum is S. In \mathbf F_2[X], put

The required answer is [X^S]P(X)^n. If n=\sum_{i=0}^{\ell-1}n_i2^i, then

Proof. In the expansion of P(X)^n, each ordered choice contributes exactly one monomial, whose exponent is the sum. In characteristic 2, Lemma 2.10 gives P(X)^{2^i}=P(X^{2^i}), yielding (2.24). Proposition 2.13 gives another interpretation: each binary 1 position of n is assigned to exactly one occurrence count; multinomial coefficients involving carries are all even. \square

Construction 2.20 (Carry states). Let A=\max_j A_j, let n>0, and let \ell be the binary length of n. Denote the ith binary digit of S by s_i. Define f_i(c)\in\mathbf F_2 to be the parity of the number of choices satisfying

among the lowest i positions, where j_t is chosen only when n_t=1. Initially f_0(0)=1 and all other entries are zero. If n_i=0, take only a=0 at each state; if n_i=1, take a=A_1,\ldots,A_k successively. Only when c+a\equiv s_i\pmod2 perform the update

The final answer is f_\ell(\lfloor S/2^\ell\rfloor); entries outside the array range are treated as zero.

Proof. Add n_i2^ia to the left-hand side of (2.25) and use S\bmod2^{i+1}=(S\bmod2^i)+s_i2^i to obtain the update. Once the parity condition is satisfied, the new carry is nonnegative. If 0\le c\le A and 0\le a\le A, then 0\le(c+a-s_i)/2\le A, so all reachable states lie in [0,A]. After all digits of n have been processed, the full sum equals S if and only if the remaining carry is \lfloor S/2^\ell\rfloor. \square

Remark. For n=0, there is exactly one empty sequence, so the answer is [S=0]. A direct implementation uses O((k+1)(A+1)\ell) Boolean operations and O(A+1) space. The states may be stored as a bitset: when n_i=1, first take the exclusive or of the shifts by all A_j, then retain and compress the odd or even positions according to s_i. The intermediate indices extend at most to 2A. If a machine word contains w bits, the shifts and exclusive-or operations use O(k\lceil(A+1)/w\rceil\ell) word operations; ordinary compression requires an additional O((A+1)\ell) time.

Chapter III — Cyclic Subgroups and Integer Factorization

§1 Orders and Primitive Roots

Definition 3.1. Let G be a finite group and g\in G. The order of g is the least positive integer d such that g^d=1, denoted by \operatorname{ord}(g). The subgroup consisting of all integer powers of g is denoted by \langle g\rangle. If G=\langle g\rangle, then G is called cyclic and g a generator. A primitive root modulo m is a generator of U_m.

Proposition 3.2. Every element of a finite group has an order. If its order is d, then

Moreover, d\mid|G|.

Proof. Finiteness implies that two terms of 1,g,g^2,\ldots are equal; cancellation gives a positive power equal to 1. Write any integer t as qd+r with 0\le r<d. Then g^t=1 if and only if r=0, which gives (3.1). The left cosets of \langle g\rangle are either equal or disjoint, and left multiplication gives a bijection between each coset and the subgroup. Thus these sets of size d partition G, and d\mid|G|. \square

Lemma 3.3. A nonzero polynomial of degree r over a field has at most r roots. If elements u,v of an abelian group have relatively prime orders a,b, respectively, then uv has order ab.

Proof. If c is a root of a polynomial, division by X-c has zero remainder. Induction on the degree gives the root bound. For the second assertion, (uv)^{ab}=1. If (uv)^t=1, then u^t=v^{-t} belongs to both \langle u\rangle and \langle v\rangle; its order divides both a,b, and hence is 1. Therefore a\mid t and b\mid t, so ab\mid t. \square

Theorem 3.4. Every finite subgroup of the multiplicative group of a field is cyclic. In particular, a primitive root modulo p exists for every prime p.

Proof. Let G be the finite subgroup, and let e be the least common multiple of all element orders. For each q^\alpha\Vert e, choose an element whose order contains the full factor q^\alpha; raising it to a suitable integer power gives an element of order exactly q^\alpha. By Lemma 3.3, multiplying these elements of pairwise relatively prime orders gives an element of order e, so e\le|G|. On the other hand, every x\in G satisfies x^e=1, and the root bound over a field gives |G|\le e. Thus e=|G|, and the constructed element generates G. \square

Remark. In the preceding proof, if g has order d, then g^t has order d/\gcd(d,t); this follows directly by finding the least positive solution of d\mid ts. The unit group for a composite modulus need not be cyclic. For example, all four units modulo 8 have square 1, so there is no element of order 4.

Proposition 3.5. Let g\in G. Suppose g^D=1 and the prime divisors of D are known. Then \operatorname{ord}(g)=D if and only if

Proof. The order d divides D. If d<D, then D/d>1; choose a prime q\mid D/d. It follows that d\mid D/q, so g^{D/q}=1. The converse follows immediately from d=D. \square

Construction 3.6. Given a unit g\in U_m and the prime factorization of \varphi(m), start with D=\varphi(m). For each prime q\mid D, as long as g^{D/q}=1, put D\leftarrow D/q. The final D is the order of g.

Proof. Each change preserves the condition that the true order divides D. At termination, every test in (3.2) holds, so D is the true order. \square

§2 Discrete Logarithms and Block Search

Definition 3.7. Given g,h in a finite group G, the discrete logarithm problem asks for x\in\mathbf N such that g^x=h. It is solvable if and only if h\in\langle g\rangle; when solvable, all integer solutions form one class modulo \operatorname{ord}(g).

Construction 3.8 (BSGS). Suppose \operatorname{ord}(g)\le N is known, with N\ge1, and put B=\lceil\sqrt N\rceil. Store

retaining only the smallest j when a key occurs more than once. For i=0,1,\ldots,\lfloor(N-1)/B\rfloor in this order, look up

If a value j is found and iB+j<N, return x=iB+j. If the search ends without finding one, there is no solution.

Proof. Every 0\le x<N has a unique expression iB+j with 0\le j<B. The equality g^x=h is equivalent to g^j=h(g^{-B})^i. Whenever the problem is solvable, it has a solution 0\le x<\operatorname{ord}(g)\le N, so the search is exhaustive. Blocks are examined in increasing order of i, and the smallest j is returned within each block. Thus the first valid candidate is the least nonnegative solution. \square

Remark. With a hash table supporting expected constant-time lookup, the number of group operations and the expected lookup cost are both O(\sqrt N); space is O(\sqrt N), with one additional inversion. If a comparison tree or binary search on a sorted table is used, the corresponding logarithmic factors must be included. For units modulo m, one may take N=m without first factoring m; a smaller known upper bound for the order permits a smaller table. If the target is not a unit and the base is a unit, the equation has no solution.

Construction 3.9 (A table for bounded powers of a fixed base). Suppose powers of a fixed g are required repeatedly, with all exponents satisfying 0\le e<K. Put B=\lceil\sqrt K\rceil, and precompute g^j and g^{iB} for 0\le j<B and 0\le i\le\lfloor(K-1)/B\rfloor, respectively. For a query, write e=iB+j and return the product of the two table entries.

Proof. Associativity gives g^{iB+j}=g^{iB}g^j. The two tables are constructed by successive multiplication by g and g^B, respectively. Preprocessing and space are O(\sqrt K), and each query requires one multiplication. \square

§3 Reduction for a Nonunit Base

Construction 3.10 (exBSGS). Given m>0 and integers a,b, find the least nonnegative solution of a^x\equiv b\pmod m. Maintain an integer t, a positive integer M, and normalized representatives B,C modulo M, initially

At each iteration, perform the following steps in order.

1. If B=C, return t.
2. Put g=\gcd(a,M). If g=1, solve a^y\equiv BC^{-1}\pmod M by BSGS; return t+y if a solution exists, and otherwise report that there is no solution.
3. If g>1 and g\nmid B, report that there is no solution.
4. Otherwise, update simultaneously

and continue. The right-hand sides of (3.6) use the variables before the update; finally normalize B' modulo M'.

Proof. Maintain three induction invariants: all x<t have been excluded; C is a unit modulo M; and, for every x\ge t, the original equation is equivalent to

They plainly hold initially. The first step tests exactly x=t. If it fails, only x-t\ge1 remains to be considered. If g\nmid B, the left-hand side is divisible by g and the right-hand side is not, so there is no solution. If g\mid B, cancel one factor g by (1.12) to obtain

which is precisely the updated condition. Since \gcd(a/g,M/g)=1 and the old C remains a unit modulo M/g, the new C' is also a unit. If g=1, cancel C and apply BSGS; its least solution plus t is the least solution of the original equation. Each reduction at least halves M, so the procedure terminates. \square

Remark. When m=1, the initial values already satisfy B=C=0, giving x=0. When a=0, first test exponent 0 according to the convention 0^0=1; thereafter, x=1 is possible only when the target remainder is 0. At each iteration, multiply the accumulated coefficient by the current a/g. There are at most O(\log(m+1)) reductions. With the direct Euclidean cost, this stage uses O(\log^2(m+1)) integer divisions, followed by O(\sqrt m) modular multiplications and expected hash-table operations for BSGS.

§4 Two Exponential Equations

Example 3.11 (P4884). Given a prime m\ge7 and 0<k<m, find the least positive integer n such that the decimal repunit satisfies

The integer identity 9(1+10+\cdots+10^{n-1})=10^n-1, together with invertibility of 9 modulo m, makes the original condition equivalent to

Since m\ne2,5, the base 10 is a unit and BSGS applies. There is no solution if the target is zero, or if it does not belong to \langle10\rangle.

Proof. Multiplication by the unit 9 gives the equivalence of the two congruences. Moreover, 0<k<m and m\ne3 imply 9k+1\not\equiv1\pmod m. Hence the least nonnegative solution of (3.10), if it exists, is already positive. \square

Remark. For more general inputs, if the target is 1, the least nonnegative solution is 0 and the least positive solution is the order of the base. Alternatively, write n=t+1 and find the least nonnegative t satisfying 10^t\equiv(9k+1)10^{-1}. When the modulus has a common factor with 9, multiplication of (3.9) by 9 gives a necessary condition; test the resulting candidates in the original congruence.

Example 3.12 (CF1106F). Let p=998244353. Given b_1,\ldots,b_k, suppose a sequence in \mathbf F_p^\times satisfies

Given n>k and m\in\mathbf F_p^\times, find f_k for which f_n=m. Choose the primitive root g=3 and write f_i=g^{z_i}. The exponents belong to \mathbf Z/(p-1)\mathbf Z and satisfy

Let T be the k\times k matrix whose first row is (b_1,\ldots,b_k), whose subdiagonal entries are 1, and whose remaining entries are 0. Then

Put A=(T^{n-k})_{1,1} and use BSGS to find B such that g^B=m. The problem becomes

Let d=\gcd(A,p-1). If d\nmid B, there is no solution. Otherwise, obtain a solution z_0 by Theorem 1.18; all possible values are

Proof. A primitive root makes z\mapsto g^z an isomorphism from the additive group \mathbf Z/(p-1)\mathbf Z to the multiplicative group \mathbf F_p^\times. Thus products correspond to sums, and powers to integer multiples. The matrix implements (3.12) at each step; among the initial exponents, only z_k can be nonzero, yielding (3.14). The solution classes of the linear congruence are distinct, and injectivity of the primitive-root parametrization makes the d values in (3.15) distinct as well. \square

Remark. Primality of the given constant is used as a known property of the problem; p-1=2^{23}\cdot7\cdot17. Binary exponentiation verifies

none of which equals 1, so Proposition 3.5 proves that 3 is a primitive root. Matrix operations are performed modulo the composite integer p-1, and (3.14) is solved by Theorem 1.18. Naive matrix multiplication and BSGS together use O(k^3\log(n+1)+\sqrt p) basic arithmetic and expected lookup operations. If A=0, the solvability condition is exactly B=0, and all p-1 nonzero initial values are then possible.

§5 Primality Testing

Lemma 3.13. In \mathbf F_p, the only possible solutions of x^2=1 are x=1 and x=-1.

Proof. One has (x-1)(x+1)=0, and a field has no zero divisors. For p=2, the two expressions denote the same element. \square

Definition 3.14. Let n>2 be odd, and write uniquely n-1=2^s d, where d is odd and s\ge1. For an integer a, if

then n is said to pass the strong pseudoprime test to base a.

Proposition 3.15. If n is an odd prime and n\nmid a, then it passes the test (3.17).

Proof. Fermat's theorem gives a^{2^sd}=1. Consider the chain of squares a^d,a^{2d},\ldots,a^{2^sd}. If the first term is not 1, take the first position where 1 occurs. The preceding term is a square root of 1 different from 1, and hence is -1 by Lemma 3.13; its index is less than s. \square

Construction 3.16 (Miller–Rabin). For the integer n to be tested, first reject n<2, accept n=2, and reject even integers greater than 2. For each chosen base a, normalize its representative modulo n; if the selection rule ensures 1\le a<n, it can be tested directly. If 1<\gcd(a,n)<n, a witness to compositeness has been found. Otherwise, compute v=a^d\bmod n. If v=1 or v=-1, the round passes. In all other cases, square at most s-1 more times; pass if -1 occurs, and declare the input composite if it never occurs.

Proof. This procedure tests exactly (3.17). Proposition 3.15 shows that a valid base never causes a prime to be declared composite. If repeated squaring produces 1 from an element other than \pm1, Lemma 3.13 itself supplies a witness to compositeness. \square

Remark. One round requires O(\log n) modular multiplications. Rabin's counting theorem for strong pseudoprimes states that, for every odd composite n, at most one quarter of the bases in 1,\ldots,n-1 pass (3.17); see M. O. Rabin, Probabilistic Algorithm for Testing Primality, 1980, Theorem 1. Consequently, for t independently and uniformly selected bases, the probability that a fixed composite passes all rounds is at most 4^{-t}. This probability is calculated for a fixed composite and random bases.

A further classical result gives a deterministic range: for 2\le n<2^{64}, testing the bases

in succession suffices to determine primality. This follows because the least integer that is a strong pseudoprime to all twelve prime bases is 318665857834031151167461>2^{64}; see J. Sorenson and J. Webster, Strong Pseudoprimes to Twelve Prime Bases, Theorem 1.1. An implementation may first handle inputs at most 37, after which all these bases lie in [1,n-1]; alternatively, reduce the bases and skip zero remainders. Modular multiplication for unsigned 64-bit inputs uses a wider product representation or a dedicated modular-multiplication implementation.

§6 Extracting Factors from Collisions

Construction 3.17 (One Pollard–Rho attempt). For an odd composite n, choose an integer c and an initial value x_0, and put, in \mathbf Z/n\mathbf Z,

Starting with x=y=x_0, repeatedly put x\leftarrow F(x) and y\leftarrow F(F(y)), and compute

If d=1, continue; if 1<d<n, return d; if d=n, the attempt has failed, so choose new parameters and retry.

Proposition 3.18. Every value returned by Construction 3.17 is a nontrivial factor of n. If a prime divisor p\mid n satisfies x\equiv y\pmod p, then p\mid d.

Proof. The first assertion follows from the definition of the greatest common divisor and the return condition. The second holds because p divides both x-y and n. \square

Remark. The nontrivial factor returned may be prime or composite. The algorithm iterates modulo the known integer n to seek collisions modulo an unknown prime divisor p. Since F is an integer polynomial, reduction of its orbit modulo n to modulo p gives an orbit of the same polynomial. The slow and fast pointers detect cycles in the orbit of a map on a finite set. If they coincide modulo n, then (3.20) gives n, and the attempt is treated as unsuccessful.

Construction 3.19 (Recursive factorization). For n>0, return the empty family if n=1. If the primality test confirms that n is prime, return n. Otherwise, extract a factor d, factor d and n/d separately, and combine their prime factors with multiplicities. Factors 2 may first be removed directly.

Proof. Provided each extracted factor satisfies 1<d<n, both subproblems are strictly smaller. Inductively, the product of the returned factors is n. When all leaves are indeed prime, unique factorization ensures that the output is correct and that the multiplicities are unique. \square

Remark. In a model approximating the orbit modulo p by a random mapping, the birthday-collision heuristic predicts a collision after about O(\sqrt p) iterations. If p is the least prime divisor of the composite integer n, then p\le\sqrt n, giving a heuristic iteration scale of O(n^{1/4}). This approximation replaces the actual quadratic-polynomial family by uniformly random mappings. The integer divisions and bit costs of the greatest common divisor computed at each iteration are counted separately. Degenerate parameters lead to new choices of c,x_0; one may also impose an iteration limit on each attempt.

With a Miller–Rabin base set having a deterministic input range, the recursive stopping criterion is deterministic within that range. With random bases, errors in primality testing can propagate to the final factorization. Extracted factors can be checked by exact division; verification of a complete prime factorization checks both the product and the primality of each factor.

Construction 3.20 (Batch greatest common divisors). Suppose a collision search produces consecutive differences \Delta_1,\ldots,\Delta_b. First compute

If d=1, the block yields no factor, so proceed to the next block. If 1<d<n, return d directly. If d=n, return to the start of this block, regenerate the \Delta_j, and compute d_j=\gcd(\Delta_j,n) one at a time; return as soon as 1<d_j<n. If the entire block yields no nontrivial factor, the current parameters have failed, so restart.

Proof. An integer product is relatively prime to n if and only if every factor is relatively prime to n: the forward implication follows because each factor divides the product, and the reverse follows from the prime-divisor criterion. Thus skipping a block with d=1 loses no factor. If 1<d<n, then d itself is a valid factor, whether it comes from one difference or several. If d=n, it does not follow that every difference has greatest common divisor n with n. For example, with n=15 and differences 3,5, the reduced product is zero, while each difference yields a factor. Replaying the block individually recovers factors merged by the product. \square

Remark. Reset Q to 1 for each block. The orbit state at the start of the block and its length suffice for replay. Batching normally uses one greatest-common-divisor call per b differences, with additional calls during replay after failure. Reducing the product modulo n at each step preserves the greatest common divisor in (3.21). A zero difference makes the batch product zero and triggers replay.

Construction 3.21 (Brent's segmented orbit). Fix F(x)=x^2+c\pmod n and a maximum block length b\ge1. Put y=x_0 and r=1. One outer iteration is as follows.

1. Save x=y, then advance y by r consecutive steps.
2. Over the next r steps, produce a difference x-y after each advance. Group at most b steps into each block and test it by Construction 3.20. Save the start of each block and its actual number of steps for replay.
3. Return if a nontrivial factor is found; restart if the current parameters are declared unsuccessful. Otherwise, put r\leftarrow2r and repeat the outer iteration.

An advance means the assignment y\leftarrow F(y). Every difference and product may immediately be normalized modulo n.

Proposition 3.22. Every return from Construction 3.21 gives a valid factorization. For any p\mid n, suppose the orbit modulo p has a preperiod of length \mu and a cycle of length \lambda. If the fixed reference point x of a segment is already in the cycle and r\ge\lambda, then at least one difference produced in that segment is divisible by p.

Proof. Correctness of the returned factor follows from Construction 3.20. Once the fixed reference point x lies in the cycle, any subsequent \lambda consecutive orbit points include the residue class of x modulo p. The segment tests r\ge\lambda consecutive points, so at least one collision modulo p occurs. As r doubles, the orbit index of the reference point increases and eventually exceeds \mu; provided no restart occurs, one eventually also has r\ge\lambda. \square

Remark. Simultaneous collisions modulo several prime factors, or an orbit that closes too early modulo n, may yield d=n, in which case the search restarts. Brent's segmentation keeps the reference point unchanged throughout a segment, facilitating batch products and replay; it changes constants and the frequency of greatest-common-divisor calls, while its iteration scale is still estimated heuristically by the preceding random model. A returned factor may contain only part of a prime-factor multiplicity; recursively factoring the remaining quotient recovers the rest.

Chapter IV — Arithmetic Functions and Convolution

§1 Algebra Indexed by Multiplication of Positive Integers

Throughout this chapter, R is a commutative ring with identity. Integer coefficients are interpreted through the unique unital ring homomorphism \mathbf Z\to R; unless otherwise stated, all divisors are positive. Division in a ring means only multiplication by an inverse known to exist.

Definition 4.1. An arithmetic function over R is a map f:\mathbf N_{>0}\to R. Let \mathcal A_R denote the set of all such functions. Define pointwise addition and Dirichlet convolution on this set by

Put \varepsilon(1)=1_R and \varepsilon(n)=0_R for n>1. Pointwise multiplication is denoted separately by fg; it is a different operation from f*g.

Proposition 4.2. Under pointwise addition and convolution, \mathcal A_R is a commutative ring with identity, whose multiplicative identity is \varepsilon.

Proof. Every positive integer has only finitely many divisors, so the sums in the definition exist. Commutativity follows from (a,b)\mapsto(b,a). Associativity follows from two groupings of the same finite sum:

Distributivity follows term by term from distributivity in R. Finally, convolution with \varepsilon leaves only the term with a=1 or b=1. \square

Theorem 4.3. A function f\in\mathcal A_R is invertible under convolution if and only if f(1) is a unit of R. In that case its inverse g=f^{*-1} is unique and is determined by the recurrence

Proof. If f*g=\varepsilon, evaluation at 1 gives f(1)g(1)=1_R, proving necessity. Conversely, n/d<n in (4.3), so g can be defined in increasing order of n. Substitution into the convolution gives f*g=\varepsilon. Uniqueness of the inverse follows from associativity. \square

In the ring of integer-valued arithmetic functions, the condition is f(1)=\pm1; in the residue class ring modulo M, it is \gcd(f(1),M)=1.

§2 Multiplicativity and Local Data at Prime Powers

Definition 4.4. A function f is multiplicative if f(1)=1_R and f(ab)=f(a)f(b) for coprime positive integers a,b. The function f is completely multiplicative if this equality holds for all positive integers a,b.

When R is nonzero, the requirement f(1)=1_R excludes the zero function and makes the empty product agree with the value at 1. In a ring with zero divisors, one cannot derive this condition by arbitrarily cancelling a nonzero function value.

Proposition 4.5. A multiplicative function is uniquely determined by its values at all prime powers. More precisely, given any family c_{p,e}\in R, where p ranges over the primes and e\ge1, there exists a unique multiplicative function satisfying f(p^e)=c_{p,e}, and

It is completely multiplicative if and only if c_{p,e}=c_{p,1}^{\,e} for all p,e.

Proof. The fundamental theorem of arithmetic expresses n uniquely as a product of pairwise coprime prime powers. Thus (4.4) is necessary, and the formula directly defines a multiplicative function. Complete multiplicativity allows one to remove a factor p repeatedly. Conversely, additivity of prime exponents ensures that function values multiply when any two integers are multiplied. \square

Proposition 4.6. The pointwise product, Dirichlet convolution, and convolution inverse of multiplicative functions are multiplicative.

Proof. The assertion for pointwise products follows directly from the definition. If \gcd(m,n)=1, every divisor of mn is uniquely of the form ab, where a\mid m and b\mid n. Hence

For the inverse, prescribe g(p^e) by (4.3) for each prime p, then extend these values to a multiplicative function by Proposition 4.5. The convolution f*g is multiplicative, takes the value 1_R at 1, and takes the value 0_R at every positive power of a prime. It therefore equals \varepsilon. \square

Definition 4.7. For a prime p, the series

is called the Bell series of f at p. Here p^0=1, and the series is a formal power series.

Proposition 4.8. For any arithmetic functions f,g,

If f is multiplicative, its Bell series jointly determine f; checking prime powers suffices to verify equality of two multiplicative functions.

Proof. The divisors of p^e are precisely 1,p,\ldots,p^e, so the coefficient of X^e on either side is \sum_{j=0}^ef(p^j)g(p^{e-j}). The last assertion follows from Proposition 4.5. \square

For a general arithmetic function, the prime-power data and other values such as f(6) may be prescribed independently.

§3 Basic Functions and Inversion

Definition 4.9. For k\in\mathbf N, define

Write \operatorname{id}=\operatorname{id}_1 and \tau=\sigma_0. Define the integer-valued functions \varphi and \mu by

The empty-product convention gives \varphi(1)=\mu(1)=1. Mapping these integer values into R gives the R-valued functions used in this chapter.

The functions \mathbf1, \operatorname{id}_k, and \varepsilon are completely multiplicative; \mu is multiplicative, but in general is not completely multiplicative. By definition, \sigma_k=\operatorname{id}_k*\mathbf1. Hence \sigma_k is multiplicative, and

The two formulas in (4.10) give sums of powers of divisors and divisor counts, respectively. Finite geometric sums are valid over any R; a quotient involving p^k-1 is evaluated by exact integer division, or by multiplication by the inverse when the denominator is a unit of R.

Proposition 4.10. One has

The last expression is first interpreted in the rational numbers; its value is an integer. The function \varphi is multiplicative.

Proof. Classify a\in[1,n] by d=n/\gcd(a,n). Writing a=(n/d)b gives b\in[1,d] and \gcd(b,d)=1, so the class indexed by d has \varphi(d) elements. This proves the convolution identity. Among the integers in the interval for a prime power, those not coprime to p^e are exactly the multiples of p, of which there are p^{e-1}. When m,n are coprime, the Chinese remainder theorem gives a bijection between the invertible residue classes modulo mn and the product of the sets of invertible residue classes modulo m,n, respectively. Thus \varphi(mn)=\varphi(m)\varphi(n). The product formula now follows from the prime-power formula. \square

Proposition 4.11. One has

Proof. If n has r distinct prime factors, the nonzero terms correspond precisely to subsets of these prime factors, and their sum is (1-1)^r. For n=1, only the empty subset occurs, contributing 1. Substitute n=\gcd(a,b) to obtain the last equality. \square

Theorem 4.12 (Divisor inversion). For any f,g\in\mathcal A_R, the following two identities are equivalent:

In particular, \varphi=\operatorname{id}*\mu.

Proof. The first identity is g=f*\mathbf1. Convolve both sides with \mu, and use associativity and (4.12) to obtain the second. The reverse implication is the same. \square

Inversion does not require f,g to be multiplicative. Multiplicativity permits reconstruction from local data, whereas inversion is an invertible operation in the convolution ring; the two have different roles.

Theorem 4.13 (Finite inversion over multiples). Fix N\ge1, and let f,g:[1,N]\to R. The following two identities are equivalent:

Proof. Substitute the first identity into the right-hand side of the second, and group terms with jk=t. This gives

The two linear transformations commute, since both compositions group terms by the product jk. The reverse implication therefore holds as well. \square

A finite truncation or finite support makes the sums over multiples finite, so that (4.14) is defined over any coefficient ring.

Proposition 4.14 (Completely multiplicative weights). If t is completely multiplicative, the map T_t(f)=tf preserves addition, convolution, and the convolution identity. In particular,

Proof. When ab=n, one has t(a)t(b)=t(n), which can be taken outside the finite sum. Then take f=\mathbf1 and g=\mu. \square

No invertibility of the individual values t(n) is required. If t is merely multiplicative, a,b need not be coprime, and the first step is invalid.

Theorem 4.15 (Weighted inversion of sums with floor arguments). Let t be completely multiplicative, and let F,G:\mathbf N\to R satisfy F(0)=G(0)=0. Then

Proof. For positive integers i,j, one has \lfloor\lfloor n/i\rfloor/j\rfloor=\lfloor n/(ij)\rfloor. Indeed, for an integer q, the inequality qj\le\lfloor n/i\rfloor is equivalent to qij\le n. Consequently, the composition of the two summation operators with floor arguments is the operator with coefficients t*(\mu t). By (4.16), this composition is the identity operator. \square

Here F,G may be arbitrary sequences; multiplicativity is unnecessary. Taking prefix sums of a convolution identity naturally leads to this form of inversion.

§4 Formal Dirichlet Generating Functions

Definition 4.16. Introduce formal symbols \mathsf e_n for n\ge1, with the rule \mathsf e_a\mathsf e_b=\mathsf e_{ab}. A formal Dirichlet series is an expression

Addition is defined term by term, and multiplication by the preceding rule. One may also write n^{-s} for \mathsf e_n; in this section, s is only a formal notation indicating multiplication of indices.

The coefficient of a fixed \mathsf e_n involves only finitely many factorizations of n. Thus multiplication requires no convergence condition, and

Introduce a variable X_p for each prime p, and put \mathsf e_n=\prod_pX_p^{\nu_p(n)}. Each monomial involves only finitely many variables; in the product of two formal sums, any prescribed monomial admits only finitely many decompositions. In this precise formal sense, convolution may also be regarded as multivariate multiplication in coordinates given by prime exponents.

Theorem 4.17 (Formal Euler product). Suppose f(1)=1_R. Then f is multiplicative if and only if the following identity holds coefficient by coefficient:

Proof. The coefficient of \mathsf e_n on the right is \prod_{p\mid n}f(p^{\nu_p(n)}). Except for the prime factors of n, every factor must contribute its constant term f(1)=1_R. Thus only finitely many choices are involved in this coefficient. Compare coefficients and apply Proposition 4.5. \square

For example, the local factors are

These denominators all have constant term 1_R and hence are invertible in R[[X]]. For complex s, an analytic interpretation of these formal series requires convergence and a valid interchange of sums.

§5 Prefix Transforms on the Divisibility Order

Definition 4.18. For f:[1,N]\to R, define the divisor zeta transform and the multiple zeta transform by

The partial order here is divisibility. Prime exponents give its coordinates: a\mid b if and only if \nu_p(a)\le \nu_p(b) for every prime. The two transforms are therefore finite downward and upward accumulations, respectively, in these coordinates.

Construction 4.19. Suppose the primes not exceeding N are known. Initialize A[n]=f(n). Enumerate the primes p\le N one at a time; for each p, in the increasing order i=1,2,\ldots,\lfloor N/p\rfloor, perform

After all operations, A=Z_\downarrow f. To recover the original array, enumerate the primes in reverse order and, for each p, perform A[ip]\leftarrow A[ip]-A[i] in the order i=\lfloor N/p\rfloor,\ldots,1.

Proof. Fix p and write each index uniquely as mp^e, where p\nmid m. Increasing i performs a one-dimensional prefix sum along each chain

The inverse transform implements convolution with \mu directly through adjacent differences along the chains; the specified traversal direction propagates contributions through the prime powers in dependency order.

Construction 4.20. For the transform over multiples, again initialize A=f. For each prime p, in the order i=\lfloor N/p\rfloor,\ldots,1, perform

The result is Z_\uparrow f. For the inverse, enumerate the primes in reverse order and, for each prime, perform A[i]\leftarrow A[i]-A[ip] in increasing order of i.

Proof. Again decompose the indices into chains mp^e, now taking suffix sums along each chain. By induction on the processed primes, A[n] is the sum of f(m) over all n\mid m\le N such that every prime factor of m/n belongs to the processed set. Correctness of the inverse follows by undoing the operations one by one. \square

Each transform and its inverse perform

ring additions or subtractions and use O(N) storage. The bound for the sum of reciprocals of primes will be proved in Proposition 6.3. Enumerating all multiples of every integer directly requires \sum_{d\le N}\lfloor N/d\rfloor=O(N\log(N+1)) operations. The algorithm accumulates along the prime coordinates of the divisibility relation one coordinate at a time and applies to arbitrary input arrays.

Chapter V — Summation over Finite Hyperbolic Regions

§1 Quotient-Value Sets and Their Closure

In this chapter, summation indices range over finite sets. Floors, comparisons, and endpoint calculations are performed in the integers; weights may lie in a commutative ring R with identity. For an arithmetic function f, write

Definition 5.1. For n\ge1, define the quotient-value set

Also put D_0=\varnothing. Quotient-value sets consist of positive integers; any prefix value at zero needed by an algorithm is stored separately.

Lemma 5.2. For a\in\mathbf N and b,c\ge1,

If m\in D_n, then D_m\subseteq D_n.

Proof. For any q\in\mathbf N, the condition q\le\lfloor\lfloor a/b\rfloor/c\rfloor is equivalent to qbc\le a, and hence to q\le\lfloor a/(bc)\rfloor. The two sides are therefore equal. If m=\lfloor n/b\rfloor and 1\le c\le m, then bc\le n, so \lfloor m/c\rfloor=\lfloor n/(bc)\rfloor\in D_n. \square

The latter conclusion keeps the states arising in recursive prefix-sum computations within a single finite set. It is the basis for organizing states in the sublinear sieves considered later.

Proposition 5.3. Put u=\lfloor\sqrt n\rfloor. Then

Proof. If i>u, then \lfloor n/i\rfloor\le u, since n<(u+1)^2. The remaining quotient values already appear in the second set. \square

Proposition 5.4. For q\in D_n, all indices satisfying \lfloor n/i\rfloor=q form the interval

If an unprocessed interval begins at l\le n, the right endpoint of the quotient block containing l is given by

Proof. The equality \lfloor n/i\rfloor=q is equivalent to q\le n/i<q+1, that is, n/(q+1)<i\le n/q. Taking integer endpoints gives the result. \square

Starting with l=1, successively compute q,r, process [l,r], and put l=r+1. This visits all maximal quotient blocks exactly once. If the upper limit is L\le n, simply truncate r to \min(r,L).

§2 Decomposition into Quotient Blocks

Theorem 5.5. Let f:[1,n]\to R and g:D_n\to R. Then

where the sum on the right runs over maximal quotient blocks. If each evaluation of S_f and g costs O(1), the summation requires O(\sqrt n) arithmetic operations and O(1) additional storage cells. Preprocessing, such as the construction of prefix tables, and its storage must be counted separately.

Proof. The argument of g is constant on each block, so its value can be taken outside the sum over that block. The blocks are disjoint and cover [1,n], and their number is bounded by Proposition 5.3. \square

Proposition 5.6. For n,i\ge1, one has \lceil n/i\rceil=\lfloor(n-1)/i\rfloor+1. Consequently, for n\ge1,

If a summand depends simultaneously on \lfloor a_1/i\rfloor,\ldots,\lfloor a_t/i\rfloor, a common refinement is obtained by computing the right endpoints of the respective quotient blocks at each step and taking their minimum.

Proof. Euclidean division verifies the first equality directly; the second separates the term i=n. In a common refinement, the dependence of the summand on these quotient values remains unchanged as long as none of them changes. The first change of any quotient value occurs immediately after the current common block. \square

Use (5.6) for positive quotient values. If a_j<i, the corresponding quotient is zero and remains zero thereafter; its right endpoint is taken to be the overall upper limit of the sum. For fixed t, the number of common blocks is at most the sum of the numbers of blocks in the individual partitions, and hence is O(\sum_j\sqrt{a_j}+t).

§3 The Hyperbola Identity

Proposition 5.7. For f,g\in\mathcal A_R and n\in\mathbf N,

Proof. The first equality assigns each ordered pair (a,b) with product k\le n to the term indexed by k; the second fixes the first coordinate. All sets are finite, so the order of summation may be interchanged. \square

Theorem 5.8 (Dirichlet hyperbola identity). Let u,v\in\mathbf N satisfy

Then

Proof. In the region H_n=\{(a,b)\in\mathbf N_{>0}^2:ab\le n\}, the inequalities a\ge u+1 and b\ge v+1 cannot both hold. Thus the two subregions a\le u and b\le v cover H_n. Their intersection is exactly [1,u]\times[1,v], since uv\le n. Inclusion–exclusion for these two finite sets gives (5.11). \square

Corollary 5.9. One may take u=v=\lfloor\sqrt n\rfloor. Alternatively, choose any 1\le u\le n and put v=\lfloor n/u\rfloor. If individual function values and prefix values are available in constant time, evaluation of (5.11) costs O(u+v), or O(\sqrt n) for the symmetric choice.

Proof. The symmetric choice plainly satisfies (5.10). For the other choice, uv\le n<u(v+1)<(u+1)(v+1). Count the terms to obtain the bound. \square

Quotient blocks and the hyperbola decomposition partition the same finite product region in different ways. The former uses interval sums of one function and the prefix function of the other; the latter uses prefix values on both sides. If the evaluation costs on the two sides differ, an asymmetric choice of u,v may give a more suitable actual cost.

§4 Summation Examples Involving Divisors and Remainders

Example 5.10 (Mathematical formulation of P3935). Compute the prefix sum of the divisor-counting function, D(n)=\sum_{k=1}^n\tau(k). Since \tau=\mathbf1*\mathbf1,

Proof. The first sum counts ordered pairs (d,e) satisfying de\le n; the second equality is the symmetric case of Theorem 5.8. \square

The two expressions can be evaluated using quotient blocks or direct enumeration up to the square root, respectively. Both take O(\sqrt n) time and O(1) extra space. If the required sum is over [L,U], the answer is D(U)-D(L-1). The bound in the original slides reaches 1.6\times10^{14}; use an integer type that accommodates the result and intermediate products.

Example 5.11 (Mathematical formulation of P2261). For n,k\in\mathbf N, compute

For k=0, the result is zero. For k>0, a block [l,r]\subseteq[1,\min(n,k)] on which the quotient q is constant requires subtraction of

Proof. Euclidean division gives k\bmod i=k-i\lfloor k/i\rfloor; the quotient is zero when i>k. Within a block, apply the formula for the sum of an arithmetic progression. \square

The complexity is O(1+\min(n,\sqrt k)), with O(1) extra space. Division by 2 in (5.14) is exact integer division. For an even modulus, first cancel the integer factor, then reduce modulo that modulus.

Example 5.12 (Mathematical formulation of ARC068E). Given N integer intervals [l_j,r_j]\subseteq[1,M], determine, for every d\in[1,M], how many of them contain at least one positive multiple of d.

The contribution of the j-th interval to the answer is

Proof. There are \lfloor r_j/d\rfloor positive multiples not exceeding r_j, exactly \lfloor(l_j-1)/d\rfloor of which lie to the left of the interval. Their difference counts the multiples in the interval; the contribution records only whether this number is positive. \square

For each interval, form common blocks for the two quotients. If a common block [a,b] satisfies (5.15), add one to \Delta[a] and subtract one from \Delta[b+1] in a difference array for the answers. After processing all intervals, take prefix sums. This takes O(N\sqrt M+M) time and O(M) space. For any fixed d, each interval contributes at most once, even if it contains several multiples.

Construction 5.13 (Decomposition by length for the same example). Let \ell_j=r_j-l_j+1 be the interval length. Process d=1,2,\ldots,M in increasing order. Insert all intervals satisfying \ell_j<d into an array supporting addition of one on an interval and point queries, and let c_d(x) denote the current coverage count at each position. Then

Proof. Any d consecutive integers contain exactly one multiple of d, so an interval of length at least d necessarily contributes. An interval of length less than d contains at most one multiple. Hence summing coverage counts over the multiples counts every contributing short interval exactly once. \square

Bucket the intervals by length, and use a Fenwick tree to maintain prefix sums of the difference array. Each interval is inserted once, and each insertion or query costs O(\log(M+1)). Since \sum_{d\le M}\lfloor M/d\rfloor=O(M\log(M+1)), the total complexity is O(N\log(M+1)+M\log^2(M+1)), with O(N+M) space. This gives a definite bound in addition to the quotient-block solution in the slides and is suitable when both N and M are large.

Example 5.14 (Mathematical formulation of P2260). For n,m\ge1, compute

Put h=\min(n,m). Separating the diagonal from all ordered pairs gives

On a block [l,r] where both q_n=\lfloor n/i\rfloor and q_m=\lfloor m/i\rfloor are constant, the diagonal contributes

Proof. Equation (5.18) deletes i=j from the rectangular index set. Equation (5.19) expands (n-iq_n)(m-iq_m) and sums term by term. \square

Using \sum_{i=1}^xi=x(x+1)/2 and \sum_{i=1}^xi^2=x(x+1)(2x+1)/6, each block requires a constant number of arithmetic operations. Thus the total complexity is O(\sqrt n+\sqrt m), with O(1) extra space; this bound does not require n,m to have the same order of magnitude. The original problem uses modulus 19940417. One may cancel 2,6 exactly among the integer factors before reduction, or use modular inverses after verifying invertibility.

§5 Combining Inversion with Quotient Blocks

Example 5.15 (Mathematical formulation of P6156). Fix a positive integer k. For several values n\le N, compute

Here \mu^2 denotes the pointwise square. Define

Then

Proof. First classify by d=\gcd(i,j) and write i=da,j=db. By (4.12),

Put t=de and interchange the finite sums. Since d^{k+1}e^k=t^kd, the coefficient of F_k(\lfloor n/t\rfloor) is t^k\sum_{d\mid t}d\mu^2(d)\mu(t/d)=h(t). \square

Proposition 5.16. Both g,h are multiplicative, and

If the twice-iterated prefix sums are defined by

then

Proof. The function \operatorname{id}\,\mu^2 is multiplicative, so its convolution with \mu is multiplicative. Its pointwise product with the completely multiplicative function \operatorname{id}_k gives the multiplicativity of h. At a prime p, the Bell series of \operatorname{id}\,\mu^2 is 1+pX. Multiplication by the Bell series 1-X of \mu gives 1+(p-1)X-pX^2, proving (5.24).

For (5.26), classify by s=a+b. There are s-1 ordered pairs when 2\le s\le x+1, and 2x-s+1 when x+2\le s\le2x. On the other hand, the coefficient of s^k in T_k(2x)-2T_k(x) is s-1 for 1\le s\le x and 2x-s+1 for x+1\le s\le2x, agreeing with these counts. \square

Equation (5.26) gives preprocessing using only addition, subtraction, and multiplication: compute s^k for s\le2N, accumulate to obtain U_k,T_k, and then compute F_k(x) for x\le N. This procedure applies in the ring modulo 2^{32}.

The linear sieve of Chapter VI can use (5.24) to compute g in O(N) operations and obtain the prefix table H(x)=\sum_{t\le x}h(t) of h. Equation (5.22) then becomes

requiring O(\sqrt n) time per query. If each s is exponentiated separately by binary exponentiation, preprocessing takes O(N\log(k+1)) time. Using the complete multiplicativity of s\mapsto s^k, one may instead compute p^k only for primes p\le2N, then multiply at every transition of the linear sieve. This gives O(N+\pi(2N)\log(k+1)) for that part. All tables together occupy O(N) space.

Chapter VI — Prime Sieves on Finite Intervals

§1 Eliminating Composite Numbers

Lemma 6.1. Every composite integer x\ge2 has a prime factor not exceeding \sqrt x.

Proof. Write x=ab with 1<a\le b<x. Then a\le\sqrt x, and any prime factor of a has the required property. \square

Construction 6.2 (Sieve of Eratosthenes). For N\in\mathbf N, create a Boolean array P[0],\ldots,P[N]. Set entries 0,1 to false wherever those indices exist, and all other entries to true. Scan p=2,3,\ldots,\lfloor\sqrt N\rfloor: if P[p] is true, set the entries at

to false. At the end, the indices whose entries remain true are exactly the primes not exceeding N.

Proof. Every marked integer has a factor p and is at least p^2, so it is composite. Conversely, let p be the least prime factor of a composite integer x. By Lemma 6.1, p\le\sqrt x\le\sqrt N, and x/p\ge p. Hence x is marked when p is processed. Primes are never marked, so the flag at p is true when the scan reaches it. \square

It suffices to start at p^2: a composite multiple less than p^2, written as ap, has a<p and has already been processed by a prime factor of a. In an implementation, the loop condition can be written as p\le N/p to avoid overflow from computing p^2 first.

Proposition 6.3. For x\ge2,

Proof. Only an upper bound is needed. Put s=1+1/\log x>1. For p\le x, one has p^{s-1}\le e, whence

The second inequality follows from -\log(1-y)\ge y for 0<y<1. Expanding the finite product of convergent geometric series by unique factorization gives a subseries of the positive-term series \sum_{n\ge1}n^{-s}. Thus the product is at most

Substitution into (6.3) proves the result. A convergent real series with s>1 is used here only to estimate complexity; its convergence has already been ensured by the integral bound. \square

Corollary 6.4. Construction 6.2 takes O(N\log\log(N+2)) time and O(N) space.

Proof. Initialization and the final scan take O(N) time. Each prime p causes at most \lfloor N/p\rfloor markings, so Proposition 6.3 gives the total bound. \square

§2 Sieving and Factorization on an Interval

Construction 6.5 (Segmented prime sieve). Given 1\le L\le U and the primes not exceeding \lfloor\sqrt U\rfloor, initialize a Boolean array on [L,U]. For each prime p in the list, start at

and mark all multiples of p not exceeding U. If 1\in[L,U], also mark 1 as nonprime. The remaining elements are exactly the primes in the interval.

Proof. The lower bound p^2 prevents removal of the prime p itself. Every composite integer x in the interval has a least prime factor p\le\sqrt U and satisfies x\ge p^2, so it is necessarily marked. The other elements cannot be removed incorrectly by the construction. \square

Put W=U-L+1 and B=\lfloor\sqrt U\rfloor. If the prime list is already available, the scan takes

with O(W) extra space. Constructing the base prime list additionally requires O(B\log\log(B+2)) time and O(B) space. The term \pi(B) in (6.6) counts the scan of the base prime list.

Construction 6.6 (Complete factorization of the integers in an interval). For every x\in[L,U], store a residual value a_x=x and an initially empty factor list. For each prime p\le B, visit all its multiples x in the interval, starting at \lceil L/p\rceil p. Repeatedly set a_x\leftarrow a_x/p until p\nmid a_x, recording the number of divisions. After all primes have been processed, if a_x>1, append a_x to the factor list as a prime factor of exponent one.

Proof. Every division is exact integer division. When a prime factor p is processed, its entire exponent in x is removed; earlier removal of other prime factors does not change that exponent. If the final a_x>1 were still composite, it would have a prime factor not exceeding \sqrt{a_x}\le\sqrt U, contradicting the removal of all such prime factors. Nor can the square of the residual prime factor divide x: that prime would then not exceed \sqrt U, again contradicting its having been completely removed. \square

Starting at the first multiple in the interval also records the factor p when x=p. For x=1, no division is performed and the factor list is empty.

For each x, the number of divisions is at most \lfloor\log_2x\rfloor. A direct uniform upper bound is therefore

integer operations, in addition to preprocessing the base prime list. Space for outputting all prime factors and exponents must be included. A coarse bound is O(W\log(U+1)); processing subsequent tasks for each integer immediately can reduce the amount retained for later use.

§3 Generating Each Composite Integer Uniquely

Definition 6.7. For n\ge2, let \operatorname{lpf}(n) denote the least prime factor of n. The canonical generating pair of a composite integer x is

Then p\le\operatorname{lpf}(i), with i\ge2.

Construction 6.8 (Linear sieve). Initialize \operatorname{lpf}[2..N]=0 and an empty prime list. Scan i=2,3,\ldots,N: if \operatorname{lpf}[i]=0, set \operatorname{lpf}[i]=i and append i to the prime list. Then enumerate the primes p in the list in increasing order. As long as

set \operatorname{lpf}[ip]=p. Terminate the current inner loop when the conditions cease to hold.

Theorem 6.9. Construction 6.8 produces the correct least-prime-factor table and generates each composite integer exactly once. Its time and space requirements are both O(N).

Proof. Any generating pair satisfies p\le\operatorname{lpf}(i), so the least prime factor of ip is exactly p; the assigned value is correct. Conversely, for a composite integer x, take its canonical generating pair (i,p). Since i<x, the index i is processed before the scan reaches x. Moreover, p\le\operatorname{lpf}(i)\le i, so p is already in the prime list and satisfies the conditions of the inner loop. Thus x must be generated. Any pair generating x must have p=\operatorname{lpf}(x) and hence i=x/p, so the generation is unique.

Consequently, an index i not yet generated is precisely a prime. The number of successful inner-loop iterations equals the number of composite integers, with at most one additional termination test per outer index. Initialization and prime insertion are also linear, giving total time O(N). \square

A common implementation first generates ip and then exits immediately if p\mid i. Since the primes are in increasing order and the first prime dividing i is \operatorname{lpf}(i), this is equivalent to (6.9).

§4 Propagating Function Values Along Canonical Factorizations

Construction 6.10. For each i\ge2, maintain the unique factorization

and store c_i,e_i and p^{e_i}. For a prime i=p, take c_i=1,e_i=1. A transition i\mapsto ip in the linear sieve has only two cases:

Proposition 6.11. Let f be multiplicative. Suppose that establishing the local state at each prime costs C_p, and that the local state for f(p^{e+1}) can be obtained from that for f(p^e) in O(1) operations. Then the linear sieve can use

to compute all f(1),\ldots,f(N) at total cost O(N+\sum_{p\le N}C_p) and with O(N) space, provided that each local state occupies O(1) storage cells.

Proof. The two factors in the canonical factorization are coprime, so (6.12) is valid. The first type of transition uses only the local value already established at the prime. The second extends an existing local state, and the function value at c_i<i is already known. Each composite integer requires one constant-time step; adding the prime initialization costs gives the stated bound. \square

Retaining f(c_i) and the local prime-power state computes f(i) by multiplication and applies over any coefficient ring. The total cost is the sum of canonical generation and local-state computation costs.

Example 6.12 (\mu and \varphi). Initialize \mu(1)=\varphi(1)=1. At primes, put \mu(p)=-1 and \varphi(p)=p-1. In a transition of the linear sieve,

Proof. If p\nmid i, use multiplicativity. If p\mid i, then ip has the square factor p^2, so \mu(ip)=0. At the same time, increasing p^e to p^{e+1} changes the local value of \varphi from (p-1)p^{e-1} to (p-1)p^e. \square

The condition in the second branch is p\mid i, where i is the index before multiplication by the prime p.

Example 6.13 (Sums of powers of divisors). Fix k\in\mathbf N. For the local prime power in (6.10), maintain

These satisfy

Proof. The first two equalities respectively increase an exponent and add a summand; the third uses multiplicativity for the coprime factorization. \square

For k=0, it suffices to store e_i, and the formula becomes \tau(i)=\tau(c_i)(e_i+1). For k=1, store p^{e_i} and its geometric sum. Both functions can be computed in O(N) time. For an arbitrary input k, the prime initialization values p^k can be computed by binary exponentiation; Proposition 6.11 gives total cost O(N+\pi(N)\log(k+1)). Over the integers, one may also use the exact quotients \tau(i)/(e_i+1) or \sigma_1(i)/B_{p,e_i}; over a general coefficient ring, use the multiplicative recurrence (6.15).

Construction 6.14 (General form with prescribed prime-power values). Suppose all values f(p^e) for p^e\le N have been supplied and can be queried in constant time. Maintain the data in (6.10) and compute f(i) directly by (6.12). The additional time and space are both O(N).

Proof. Each step of the canonical factorization adds only one prime factor, and the local value can be read directly from the table. The prime powers are distinct integers in [1,N], so the local table contains at most N-1 entries. The cost of constructing it has been separated from this construction. \square

The enumeration by cp^e used in the slides can also be expressed through this canonical factorization; its total cost is the sum of integer generation and prime-power data preparation costs.

§5 An Application of Factorization on a Short Interval

Example 6.15 (Mathematical formulation of the strengthened Haitang and Math problem). Given 1\le n\le10^{12}, determine the number of solutions of

in positive integers m, where S(m) is the sum of the decimal digits of m.

Within the given range, S(m)\le108: a positive integer below 10^{12} has at most twelve digits, and the endpoint 10^{12} has digit sum 1. Since a remainder is less than the divisor, every solution also satisfies 1\le S(m)<m.

Proposition 6.16. Put B=\min(108,n-1). The solutions of (6.16) are in bijection with the accepted candidates in the following procedure: enumerate s\in[1,B], enumerate the positive divisors m of n-s, and retain those satisfying

as candidates.

Proof. If m solves the original problem, put s=S(m). Then n-s is a positive multiple of m, and s<m\le n, so s\le n-1. Conversely, if m\mid n-s and 0<s<m, the equality n=(n-s)+s is the Euclidean division decomposition with divisor m and remainder s. Together with S(m)=s, it gives the original equality. Moreover, m\le n-s<n, so the range condition holds automatically. Each m corresponds to the unique value s=S(m) and is not counted again for a different s. \square

It therefore suffices to factor all integers in [n-B,n-1] simultaneously. There are at most 108 positive integers in this interval; Construction 6.6 with a list of primes not exceeding 10^6 gives their factorizations. If n=1, then B=0, the candidate set is empty, and the answer is zero.

Once the factorization x=\prod_{j=1}^{r}p_j^{e_j} is known, enumerate its divisors through the Cartesian product map

Unique factorization ensures that this is a bijection. Successively adjoining each prime power enumerates all divisors in O(\tau(x)) multiplications. Computing each digit sum takes O(\log(n+1)) integer divisions; in the original range there are at most twelve digits, with the endpoint handled separately.

With the base prime list already preprocessed, a sufficient upper bound for one query is

The three terms respectively count the scan of the base primes, divisions during factorization of the interval, and checks of all candidate divisors. Divisor enumeration is charged according to the sum of divisor counts over the entire short interval. Multiple queries may share the prime list up to 10^6; the residual values on each interval and the candidate checks are handled separately for each query.

Chapter VII — Convolution and Division on Quotient-Value Sets

§1 Finite Quotient-Value Representations

Definition 7.1. Fix N\ge1 and put

Write the prefix sum of an arithmetic function f as F(x)=\sum_{1\le n\le x}f(n), and set F(0)=0. The indexed finite family (F(v))_{v\in\mathcal Q_N} is called the quotient-value representation of f at N, corresponding to the block sieve in the lecture notes. Here \mathcal Q_N is precisely D_N from Chapter V. The indexing retains the correspondence between arguments and values for lookup in the operations below.

Proposition 7.2. There is a disjoint decomposition

In particular, |\mathcal Q_N|=s+K\le2s. If v\in\mathcal Q_N and 1\le d\le v, then \lfloor v/d\rfloor\in\mathcal Q_N.

Proof. For 1\le t\le s, take i=\lfloor N/t\rfloor\ge t. From ti\le N<t(i+1) we obtain t\le N/i<t+t/i\le t+1, so t is a quotient value. The values greater than s are precisely those obtained from i\le K; in this range, successive values of N/i differ by at least 1, so their floors are distinct. Finally, if v=\lfloor N/i\rfloor, then

Whenever the right-hand side is positive, it still belongs to \mathcal Q_N. \square

A small quotient value can be addressed directly by v, and a large one by \lfloor N/v\rfloor. Thus the quotient-value table can be represented by two arrays of length O(\sqrt N), supporting worst-case constant-time queries.

Proposition 7.3. For an integer 0\le t\le N, write r_N(t)=|\mathcal Q_N\cap[1,t]|. Then

Proof. In the second case, a large quotient value satisfies \lfloor N/i\rfloor\le t if and only if i>N/(t+1). Count such indices in [1,K] and add all s small quotient values. \square

§2 Prefix Representations of Convolution

Proposition 7.4. Let h=f*g, with respective prefix sums H,F,G. For x\ge1 and u=\lfloor\sqrt x\rfloor,

Proof. Write H(x) as \sum_{ab\le x}f(a)g(b). If a,b>u, then ab\ge(u+1)^2>x. Every admissible pair therefore satisfies a\le u or b\le u, and the part counted twice is exactly [1,u]^2. \square

Construction 7.5. Given only the values of F,G on \mathcal Q_N, the quotient-value representation of H can be computed using O(N^{3/4}) ring operations and O(\sqrt N) storage cells.

Proof. Since [1,s]\subseteq\mathcal Q_N, successive differences give all the required values f(a),g(a). Apply (7.5) for every x\in\mathcal Q_N; Proposition 7.2 ensures that the remaining prefix queries are available. The total operation count is bounded by a constant multiple of

Estimate each sum by comparison with an integral. \square

Proposition 7.6. Suppose, in addition, that h(1),\ldots,h(B) and their prefix sums can be computed at total cost P(B), where \sqrt N\le B\le N. The remaining quotient values can then be computed by (7.5), at total cost

Proof. Quotient values greater than B have the form \lfloor N/i\rfloor, with i<N/B. Their total evaluation cost is at most

The small quotient values have already been obtained by preprocessing. \square

Remark. If a linear-time method is also available for precomputing the first B values and their prefix sums, so that P(B)=O(B), taking B\asymp N^{2/3} gives O(N^{2/3}). Construction 7.5 applies to two given quotient-value tables; with this additional preprocessing method, formula (7.7) applies.

§3 Division for Convolution

Theorem 7.7. Suppose h=f*g and g(1) is a unit in the coefficient ring. Then F is uniquely determined by

If [l,r]\subseteq[2,x] is a maximal interval on which the quotient \lfloor x/b\rfloor is constant, its contribution to the sum is

Proof. Summing over the second factor gives H(x)=\sum_{b\le x}g(b)F(\lfloor x/b\rfloor). Separate the term b=1 and multiply by g(1)^{-1} to obtain (7.8). Every unknown prefix on its right-hand side has argument at most \lfloor x/2\rfloor, so the values of F are uniquely determined in increasing order. Formula (7.9) follows by factoring out the common value on a constant interval. \square

Construction 7.8. Suppose G(x),H(x) can be evaluated with a constant number of elementary operations, and F(1),\ldots,F(B) can be precomputed in linear time. Memoized evaluation of (7.8)–(7.9) on quotient values greater than B is called the Du Jiao sieve. For \sqrt N\le B\le N, its time and space bounds are, respectively,

Proof. The recursive argument strictly decreases, so the computation terminates; Proposition 7.2 ensures that a computation starting at N visits only its quotient values. Each large state is evaluated once, processing at most O(\sqrt x) constant intervals. Sum over these states using the estimate from Proposition 7.6. Small states are looked up in the precomputed table, and there are at most N/B large states. \square

Taking B=\lceil N^{2/3}\rceil, capped at N, gives the bound O(N^{2/3}). For evaluations of G,H that take more than constant time, multiply the operation count by the corresponding evaluation costs.

Example 7.9. Write M(x)=\sum_{n\le x}\mu(n) and \Phi(x)=\sum_{n\le x}\varphi(n). From \mu*\mathbf1=\varepsilon and \varphi*\mathbf1=\operatorname{id}, for x\ge1 we have

Both expressions can be evaluated by Construction 7.8. The rational expression here is first evaluated as an integer, then mapped into the coefficient ring; this evaluation method also applies to even moduli.

§4 Completely Multiplicative Weights

Proposition 7.10. Let w be completely multiplicative, and let wf denote pointwise multiplication. Then

Proof. For every d\mid n, we have w(d)w(n/d)=w(n), which can therefore be factored out of the convolution sum. \square

Example 7.11. Put

Then

and

Proof. Substitute \gcd(i,j)=\sum_{d\mid i,\ d\mid j}\varphi(d), interchange the finite sums, and put i=du,j=dv to obtain (7.13). Applying Proposition 7.10 with w(n)=n^2 gives

Theorem 7.7 and \sum_{i\le x}i^3=T(x)^2 then give (7.14). \square

In (7.14), an interval weight is a difference of sums of squares; in (7.13), it is F(r)-F(l-1). Compute F by the Du Jiao sieve and then perform one quotient-block summation to obtain the computation required by P3768 in the lecture notes. Modulo a prime p>3, both 2,6 are invertible; for other moduli, the integer identities still permit cancellation before reduction.

Chapter VIII — Prime Weights and Decomposition by the Least Prime Factor

§1 A Finite Process for Removing Composite Numbers

Definition 8.1. Let w be completely multiplicative, and let p_1<p_2<\cdots be the sequence of primes. Put

where \operatorname{lpf}(n) denotes the least prime factor. Set A_0(v)=\sum_{n=2}^{v}w(n) and write P_w(t)=\sum_{p\le t}w(p).

Proposition 8.2. If p=p_j and v\ge p^2, then

If v<p^2, then A_j(v)=A_{j-1}(v).

Proof. The difference between the two stages consists exactly of the composite numbers whose least prime factor is p. They are uniquely expressible as pm, where m\ge p and \operatorname{lpf}(m)\ge p. Evaluating (8.1) at \lfloor v/p\rfloor\ge p and subtracting the smaller primes already retained enumerates precisely these m. Complete multiplicativity gives w(pm)=w(p)w(m), including when p\mid m. If v<p^2, there is no composite number to remove. \square

Construction 8.3. Initialize A(v)=\sum_{n=2}^{v}w(n) on \mathcal Q_N. Process primes p\le\sqrt N in increasing order; for each prime, apply (8.2) to the quotient values satisfying v\ge p^2 in decreasing order. On completion, A(v)=P_w(v).

Proof. Within a round, \lfloor v/p\rfloor<v. Thus, when updating in decreasing order of v, the smaller argument on the right still holds its value from the preceding round. Proposition 7.2 ensures that this argument is present in the table. Induction using Proposition 8.2 establishes the invariant at every stage. At completion, every composite n\le v\le N has a prime factor at most \sqrt N and has therefore been removed; primes have been retained throughout. \square

Obtain P_w(p-1) from a separate prefix table for the small primes already enumerated. Initialize by subtracting w(1)=1 from the full prefix sum, so that the candidate set starts at 2.

§2 Operation Counts

Lemma 8.4. Let \pi(x) denote the number of primes not exceeding x. Then \pi(x)=O(x/\log x) for x\ge2.

Proof. For an integer m\ge2, the product of all primes in (m,2m] divides \binom{2m}{m}<4^m. Taking logarithms bounds the number of primes in this interval by O(m/\log m). Partition (\sqrt x,x] into dyadic intervals: each denominator is \Omega(\log x), and the numerators sum to O(x). There are at most \sqrt x primes not exceeding \sqrt x. Combining the bounds proves the assertion; rounding affects only absolute constants. \square

Proposition 8.5. The number of updates in Construction 8.3 is

Proof. Put L=N^{1/4}. The contribution from p\le L is at most 2\sqrt N\pi(L)=O(N^{3/4}/\log N). For p>L, all required quotient values exceed \sqrt N, and their number is at most N/p^2. Grouping the primes into intervals (2^jL,2^{j+1}L] and applying Lemma 8.4 gives

Add the two contributions. \square

The cost of generating the small primes and computing all initial values and local weights must be added to this bound. If these costs are O(\sqrt N) and weight evaluation and ring operations are elementary operations, the overall bound remains (8.3), with O(\sqrt N) space.

Example 8.6. For fixed k\ge0, take w(n)=n^k. The initial value is the power sum \sum_{n\le v}n^k-1. Since k is fixed, it can be evaluated by an integer-valued polynomial of fixed degree with rational coefficients, performing exact cancellation in the integers before mapping into the coefficient ring; alternatively, interpolation may be used over a suitable coefficient field. Construction 8.3 yields all values \sum_{p\le v}p^k. The expression required by P5493 in the lecture notes,

then requires only O(\sqrt N) table queries. If f(p)=\sum_{k=0}^{d}c_kp^k, compute the d+1 tables of sums of powers of primes separately and combine them linearly to obtain P_f(v)=\sum_{p\le v}f(p); for fixed d, only the constant changes. Polynomial interpolation is performed in a coefficient field or residue class ring in which all the denominators involved are invertible.

§3 Reconstructing a Multiplicative Function

Definition 8.7. Let f be multiplicative and define

with the convention p_0=1. If x\le p_j, then S_f(x,j)=0.

Theorem 8.8. For x>p_j,

In particular, \sum_{n\le N}f(n)=1+S_f(N,0).

Proof. The first term gives the contribution of primes. Every composite n has a unique least prime factor p_k>p_j, with p_k^2\le n. If n=p_k^t, where t\ge2, it occurs exactly in the second summand for e=t-1. Otherwise, it is uniquely expressible as p_k^eu, where e\ge1, u>1, and \operatorname{lpf}(u)>p_k; its contribution is f(p_k^e)f(u) by multiplicativity. In this case p_k^{e+1}<n\le x, so the first summand enumerates it exactly once. The cases are pairwise disjoint. \square

Given the prime-weight table, (8.6) defines a finite recursion with increasing least prime factors, constituting the reconstruction stage of the Min_25 method. The integer argument strictly decreases, while the lower bound on prime factors strictly increases; the unit contribution 1 is added only at the outermost level.

Remark. The total cost of the Min_25 method consists of prime-table construction and recursive reconstruction. Proposition 8.5 counts the former; §5 counts the latter through its recursive states and prime-power branches. Chapter IX gives another reconstruction procedure and its operation bound under the same polynomial condition at primes.

§4 The Inverse Process in the Completely Multiplicative Case

Theorem 8.9. Let g be completely multiplicative. Suppose the values of P_g on \mathcal Q_N and the weights at primes p\le\sqrt N are known. Set A(v)=P_g(v). Process the small primes p in decreasing order; for each prime, process quotient values satisfying v\ge p^2 in increasing order, applying

On completion, 1+A(v)=\sum_{n\le v}g(n). The number of updates is the same as in (8.3).

Proof. Fix a round p=p_j. By Proposition 8.2, restoring A_{j-1}(v) requires A_{j-1}(\lfloor v/p\rfloor). With increasing v, the smaller argument has either already been restored in this round or was unchanged because it is less than p^2. Thus (8.7) undoes the round entry by entry. Undoing all rounds in decreasing prime order finally restores A_0(v)=\sum_{2\le n\le v}g(n). \square

This process inverts the triangular transformation by reversing the loop order. Formula (8.7) uses only addition, subtraction, and multiplication, so it also applies when g(p)=0.

§5 An Independent Count of the Recursion Tree

Proposition 8.10. Let P^+(d) denote the greatest prime factor of d>1, and put

Starting at S_f(N,0), evaluate (8.6) directly by depth-first recursion without merging equal states. If local function values and prime-prefix queries have constant cost, the recursive stage uses O(1+C(N)) operations. Moreover, C(N)=o(N); a coarse explicit upper bound is O(N/\sqrt{\log N}).

Proof. A recursive path selects complete prime powers in strictly increasing prime order; denote their product by d. If the last selected power is p^e, the loop condition is exactly dp\le N, and p=P^+(d). Conversely, for an integer satisfying dP^+(d)\le N, every initial product of its prime-power factors also satisfies the corresponding loop condition. Unique factorization therefore gives a unique recursive path. Thus the nonroot calls are indexed precisely by the set in (8.8). Each call computes its prime contribution at constant cost; every successful exponent iteration produces a child call. Every prime branch entered contains at least one successful exponent iteration, so its final exponent-loop check is amortized against a child call. Each function call has only a constant number of additional outer-loop termination checks, proving the operation bound.

Fix Y\ge2. If P^+(d)>Y, then d\le N/Y; if P^+(d)\le Y, then d uses at most \pi(Y) primes, each with exponent at most \lfloor\log_2N\rfloor. Consequently,

First let N\to\infty, then let the fixed Y\to\infty, to obtain C(N)=o(N). Taking Y=\sqrt{\log N} and using \pi(Y)\le Y, the logarithm of the second term is O(\sqrt{\log N}\log\log N)=o(\log N), so that term is eventually less than N^{1/2}; the first term is N/\sqrt{\log N}. \square

This gives a sublinear bound for the recursion. The next chapter gives another complete reconstruction procedure with a fixed-power operation bound.

Chapter IX — Sparse Convolution Supported on Powerful Numbers

§1 Local Differences

Definition 9.1. A positive integer d is called powerful if every prime p dividing d satisfies p^2\mid d. Denote the set of such integers by \mathcal P, with the convention 1\in\mathcal P. These are the powerful numbers of the lecture notes.

Theorem 9.2. Let f,g be multiplicative functions satisfying f(p)=g(p) for every prime p. There is a unique multiplicative function h such that f=g*h, and h(d)\ne0 implies d\in\mathcal P.

Proof. Since g(1)=1, its convolution inverse exists and is multiplicative. Thus h=g^{-1}*f exists, is unique, and is multiplicative. Comparing the two sides at a prime gives f(p)=g(p)+h(p), so h(p)=0. If a prime factor p occurs in d with exponent exactly one, then h(d)=h(p)h(d/p)=0. \square

Proposition 9.3. Given the local values f(p^e),g(p^e), the values of h can be computed successively from

If, in particular, we choose the completely multiplicative function g(p^e)=f(p)^e, this simplifies to

Proof. Formula (9.1) follows by separating the term i=0 in the convolution identity at a prime power. For (9.2), direct substitution gives

The intermediate terms cancel in pairs, and uniqueness proves the assertion. \square

The telescoping sum uses only addition, subtraction, and multiplication, so it applies when f(p)=0, in coefficient rings with zero divisors, and for composite moduli.

§2 Counting and Enumerating the Support

Lemma 9.4. Every d\in\mathcal P has a unique expression

Consequently, |\mathcal P\cap[1,X]|=O(\sqrt X), and

Proof. For each positive prime exponent e in d, if e is even, give that prime exponent e/2 in a and 0 in b; if e is odd, then e\ge3, and the corresponding exponents are (e-3)/2 and 1. Parity at each prime also proves uniqueness. Hence

Finally, group the terms by 2^j\le d<2^{j+1}. Each group's weighted sum is O(1), and there are O(\log(2X)) groups. \square

Construction 9.5. First list the primes not exceeding \sqrt N. A recursive state (d,j,z) records a powerful number d formed from chosen primes not exceeding p_j, together with its weight z=h(d). Starting from (1,0,1), first output (d,z). For every k>j satisfying dp_k^2\le N, successively take all e\ge2 with dp_k^e\le N, and enter

Proof. Along each path, the primes strictly increase and all exponents are at least 2, so every output is powerful. Conversely, the factorization of a powerful number into increasing primes with their unique exponents determines its unique generation path. The parent state's integer is coprime to the newly added prime power, so multiplicativity ensures the correct weight. Each edge produces a new node; prime enumeration at each node adds only a constant number of final boundary checks. Thus, when local values are available in constant time, the enumeration cost is of the same order as the number of outputs. \square

In an implementation, test p\le\lfloor\sqrt{N/d}\rfloor before entering a branch, and use division to test whether another multiplication by p is possible, avoiding overflow before the test. Even when z=0, correctness does not depend on skipping the branch; if zero in the coefficient ring can be recognized in constant time, all descendants may be pruned.

§3 A Single Prefix Sum

Theorem 9.6. Under the hypotheses of Theorem 9.2, write G(x)=\sum_{n\le x}g(n). Then

If every required value of G and every local value h(p^e) is available in constant time, then, once the small primes are known, F(N) can be computed with O(\sqrt N) elementary operations.

Proof. Interchanging the finite sums in f=g*h gives F(N)=\sum_{d\le N}h(d)G(\lfloor N/d\rfloor). Remove the terms outside the support by Theorem 9.2, then apply Lemma 9.4 and Construction 9.5. \square

Remark. Choosing the auxiliary function g has two parts: obtaining sparse support through g(p)=f(p), and determining how to compute its prefix sums G. For a completely multiplicative choice, Theorem 8.9 recovers the prefix sums; when a simple convolution relation is available, the methods of Chapter VII apply.

Example 9.7. Fix an integer c\ge0, and let the multiplicative function q satisfy

Take g=\mathbf1. By (9.2), h(p^{2j+1})=0, whereas

Define the multiplicative function J_c by the local values J_c(p^j)=p^{cj}-p^{c(j-1)}, and set J_c(1)=1. Then

Proof. The function h can be nonzero only at squares, and h(a^2)=J_c(a). Substitute this fact and G(x)=x into (9.6). \square

The values of J_c can be computed by a linear sieve using J_c(p)=p^c-1 and

Thus, for fixed c and with the cost of computing p^c included, the time and space required by (9.8) are both O(\sqrt N). Only addition, subtraction, and multiplication occur, so the computation can be performed directly in \mathbf Z/2^{32}\mathbf Z. This gives the auxiliary function for the PN example in the lecture notes.

§4 Recovering All Quotient Values Simultaneously

Lemma 9.8. Fix 1\le d\le N, and put M=\lfloor N/d\rfloor. Then

For q\in\mathcal Q_M, the corresponding preimage is exactly

Proof. For v=\lfloor N/i\rfloor, we have \lfloor v/d\rfloor=\lfloor M/i\rfloor; the nonzero values correspond to i\le M. This proves the equality of image sets. The description of the preimage follows directly from integer division. \square

Construction 9.9. Suppose \mathcal Q_N has been indexed in increasing order as v_1<\cdots<v_L, and G(v) is known. Initialize a zero array \Delta of length L+1. Enumerate each (d,h(d)) and put M=\lfloor N/d\rfloor; then enumerate q\in\mathcal Q_M and set

Perform \Delta_a\leftarrow\Delta_a+c and \Delta_{b+1}\leftarrow\Delta_{b+1}-c. Finally, take one ordinary prefix sum, setting

Formula (7.4) computes ranks in constant time; (7.2) enumerates each \mathcal Q_M without repetition. Proposition 7.2 guarantees q\in\mathcal Q_N, so every G(q) is available for lookup.

Theorem 9.10. Construction 9.9 correctly recovers the entire quotient-value representation of F. When local weights and queries have constant cost, its time, excluding the construction of G, is O(\sqrt N\log(2N)), and its space is O(\sqrt N).

Proof. For fixed d, (9.10) partitions all v\ge d into disjoint intervals of consecutive indices according to equal quotients. The difference-array updates add exactly h(d)G(\lfloor v/d\rfloor) to every such v. Summing over all d gives correctness by Theorem 9.6. The number of interval updates is

The final prefix accumulation and the enumeration of powerful numbers both take O(\sqrt N). Only the quotient-value table, the difference array, the small-prime table, and the recursion stack need to be retained. \square

Applying Theorem 9.6 independently to each quotient value yields only the general bound O(N^{3/4}) from (7.6). The improvement in (9.13) comes from accumulating identical contributions over intervals.

§5 A Complete Reconstruction Procedure

Theorem 9.11. Let f be multiplicative, and suppose f(p)=\sum_{k=0}^{D}c_kp^k at primes. Assume that D is fixed, the required local prime-power values are available in constant time, and sums of fixed powers can be evaluated in constant time. Then the quotient-value representation of f can be computed in

elementary operations, using O((D+1)\sqrt N) space.

Proof. Perform the following finite constructions in order.

1. For k=0,\ldots,D, use Construction 8.3 to compute the prime power sums, then combine them linearly into P_f(v).
2. Define the completely multiplicative function g by g(p)=f(p); thus P_g=P_f. Obtain G(v) by Theorem 8.9.
3. Define h by (9.2), and enumerate its powerful-number support by Construction 9.5.
4. Obtain all F(v) by Construction 9.9.

The preceding propositions establish the correctness of each stage. Adding the costs in Proposition 8.5 and Theorem 9.10 gives the result. \square

Thus the polynomial condition at primes is used only to construct the prime-weight table; deviations at higher prime powers are carried separately by h. This explains why the prime-weight sieve, its inverse process, and the PN correction can be combined.

Chapter X — Divisibility Matrices and Sums in Two Variables

§1 The Matrix of the Divisibility Relation

Definition 10.1. Fix n\ge1, and let E be the n\times n integer matrix

It is lower triangular with every diagonal entry equal to 1. For an arithmetic function a, put b=\mathbf1*a.

Proposition 10.2. We have

In particular, its determinant is \prod_{d=1}^{n}a(d).

Proof. The (i,j) entry on the right is \sum_{d\mid i,\ d\mid j}a(d)=b(\gcd(i,j)). Since \det E=1, taking determinants of the matrix product proves the remaining assertion. \square

Proposition 10.3. The matrix E is invertible over the integers, with

Proof. The (i,j) entry of the product of the two matrices is [j\mid i]\sum_{d\mid i/j}\mu(d)=[i=j]. \square

Elimination on the divisibility matrix is therefore a finite linear representation of Möbius inversion. Formula (10.3) directly gives an integer elimination transformation that can be applied entrywise in the coefficient ring.

§2 Matrices with Integer Weights

Example 10.4. Let A_{ij}=ij\gcd(i,j). Define the lower triangular integer matrix

Then

Proof. The (i,j) entry of the product is

Every diagonal entry of L equals 1, so the determinant identity follows. \square

This factorization holds entrywise over the integers. Left multiplication by L^{-1} gives the upper triangular matrix

which is exactly the elimination result used in P6384 in the lecture notes. To implement it solely by elementary row operations, process i in increasing order and, for k=2,\ldots,\lfloor n/i\rfloor, subtract k times row i from row ki. The multiplier k corresponds to the original matrix's row weight i. The notes' method of first dividing out row and column weights and then eliminating is also valid over the rationals; (10.4) is itself an integer identity and can be mapped directly into a residue class ring.

§3 The Multiplicative Relation between GCD and LCM

Theorem 10.5. Suppose f(1)=1. The function f is multiplicative if and only if, for all positive integers a,b,

Proof. If f is multiplicative, factor a,b into primes. For each prime p, the local factors on the two sides are f(p^u)f(p^v) and f(p^{\min(u,v)})f(p^{\max(u,v)}), respectively, and are equal. Conversely, for coprime a,b, the right-hand side is f(1)f(ab)=f(ab), proving multiplicativity. \square

This assertion requires neither complete multiplicativity nor nonzero function values.

Example 10.6. Let \omega(n) denote the number of distinct prime factors of n, with \omega(1)=0. For t\in\{0,1\}, put

Then the double sum in P12487 in the lecture notes is

Proof. The function f_t is multiplicative. By Theorem 10.5, rewrite every summand as f_t(i)f_t(j) and separate the two finite sums. \square

§4 Eliminating Squarefree Factors

Lemma 10.7. For every positive integer n,

Proof. The right-hand side of the first identity counts squarefree divisors, with each distinct prime factor independently included or omitted. Both sides of the second identity are multiplicative. At p^e, the right-hand side is 1 for e=0,1, and 1+\mu(p)=0 for e\ge2, agreeing with the left-hand side. \square

Definition 10.8. Write

In the cases needed here, S_0(x)=x and S_1(x)=x(x+1)/2.

Proposition 10.9. We have

Proof. Substitute the two identities in Lemma 10.7 successively and put d=k^2a, obtaining

The inner sum, enumerated by a, is exactly D_t(\lfloor N/k^2\rfloor). All sums are finite, so interchanging their order requires no convergence assumptions. \square

Construction 10.10. Put r=\lfloor\sqrt x\rfloor. Then

For t=0,1, this gives an exact algorithm using O(\sqrt x) elementary operations.

Proof. Apply Proposition 7.4 symmetrically, taking both functions to be a\mapsto a^t. \square

First compute \mu(k) by a linear sieve, then apply (10.11) term by term. Formula (10.10) can thereby be evaluated with O(\sqrt N\log(2N)) elementary operations, because

This gives a complete baseline implementation against which to check more advanced algorithms. The largest parameters in the lecture notes call for an improvement to the inner lattice-point sum; the next chapter develops this geometric step.

Chapter XI — The Lattice Boundary of a Hyperbola and Weighted Summation

§1 From Row-by-Row Enumeration to Rational Directions

This chapter supplies the geometric algorithm omitted from the final example in the original lecture notes. Fix N\ge1, and write

The variables x,y in the sum are positive integers. The hyperbola method in the preceding chapter computes \lfloor N/y\rfloor one row at a time. The alternative representation in this chapter combines a consecutive group of rows into a rectangle and a triangle with a rationally sloped hypotenuse, and then evaluates their weights together.

We use the subdivision by rational directions developed in the original papers, together with a polygonal chain of lattice points lying strictly above the hyperbola. All update formulas and an executable program are given here. The sources of the method are Section 4 of Zhu Zhenting's Some Special Problems in Summing Arithmetic Functions and Richard Sladkey's original paper.

Definition 11.1. A point (X,Y)\in\mathbf Z_{>0}^2 is an exterior point if XY>N. A lattice point with positive coordinates satisfying XY\le N is an interior point; lattice points on the hyperbola are included among the interior points. A primitive direction pointing downward and to the right is written as

The nonnegative pair (p,q) is used for storage; the second component is negated only when applying the direction to a point.

Definition 11.2. Two directions L=(p,q) and R=(r,s) are adjacent directions if

They satisfy q/p>s/r. Their mediant direction is defined as M=L+R=(p+r,q+s).

Proposition 11.3. Equation (11.3) implies that M is primitive and that both (L,M) and (M,R) are pairs of adjacent directions. If the absolute slope of a primitive direction (a,b) lies strictly between s/r and q/p, there are unique positive integers u,v such that

In particular, a\ge p+r and b\ge q+s.

Proof. Both new determinants equal the original determinant. If an integer divides both p+r and q+s, it divides q(p+r)-p(q+s)=1. Since their determinant is -1, the two directions form a basis of \mathbf Z^2; solving the linear equations gives integer coordinates. The strict inequalities between the slopes make both coordinates positive. \square

Repeated insertion of mediant directions gives the Stern–Brocot search among rational numbers, with comparisons performed on integer directions.

§2 Four Moments of a Triangle

Definition 11.4. For a primitive direction (p,q), define the lattice-point set in a closed triangle by

When (p,q)=(1,0), the set degenerates to the two points \{(0,0),(1,0)\}. Define four moments:

The full record for a direction is (p,q,C,U,V,W). The initial records are

Lemma 11.5. Suppose that L=(p,q) and R=(r,s) satisfy (11.3). The lattice-point set \mathcal T(p+r,q+s) can be assembled from three parts: \mathcal T(p,q); the translate of \mathcal T(r,s) by (p,-q); and the rectangle

The first two parts overlap exactly once, at (p,-q); the rectangle is disjoint from both.

Proof. The triangle between the two smaller hypotenuses and the larger hypotenuse has area 1/2. Since the determinant of the two directions is \pm1, a linear transformation invertible over the integers maps this small triangle to the triangle with vertices (0,0),(1,0),(1,1). The latter has no lattice points other than its three vertices. Replacing the two smaller hypotenuses by the larger one therefore adds no lattice points. The remaining parts follow by comparing their boundaries row by row or column by column. \square

Proposition 11.6. Denote the moments in the left and right records by (C_L,U_L,V_L,W_L) and (C_R,U_R,V_R,W_R), respectively, and put

The moments of the mediant direction are

Proof. On the translated right triangle, the four functions being summed become, respectively,

The four moments of the rectangle are ps,sA,pB,AB. Finally, for the repeated point (p,-q), subtract its four contributions 1,p,-q,-pq. \square

The formulas also apply to R=(1,0), because Definition 11.4 assigns a two-point set to the degenerate triangle.

Remark. Column-by-column summation explains why the weighted problem needs an additional quantity. If we also define

then \gcd(p,q)=1 gives

Indeed, for fixed 1\le u<p, the vertical coordinate ranges from -q to -\lfloor qu/p\rfloor-1, whose sum is

The same boundary gives both D_0 and D_1, with each edge's contribution computed from a constant-dimensional moment record.

§4 Boundary Search, the Stack, and the Stopping Criterion

Construction 11.8. Start at

This is an exterior point: writing N=km+r, where k=\lfloor N/m\rfloor\ge m and 0\le r<m, gives k(m+1)-N=k-r>0.

Initially the stack contains the directions (1,0),(1,1) from bottom to top, each with the moments in (11.7). The current accumulators are (Z_0,Z_1)=(0,0). Repeat the following procedure.

1. Pop the top direction L. As long as P+L remains an exterior point, add (11.14) and (11.15) to the accumulators and set P\leftarrow P+L. Repeat until the product of the next candidate's coordinates is at most N.
2. If the current value satisfies Y^3\le N, end the geometric procedure.
3. Let R be the current top of the stack. As long as the product of the coordinates of P+R is at most N, set L\leftarrow R, pop it, and take the new top as R.
4. Set M=L+R and compute its record using (11.10). If P+M is an exterior point, set R\leftarrow M, push it onto the stack, and repeat this step.
5. If the product of the coordinates of P+M is at most N, test

If the condition holds, end this rational-direction search and return to step 1; otherwise, set L\leftarrow M and return to step 4.

Here P+(p,q) always means (X+p,Y-q). All geometric tests use integer arithmetic. Search candidates may have Y-q\le0; such a candidate also enters the branch where the coordinate product is at most N, but it is not called an interior point under Definition 11.1. Actual advances accept only candidates whose coordinate product is greater than N. Since X+p>0, every endpoint of an edge actually emitted has a positive vertical coordinate.

Lemma 11.9. Equation (11.17) is a justified exclusion criterion. If the product of the coordinates of P+M is at most N and (11.17) holds, then starting from P+M and advancing along R, or along any direction downward and to the right that is steeper than R, cannot reach an exterior point with positive coordinates. The starting point may have a nonpositive vertical coordinate.

Proof. Put x_0=X+p_M and f(x)=N/x. Equation (11.17) is equivalent to

The function f is convex, so for x\ge x_0 it lies above its tangent at x_0. The starting point is no higher than f(x_0), and the slope of the ray under consideration is at most the tangent slope. The entire ray therefore lies no higher than the curve. \square

Proposition 11.10. Every edge emitted by Construction 11.8 satisfies the hypotheses of Proposition 11.7, and the processed row intervals are pairwise disjoint. At termination, the processed rows are exactly Y,Y+1,\ldots,m.

Proof. First establish the convexity fact repeatedly used in the search. If P=(X,Y) is an exterior point and Q=(X+a,Y-b) satisfies a>0 with coordinate product at most N, then

The first inequality uses Y>N/X and Y-b\le N/(X+a); the second uses convexity. Consequently, starting from Q and advancing along (a,-b) or any steeper direction downward and to the right always stays no higher than the hyperbola. This conclusion also permits a nonpositive vertical coordinate for Q.

The search invariant is as follows: L,R are adjacent, the coordinate product of P+L is at most N, and P+R is an exterior point; viewed from P, all lattice points with positive coordinates in directions no shallower than L are interior points. If the mediant candidate P+M has coordinate product at most N, any integer direction between L and M can be written as uL+vM, where u,v\ge1. Starting from P+M, the remaining direction uL+(v-1)M is no shallower than M, so the convexity fact above shows that this point also lies below the curve. The same holds for subsequent points on the ray in direction M. Thus, when setting L\leftarrow M, the entire newly enlarged region of steep directions has been excluded. The same argument applies when step 3 pops a right-hand candidate whose coordinate product is at most N. If P+M is an exterior point, only R is replaced by M, while the old, shallower direction is retained in the stack. The invariant is therefore preserved.

When (11.17) holds, any integer direction lying strictly below the ray in direction R can, by the integer-basis property, be written as \alpha L+\beta R with \alpha\ge1. If \beta\le0, it belongs to the already excluded region of steep directions. If \beta\ge1, it equals

By Lemma 11.9, starting from P+M and following the remaining direction, which is no shallower than R, cannot encounter an exterior point. Thus there are no exterior lattice points strictly below the entire rightward ray in direction R. Since R is primitive, one unit segment in this direction has no other lattice points on its edge; all points of the translated triangle except its two endpoints are therefore interior points.

Initialization and repeated advances must also be checked. Initially P=(k,m+1), where k=\lfloor N/m\rfloor\ge m. Integer points strictly below the initial ray of slope -1 have the form (k+u,y), where u\ge0 and y\le m-u. Whenever both coordinates are positive,

The initial edge therefore has the required property. During repeated advances along a direction already established, lattice points farther to the right and strictly below the line still belong to the previously excluded set. Once the next point on the ray has coordinate product at most N, the convexity fact above also excludes all subsequent points on the ray, restoring the full invariant for the next search.

For every edge actually emitted, Proposition 11.7 proves that the triangle and the rectangle to its left count exactly all the relevant interior points; the exterior right endpoint rules out omissions to the right. Each edge changes the row index from Y to Y-q and counts only [Y-q,Y-1]. Induction from the initial row index m+1 gives the disjointness and coverage of the row intervals. \square

Stack boundary conditions. The horizontal sentinel (1,0) at the bottom of the stack is never popped by step 3: the current point satisfies XY>N, so (X+1)Y>N. When R is the sentinel, (11.17) becomes 0\ge N, which is impossible. A search can therefore end only after at least one nonhorizontal direction has been pushed. On returning to step 1, the top to be popped is always above the sentinel. Thus neither access to the top of the stack can encounter an empty stack.

Termination. If R is the horizontal sentinel, adjacency forces q_L=1. Before the cutoff condition is met, Y\ge2. Repeatedly setting L\leftarrow L+R keeps the vertical difference equal to 1 while increasing the horizontal difference, and eventually produces a candidate whose coordinate product is greater than N, pushing a nonhorizontal direction. Thereafter, each exterior candidate sets R\leftarrow L+R and strictly increases its vertical difference; the exterior-point condition requires this difference to be less than Y, so only finitely many such updates can occur. With a fixed nonhorizontal R, consecutive candidates on the non-exterior branch strictly increase the horizontal difference and eventually satisfy (11.17). Each direction search is therefore finite. At the end of a search, the direction on top of the stack leads to an exterior point, so the next iteration successfully advances at least once, decreasing the positive integer Y by at least one. The entire geometric procedure is therefore finite as well.

Remark. The stack retains adjacent directions obtained in earlier searches for reuse in the next search; (11.17) determines where each search stops.

Theorem 11.11. After the geometric procedure ends, add the remaining rows:

Then

Proof. Together, the geometric procedure and the remaining rows count exactly all interior points with y\le m. Exchanging the horizontal and vertical coordinates gives the part with x\le m. These two parts cover xy\le N and intersect in [1,m]^2. Summing the respective weights over this square gives the subtracted terms. \square

§5 An Executable Integer Construction

The following Python program implements the procedure above using arbitrary-precision integers. The input is N\ge1, and the output is the pair of exact integers (D_0(N),D_1(N)). The next section analyzes its complexity by counting elementary arithmetic operations.

from math import isqrt

def merge(L, R):
    p, q, C, U, V, W = L
    r, s, D, E, F, G = R
    assert q * r - p * s == 1
    A = p * (p - 1) // 2
    B = -s * (2 * q + s + 1) // 2
    return (p + r, q + s,
            C + D + p * s - 1,
            U + E + p * D + s * A - p,
            V + F - q * D + p * B + q,
            W + G + p * F - q * E - p * q * D + A * B + p * q)

def edge_weight(X, Y, R):
    p, q, C, U, V, W = R
    rectangle = X * (X - 1) // 2 * (q * (2 * Y - q - 1) // 2)
    return ((X - 1) * q + C - 2,
            rectangle + W + X * V + Y * U + X * Y * C
            - X * Y - (X + p) * (Y - q))

def hyperbola_moments(N):
    m = isqrt(N)
    X, Y = N // m, m + 1
    stack = [(1, 0, 2, 1, 0, 0), (1, 1, 3, 1, -2, -1)]
    Z0 = Z1 = 0

    def outside(R):
        return (X + R[0]) * (Y - R[1]) > N

    while True:
        L = stack.pop()
        while outside(L):
            E0, E1 = edge_weight(X, Y, L)
            Z0 += E0
            Z1 += E1
            X += L[0]
            Y -= L[1]
        if Y ** 3 <= N:
            break

        R = stack[-1]
        while not outside(R):
            L = R
            stack.pop()
            R = stack[-1]

        while True:
            M = merge(L, R)
            if outside(M):
                R = M
                stack.append(R)
            elif (X + M[0]) ** 2 * R[1] >= N * R[0]:
                break
            else:
                L = M

    for y in range(1, Y):
        q = N // y
        Z0 += q
        Z1 += y * q * (q + 1) // 2

    return 2 * Z0 - m * m, 2 * Z1 - (m * (m + 1) // 2) ** 2

Example 11.12. For N=10, the program returns (27,170). Direct classification by the product also gives

Computing both moments simultaneously increases each direction record to a constant size and does not change the geometric branches.

§6 Operation Counts and Region Counts

Let B=\lfloor N^{1/3}\rfloor. The construction switches to direct enumeration when Y\le B, so processing the final rows requires O(B) operations. The geometric part consists of edge advances, direction searches, and stack operations.

Let E(N) be the number of advances in step 1 and J(N) the number of mediant directions formed in step 4. The program's total operation count is

Each advance and merge uses only a constant number of operations. Every push after initialization corresponds to a mediant-direction merge, and each record is popped at most once, so all stack operations are also controlled by J(N)+2. Each advance along the same direction is counted separately in E(N). Since every step decreases Y by at least one, E(N)\le m.

We next estimate the number of regions produced by subdivision in rational directions in tangent-region algorithms.

In a subdivision into regions by rational directions, nontrivial processing can be assigned to adjacent fractions a/b>c/d>0, where ad-bc=1. The two corresponding tangency points have horizontal coordinates \sqrt{Nb/a} and \sqrt{Nd/c}. A region with nonzero width in both integer basis directions must satisfy

This is because the horizontal span of the two unit basis directions is b+d.

Lemma 11.13. The number of positive-integer quadruples satisfying (11.21) is

Proof. Put z=bc. Since ad=z+1,

Hence 4abc^2(b+d)^2\le N. After substituting c=z/b and d=(z+1)/a, the left-hand side is at least 8z^3, so z=O(N^{1/3}). For a given z, choosing b\mid z and a\mid z+1 uniquely determines c,d. The required number is therefore at most

Let d_4(n) be the number of ways to express n as an ordered product of four positive integers. Comparing prime-power values using

Lemma 11.13 bounds the number of regions associated with adjacent fractions. The program's time and space are described by (11.20) and O(J(N)+1), respectively. Maintaining four moments increases the cost of the same geometric procedure by only a constant factor.

Source. Printed pages 108–109 of Zhu Zhenting's original paper give the region count involving \log^3 N; page 110 discusses a geometric comparison with the polygonal-chain method.

When the answer is required modulo an integer, additions and multiplications involving C,U,V,W,Z_0,Z_1 may be reduced modulo it. Coordinates, order comparisons, divisibility tests, and (11.17) must always be evaluated over the integers. Every factor of one-half in the formulas comes from an integer summation formula and should either be canceled over the integers first or interpreted using the inverse of 2 in a residue ring where it is invertible.

§7 Returning to the Final Möbius Sum

The preceding chapter established

This chapter makes D_0,D_1 explicit executable subroutines. An implementation should cache results for repeated arguments \lfloor N/k^2\rfloor. For small arguments, prefix sums of n^t\tau(n) can be precomputed directly, avoiding geometric searches.

Calling the program separately for each term gives total cost

where T is given by (11.20), and the initial O(\sqrt N) includes the Möbius sieve. Using the ordinary hyperbola subroutine from the preceding chapter gives total bound O(\sqrt N\log(N+1)).

Corollary 11.14. Suppose that an exact hyperbola subroutine uses O(x^{1/3}\log^3(x+1)) arithmetic operations per call. Then term-by-term evaluation of (11.24) has total cost O(\sqrt N), including the time for a linear sieve of the Möbius function.

Proof. Group the indices k by \sqrt N/2^{j+1}<k\le\sqrt N/2^j. Such a group has at most O(\sqrt N/2^j) elements, and every argument is at most 2^{2j+2}. The assumed subroutine bound within this group is O(2^{2j/3}(j+1)^3), so the total cost over all groups is at most

The linear sieve for the Möbius function satisfies the same time bound. \square

Editorial Notes

This text reconstructs Su Ziran (minstdfx), Cirno's Perfect Number Theory Classroom: Summer-Term Course on Algorithmic Competitions, file version 20260619-2 ver., comprising 128 pages. Its chapter organization and proof style follow the Bourbaki-style reconstruction of the introductory data-structures supplement, 数据结构入门补充材料_布尔巴基式重构.md. Informal explanations, diagrams, and code from the original notes have been recast as explicit input conditions, states, transformations, and correctness propositions, reordered by mathematical dependence. To develop algorithms only sketched in the original, the text adds local recurrences, the inverse process of the prime sieve, interval accumulation for sparse convolution, and a geometric construction for lattice-point summation.

Correspondence with the Original

All page numbers below refer to physical pages of the PDF.

Corrections and Additions

1. The greatest common divisor is defined as the nonnegative generator of an integer ideal, including the all-zero family. The least common multiple, in the range used here, is the least positive common multiple. Primes, composite numbers, and units are distinguished; 0,\pm1 are not classified as composite. Divisor pairing is restricted to positive integers.
2. Cancelling a common factor in a congruence also changes the modulus. ExCRT merging maintains the complete solution class. GCD Table checks candidates against the GCD equalities, treating a bounded starting position and an entire segment contained in the table separately.
3. Euler exponent reduction without coprimality is proved componentwise for the full prime-power factors, including the condition that the exponent has reached the threshold. Zero exponents, modulus one, and noninvertible denominators are treated separately. The ex-Lucas construction distinguishes a factorial with all factors of p removed from a product formed by omitting the multiples of p.
4. The exBSGS construction retains the coefficient accumulated by multiplying by a/g in each round and first tests the current smallest exponent. A single invariant proves the stripping step, the change of modulus, minimality of the solution, and termination.
5. The modulus of the exponent ring for a discrete logarithm is the group order; in CF1106F it is p-1, and the linear congruence is solved for this general modulus. Primality testing is presented with one-sided correctness, the random error theorem, and finite-range base sets. Pollard–Rho returns a factor when 1<d<n; failed batches are replayed or the parameters are changed.
6. The constant-one function on pages 65 and 67 is consistently denoted by \mathbf1, and the identity function by \operatorname{id}. The local factors in the divisor-sum formula are corrected to geometric sums; \prod_p(\nu_p(n)+1) counts divisors. The domain counted by \varphi is specified as a finite complete residue system, and \mu(1)=1 is stated explicitly.
7. Inversion applies to general arithmetic functions; the corresponding weighted identities use complete multiplicativity. Sums over multiples are finitely truncated. Generating functions are constructed as formal series, and convergence of real series is proved where they are used.
8. With the original definition of the double prefix sum on page 74, the correct identity is F_k(x)=T_k(2x)-2T_k(x). The text proves it by counting the coefficient corresponding to each possible sum of an ordered pair, restoring the missing factor 2. Preprocessing modulo 2^{32} uses no division.
9. Prime detection on an interval and complete factorization on an interval use different starting points. In the Haitang example, the candidate conditions are S(m)=s, m>s, and m\mid n-s. The entire short interval must be factored, and the checking cost depends on the total number of divisors of all its elements.
10. The linear-sieve branch on page 85 is corrected to ask whether the newly multiplied prime p divides the old index i, with the result stored at ip. Computations over general coefficient rings use the coprime cofactor and local data; their complexity includes the cost of obtaining prime-power data.
11. Given two quotient-value tables, direct hyperbola evaluation of a general convolution uses O(N^{3/4}) operations; when the smaller prefix values admit linear preprocessing, truncation gives O(N^{2/3}). The Du Jiao sieve uses a convolution factor with invertible first value and quickly evaluable prefixes.
12. The prime-weight sieve updates in increasing prime order and decreasing quotient-value order. Recovering completely multiplicative prefixes reverses both orders through the inverse triangular transformation. The cost of Min_25 recursion is counted by its nodes.
13. The auxiliary function actually recovered in the PN example is h; the end of page 114 incorrectly labels it g. For the completely multiplicative choice, the local difference formula is proved by a telescoping identity. Accumulating equal contributions over intervals of ranks recovers all quotient values in O(\sqrt N\log N) time.
14. The weighted GCD matrix is factored through a unit lower triangular integer matrix, including for composite moduli. The row multipliers used in elimination, the diagonal, and the determinant all follow from the same factorization.
15. The last problem's double sum is first reduced to the square of one prefix sum, then subjected to squarefree inversion. Chapter XI develops the inner lattice-point algorithm and calculates the total cost of the outer Möbius sum.

Problem Lists from the Original Notes

Page 125 of the original PDF provides the following lists for further practice:

• gza's Recreational Number Theory 1
• gza's Recreational Number Theory 2
• Leasier's Number Theory Problem List
• CF2081G2.

约定

本文以通常的集合论、整数与有理数的构造为基础。有限集合、映射、等价关系、有限序列、归纳法及有限递归定理视为已建立。关于实数的对数、积分及初等渐近估计，仅在复杂度证明中使用。记

当 l>r 时，[l,r]=\varnothing。空和为 0，空积为 1。[P] 表示命题 P 的示性值：成立时为 1，否则为 0。未特别限定的约数均为正约数，素数均为正素数；模数总为正整数，余数总取最小非负代表。p^e\parallel n 表示 p^e\mid n 且 p^{e+1}\nmid n。整数平方根及其他根式边界均取其明确标出的整数取整值。未标底数的 \log 为自然对数；在渐近操作数中换用任意固定底数只改变常数。

环均为含幺交换环；群、单位与剩余类的具体用途在正文中定义。允许模数为 1 的退化情形。一般系数环记为 R，整数通过自然的含幺同态映入 R。环中的除法使用单位的逆元。全文区分普通乘法、逐点乘法与 Dirichlet 卷积；后者统一记为 *。形式级数的等式均指逐系数相等，除非另外给出收敛条件。

计算约定。 时间界计基本算术操作次数，存储界计整数或环元素个数。若整数、地址、模数及中间结果可用常数个机器字精确表示，则相应的加减、乘除、比较、取整和随机存取按基本操作计；任意精度整数的位复杂度另计。环元素的运算与存储代价在各构造中说明。快速幂、求逆、模数分解、局部函数求值和预处理的成本计入各算法的复杂度。

散列查找按期望代价计；随机算法的误差概率与轨道模型在相应章节说明。

各章的定义、引理、命题、定理、推论、构造与例共用“章号.序号”；公式另行编号。章节按依赖组织：先建立整除与剩余类的结构，再建立算术函数的卷积，最后讨论有限表示、筛法与格点求和。算法名称保留在相应构造处，以便与原讲义对照；原讲义的页码对应、校订及补充来源列于末尾。

第一章　整除与整数理想

§1　整除与带余除法

定义 1.1。 对 a,b\in\mathbf Z，若存在 q\in\mathbf Z 使 b=aq，则称 a 整除 b，记为 a\mid b。对 n>0，记

由定义，任意整数整除 0，而 0\mid b 当且仅当 b=0。对 n>0，因子配对 d\mapsto n/d 是 D(n) 上的对合。

命题 1.2。 对整数 a,b,c，有如下性质：整除具有传递性；若 a\mid b 且 a\mid c，则 a\mid ub+vc 对任意 u,v\in\mathbf Z 成立；若 a\mid b 且 b\mid a，则 |a|=|b|。对 t\ne0，

证明。 前两项由代入整除的定义得到。若 a,b 均非零，由 b=aq、a=br 得 qr=1，故 q,r\in\{1,-1\}；有一个为零时二者均为零。式 (1.2) 的反向使用整数乘法的消去律。\square

定理 1.3（带余除法）。 设 m>0，a,d\in\mathbf Z。存在唯一一对整数 (q,r) 满足

证明。 取 q=\lfloor(a-d)/m\rfloor、r=a-qm 即得存在性。若另有 (q',r')，则 m(q-q')=r'-r，而 |r'-r|<m；故两边均为零。\square

当 d=0 时，记所得 r 为 a\bmod m。任意连续 m 个整数的余数互异，并恰好遍历 0,\ldots,m-1：其中两个数之差的绝对值小于 m，故同余只能在二者相等时成立。

§2　整数理想与最大公约数

定义 1.4。 整数环的理想为满足下列条件的非空子集 I\subseteq\mathbf Z：若 x,y\in I，则 x-y\in I；若 x\in I、t\in\mathbf Z，则 tx\in I。对有限整数族 a_1,\ldots,a_r，置

它是包含各 a_i 的最小理想。空族生成的理想约定为 \{0\}。

定理 1.5。 每个整数理想唯一地写成 d\mathbf Z，其中 d\in\mathbf N。

证明。 若 I=\{0\}，取 d=0。否则 I 含有正整数，取其中最小者 d。对任意 a\in I，带余除法给出 a=qd+r、0\le r<d。因 r\in I，最小性迫使 r=0，故 I\subseteq d\mathbf Z；反向包含由理想性质得到。若 d\mathbf Z=e\mathbf Z，则二者互相整除；由 d,e\ge0 得 d=e。\square

定义 1.6。 有限整数族的最大公约数为其生成理想的非负生成元，记为 \gcd(a_1,\ldots,a_r)。两个整数互素，意指其最大公约数为 1。特别地，\gcd(0,0)=0。

定理 1.7（Bézout）。 设 d=\gcd(a_1,\ldots,a_r)。则 d\mid a_i 对所有 i 成立，且存在整数 u_i 满足

任意公约数均整除 d。因此，只要该族不全为零，d 就是通常大小次序中的最大正公约数；方程 \sum_i a_i x_i=c 有整数解，当且仅当 d\mid c。

证明。 等式 (a_1,\ldots,a_r)=d\mathbf Z 同时给出前两项。公约数整除每个整数线性组合，故整除 d。最后一项等价于 c\in d\mathbf Z，包括 d=0 的情形。\square

命题 1.8。 对整数 a,b,q，有

若 \gcd(a,b)=1 且 b\mid ac，则 b\mid c。若另有 a\mid c，则 ab\mid c。

证明。 第一式由 (a,b)=(b,a-qb) 得到，第二式由生成理想整体乘以 t 得到。取 ua+vb=1，乘以 c 后可知 b\mid c。若 c=at，则 b\mid at 蕴含 b\mid t，故 ab\mid c。\square

构造 1.9（扩展 Euclid 算法）。 先对 |a|,|b| 连续实施带余除法，直到余数为零。最后一个非零余数为 \gcd(a,b)；二者初始均为零时直接返回 0。若递归已得到

则返回 x=y'、y=x'-qy'，得到 ax+by=d。最后按原输入的符号调整系数。

证明。 式 (1.6) 保持最大公约数，严格递减的非负余数保证终止。回代恒等式给出所述系数。设相邻正余数为 u>v；若 v\le u/2 已经减半，否则下一余数为 u-v<u/2。故每两步至少减半，除法次数为 O(1+\log(\min(|a|,|b|)+1))。\square

注。 这里计整数除法次数；系数大小、整数乘法及除法的位代价不包含在此界中。对多个整数逐次求最大公约数，所得结果与分组次序无关，因为其生成理想相同。

§3　素数与唯一分解

定义 1.10。 素数为恰有两个正约数的正整数。大于 1 而非素数的整数称为合数。对素数 p 和正整数 n，记

引理 1.11。 若素数 p 整除 ab，则 p\mid a 或 p\mid b。

证明。 若 p\nmid a，则 \gcd(p,a) 是 p 的正约数且不等于 p，故等于 1。应用命题 1.8。\square

定理 1.12（唯一分解）。 每个 n>0 唯一地写成

其中仅有限个指数非零，且 n=1 对应空积。

证明。 对 n 归纳。若 n>1 非素数，取 n=uv、1<u,v<n，分别分解 u,v，得到存在性。若有两种分解，左边任一素因子由引理 1.11 整除右边某一素因子，故二者相等；约去这一因子并归纳，得到唯一性。\square

定义 1.13。 对非空有限正整数族，最小公倍数为其公共正倍数中的最小者，记为 \operatorname{lcm}(a_1,\ldots,a_r)。

命题 1.14。 对正整数 a,b，

相应的最小值、最大值公式适用于任意非空有限正整数族。

证明。 整除的条件直接来自唯一分解。公约数的每个指数至多为逐项最小值，且该最小值确实给出公约数；公倍数同理。最后使用 \min(s,t)+\max(s,t)=s+t。\square

推论 1.15。 若正整数 a_1,\ldots,a_r 两两互素，且其乘积为某正整数的 k 次幂，其中 k>0，则每个 a_i 都是正整数的 k 次幂。

证明。 任一素数至多出现在一个 a_i 中，其在总乘积中的指数为 k 的倍数，故在该 a_i 中的指数也是 k 的倍数。\square

§4　同余与线性方程

定义 1.16。 对 m>0，称 a,b 模 m 同余，若 m\mid a-b，记为 a\equiv b\pmod m。整数 a 的等价类记为 [a]_m。这些类组成的集合记为 \mathbf Z/m\mathbf Z，其加法、乘法定义为

命题 1.17。 同余是等价关系，式 (1.11) 与代表元的选择无关，并使 \mathbf Z/m\mathbf Z 成为含幺交换环。若正整数 d 满足 d\mid a,b,m，则

对任意整数 c，令 g=\gcd(c,m)，则

证明。 等价关系及环运算由整除对加法、乘法的稳定性得到。式 (1.12) 是整数等式 a-b=tm 的约分。写 c=gc'、m=gm'，则 \gcd(c',m')=1；命题 1.8 给出 m'\mid c'(a-b) 当且仅当 m'\mid a-b。\square

注。 同余式中约去非单位因子通常必须同时改变模数。模 1 的环只有一个元素，[0]_1=[1]_1；所有整数模 1 同余。

定理 1.18（线性同余）。 设 m>0，d=\gcd(a,m)。方程

有解当且仅当 d\mid b。若 ua+vm=d 且 d\mid b，则全部整数解为

因此恰有 d 个模 m 的解类。

证明。 原方程等价于 ax+my=b 有整数解，故由定理 1.7 得到可解条件与一个特解。两个解之差满足 a(x-x_0)\equiv0\pmod m，再用式 (1.13)。长度为 m 的完整剩余系中，间隔 m/d 的类恰有 d 个。\square

§5　同余条件的合并

定理 1.19（中国剩余定理）。 设正整数 m_1,\ldots,m_r 两两互素，M=\prod_i m_i。映射

为环同构。特别地，任意给定 a_i 的方程组 x\equiv a_i\pmod{m_i} 恰有一个模 M 的解类。

证明。 置 M_i=M/m_i，由 Bézout 选取 u_i 使 u_iM_i\equiv1\pmod{m_i}。则

满足所有方程：在第 i 个分量中，第 i 项给出 a_i，其余项均为零。若两个整数满足同一组条件，每个 m_i 均整除其差；由两两互素性，M 整除其差。环运算显然保持，故得同构。模数 1 的分量不施加条件，可直接略去。\square

定理 1.20（两个任意模数的合并）。 设 M,N>0，d=\gcd(M,N)。方程组

有解当且仅当 d\mid B-A。可解时，令 uM+vN=d，取

则全部解为 x\equiv A+Mt_0\pmod L。

证明。 写 x=A+Mt，问题化为 Mt\equiv B-A\pmod N。定理 1.18 给出可解条件及 t 模 N/d 的唯一类。改变 t 一个周期使 x 改变 L，且 L=\operatorname{lcm}(M,N)。\square

构造 1.21（ExCRT）。 以 (A,M)=(0,1) 开始，逐一把当前类 x\equiv A\pmod M 与下一个同余条件按定理 1.20 合并。每次将 A 规范到 [0,M-1]。若某步不兼容，方程组无解；否则最终类即全部整数解。

证明。 归纳不变量为“当前剩余类恰为已处理方程的解集”。定理 1.20 保持此不变量；空方程组的解集为所有整数，与初始状态一致。\square

注。 合并后的模数是已经处理模数的最小公倍数。计算时先除以最大公约数再乘，并采用能容纳所得乘积的整数表示。若所有中间模数不超过 L，r 次合并用 O(r\log(L+1)) 次整数除法；任意精度整数的位代价另计。

§6　剩余类上的两个存在性问题

例 1.22（小凯的疑惑）。 设 a,b>1 且 \gcd(a,b)=1。不能写成 ax+by、x,y\in\mathbf N 的非负整数中的最大者为

证明。 0,b,\ldots,(a-1)b 模 a 的余数互异，故每一类恰含其中一项 jb。在该类内，可表示的非负整数恰为 jb+ta、t\in\mathbf N：若 n=ax+by，写 y=j+sa、0\le j<a，则 n=jb+a(x+sb)\ge jb；反向显然。因此该类中小于 jb 的整数均不能表示，最大的候选为 jb-a。取所有类的最大值，得到 (a-1)b-a。因 a,b>1 且互素，此值非负，故确实属于所求集合。\square

注。 若 a=1 或 b=1，所有非负整数均可表示，没有“最大不能表示的非负整数”；式 (1.20) 此时给出 -1，只能作为另行规定的哨兵值。若 a,b 不互素，则有无穷多个非负整数不能表示。

例 1.23（CF338D，GCD Table）。 给定正整数 a_1,\ldots,a_k 与上界 n,Y，求是否存在

置 L=\operatorname{lcm}(a_1,\ldots,a_k)。先检查 L\le n，再用 ExCRT 解

若无解则原问题无解。否则其全部解为 y\equiv y_0\pmod L。取最小正代表 y_*=y_0（若 y_0>0），或 y_*=L（若规范代表 y_0=0）。原问题可解当且仅当

证明。 任意原解均满足 L\mid x 与式 (1.22)。因 a_i 同时整除 L 与 y+i-1，且 L\mid x，有

故原解存在时可把 x 换成 L。同时，给 y 加减 L 不改变所有这些最大公约数；因此最小正代表已经足以判定边界与等式。反向由 (x,y)=(L,y_*) 验证。\square

注。 对 n\times m 表内长度为 k 的连续片段，边界条件为 y+k-1\le m，即取 Y=m-k+1；若 Y<1 则无解。课件简述中的条件 y\le m 对应 Y=m。求 L 时一旦确定 L>n 即可停止；求得候选后，以式 (1.23) 的最大公约数等式作最终检验。

第二章　剩余类环中的有限运算

§1　单位、逆元与幂

定义 2.1。 含幺交换环 R 中的单位为存在乘法逆元的元素。单位全体记为 R^\times，以环乘法组成交换群。对 m>0，置

对任意幺半群中的元素 a，约定 a^0=1，包括零环及环中 a=0 的情形。

命题 2.2。 [a]_m 是单位当且仅当 \gcd(a,m)=1；存在时其逆元唯一，且扩展 Euclid 算法所得 ua+vm=1 中的 [u]_m 即为逆元。

证明。 可逆性等价于 ax\equiv1\pmod m 可解，应用定理 1.18。若 ax=ay=1，则 x=x(ay)=(xa)y=y。\square

注。 U_1 为单元素群，\varphi(1)=1。素数 p 的所有非零剩余类都是单位，故 \mathbf F_p=\mathbf Z/p\mathbf Z 是域。一般合数模数下的非零元素可能不可逆。

构造 2.3（二进制求幂）。 对幺半群元素 a 与 e\in\mathbf N，维护 (z,b,t)，初值为 (1,a,e)。当 t>0，若 t 为奇数则令 z\leftarrow zb，随后令 b\leftarrow b^2、t\leftarrow\lfloor t/2\rfloor。返回 z。

证明。 不变量为 zb^t=a^e。写 t=2q+\epsilon，其中 \epsilon\in\{0,1\}，一次更新后的乘积为 zb^\epsilon(b^2)^q，与更新前相同。终止时 t=0，故 z=a^e。所用乘法次数为 O(1+\log(e+1))。\square

注。 算法只要求结合律，因而也适用于矩阵。模乘采用精确整数乘积后约简，或采用等价的无溢出实现。指数 0 的结果是单位元；模 1 时其规范代表仍为 0。

构造 2.4（批量逆元）。 设 a_1,\ldots,a_r 均为 R 的单位。置 P_0=1、P_i=P_{i-1}a_i。只求一次 P_r^{-1}，随后从 i=r 到 1 依次计算

证明。 两式分别乘以 a_i、P_{i-1} 即得单位元。此过程只用一次求逆和 O(r) 次乘法。\square

命题 2.5。 对素数 p 与 1<i<p，令 q=\lfloor p/i\rfloor、r=p\bmod i，则

因此可由 1^{-1}=1 递推得到 1,\ldots,p-1 的逆元。

证明。 p=qi+r，故 qi\equiv-r\pmod p。这里 1\le r<i<p，因 p 为素数；乘以 i^{-1}r^{-1} 即得结果。\square

注。 此递推利用素数模数下非零余数可逆的性质。对合数模数，即使 i 是单位，余数 m\bmod i 仍可能不是单位，例如 m=8,i=3 时余数为 2。

§2　单位群与指数约简

定理 2.6（Euler）。 若 \gcd(a,m)=1，则

特别地，对素数 p 与 p\nmid a，有 a^{p-1}\equiv1\pmod p，且 a^{-1}\equiv a^{p-2}\pmod p。

证明。 乘以 [a]_m 给出有限集合 U_m 的置换。把其中所有元素相乘，得

右侧公共乘积是单位，可以约去。素数情形有 \varphi(p)=p-1，再乘以 a^{-1}。\square

命题 2.7。 若 m=\prod_{p\mid m}p^{\alpha_p}，则

空积给出 \varphi(1)=1。

证明。 模 p^\alpha 的单位恰为不被 p 整除的剩余类，共有 p^\alpha-p^{\alpha-1} 个。中国剩余定理把模 m 的单位与各素数幂分量的单位逐一对应，故基数相乘。\square

定理 2.8（非单位的指数约简）。 设 m>0、a\in\mathbf Z、e\in\mathbf N。若 \gcd(a,m)=1，可把指数替换为 e\bmod\varphi(m)。对任意 a，当 e\ge\varphi(m) 时有

当 e<\varphi(m) 时保留原指数即可。

证明。 m=1 时结论自动成立。以下逐一考察 p^\alpha\mid m 的完整素数幂分量。若 p\nmid a，因 \varphi(p^\alpha)\mid\varphi(m)，由 Euler 定理可按 \varphi(m) 的倍数改变指数。若 p\mid a，则 a^t 在 t\ge\alpha 时为模 p^\alpha 的零元。又

最后的不等式可由 2^{\alpha-1}\ge\alpha 的归纳证明得到。因此式 (2.7) 两边在此分量中均为零。各分量结论由中国剩余定理合并。单位情形直接应用定理 2.6。\square

注。 非互素情形按定理 2.8 的阈值选择指数。例如 2^2\equiv0\pmod4，但 2^{2\bmod\varphi(4)}=2^0\equiv1\pmod4。递归计算幂塔时，同时保存指数模 \varphi(m) 的余数及其是否达到 \varphi(m)；指数为零时按 a^0=1 处理。

§3　组合数的局部计算

定义 2.9。 对 n\in\mathbf N、k\in\mathbf Z，\binom nk 为 n 元集合中 k 元子集的个数；当 k<0 或 k>n 时置为 0。因而在整数多项式环中

在 0\le k\le n 时，由对子集元素及其补集分别排序，得到整数恒等式 n!=\binom nk k!(n-k)!。

引理 2.10。 在特征为素数 p 的交换环中，

证明。 对 0<j<p，整数恒等式 j!\binom pj=p(p-1)\cdots(p-j+1) 表明 p\mid j!\binom pj。因 p\nmid j!，引理 1.11 给出 p\mid\binom pj，故二项式展开的中间项均消失。\square

定理 2.11（Lucas）。 设 p 为素数，n=\sum_i n_ip^i、k=\sum_i k_ip^i 为 p 进制展开，其中 0\le n_i,k_i<p。则

证明。 在 \mathbf F_p[X] 中反复应用引理 2.10，得到

右侧每个因子选择的指数形如 j_ip^i，且 0\le j_i\le n_i<p。这些指数相加没有进位，故得到 X^k 的唯一可能是逐位取 j_i=k_i；若某位 k_i>n_i，则该系数为零。比较系数即得结论。\square

注。 可用 O(p) 时间与空间预处理 0,\ldots,p-1 的阶乘和逆阶乘，每次求值用 O(1+\log_p(n+1)) 个数位操作，适用于能够存储长度 p 的表的情形。对 n\ge p，n! 模 p 为零，无乘法逆元。

定理 2.12（Legendre）。 对素数 p 与 n\in\mathbf N，记 s_p(n) 为 n 的 p 进制数位和。则

这里 0!=1，故两边在 n=0 时均为零。

证明。 每个 1\le t\le n 为总指数贡献 \nu_p(t)，即对每个满足 p^j\mid t 的 j\ge1 贡献一次。交换这两个有限求和次序，得第一式。若 n=\sum_i n_ip^i，则

即第二式。\square

命题 2.13。 设 t_1+\cdots+t_r=n，其中 t_i\in\mathbf N。多项式系数

为奇数，当且仅当这些 t_i 的二进制加法不发生进位，即它们的二进制 1 位两两不交且并集为 n 的 1 位。

证明。 由定理 2.12，其 2 进赋值等于 \sum_i s_2(t_i)-s_2(n)。在每一位把两份 2^j 合并成一份 2^{j+1}，数位总和恰减少 1；所有进位完成后数位总和变为 s_2(n)。因此差为零当且仅当根本无需进位。\square

§4　素数幂模数与去因子阶乘

定义 2.14。 固定素数 p，定义正整数

两者在 n=0 时均为 1，且均与 p 互素。

命题 2.15。 有整数恒等式

证明。 把 n! 中的因子分为不被 p 整除者与 p,2p,\ldots,\lfloor n/p\rfloor p。后一部分去掉每个因子的一个 p 后成为 \lfloor n/p\rfloor!；再去掉其中全部 p 因子，恰得 U_p(\lfloor n/p\rfloor)。\square

构造 2.16。 设 q=p^\alpha、\alpha\ge1。预处理

若 n=vq+r、0\le r<q，则

把此式代入式 (2.17)，可递归计算 U_p(n)\bmod q。

证明。 对每个完整长度为 q 的区间，所有不被 p 整除的因子模 q 与第一块一一对应；末尾残块对应 1,\ldots,r。递归参数每次除以 p，故最终到达 0。\square

定理 2.17（组合数模素数幂）。 设 0\le k\le n，q=p^\alpha，置

若 e\ge\alpha，则 \binom nk\equiv0\pmod q；否则

证明。 由组合数的整数性，e\ge0。把三个阶乘分别分离出 p 的幂次，在整数中约去公共 p 因子，得到

等式左侧两个 U_p 因子均与 q 互素，故在模 q 的环中可求逆，得到式 (2.21)。若 e\ge\alpha，右侧为零。\square

构造 2.18（一般模数下的组合数）。 对 M>1，先取得分解 M=\prod_i p_i^{\alpha_i}，用定理 2.17 计算每个素数幂分量，再由中国剩余定理合并；M=1 时直接返回规范余数 0。此构造通常称为 ex-Lucas。

注。 素数分解的成本另计。对一个分量 q，表 (2.18) 用 O(q) 时间及空间；按每层二进制求幂的直接实现，一次去因子阶乘计算用 O((1+\log_p(n+1))(1+\log(n+1))) 次模乘。长度 q 的表超出存储容量时，可另选计算 G_p 的方法。组合数先在整数中分离素数幂，再对分母中的单位求逆。

§5　奇偶计数与二进制递推

例 2.19（gym102586E）。 给定互异的非负整数 A_1,\ldots,A_k，求长度为 n、每项从这 k 个值中选择、总和为 S 的有序序列个数模 2。在 \mathbf F_2[X] 中置

所求答案为 [X^S]P(X)^n。若 n=\sum_{i=0}^{\ell-1}n_i2^i，则

证明。 展开 P(X)^n，每个有序选择恰贡献一个单项式，指数为总和。特征 2 下由引理 2.10，P(X)^{2^i}=P(X^{2^i})，于是得到式 (2.24)。它也可由命题 2.13 解释：每个 n 的二进制 1 位恰分配给一个出现次数，发生进位的多项式系数全部为偶数。\square

构造 2.20（进位状态）。 设 A=\max_j A_j，n>0，\ell 为 n 的二进制位数。令 s_i 为 S 的第 i 个二进制位。定义 f_i(c)\in\mathbf F_2 为满足

的低 i 位选择的奇偶数，其中仅在 n_t=1 时选择 j_t。初值 f_0(0)=1，其余为零。若 n_i=0，对每个状态仅取 a=0；若 n_i=1，则依次取 a=A_1,\ldots,A_k。仅当 c+a\equiv s_i\pmod2 时执行

最终答案为 f_\ell(\lfloor S/2^\ell\rfloor)；数组范围外的项视为零。

证明。 把 n_i2^ia 加到式 (2.25) 的左侧，并用 S\bmod2^{i+1}=(S\bmod2^i)+s_i2^i，恰得更新式。通过奇偶条件后，新的进位非负。若 0\le c\le A 且 0\le a\le A，则 0\le(c+a-s_i)/2\le A，故所有可达状态均在 [0,A] 中。处理完 n 的全部位后，完整和等于 S 当且仅当剩余进位为 \lfloor S/2^\ell\rfloor。\square

注。 n=0 时空序列仅有一种，答案为 [S=0]。直接实现用 O((k+1)(A+1)\ell) 次布尔运算、O(A+1) 空间。可把状态写成位集：在 n_i=1 时，先对各 A_j 的平移取异或，再按 s_i 保留奇位或偶位并压缩；中间下标范围至多 2A。若一个机器字含 w 位，移位与异或阶段用 O(k\lceil(A+1)/w\rceil\ell) 次字操作，普通压缩阶段另用 O((A+1)\ell) 时间。

第三章　循环子群与整数分解

§1　阶与原根

定义 3.1。 设 G 为有限群、g\in G。g 的阶为满足 g^d=1 的最小正整数 d，记为 \operatorname{ord}(g)。由 g 的全部整数次幂组成的子群记为 \langle g\rangle；若 G=\langle g\rangle，则称 G 为循环群，g 为其生成元。模 m 的原根为 U_m 的生成元。

命题 3.2。 有限群元素的阶存在。若其阶为 d，则

且 d\mid|G|。

证明。 有限性使 1,g,g^2,\ldots 中存在相等的两项，消去后得到一个正指数的幂为 1。任一整数 t 写成 qd+r、0\le r<d，则 g^t=1 当且仅当 r=0，由此得到式 (3.1)。子群 \langle g\rangle 的左陪集两两相同或不交，且每个陪集经左乘映射与该子群等势；因此这些大小为 d 的集合划分 G，得 d\mid|G|。\square

引理 3.3。 域上次数为 r 的非零多项式至多有 r 个根。若交换群中元素 u,v 的阶分别为互素的 a,b，则 uv 的阶为 ab。

证明。 若多项式以 c 为根，则可除以 X-c，所得余数为零；对次数归纳即得根数界。对第二项，(uv)^{ab}=1。若 (uv)^t=1，则 u^t=v^{-t} 同时属于 \langle u\rangle 与 \langle v\rangle；其阶同时整除 a,b，因而为 1。于是 a\mid t 且 b\mid t，故 ab\mid t。\square

定理 3.4。 域的乘法群的每个有限子群均为循环群。特别地，对每个素数 p，模 p 的原根存在。

证明。 设有限子群为 G，令 e 为所有元素阶的最小公倍数。对每个 q^\alpha\Vert e，选择一个阶含完整因子 q^\alpha 的元素；把它提升到适当的整数次幂，可得到恰为 q^\alpha 阶的元素。由引理 3.3，把这些阶两两互素的元素相乘，得到一个阶为 e 的元素，故 e\le|G|。另一方面，所有 x\in G 都满足 x^e=1，域上的根数界给出 |G|\le e。因此 e=|G|，所构造元素生成整个 G。\square

注。 上述证明中，若 g 的阶为 d，则 g^t 的阶为 d/\gcd(d,t)，这直接由 d\mid ts 的最小正解得到。模合数的单位群不一定循环，例如模 8 的四个单位的平方均为 1，所以没有阶为 4 的元素。

命题 3.5。 设 g\in G，已知 g^D=1，并已知 D 的素因子。则 \operatorname{ord}(g)=D 当且仅当

证明。 阶 d 整除 D。若 d<D，则 D/d>1，取素数 q\mid D/d，便有 d\mid D/q，所以 g^{D/q}=1。反向由 d=D 立即得到。\square

构造 3.6。 给定单位 g\in U_m 及 \varphi(m) 的素因数分解，从 D=\varphi(m) 开始：对素数 q\mid D，只要 g^{D/q}=1 就令 D\leftarrow D/q。最后所得 D 为 g 的阶。

证明。 每次修改仍保持真实阶整除 D。停止时式 (3.2) 的全部检验成立，故 D 就是真实阶。\square

§2　离散对数与分块查找

定义 3.7。 给定有限群 G 中的 g,h，离散对数问题为求 x\in\mathbf N 使 g^x=h。可解当且仅当 h\in\langle g\rangle；可解时全部整数解构成模 \operatorname{ord}(g) 的一个类。

构造 3.8（BSGS）。 设已知 \operatorname{ord}(g)\le N，N\ge1，令 B=\lceil\sqrt N\rceil。存储

同一键出现多次时只保留最小 j。按 i=0,1,\ldots,\lfloor(N-1)/B\rfloor 的顺序，查询

若查询到 j 且 iB+j<N，返回 x=iB+j；若遍历结束仍未找到，则无解。

证明。 每个 0\le x<N 唯一写成 iB+j、0\le j<B。等式 g^x=h 等价于 g^j=h(g^{-B})^i。任何可解问题都存在 0\le x<\operatorname{ord}(g)\le N 的解，故查找完备。块按 i 递增，同一块返回最小的 j，因而第一个有效候选为最小非负解。\square

注。 使用期望常数时间查找的散列表时，所用群运算与期望查找次数均为 O(\sqrt N)，空间为 O(\sqrt N)，另加一次求逆。若使用比较树或排序后的二分查找，应另计相应对数因子。模 m 的单位情形可直接取 N=m，无需先分解 m；若知道更小的阶上界，可相应减小表。目标不是单位时，单位底数的方程无解。

构造 3.9（固定底数的有界指数表）。 若需多次求固定 g 的幂，且所有指数均满足 0\le e<K，令 B=\lceil\sqrt K\rceil，分别预处理 g^j、g^{iB}，其中 0\le j<B、0\le i\le\lfloor(K-1)/B\rfloor。查询时写 e=iB+j，返回两表元素之积。

证明。 幂的结合律给出 g^{iB+j}=g^{iB}g^j。两表可分别连续乘以 g 与 g^B 建成，预处理及空间为 O(\sqrt K)，单次查询只需一次乘法。\square

§3　非单位底数的剥离

构造 3.10（exBSGS）。 给定 m>0、整数 a,b，求 a^x\equiv b\pmod m 的最小非负解。维护整数 t、正整数 M 及模 M 的规范代表 B,C，初值为

每轮依次作如下处理。

1. 若 B=C，返回 t。
2. 令 g=\gcd(a,M)。若 g=1，以 BSGS 解 a^y\equiv BC^{-1}\pmod M；有解则返回 t+y，无解则报告无解。
3. 若 g>1 且 g\nmid B，报告无解。
4. 否则同步更新

并继续。式 (3.6) 右侧均使用更新前的变量，最后把 B' 规范到模 M' 的余数。

证明。 归纳保持三个不变量：所有 x<t 已被排除；C 为模 M 的单位；对任意 x\ge t，原方程等价于

初值显然满足。第一步恰检验 x=t。若其失败，则只需考虑 x-t\ge1。若 g\nmid B，左侧被 g 整除而右侧不能，故无解。若 g\mid B，按式 (1.12) 约去一个 g，得到

正是更新后的条件。因 \gcd(a/g,M/g)=1，且旧 C 在模 M/g 下仍为单位，新 C' 也是单位。若 g=1，则可约去 C，交给 BSGS；其最小解加上 t 就是原方程的最小解。每次剥离都使 M 至少减半，因此过程终止。\square

注。 m=1 时初值即有 B=C=0，返回 x=0。a=0 时，按 0^0=1 的约定先检验指数 0；其后只有目标余数为 0 才能取 x=1。累计系数逐轮乘入当前的 a/g。至多 O(\log(m+1)) 次剥离；按直接 Euclid 成本计，剥离阶段用 O(\log^2(m+1)) 次整数除法，末尾 BSGS 用 O(\sqrt m) 次模乘与期望散列操作。

§4　两个指数方程

例 3.11（P4884）。 给定素数 m\ge7 与 0<k<m，求最小正整数 n 使十进制循环单位数满足

由整数恒等式 9(1+10+\cdots+10^{n-1})=10^n-1，且 9 在模 m 下可逆，原条件等价于

因为 m\ne2,5，底数 10 为单位，可用 BSGS。若目标为零，则无解；目标不在 \langle10\rangle 中时也无解。

证明。 两个同余条件的等价性由乘以单位 9 得到。另因 0<k<m 且 m\ne3，有 9k+1\not\equiv1\pmod m，故式 (3.10) 的最小非负解一旦存在就已经为正数。\square

注。 对更一般的输入，目标为 1 时最小非负解为 0，最小正解为底数的阶。也可改写 n=t+1，对 10^t\equiv(9k+1)10^{-1} 求最小非负 t。模数与 9 不互素时，把式 (3.9) 乘以 9 得到必要条件，候选解再由原式检验。

例 3.12（CF1106F）。 设 p=998244353，给定 b_1,\ldots,b_k，序列在 \mathbf F_p^\times 中满足

给定 n>k 与 m\in\mathbf F_p^\times，求满足 f_n=m 的 f_k。取原根 g=3，令 f_i=g^{z_i}，则指数属于 \mathbf Z/(p-1)\mathbf Z，并满足

令 T 为 k\times k 矩阵，其第一行为 (b_1,\ldots,b_k)，其次对角线为 1，其余为 0。则

置 A=(T^{n-k})_{1,1}，用 BSGS 求 B 使 g^B=m。问题即化为

令 d=\gcd(A,p-1)。若 d\nmid B，则无解；否则由定理 1.18 求一个解 z_0，全部可能值为

证明。 原根使 z\mapsto g^z 成为从加法群 \mathbf Z/(p-1)\mathbf Z 到乘法群 \mathbf F_p^\times 的同构，故乘积对应和、乘方对应整数倍。矩阵逐次实现式 (3.12)，初始指数中只有 z_k 可能非零，于是得到式 (3.14)。线性同余的解类互异，原根的单射性使式 (3.15) 的 d 个值亦互异。\square

注。 本题的素数性质作为给定常数的已知性质使用；p-1=2^{23}\cdot7\cdot17。由二进制求幂可核验

三者均不为 1，故命题 3.5 证明 3 为原根。矩阵运算的模数为合数 p-1，式 (3.14) 按定理 1.18 求解。朴素矩阵乘法与 BSGS 合计用 O(k^3\log(n+1)+\sqrt p) 次基本算术及期望查找操作。若 A=0，可解条件恰为 B=0，此时全部 p-1 个非零初值都可行。

§5　素性测试

引理 3.13。 在 \mathbf F_p 中，方程 x^2=1 的解只能为 x=1 或 x=-1。

证明。 (x-1)(x+1)=0，而域中没有零因子。对 p=2，这两个写法表示同一个元素。\square

定义 3.14。 设奇数 n>2，唯一地写成 n-1=2^s d，其中 d 为奇数、s\ge1。对整数 a，若

则称 n 通过底数 a 的强伪素数检验。

命题 3.15。 若 n 为奇素数且 n\nmid a，则它通过式 (3.17) 的检验。

证明。 Fermat 定理给出 a^{2^sd}=1。考察平方链 a^d,a^{2d},\ldots,a^{2^sd}。若首项不为 1，取第一次变成 1 的位置；其前一项是 1 的一个不等于 1 的平方根，由引理 3.13 必为 -1，而其下标小于 s。\square

构造 3.16（Miller–Rabin）。 对待测整数 n，先排除 n<2，接受 n=2，排除大于 2 的偶数。对每个选定的底数 a，先规范其模 n 的代表；若选择规则保证 1\le a<n，可直接检验。若 1<\gcd(a,n)<n，则已找到合数证据。否则计算 v=a^d\bmod n：若 v=1 或 v=-1 则本轮通过；其余情况下至多再平方 s-1 次，出现 -1 则通过，始终未出现则判定为合数。

证明。 该流程精确检验式 (3.17)。命题 3.15 表明，合法底数不会把素数误判为合数。若平方途中出现由非 \pm1 元素平方得到的 1，引理 3.13 本身也已给出合数证据。\square

注。 一轮用 O(\log n) 次模乘。Rabin 的强伪素数计数定理断言：对任意奇合数 n，在 1,\ldots,n-1 中通过式 (3.17) 的底数至多占四分之一，见 M. O. Rabin，Probabilistic Algorithm for Testing Primality，1980，定理 1。因此，独立均匀选取 t 个底数时，固定合数全部通过的概率至多为 4^{-t}。该概率针对固定合数与随机底数计算。

另有经典的确定性范围结果：对 2\le n<2^{64}，依次检验底数

足以确定素性。其依据是前十二个素数底数共同的最小强伪素数为 318665857834031151167461>2^{64}，见 J. Sorenson 与 J. Webster，Strong Pseudoprimes to Twelve Prime Bases，定理 1.1。实现中先处理小于等于 37 的输入，之后这些底数均在 [1,n-1] 内；或者对底数取模并跳过零余数。无符号 64 位输入的模乘采用更宽的乘积表示或专门的模乘实现。

§6　从碰撞提取因子

构造 3.17（Pollard–Rho 的单次尝试）。 对奇合数 n，选取整数 c 与初值 x_0，在 \mathbf Z/n\mathbf Z 中置

从 x=y=x_0 开始，反复令 x\leftarrow F(x)、y\leftarrow F(F(y))，并计算

若 d=1，继续迭代；若 1<d<n，返回 d；若 d=n，本次尝试失败，换取参数重试。

命题 3.18。 构造 3.17 返回的值一定为 n 的非平凡因子。若某个素因子 p\mid n 满足 x\equiv y\pmod p，则 p\mid d。

证明。 第一项来自最大公约数的定义及返回条件。第二项因为 p 同时整除 x-y 与 n。\square

注。 返回的非平凡因子可以是素数或合数。算法在已知模数 n 下迭代，寻找模未知素因子 p 的碰撞。由于 F 为整数多项式，模 n 的轨道映射到模 p 后仍是同一多项式的轨道。快慢指针用于寻找有限函数轨道的循环；二者模 n 相等时，式 (3.20) 给出 n，本次尝试按失败处理。

构造 3.19（递归分解）。 对 n>0，若 n=1 则返回空族；若素性判定确认 n 为素数，则返回 n；否则提取因子 d，分别分解 d 与 n/d，合并所得素因子及其重数。偶因子 2 可先直接取出。

证明。 只要每次取出的因子满足 1<d<n，两个子问题都严格缩小。归纳所得因子乘积为 n；当所有叶子确实为素数时，唯一分解定理保证输出正确且重数唯一。\square

注。 在把轨道模 p 的行为近似为随机映射的模型下，生日碰撞启发式预测约 O(\sqrt p) 次迭代出现碰撞。若 p 为合数 n 的最小素因子，则 p\le\sqrt n，得到 O(n^{1/4}) 的启发式迭代量级。这一近似以均匀随机映射代替实际的二次多项式族。每次迭代中的最大公约数计算另计整数除法及位代价。参数退化时重新选择 c,x_0，亦可设定单次尝试的迭代上限。

采用有确定输入范围的 Miller–Rabin 底数集时，递归停止条件在该范围内是确定的；采用随机底数时，素性测试的误判可能传递到最终分解。提取的因子可用精确整除核验，完整素因数分解的核验包括因子乘积与各因子的素性。

构造 3.20（批量最大公约数）。 设一次碰撞搜索连续产生差值 \Delta_1,\ldots,\Delta_b。先计算

若 d=1，整块没有提供因子，继续处理下一块；若 1<d<n，直接返回 d；若 d=n，回到这块的起点，重新产生各 \Delta_j，逐个计算 d_j=\gcd(\Delta_j,n)，一旦 1<d_j<n 就返回。若整块均未给出非平凡因子，则本次参数失效，重新开始。

证明。 一个整数乘积与 n 互素，当且仅当每个因子与 n 互素：正向因为每个因子整除乘积，反向由素因子判别得到。因此 d=1 时跳过整块不遗漏因子。1<d<n 时，无论 d 来自一个还是多个差值，它本身都是合法因子。d=n 则不能断定每个差值都与 n 有最大公约数 n；例如 n=15、差值为 3,5 时，积为零而各自可分解。逐项回放可恢复这些被乘积合并的因子。\square

注。 每块的 Q 从 1 重新开始；保存块首的轨道状态与块长即可回放。批处理通常每 b 个差值调用一次最大公约数，失败回放时增加调用。积在每一步模 n 约简，保持式 (3.21) 的最大公约数。某个差值为零时，批量积为零，进入失败回放。

构造 3.21（Brent 分段轨道）。 固定 F(x)=x^2+c\pmod n 与块长上限 b\ge1。令 y=x_0、r=1。一次外层过程如下。

1. 保存 x=y，再让 y 连续前进 r 步。
2. 在随后的 r 步中，每前进一步就产生差值 x-y；每至多 b 步组成一块，用构造 3.20 检验。块的起点、实际步数都须保存，以便回放。
3. 若得到非平凡因子则返回，若判定本次参数失效则重启；否则令 r\leftarrow2r，重复外层过程。

这里“前进一步”指赋值 y\leftarrow F(y)。所有减法及乘积均可立即规范到模 n 的余数。

命题 3.22。 构造 3.21 的每次返回均为有效分解。对任意 p\mid n，若模 p 的轨道有长度为 \mu 的前导部分及长度为 \lambda 的循环部分，则当某一段的固定点 x 已落在循环中且 r\ge\lambda 时，该段产生的差值至少有一个被 p 整除。

证明。 返回正确性由构造 3.20 得到。固定点 x 已处于循环时，后续任何连续 \lambda 个轨道点都会经过 x 的模 p 剩余类。该段检测连续 r\ge\lambda 个点，故至少出现一次模 p 碰撞。随着 r 倍增，固定点的轨道下标不断增加，最终越过 \mu；只要尚未重启，也最终满足 r\ge\lambda。\square

注。 多个素因子的碰撞同时发生，或模 n 的轨道过早闭合时，可能得到 d=n，此时重启。Brent 分段使固定点在一整段内保持不变，便于批量乘积与回放；它调整运算常数与最大公约数调用频率，迭代量级仍按前述随机模型作启发式估计。返回因子可能只含某个素因子的部分重数，递归分解剩余商即可恢复其余重数。

第四章　算术函数与卷积

§1　以正整数乘法为下标的代数

本章固定含幺交换环 R。整数系数通过唯一的含幺环同态 \mathbf Z\to R 解释；除非另行指出，所有约数均为正约数。环中的除法仅表示乘以已知存在的逆元。

定义 4.1。 R 上的算术函数为映射 f:\mathbf N_{>0}\to R。全体此类函数组成集合 \mathcal A_R。在此集合上定义逐点加法及 Dirichlet 卷积：

记 \varepsilon(1)=1_R，\varepsilon(n)=0_R（n>1）。逐点乘法另记为 fg；它与 f*g 为不同运算。

命题 4.2。 \mathcal A_R 在逐点加法与卷积下为含幺交换环，其乘法单位元为 \varepsilon。

证明。 每个正整数仅有有限多个约数，故定义中的求和存在。交换律来自 (a,b)\mapsto(b,a)。结合律来自同一有限和的两种分组：

分配律由 R 的分配律逐项得到。最后，\varepsilon 使卷积中仅留下 a=1 或 b=1 的一项。\square

定理 4.3。 f\in\mathcal A_R 在卷积下可逆，当且仅当 f(1) 为 R 的单位。此时逆元 g=f^{*-1} 唯一，且由下列递推确定：

证明。 若 f*g=\varepsilon，在 1 处求值得 f(1)g(1)=1_R，故条件必要。反之，式 (4.3) 中 n/d<n，可按 n 递增定义 g。代回卷积即得 f*g=\varepsilon。逆元唯一性由结合律得出。\square

在整数值算术函数环中，条件为 f(1)=\pm1；在模 M 的剩余类环中，条件为 \gcd(f(1),M)=1。

§2　积性与素数幂上的局部数据

定义 4.4。 若 f(1)=1_R，且对互素的正整数 a,b 有 f(ab)=f(a)f(b)，则称 f 为积性函数。若上式对任意正整数 a,b 成立，则称 f 为完全积性函数。

当 R 非零时，要求 f(1)=1_R 排除了零函数，并使空乘积与 1 的取值一致。在有零因子的环中，不能从某个非零函数值随意约去因子来补出这一条件。

命题 4.5。 一个积性函数由所有素数幂上的值唯一确定。更确切地，任给元素族 c_{p,e}\in R，其中 p 遍历素数、e\ge1，存在唯一积性函数满足 f(p^e)=c_{p,e}，并有

它完全积性，当且仅当 c_{p,e}=c_{p,1}^{\,e} 对所有 p,e 成立。

证明。 算术基本定理把 n 唯一分解为两两互素的素数幂，故式 (4.4) 为必要条件，且该式直接满足积性。完全积性时可反复拆出一个 p；反之，素因子指数的可加性保证任意两整数相乘时函数值相乘。\square

命题 4.6。 积性函数的逐点积、Dirichlet 卷积及卷积逆元仍为积性函数。

证明。 逐点积的情形直接由定义得到。若 \gcd(m,n)=1，则 mn 的每个约数唯一写成 ab，其中 a\mid m、b\mid n。因此

为证明逆元的情形，对每个素数 p 按式 (4.3) 规定 g(p^e)，再由命题 4.5 将它们延拓为积性函数。卷积 f*g 积性，并在 1 处为 1_R、在每个正次素数幂处为 0_R，故等于 \varepsilon。\square

定义 4.7。 对素数 p，称

为 f 在 p 处的 Bell 级数。这里 p^0=1，且级数是形式幂级数。

命题 4.8。 对任意算术函数 f,g 均有

若 f 积性，则全部 Bell 级数共同确定 f；仅考察素数幂足以验证两个积性函数相等。

证明。 p^e 的约数恰为 1,p,\ldots,p^e，故两边 X^e 的系数均为 \sum_{j=0}^ef(p^j)g(p^{e-j})。后一断言由命题 4.5 得出。\square

对于一般算术函数，素数幂上的数据与诸如 f(6) 的其余取值可独立指定。

§3　基本函数与反演

定义 4.9。 对 k\in\mathbf N，定义

写 \operatorname{id}=\operatorname{id}_1、\tau=\sigma_0。整数值函数 \varphi 与 \mu 定义为

空乘积约定给出 \varphi(1)=\mu(1)=1。将这些整数值映入 R，便得到本章所用的 R 值函数。

命题 4.11。 有

证明。 若 n 有 r 个不同素因子，非零项恰对应这些素因子的子集，故和为 (1-1)^r。n=1 时只有空子集，其贡献为 1。最后一式以 n=\gcd(a,b) 代入。\square

定理 4.12（约数反演）。 对任意 f,g\in\mathcal A_R，下列两式等价：

特别地，\varphi=\operatorname{id}*\mu。

证明。 第一式为 g=f*\mathbf1。两边与 \mu 卷积，并使用结合律及式 (4.12)，即得第二式；反向推导相同。\square

反演不要求 f,g 积性。积性使函数可由局部数据重建，反演则是卷积环中的可逆运算；二者承担不同作用。

定理 4.13（有限倍数反演）。 固定 N\ge1，设 f,g:[1,N]\to R。下列两式等价：

证明。 把第一式代入第二式的右边，并把 jk=t 的项合并，得到

两种线性变换交换，因为复合时均按乘积 jk 合并；于是反向等价亦成立。\square

有限截断或有限支撑使倍数求和成为有限和，从而使式 (4.14) 在任意系数环上有定义。

命题 4.14（完全积性加权）。 若 t 完全积性，则映射 T_t(f)=tf 保持加法、卷积与卷积单位元。特别地，

证明。 在 ab=n 时，t(a)t(b)=t(n)，可从有限和中提出。再取 f=\mathbf1、g=\mu。\square

这里不要求每个 t(n) 可逆。若 t 仅积性，a,b 可能不互素，第一步便不成立。

定理 4.15（取整求和的加权反演）。 设 t 完全积性，F,G:\mathbf N\to R 且 F(0)=G(0)=0。则

证明。 对正整数 i,j，有 \lfloor\lfloor n/i\rfloor/j\rfloor=\lfloor n/(ij)\rfloor；这可直接由整数 q 满足 qj\le\lfloor n/i\rfloor 当且仅当 qij\le n 得出。因此两种取整求和算子的复合等于以 t*(\mu t) 为系数的算子。式 (4.16) 说明此复合为恒等算子。\square

例如，各个局部因子为

这里的分母均具有常数项 1_R，所以在 R[[X]] 中可逆。取复数 s 时，这些形式级数的解析解释以收敛和交换求和的合法性为条件。

§5　除数偏序上的前缀变换

定义 4.18。 对 f:[1,N]\to R，定义除数 zeta 变换与倍数 zeta 变换为

这里的偏序为整除关系。素因子指数给出其坐标：a\mid b 当且仅当 \nu_p(a)\le \nu_p(b) 对所有素数成立。由此，两种变换分别是这些坐标上向下与向上的有限累加。

构造 4.19。 已知不大于 N 的素数表。初始化 A[n]=f(n)。逐个枚举素数 p\le N；对每个 p，按 i=1,2,\ldots,\lfloor N/p\rfloor 的递增次序执行

全部操作完成后，A=Z_\downarrow f。恢复原数组时，倒序枚举上述素数，并对每个 p 按 i=\lfloor N/p\rfloor,\ldots,1 执行 A[ip]\leftarrow A[ip]-A[i]。

证明。 固定 p，将下标唯一写成 mp^e，其中 p\nmid m。递增的 i 使链

逆变换通过链上的相邻差分直接实现与 \mu 的卷积；规定的遍历方向使各素数幂按依赖次序传播。

构造 4.20。 对倍数变换，仍初始化 A=f；对每个素数 p，按 i=\lfloor N/p\rfloor,\ldots,1 执行

其结果为 Z_\uparrow f。逆变换倒序枚举素数，对每个素数按 i 递增执行 A[i]\leftarrow A[i]-A[ip]。

证明。 仍分解为 mp^e 的链；此时在每条链上执行后缀和。对已处理素数归纳，A[n] 为所有满足 n\mid m\le N 且 m/n 的素因子均在已处理集合中的 f(m) 之和。逆变换的正确性由逐项撤销得到。\square

两种变换及其逆变换均执行

次环加减法，存储量为 O(N)；素数倒数和的上界将在命题 6.3 中证明。朴素地枚举每个数的全部倍数则需 \sum_{d\le N}\lfloor N/d\rfloor=O(N\log(N+1)) 次操作。上述算法按素数坐标逐维累加整除关系，适用于任意输入数组。

第五章　有限双曲线上的求和

§1　商值集合及其闭包

本章将求和指标限制在有限集合中。取整、比较及端点计算在整数中进行；权值可取于含幺交换环 R。对算术函数 f，记

定义 5.1。 对 n\ge1，定义商值集合

另置 D_0=\varnothing。商值集合为正整数集合；在算法中需要的零前缀值单独保存。

引理 5.2。 对 a\in\mathbf N、b,c\ge1，有

若 m\in D_n，则 D_m\subseteq D_n。

证明。 对任意 q\in\mathbf N，条件 q\le\lfloor\lfloor a/b\rfloor/c\rfloor 等价于 qbc\le a，也等价于 q\le\lfloor a/(bc)\rfloor；两边因此相等。若 m=\lfloor n/b\rfloor 且 1\le c\le m，则 bc\le n，所以 \lfloor m/c\rfloor=\lfloor n/(bc)\rfloor\in D_n。\square

后一结论使递归求前缀和时产生的状态留在同一个有限集合内。它是后续亚线性筛组织状态的依据。

命题 5.3。 令 u=\lfloor\sqrt n\rfloor，则

证明。 若 i>u，则 \lfloor n/i\rfloor\le u，因为 n<(u+1)^2。其余商值已列在第二个集合中。\square

命题 5.4。 对 q\in D_n，所有使 \lfloor n/i\rfloor=q 的下标构成区间

若某一尚未处理的区间从 l\le n 开始，则包含 l 的商值块右端点为

证明。 \lfloor n/i\rfloor=q 等价于 q\le n/i<q+1，即 n/(q+1)<i\le n/q。取整数端点即得结论。\square

从 l=1 开始，依次求 q,r，处理 [l,r]，再令 l=r+1，恰好遍历全部最大商值块。若上限为 L\le n，将 r 截断为 \min(r,L) 即可。

§2　按商值分块

定理 5.5。 设 f:[1,n]\to R，g:D_n\to R。则

右侧遍历最大商值块。若 S_f 与 g 的每次求值均为 O(1)，则求和需 O(\sqrt n) 次算术操作与 O(1) 个额外存储单元；前缀表等预处理及其存储量须另计。

证明。 每个块上 g 的自变量不变，故可从块内求和中提出。全部块不交且覆盖 [1,n]，而块数由命题 5.3 控制。\square

命题 5.6。 对 n,i\ge1，有 \lceil n/i\rceil=\lfloor(n-1)/i\rfloor+1。因此，对 n\ge1，

若一个被加项同时依赖 \lfloor a_1/i\rfloor,\ldots,\lfloor a_t/i\rfloor，则每一步分别求这些商值块的右端点，并取最小值，便得到公共细分。

证明。 第一式由带余除法直接验证，第二式将 i=n 单独取出。公共细分中，只要所有商值尚未改变，被加项对这些商值的依赖便保持不变；任一商值首次改变的位置正是当前公共块之后。\square

对正商值使用式 (5.6)。若 a_j<i，该商值为零并在之后恒为零，其右端点取为本次求和的总上限。当 t 固定时，公共块数至多为各个分割的块数之和，因而为 O(\sum_j\sqrt{a_j}+t)。

§3　双曲线恒等式

命题 5.7。 对 f,g\in\mathcal A_R 与 n\in\mathbf N，有

证明。 第一式将每个乘积为 k\le n 的有序对 (a,b) 归入第 k 项；第二式固定第一坐标。所有集合有限，故可交换求和次序。\square

定理 5.8（Dirichlet 双曲线恒等式）。 设 u,v\in\mathbf N 满足

则

证明。 在区域 H_n=\{(a,b)\in\mathbf N_{>0}^2:ab\le n\} 中，不可能同时有 a\ge u+1 与 b\ge v+1。故 a\le u 与 b\le v 两个子区域覆盖 H_n。二者的交集恰为 [1,u]\times[1,v]，因为 uv\le n。对这两个有限集合应用容斥，得到式 (5.11)。\square

推论 5.9。 可取 u=v=\lfloor\sqrt n\rfloor；也可任选 1\le u\le n，令 v=\lfloor n/u\rfloor。若单点函数值及前缀值均可在常数时间取得，则式 (5.11) 的求值代价为 O(u+v)，对称选取时为 O(\sqrt n)。

证明。 对称选取显然满足式 (5.10)。后一选取有 uv\le n<u(v+1)<(u+1)(v+1)。逐项计数即可。\square

分块与双曲线分割使用同一有限乘积区域的不同分解。前者需要一个函数的区间和及另一个前缀函数；后者同时使用两侧的前缀值。若两侧求值成本不同，非对称的 u,v 可以给出更合适的实际代价。

§4　除数与余数的求和实例

例 5.10（P3935 的数学形式）。 求约数个数前缀和 D(n)=\sum_{k=1}^n\tau(k)。由 \tau=\mathbf1*\mathbf1，

证明。 第一式数有序对 (d,e) 满足 de\le n；第二式是定理 5.8 的对称特例。\square

两种式子分别可用商值分块或直接枚举至平方根求值，复杂度均为 O(\sqrt n)，额外空间为 O(1)。若问区间 [L,U] 上的约数个数总和，答案为 D(U)-D(L-1)。原课件给出的规模达到 1.6\times10^{14}；实现使用能够容纳求和结果与中间乘积的整数类型。

例 5.11（P2261 的数学形式）。 对 n,k\in\mathbf N，求

当 k=0 时结果为零。对 k>0，在商值 q 恒定的块 [l,r]\subseteq[1,\min(n,k)] 上，需要减去

证明。 带余除法给出 k\bmod i=k-i\lfloor k/i\rfloor；i>k 时商为零。块内使用等差数列求和公式。\square

复杂度为 O(1+\min(n,\sqrt k))，额外空间为 O(1)。式 (5.14) 中的除以 2 是整数精确除法；在偶模数下，先消去整数因子，再取模。

例 5.12（ARC068E 的数学形式）。 给定 N 个整数区间 [l_j,r_j]\subseteq[1,M]。对每个 d\in[1,M]，求其中包含至少一个 d 的正倍数的区间数。

第 j 个区间对答案的贡献为

证明。 不大于 r_j 的正倍数有 \lfloor r_j/d\rfloor 个，其中恰有 \lfloor(l_j-1)/d\rfloor 个位于区间左边。二者之差为区间内倍数的数目，贡献只取其是否为正。\square

对每个区间同时分块两个商。若某一公共块 [a,b] 满足式 (5.15)，便在答案的差分数组中令 \Delta[a] 加一、\Delta[b+1] 减一。所有区间处理完后取前缀和。算法需 O(N\sqrt M+M) 时间与 O(M) 空间；每个区间在一个 d 上至多贡献一次，即使它包含多个倍数亦然。

构造 5.13（同一实例的长度分解）。 记区间长度为 \ell_j=r_j-l_j+1。按 d=1,2,\ldots,M 递增处理，将所有满足 \ell_j<d 的区间加入一个支持区间加一、单点查询的数组，记当前位置覆盖数为 c_d(x)。则

证明。 任意连续 d 个整数中恰含一个 d 的倍数，故长度至少为 d 的区间必有贡献。长度小于 d 的区间至多含一个倍数，故对各倍数累加覆盖数时，每个有贡献的短区间恰被统计一次。\square

按长度分桶，利用树状数组维护差分的前缀和，每个区间只加入一次，单次加入或查询需 O(\log(M+1))。由于 \sum_{d\le M}\lfloor M/d\rfloor=O(M\log(M+1))，总复杂度为 O(N\log(M+1)+M\log^2(M+1))，空间为 O(N+M)。这给出课件分块解法之外的一个确定上界，适合其 N 与 M 均较大的参数范围。

例 5.14（P2260 的数学形式）。 对 n,m\ge1，求

令 h=\min(n,m)。全部有序对与对角线分离后，得到

在 q_n=\lfloor n/i\rfloor 与 q_m=\lfloor m/i\rfloor 同时恒定的块 [l,r] 上，对角线的贡献为

证明。 式 (5.18) 从矩形指标集删去 i=j；式 (5.19) 展开 (n-iq_n)(m-iq_m) 后逐项求和。\square

使用 \sum_{i=1}^xi=x(x+1)/2 与 \sum_{i=1}^xi^2=x(x+1)(2x+1)/6，每个块需常数次算术操作。因此总复杂度为 O(\sqrt n+\sqrt m)，额外空间为 O(1)；这一界无需假设 n,m 同阶。原题要求模 19940417，可先在整数因子中精确消去 2,6 再取模，也可在验证其可逆后使用模逆元。

§5　反演与分块的组合

例 5.15（P6156 的数学形式）。 固定正整数 k。对若干 n\le N，求

这里 \mu^2 表示逐点平方。定义

则

证明。 先按 d=\gcd(i,j) 分类，写 i=da,j=db。由式 (4.12)，有

令 t=de 并交换有限求和次序，d^{k+1}e^k=t^kd，故 F_k(\lfloor n/t\rfloor) 的系数为 t^k\sum_{d\mid t}d\mu^2(d)\mu(t/d)=h(t)。\square

命题 5.16。 g,h 均积性，并有

若定义两重前缀和

则

证明。 \operatorname{id}\,\mu^2 积性，因此其与 \mu 的卷积积性；再与完全积性的 \operatorname{id}_k 逐点相乘，得到 h 积性。在素数 p 处，\operatorname{id}\,\mu^2 的 Bell 级数为 1+pX，与 \mu 的 1-X 相乘，得 1+(p-1)X-pX^2，即式 (5.24)。

对式 (5.26)，按 s=a+b 分类。2\le s\le x+1 时有 s-1 个有序对；x+2\le s\le2x 时有 2x-s+1 个。另一方面，T_k(2x)-2T_k(x) 中 s^k 的系数在 1\le s\le x 时为 s-1，在 x+1\le s\le2x 时为 2x-s+1，与上述计数一致。\square

式 (5.26) 给出仅用加法、减法和乘法的预处理方法：计算 s^k（s\le2N），累加得到 U_k,T_k，再求出 F_k(x)（x\le N）。这一过程适用于模 2^{32} 的环。

第六章的线性筛可利用式 (5.24) 在 O(N) 次操作内求出 g，并得到 h 的前缀表 H(x)=\sum_{t\le x}h(t)。式 (5.22) 随后写成

每次询问需 O(\sqrt n) 时间。若对每个 s 独立快速幂，预处理时间为 O(N\log(k+1))。进一步利用 s\mapsto s^k 的完全积性，只对每个素数 p\le2N 计算 p^k，再在线性筛的每次转移中相乘，可将该部分写为 O(N+\pi(2N)\log(k+1))。全部表的空间为 O(N)。

第六章　有限区间上的素数筛

§1　消去合数

引理 6.1。 每个合数 x\ge2 都有不超过 \sqrt x 的素因子。

证明。 写 x=ab，其中 1<a\le b<x，则 a\le\sqrt x；a 的任一素因子满足所求条件。\square

构造 6.2（Eratosthenes 筛）。 对 N\in\mathbf N，建立布尔数组 P[0],\ldots,P[N]。将存在的下标 0,1 置为假，其余置为真。按 p=2,3,\ldots,\lfloor\sqrt N\rfloor 扫描：若 P[p] 为真，就把

处的值置为假。最后所有仍为真的下标恰为不大于 N 的素数。

证明。 被标记的数有因子 p 且至少为 p^2，故是合数。反之，取合数 x 的最小素因子 p。由引理 6.1，p\le\sqrt x\le\sqrt N，且 x/p\ge p，所以 x 在处理 p 时被标记。素数从未被标记，因此扫描到 p 时其标志为真。\square

从 p^2 开始足够：小于 p^2 的合数倍数 ap 满足 a<p，已经由 a 的某个素因子处理。实现中的循环条件可写为 p\le N/p，以避免先计算 p^2 引起溢出。

命题 6.3。 当 x\ge2 时，

证明。 这里只需上界。置 s=1+1/\log x>1。对 p\le x 有 p^{s-1}\le e，故

第二个不等式来自 -\log(1-y)\ge y（0<y<1）。有限个收敛几何级数之积，按唯一分解展开，是正项级数 \sum_{n\ge1}n^{-s} 的子和。因此该乘积不超过

代回式 (6.3) 得结论。这里仅为复杂度估计使用 s>1 的实数收敛级数；其收敛已经由积分上界保证。\square

推论 6.4。 构造 6.2 的时间为 O(N\log\log(N+2))，空间为 O(N)。

证明。 初始化与最终扫描需 O(N)。每个素数 p 至多标记 \lfloor N/p\rfloor 次，故由命题 6.3 得到总上界。\square

§2　区间筛与区间分解

构造 6.5（区间素数筛）。 给定 1\le L\le U 及不大于 \lfloor\sqrt U\rfloor 的素数表，在 [L,U] 上初始化布尔数组。对每个表中素数 p，从

开始标记所有不大于 U 的 p 的倍数；若 1\in[L,U]，另外将 1 标为非素数。剩余元素恰为该区间中的素数。

证明。 p^2 的下界使素数 p 自身不会被删除。任取区间中的合数 x，其最小素因子 p\le\sqrt U 且 x\ge p^2，故一定被标记。其余元素不可能因构造中的操作被误删。\square

令 W=U-L+1、B=\lfloor\sqrt U\rfloor。在素数表已经存在时，扫描时间为

额外空间为 O(W)；构造基础素数表的代价需另加 O(B\log\log(B+2)) 时间与 O(B) 空间。式 (6.6) 中的 \pi(B) 项计入基础素数表的扫描。

构造 6.6（区间整数的完全分解）。 对每个 x\in[L,U] 保存残余值 a_x=x 与初始为空的因子表。对每个素数 p\le B，从 \lceil L/p\rceil p 开始访问区间内的所有倍数 x；不断令 a_x\leftarrow a_x/p，直到 p\nmid a_x，并记录所除次数。全部素数处理后，若 a_x>1，则把 a_x 作为指数为一的素因子加入因子表。

证明。 每次除法均为整数精确除法。一个素因子 p 被处理时，其在 x 中的全部次数被除去；此前除去其他素因子不改变该次数。最终若 a_x>1 仍为合数，它必有不超过 \sqrt{a_x}\le\sqrt U 的素因子，与所有这类素数已被除尽矛盾。残余素因子的平方也不可能整除 x，否则该素因子不超过 \sqrt U，仍会与已被除尽矛盾。\square

从区间中的第一个倍数开始，使 x=p 时的因子 p 也得到记录。对 x=1 不执行除法，其因子表为空。

对每个 x，除法次数至多为 \lfloor\log_2x\rfloor，故一个直接且统一的上界为

次整数操作，另加基础素数表的预处理。输出全部质因子及指数所需空间应计入；粗略上界为 O(W\log(U+1))，若即时处理每个数的后续任务则可减少长期保存量。

§3　唯一生成合数

定义 6.7。 对 n\ge2，记 \operatorname{lpf}(n) 为 n 的最小素因子。一个合数 x 的规范生成对为

此时 p\le\operatorname{lpf}(i)；其中 i\ge2。

构造 6.8（线性筛）。 初始化 \operatorname{lpf}[2..N]=0 与空素数表。按 i=2,3,\ldots,N 扫描：若 \operatorname{lpf}[i]=0，则置 \operatorname{lpf}[i]=i 并将 i 追加到素数表。随后按升序枚举表中素数 p，只要

就置 \operatorname{lpf}[ip]=p。条件不再成立时终止本次内层循环。

定理 6.9。 构造 6.8 给出正确的最小素因子表，每个合数恰被生成一次。时间与空间均为 O(N)。

证明。 任一生成对满足 p\le\operatorname{lpf}(i)，因此 ip 的最小素因子正是 p；所写值正确。反之，给定合数 x，取其规范生成对 (i,p)。由于 i<x，扫描到 x 之前已处理 i；且 p\le\operatorname{lpf}(i)\le i，所以 p 已进入素数表，且满足内层循环条件，故 x 必被生成。任何生成 x 的对都必须具有 p=\operatorname{lpf}(x)，从而 i=x/p，所以生成唯一。

尚未被生成的 i 因而恰为素数。内层循环的有效迭代次数等于合数个数，另有每个外层下标至多一次终止条件检查；初始化与素数追加亦为线性，故总时间为 O(N)。\square

常见实现写成先生成 ip，然后在 p\mid i 时立即退出。因素数升序且首次整除 i 的素数恰为 \operatorname{lpf}(i)，该写法与式 (6.9) 等价。

§4　随规范分解传递函数值

构造 6.10。 对每个 i\ge2，维护唯一分解

并保存 c_i,e_i 与 p^{e_i}。对素数 i=p，取 c_i=1,e_i=1。在线性筛的转移 i\mapsto ip 中，只有两种情况：

命题 6.11。 设 f 积性。若在每个素数处建立局部状态的代价为 C_p，且由 f(p^e) 的局部状态可用 O(1) 次操作得到 f(p^{e+1}) 的局部状态，则可在线性筛中用

计算全部 f(1),\ldots,f(N)，总代价为 O(N+\sum_{p\le N}C_p)，空间为 O(N)，前提是每个局部状态占用 O(1) 个存储单元。

证明。 规范分解中的两因子互素，故式 (6.12) 正确。第一种转移仅需使用已建立的素数局部值；第二种转移扩展已有局部状态，且 c_i<i 的函数值已知。每个合数只需一次常数时间处理，再加素数初始化成本，即得所述界。\square

保存 f(c_i) 与局部素数幂状态，以乘法计算 f(i)，适用于任意系数环。总成本为规范生成与局部状态计算两部分之和。

例 6.12（\mu 与 \varphi）。 初始化 \mu(1)=\varphi(1)=1。素数处置 \mu(p)=-1、\varphi(p)=p-1。在线性筛的转移中，

证明。 p\nmid i 时使用积性。p\mid i 时，ip 含平方因子 p^2，故 \mu(ip)=0；同时将 p^e 提升为 p^{e+1}，使 \varphi 的局部值从 (p-1)p^{e-1} 变为 (p-1)p^e。\square

第二分支的条件为 p\mid i，其中 i 是乘入素数 p 之前的下标。

例 6.13（约数幂和）。 固定 k\in\mathbf N。对式 (6.10) 的局部素数幂，维护

它们满足

证明。 前两式分别增加一个指数与一个求和项；第三式使用互素分解的积性。\square

对 k=0，只需保存 e_i，公式化为 \tau(i)=\tau(c_i)(e_i+1)。对 k=1，保存 p^{e_i} 与其几何和即可。它们均能在 O(N) 时间内求出。对任意输入 k，素数初始化中的 p^k 可用快速幂计算，总代价由命题 6.11 得到 O(N+\pi(N)\log(k+1))。在整数中亦可使用精确商 \tau(i)/(e_i+1) 或 \sigma_1(i)/B_{p,e_i}；一般系数环采用式 (6.15) 的乘法递推。

构造 6.14（预给素数幂值的通用形式）。 若全部 f(p^e)（p^e\le N）已给出，并支持常数时间查询，则维护式 (6.10) 的数据，直接以式 (6.12) 计算 f(i)。所需附加时间与空间均为 O(N)。

证明。 规范分解每次只增加一个素因子，表中局部值可直接读取。素数幂均为 [1,N] 中的不同整数，故局部表大小至多为 N-1；其建立成本则已从本构造中分离。\square

课件中按 cp^e 枚举的过程也可用上述规范分解表示；其总成本为整数生成与素数幂数据准备两部分之和。

§5　短区间分解的应用

例 6.15（Haitang and Math 加强版的数学形式）。 给定 1\le n\le10^{12}，求满足

的正整数 m 的数目，其中 S(m) 为 m 的十进制数位和。

在给定范围内，S(m)\le108：小于 10^{12} 的正整数至多有十二位，端点 10^{12} 的数位和为 1。由于余数小于除数，任何解还满足 1\le S(m)<m。

命题 6.16。 令 B=\min(108,n-1)。式 (6.16) 的解与下列候选的通过者一一对应：枚举 s\in[1,B]，枚举 n-s 的正约数 m，并保留满足

的候选。

证明。 若 m 是原问题的解，取 s=S(m)，则 n-s 为 m 的正倍数，且 s<m\le n，所以 s\le n-1。反之，若 m\mid n-s 且 0<s<m，则 n=(n-s)+s 是以 m 为除数、余数为 s 的带余分解；再由 S(m)=s 得到原等式。又 m\le n-s<n，所以范围条件自动满足。每个 m 对应唯一的 s=S(m)，不会在不同的 s 下重复计数。\square

因此只需同时分解区间 [n-B,n-1]。它至多有 108 个正整数，可用构造 6.6 和不大于 10^6 的素数表完成分解。n=1 时 B=0，候选集为空，答案为零。

若 x=\prod_{j=1}^{r}p_j^{e_j} 的分解已经得到，则以笛卡尔积

枚举约数。唯一分解定理保证这是双射。递增添加每个素数幂可在 O(\tau(x)) 次乘法内枚举全部约数；每个数位和的计算需 O(\log(n+1)) 次整数除法，原题范围内至多十二位，另加端点的处理。

在基础素数表已预处理的条件下，一次询问有下列充分上界：

该式分别计入基础素数的扫描、区间分解中的除法，以及所有候选约数的检查。约数枚举按整个短区间上的约数数目之和计费。多组询问可以共享 10^6 以内的素数表；每组的区间残余值及候选检查独立进行。

第七章　商值集合上的卷积与除法

§1　有限商值表示

定义 7.1。 固定 N\ge1，置

算术函数 f 的前缀和记为 F(x)=\sum_{1\le n\le x}f(n)，并置 F(0)=0。有标号的有限族 (F(v))_{v\in\mathcal Q_N} 称为 f 在 N 处的商值表示，亦即讲义中的块筛。此处 \mathcal Q_N 就是第五章的 D_N。标号保留自变量与函数值的对应，供后续运算查取。

命题 7.2。 有不交分解

特别地，|\mathcal Q_N|=s+K\le2s。若 v\in\mathcal Q_N 且 1\le d\le v，则 \lfloor v/d\rfloor\in\mathcal Q_N。

证明。 对 1\le t\le s，取 i=\lfloor N/t\rfloor\ge t。由 ti\le N<t(i+1) 得 t\le N/i<t+t/i\le t+1，故 t 是商值。大于 s 的商值恰由 i\le K 给出；在这一范围内，相邻的 N/i 相差至少 1，故取整后互异。最后，若 v=\lfloor N/i\rfloor，则

右式为正时仍在 \mathcal Q_N 中。\square

小商值直接以 v 为地址，大商值可用 \lfloor N/v\rfloor 为地址。因此商值表可以用两个长度 O(\sqrt N) 的数组表示，并支持最坏常数时间查询。

命题 7.3。 对整数 0\le t\le N，记 r_N(t)=|\mathcal Q_N\cap[1,t]|，则

证明。 在第二种情形中，大商值 \lfloor N/i\rfloor\le t 等价于 i>N/(t+1)。在 [1,K] 中计数，再加上全部 s 个小商值。\square

§2　卷积的前缀表示

命题 7.4。 设 h=f*g，其前缀和分别为 H,F,G。对 x\ge1 及 u=\lfloor\sqrt x\rfloor，有

证明。 将 H(x) 写为 \sum_{ab\le x}f(a)g(b)。若 a,b>u，则 ab\ge(u+1)^2>x。故每一允许的数对至少满足 a\le u 或 b\le u，而被计两次的部分恰为 [1,u]^2。\square

构造 7.5。 仅给定 F,G 在 \mathcal Q_N 上的值时，可以在 O(N^{3/4}) 次环运算和 O(\sqrt N) 个存储单元内求出 H 的商值表示。

证明。 因 [1,s]\subseteq\mathcal Q_N，可由相邻差得到所有所需的 f(a),g(a)。对每个 x\in\mathcal Q_N 应用式 (7.5)，其余前缀查询由命题 7.2 保证可用。总操作数不超过常数倍的

两个和式分别用积分比较估计。\square

命题 7.6。 若另外能以总代价 P(B) 求得 h(1),\ldots,h(B) 及其前缀和，其中 \sqrt N\le B\le N，则其余商值可用式 (7.5) 求出，总代价为

证明。 大于 B 的商值可写为 \lfloor N/i\rfloor，其中 i<N/B。其总求值代价不超过

小商值已由预处理给出。\square

注。 若另有求出前 B 项及其前缀和的线性预处理方法，即 P(B)=O(B)，取 B\asymp N^{2/3} 即得 O(N^{2/3})。给定两张商值表时可使用构造 7.5；具备上述预处理方法时可采用式 (7.7)。

§3　卷积除法

定理 7.7。 设 h=f*g 且 g(1) 为系数环中的单位。则 F 唯一满足

若 [l,r]\subseteq[2,x] 是商 \lfloor x/b\rfloor 的一个最大常值区间，则其中的和为

证明。 按第二个因子求和得 H(x)=\sum_{b\le x}g(b)F(\lfloor x/b\rfloor)。分离 b=1 并乘以 g(1)^{-1}，得到式 (7.8)。其右侧所有未知前缀的自变量至多为 \lfloor x/2\rfloor，故从小到大唯一确定 F。式 (7.9) 由常值区间内提取共同因子而得。\square

构造 7.8。 若 G(x),H(x) 能以常数次基本操作求值，且 F(1),\ldots,F(B) 能在线性时间内预处理，则对大于 B 的商值按式 (7.8)—(7.9) 记忆化求值，称为杜教筛。取 \sqrt N\le B\le N，其时间与空间分别为

证明。 递归参数严格减小，因而终止；命题 7.2 保证从 N 出发只会访问其商值。每个大状态求值一次，每次至多处理 O(\sqrt x) 个常值区间。对这些状态求和，使用命题 7.6 的同一估计。小状态直接查预处理表，大状态不超过 N/B 个。\square

取 B=\lceil N^{2/3}\rceil 并截断至 N，即得到 O(N^{2/3}) 的界。对代价高于常数的 G,H 求值，在操作数中乘入相应的求值代价。

例 7.9。 记 M(x)=\sum_{n\le x}\mu(n)，\Phi(x)=\sum_{n\le x}\varphi(n)。由 \mu*\mathbf1=\varepsilon 与 \varphi*\mathbf1=\operatorname{id}，对 x\ge1 有

两式均可用构造 7.8 求值。这里的有理式先在整数中取其整值，再映入系数环；这种求值方式也适用于偶数模数。

§4　完全积性权重

命题 7.10。 设 w 完全积性，以 wf 表示逐点乘积，则

证明。 对每个 d\mid n，w(d)w(n/d)=w(n)，故可从卷积和中提出共同因子。\square

例 7.11。 置

则

且

证明。 将 \gcd(i,j)=\sum_{d\mid i,\ d\mid j}\varphi(d) 代入，交换有限求和，再置 i=du,j=dv，得到式 (7.13)。对 w(n)=n^2 使用命题 7.10，得

再用定理 7.7 与 \sum_{i\le x}i^3=T(x)^2，即得式 (7.14)。\square

式 (7.14) 中区间权重为平方和之差，式 (7.13) 中区间权重为 F(r)-F(l-1)。先以杜教筛求 F，再作一次商值分块，得到讲义 P3768 所需的计算。若在模素数 p>3 下执行，2,6 均可逆；其他模数下仍可先用整数恒等式约分。

第八章　素数权重与按最小素因子的分解

§1　删去合数的有限过程

定义 8.1。 设 w 完全积性，p_1<p_2<\cdots 为素数序列。置

其中 \operatorname{lpf}(n) 表示最小素因子。约定 A_0(v)=\sum_{n=2}^{v}w(n)，并记 P_w(t)=\sum_{p\le t}w(p)。

命题 8.2。 若 p=p_j 且 v\ge p^2，则

若 v<p^2，则 A_j(v)=A_{j-1}(v)。

证明。 两阶段之差恰为最小素因子等于 p 的合数。它们唯一写为 pm，其中 m\ge p 且 \operatorname{lpf}(m)\ge p。式 (8.1) 在参数 \lfloor v/p\rfloor\ge p 处减去已经保留的较小素数，恰枚举这些 m。完全积性给出 w(pm)=w(p)w(m)，包括 p\mid m 的情形。若 v<p^2，不存在待删去的合数。\square

构造 8.3。 在 \mathcal Q_N 上初始化 A(v)=\sum_{n=2}^{v}w(n)。按 p\le\sqrt N 的递增次序，对满足 v\ge p^2 的商值按递减次序执行式 (8.2)。结束后 A(v)=P_w(v)。

证明。 一轮中 \lfloor v/p\rfloor<v，所以按 v 递减更新时，右端的较小自变量仍保持上一轮的值。命题 7.2 保证该自变量已在表中。归纳使用命题 8.2 得到各阶段不变量。结束时，任何 n\le v\le N 的合数均有不超过 \sqrt N 的素因子，故全部删去；素数一直保留。\square

命题 8.5。 构造 8.3 的更新次数为

证明。 令 L=N^{1/4}。p\le L 的部分至多 2\sqrt N\pi(L)=O(N^{3/4}/\log N)。对 p>L，所需商值均大于 \sqrt N，个数至多 N/p^2。再将素数按 (2^jL,2^{j+1}L] 分组，使用引理 8.4，得到

两部分相加即得。\square

此界尚须加上生成小素数、计算所有初值和局部权重的代价。若这些代价为 O(\sqrt N)，且权重求值和环运算均为基本操作，则总界仍为式 (8.3)，空间为 O(\sqrt N)。

例 8.6。 对固定 k\ge0，取 w(n)=n^k。其初值为幂和 \sum_{n\le v}n^k-1；因 k 固定，可用固定次数的有理系数整值多项式求值，先在整数中精确约分再映入系数环，或在允许的系数域中插值。构造 8.3 给出所有 \sum_{p\le v}p^k。讲义 P5493 所求

随后只需 O(\sqrt N) 次表查询。若 f(p)=\sum_{k=0}^{d}c_kp^k，分别求得 d+1 张素数幂和表并线性组合，即得到 P_f(v)=\sum_{p\le v}f(p)；固定 d 时仅改变常数。多项式插值在所涉及分母均可逆的系数域或模环中进行。

§3　积性函数的重建

定义 8.7。 设 f 积性，定义

并约定 p_0=1。若 x\le p_j，则 S_f(x,j)=0。

定理 8.8。 当 x>p_j 时，有

特别地，\sum_{n\le N}f(n)=1+S_f(N,0)。

证明。 素数贡献由第一项给出。任一合数 n 有唯一的最小素因子 p_k>p_j，且 p_k^2\le n。若 n=p_k^t，其中 t\ge2，它恰出现于 e=t-1 的第二项；否则唯一写为 p_k^eu，其中 e\ge1、u>1 且 \operatorname{lpf}(u)>p_k，由积性贡献 f(p_k^e)f(u)。此时 p_k^{e+1}<n\le x，故恰被第一项枚举一次。所有分类互不相交。\square

已知素数权重表时，式 (8.6) 给出按最小素因子递增的有限递归，构成 Min_25 方法的重建部分。递归的整数参数严格减小，素数下界严格增大；单位项 1 只在最外层添加。

注。 Min_25 方法的总成本由素数表构造与递归重建两部分组成。前者由命题 8.5 计量，后者在 §5 按递归状态和素数幂分支计数。第九章在同样的素数多项式条件下给出另一套恢复程序及其操作数界。

§4　完全积性情形的逆过程

定理 8.9。 设 g 完全积性，且已知 P_g 在 \mathcal Q_N 上的值及 p\le\sqrt N 处的权重。令 A(v)=P_g(v)。按小素数 p 的递减次序，对满足 v\ge p^2 的商值按递增次序执行

全部执行后，1+A(v)=\sum_{n\le v}g(n)。其更新次数与式 (8.3) 相同。

证明。 固定一轮 p=p_j。由命题 8.2，恢复 A_{j-1}(v) 时需要 A_{j-1}(\lfloor v/p\rfloor)。按 v 递增执行，较小自变量或者已在本轮恢复，或者因小于 p^2 本就未变。故式 (8.7) 逐项撤销这一轮。再按素数递减撤销全部轮次，最终恢复 A_0(v)=\sum_{2\le n\le v}g(n)。\square

这个过程按相反的循环次序求取三角变换的逆。式 (8.7) 仅使用加减乘法，因而同样适用于 g(p)=0。

§5　递归树的独立计数

命题 8.10。 令 P^+(d) 表示 d>1 的最大素因子，置

从 S_f(N,0) 开始，按式 (8.6) 直接深度优先求值，不合并相同状态；在局部函数值和素数前缀查询均为常数代价的条件下，递归阶段的操作数为 O(1+C(N))。进一步，有 C(N)=o(N)；一个较粗的显式上界为 O(N/\sqrt{\log N})。

证明。 一条递归路径按严格递增的素因子依次选取完整素数幂，其乘积记为 d。若最后选取的是 p^e，循环条件正是 dp\le N，而 p=P^+(d)。反之，满足 dP^+(d)\le N 的整数，其全部初始素数幂乘积也满足相应循环条件，故唯一分解给出唯一递归路径。因此非根调用恰由式 (8.8) 的集合标号。每次调用有一次常数代价的素数项计算；每个成功的指数迭代产生一个子调用。每个进入的素数分支至少有一次成功的指数迭代，所以该分支的指数循环终止检查可由其子调用摊销；每次函数调用另有常数次外层终止检查，得到操作数界。

固定 Y\ge2。若 P^+(d)>Y，则 d\le N/Y；若 P^+(d)\le Y，则 d 只使用至多 \pi(Y) 个素数，每个指数至多 \lfloor\log_2N\rfloor。因此

先令 N\to\infty，再令固定的 Y\to\infty，便得 C(N)=o(N)。若取 Y=\sqrt{\log N} 并使用 \pi(Y)\le Y，第二项的对数为 O(\sqrt{\log N}\log\log N)=o(\log N)，因而最终小于 N^{1/2}；第一项为 N/\sqrt{\log N}。\square

由此得到递归的亚线性上界。下一章给出另一套具有固定幂次操作数界的完整恢复程序。

第九章　平方富数上的稀疏卷积

§1　局部差分

定义 9.1。 正整数 d 称为平方富数，若每个整除 d 的素数 p 都满足 p^2\mid d。记其集合为 \mathcal P，并规定 1\in\mathcal P。这就是讲义中的 powerful number。

定理 9.2。 设 f,g 为积性函数且对每个素数 p 有 f(p)=g(p)。存在唯一积性函数 h 使 f=g*h，且 h(d)\ne0 蕴含 d\in\mathcal P。

证明。 因 g(1)=1，其卷积逆存在且积性，故 h=g^{-1}*f 存在、唯一并积性。在素数处比较两边得 f(p)=g(p)+h(p)，故 h(p)=0。若某个 d 恰含一次素因子 p，则 h(d)=h(p)h(d/p)=0。\square

命题 9.3。 给定局部值 f(p^e),g(p^e)，可由

依次求得 h。若特取完全积性函数 g(p^e)=f(p)^e，则简化为

证明。 式 (9.1) 是素数幂处卷积等式中分离 i=0 的结果。对式 (9.2)，直接代入可得

中间各项两两抵消，故由唯一性得到结论。\square

望远镜求和只使用加减乘法，适用于 f(p)=0、含零因子的系数环以及合数模数。

§2　支撑集合的计数与枚举

引理 9.4。 每个 d\in\mathcal P 唯一写为

因此 |\mathcal P\cap[1,X]|=O(\sqrt X)，并且

证明。 对 d 的每个正素因子指数 e，若 e 为偶数，令它在 a 中的指数为 e/2，在 b 中为 0；若 e 为奇数，则 e\ge3，相应指数为 (e-3)/2 与 1。逐素数的奇偶性也证明唯一性。于是

最后按 2^j\le d<2^{j+1} 分组，每组的加权和为 O(1)，组数为 O(\log(2X))。\square

构造 9.5。 预先列出不超过 \sqrt N 的素数。递归状态为 (d,j,z)，表示已经选择素数不超过 p_j 的平方富数 d 及其权重 z=h(d)。从 (1,0,1) 出发，先输出 (d,z)；对每个 k>j 且 dp_k^2\le N，依次取所有 e\ge2、dp_k^e\le N，进入

证明。 每条路径的素因子严格递增且指数至少为 2，故输出均为平方富数。反之，将任一平方富数按递增素因子及其唯一指数分解，即得到唯一生成路径。父状态与新增素数幂互素，故积性保证权重正确。每条边产生一个新节点；各节点枚举素数时只在末尾多作常数次边界检查，故在局部值能常数时间取得的条件下，枚举成本与输出数同阶。\square

实现中用 p\le\lfloor\sqrt{N/d}\rfloor 判断是否进入分支，用除法判断能否再乘 p，以免先发生整数溢出。即使 z=0，正确性也不依赖跳过该分支；若系数环中的零可常数时间判定，则可以剪去其全部后代。

§3　一个前缀和

定理 9.6。 在定理 9.2 的条件下，记 G(x)=\sum_{n\le x}g(n)，则

若每个所需的 G 值与局部 h(p^e) 均能常数时间取得，则已知小素数后可用 O(\sqrt N) 次基本操作求出 F(N)。

证明。 对 f=g*h 交换有限求和，得 F(N)=\sum_{d\le N}h(d)G(\lfloor N/d\rfloor)。由定理 9.2 删除支撑以外的项，再用引理 9.4 与构造 9.5。\square

注。 辅助函数 g 的选择包括两个部分：用 g(p)=f(p) 获得稀疏支撑，并确定其前缀和 G 的计算方法。完全积性的选择可用定理 8.9 恢复前缀和；已有简单卷积关系时可用第七章的方法。

例 9.7。 固定整数 c\ge0，设积性函数 q 满足

取 g=\mathbf1。由式 (9.2)，h(p^{2j+1})=0，而

定义积性函数 J_c 的局部值为 J_c(p^j)=p^{cj}-p^{c(j-1)}，并令 J_c(1)=1。则

证明。 h 只可能在平方数处非零，且 h(a^2)=J_c(a)。将这一事实与 G(x)=x 代入式 (9.6)。\square

构造 9.9。 设 \mathcal Q_N 已按递增顺序编号为 v_1<\cdots<v_L，并已知 G(v)。建立长度 L+1 的零数组 \Delta。枚举每个 (d,h(d))，令 M=\lfloor N/d\rfloor；再枚举 q\in\mathcal Q_M，置

执行 \Delta_a\leftarrow\Delta_a+c，\Delta_{b+1}\leftarrow\Delta_{b+1}-c。最后取一次普通前缀和，令

式 (7.4) 给出常数时间的排名计算；式 (7.2) 给出每个 \mathcal Q_M 的无重复枚举方法。q\in\mathcal Q_N 由命题 7.2 保证，故所有 G(q) 都可查询。

定理 9.10。 构造 9.9 正确恢复 F 的整个商值表示。在局部权重与查询均为常数代价的条件下，除构造 G 外，时间为 O(\sqrt N\log(2N))，空间为 O(\sqrt N)。

证明。 对固定 d，式 (9.10) 将所有 v\ge d 按相同商分成互不相交的连续编号段。差分更新使每个这样的 v 恰增加 h(d)G(\lfloor v/d\rfloor)。再对全部 d 求和，由定理 9.6 得到正确性。区间修改次数为

最后一次前缀累加和 PN 枚举均为 O(\sqrt N)。只需保留商值表、差分数组、小素数表和递归栈。\square

逐个商值独立使用定理 9.6，只能从式 (7.6) 得出 O(N^{3/4}) 的一般界；式 (9.13) 的改进来自对相同贡献作区间累加。

§5　一套完整的恢复程序

定理 9.11。 设 f 积性，且对素数有 f(p)=\sum_{k=0}^{D}c_kp^k。假定 D 固定、所需素数幂局部值可常数时间取得、固定次幂和可常数时间求值。则 f 的商值表示可在

次基本操作内求出，空间为 O((D+1)\sqrt N)。

证明。 依次执行以下有限构造。

1. 对 k=0,\ldots,D 用构造 8.3 求素数幂和，线性组合成 P_f(v)。
2. 定义完全积性的 g，使 g(p)=f(p)；于是 P_g=P_f。用定理 8.9 得到 G(v)。
3. 以式 (9.2) 定义 h，用构造 9.5 枚举其平方富数支撑。
4. 用构造 9.9 得到全部 F(v)。

各阶段正确性分别由上述命题保证；相加命题 8.5 与定理 9.10 的代价即得。\square

因此，素数处的多项式条件只用于建立素数权重表；高次素数幂的偏离由 h 单独承担。这解释了叶筛、逆过程与 PN 修正为什么能够组合。

第十章　整除矩阵与双变量求和

§1　整除关系的矩阵

定义 10.1。 固定 n\ge1，令 E 为 n\times n 整数矩阵，

它是对角线全为 1 的下三角矩阵。对算术函数 a，置 b=\mathbf1*a。

命题 10.2。 有

特别地，其行列式为 \prod_{d=1}^{n}a(d)。

证明。 右端 (i,j) 元素为 \sum_{d\mid i,\ d\mid j}a(d)=b(\gcd(i,j))。又因 \det E=1，对矩阵乘积取行列式即得。\square

命题 10.3。 E 在整数上可逆，且

证明。 两矩阵乘积的 (i,j) 元素为 [j\mid i]\sum_{d\mid i/j}\mu(d)=[i=j]。\square

因而对整除矩阵的消元就是 Möbius 反演的有限线性表示。式 (10.3) 直接给出整数消元变换，可在系数环中逐项使用。

§2　带整数权重的矩阵

例 10.4。 令 A_{ij}=ij\gcd(i,j)。定义整数下三角矩阵

则

证明。 乘积的 (i,j) 元素为

这个分解逐项都在整数中成立。左乘 L^{-1} 后得到上三角矩阵

正是讲义 P6384 所用的消元结果。若仅用基本行变换实现，则按 i 递增，对 k=2,\ldots,\lfloor n/i\rfloor，从第 ki 行减去第 i 行的 k 倍；倍数 k 与原矩阵的行权重 i 对应。讲义先除去行列权重再消元的办法在有理数域中同样成立；式 (10.4) 本身是整数恒等式，可直接映入模环。

§3　最大公约数与最小公倍数的乘法关系

定理 10.5。 设 f(1)=1。函数 f 积性，当且仅当对所有正整数 a,b 有

证明。 若 f 积性，将 a,b 按素因子分解。对每个素数 p，两侧局部因子分别为 f(p^u)f(p^v) 与 f(p^{\min(u,v)})f(p^{\max(u,v)})，二者相等。反之，对互素的 a,b，右端等于 f(1)f(ab)=f(ab)，即得积性。\square

这个结论不要求 f 完全积性，也不要求其值非零。

例 10.6。 令 \omega(n) 为 n 的不同素因子个数，\omega(1)=0。对 t\in\{0,1\}，置

则讲义 P12487 的双重和为

证明。 f_t 积性，故由定理 10.5 将每一项改写为 f_t(i)f_t(j)，再分离两个有限求和。\square

§4　平方自由因子的消去

引理 10.7。 对每个正整数 n，

证明。 第一式右端计数平方自由约数，每个不同素因子独立选择出现或不出现。第二式两端均积性；在 p^e 处，e=0,1 时右端为 1，e\ge2 时为 1+\mu(p)=0，恰与左端一致。\square

定义 10.8。 记

在本文所需情形中，S_0(x)=x，S_1(x)=x(x+1)/2。

命题 10.9。 有

证明。 依次代入引理 10.7 的两个等式，并置 d=k^2a，得到

内层按 a 枚举的就是 D_t(\lfloor N/k^2\rfloor)。全部求和有限，交换次序无需收敛条件。\square

构造 10.10。 令 r=\lfloor\sqrt x\rfloor，则

对 t=0,1，该式给出 O(\sqrt x) 次基本操作的精确算法。

证明。 对称应用命题 7.4，其中两个函数都取 a\mapsto a^t。\square

先线性筛出 \mu(k)，再逐项使用式 (10.11)，可在 O(\sqrt N\log(2N)) 次基本操作内计算式 (10.10)，因为

这给出用于核验高级算法的完整基准实现。对讲义的最大参数，还应改进内层格点求和；下一章展开这一几何步骤。

第十一章　双曲线的整点边界与带权求和

§1　由逐行枚举到有理方向

本章补足原讲义最后一例所略去的几何算法。固定 N\ge1，并记

求和中的 x,y 均为正整数。上一章的双曲线法逐行求取 \lfloor N/y\rfloor。本章的另一种表示把一批连续的行合成一个矩形与一个有理斜边三角形，再一次计算其权重。

采用原始论文中的有理方向细分思想，以及严格位于双曲线上方的整点折线表示。这里给出全部更新式和一份可执行程序；所用方法的来源为朱震霆《一些特殊的数论函数求和问题》第 4 节，以及 Richard Sladkey 的原始论文。

定义 11.1。 点 (X,Y)\in\mathbf Z_{>0}^2 称为外点，若 XY>N。满足 XY\le N 的正整点称为内点；双曲线上的整点属于内点。一个向右下方的本原方向记为

存储时使用非负二元组 (p,q)；向点施加该方向时才将第二个分量取负。

定义 11.2。 两个方向 L=(p,q)、R=(r,s) 称为相邻方向，若

它们满足 q/p>s/r。其中间方向定义为 M=L+R=(p+r,q+s)。

命题 11.3。 式 (11.3) 蕴含 M 本原，且 (L,M)、(M,R) 均为相邻方向。若本原方向 (a,b) 的斜率绝对值严格位于 s/r 与 q/p 之间，则存在唯一正整数 u,v，使

特别地，a\ge p+r、b\ge q+s。

证明。 两个新行列式均等于原行列式。若一个整数整除 p+r 和 q+s，则它整除 q(p+r)-p(q+s)=1。由行列式为 -1，两个方向构成 \mathbf Z^2 的一组基；解线性方程得到的坐标为整数。严格的斜率夹持条件使两个坐标均为正。\square

反复插入中间方向即得到 Stern–Brocot 有理数搜索，各次比较直接作用于整数方向。

§2　三角形的四个矩量

定义 11.4。 对本原方向 (p,q)，定义闭三角形中的整点集

当 (p,q)=(1,0) 时，该集合退化为两点 \{(0,0),(1,0)\}。定义四个矩量

一个方向的完整记录为 (p,q,C,U,V,W)。初始记录为

引理 11.5。 设 L=(p,q)、R=(r,s) 满足式 (11.3)。\mathcal T(p+r,q+s) 的整点集可由以下三部分合成：\mathcal T(p,q)；将 \mathcal T(r,s) 平移 (p,-q) 所得集合；矩形

前两部分恰在 (p,-q) 重复一次；矩形与它们不交。

证明。 两条小斜边与大斜边之间的三角形面积为 1/2。由于两方向的行列式为 \pm1，整数可逆的线性变换将该小三角形变为顶点为 (0,0),(1,0),(1,1) 的三角形。后者除三个顶点外没有整点。因此将两条小斜边替换为大斜边不增加任何整点。其余部分由逐行或逐列比较边界即得。\square

命题 11.6。 记左右记录的矩量分别为 (C_L,U_L,V_L,W_L)、(C_R,U_R,V_R,W_R)，并置

中间方向的矩量为

证明。 在平移后的右三角形上，四个被求和的函数分别变为

矩形的四个矩量分别为 ps,sA,pB,AB。最后减去重复点 (p,-q) 的四个贡献 1,p,-q,-pq。\square

该公式也适用于 R=(1,0)，因为定义 11.4 给退化三角形赋予了两点集合。

注。 为什么带权问题需要一个额外量，可以从逐列求和看出。若另外定义

则由 \gcd(p,q)=1，

这是因为固定 1\le u<p 后，纵坐标从 -q 取到 -\lfloor qu/p\rfloor-1，其和为

同一条边界可同时给出 D_0 和 D_1，每条边的贡献由常数维矩量记录求得。

§4　边界搜索、栈与停止条件

构造 11.8。 以

为起点。该点为外点，因为写 N=km+r，其中 k=\lfloor N/m\rfloor\ge m、0\le r<m，则 k(m+1)-N=k-r>0。

栈从底至顶初始为方向 (1,0),(1,1)，各附式 (11.7) 的矩量。当前累计量为 (Z_0,Z_1)=(0,0)。重复以下过程。

1. 弹出栈顶方向 L。只要 P+L 仍为外点，就把式 (11.14)、(11.15) 加入累计量，并令 P\leftarrow P+L；重复直到下一候选的坐标乘积不大于 N。
2. 若当前 Y^3\le N，结束几何过程。
3. 令 R 为当前栈顶。只要 P+R 的坐标乘积不大于 N，就令 L\leftarrow R 并弹出它，再取新的栈顶为 R。
4. 置 M=L+R，按式 (11.10) 计算记录。若 P+M 为外点，置 R\leftarrow M，将它压栈，继续本步骤。
5. 若 P+M 的坐标乘积不大于 N，检验

成立时结束此次有理方向搜索并回到步骤 1；否则置 L\leftarrow M，回到步骤 4。

这里 P+(p,q) 始终表示 (X+p,Y-q)。所有几何判断均在整数上完成。搜索候选允许 Y-q\le0；此时它也进入乘积不大于 N 的分支，但不称为定义 11.1 的内点。实际前进只接受乘积大于 N 的候选，而 X+p>0，所以真正输出的边端点始终具有正纵坐标。

引理 11.9。 式 (11.17) 是一个可证明的排除条件。若 P+M 的坐标乘积不大于 N，且式 (11.17) 成立，则从 P+M 出发，沿 R 或比 R 更陡的右下方向前进，所得正坐标点均不可能成为外点。起点纵坐标可以非正。

证明。 置 x_0=X+p_M，f(x)=N/x。式 (11.17) 等价于

函数 f 为凸函数，故在 x\ge x_0 上位于它在 x_0 的切线上方。起点不高于 f(x_0)，而所考虑射线的斜率不大于切线斜率，因此整条射线均不高于曲线。\square

命题 11.10。 构造 11.8 输出的各条边均满足命题 11.7 的条件，并且处理的行区间两两不交。结束时，已经处理的行恰为 Y,Y+1,\ldots,m。

证明。 先给出搜索中反复使用的凸性事实。若 P=(X,Y) 为外点，而 Q=(X+a,Y-b) 满足 a>0、坐标乘积不大于 N，则

第一个不等式使用 Y>N/X 与 Y-b\le N/(X+a)，第二个使用凸性。因此从 Q 沿 (a,-b) 或更陡的右下方向前进，始终不高于双曲线。这一结论也允许 Q 的纵坐标非正。

搜索的不变量如下：L,R 相邻，P+L 的坐标乘积不大于 N，P+R 为外点；从 P 看去，不比 L 平缓的全部正整点均为内点。若中间候选 P+M 的乘积不大于 N，则 L 与 M 之间的任一整数方向均可写为 uL+vM，其中 u,v\ge1。从 P+M 出发，剩余方向 uL+(v-1)M 不比 M 平缓，由前述凸性事实可知该点仍在曲线下。M 射线上的后续点也一样。因此令 L\leftarrow M 时，新扩大的整个陡方向区域都已排除。步骤 3 弹出一个乘积不大于 N 的右候选时，应用同一论证；若 P+M 为外点，则仅将 R 换为 M，旧的较平缓方向保存在栈中。不变量由此保持。

当式 (11.17) 成立时，任一严格位于 R 射线下方的整数方向，因整数基性质可写为 \alpha L+\beta R，其中 \alpha\ge1。若 \beta\le0，它属于已排除的陡方向区域；若 \beta\ge1，则它等于

由引理 11.9，从 P+M 出发沿剩余的、不比 R 平缓的方向前进不可能遇到外点。于是 R 整条右射线严格下方均无外整点。R 本原，故其一个单位线段内没有其他边上整点，平移三角形除两个端点外全为内点。

还需初始化并检查重复前进。初始 P=(k,m+1)，其中 k=\lfloor N/m\rfloor\ge m。初始斜率 -1 射线严格下方的整数点写成 (k+u,y)，其中 u\ge0、y\le m-u。只要坐标为正，

故初始边满足所需性质。沿已确定的方向重复前进，严格在线下的右侧整点仍属于原先排除的集合。一旦下一射线点的乘积不大于 N，前述凸性事实又排除射线上全部后续点，从而恢复下一轮搜索的完整不变量。

对每一条实际输出的边，命题 11.7 已证明三角形与左侧矩形恰计入全部相关内点，右端外点排除了右侧遗漏。每条边把行标从 Y 改为 Y-q，并只计入 [Y-q,Y-1]。自初始行标 m+1 归纳，即得行区间的不交性与覆盖性。\square

栈的边界条件。 栈底的水平哨兵 (1,0) 永远不会被步骤 3 弹出：当前 XY>N，所以 (X+1)Y>N。当 R 是哨兵时，式 (11.17) 为 0\ge N，不可能成立。因此一次搜索只有在压入至少一个非水平方向之后才可能结束；回到步骤 1 时，可弹出的栈顶总位于哨兵之上。这说明两处取栈顶操作不会访问空栈。

终止性。 若 R 为水平哨兵，相邻条件强制 q_L=1。未触发截止条件时 Y\ge2，反复令 L\leftarrow L+R 保持纵差为 1 而增加横差，因此最终会产生乘积大于 N 的候选，压入一个非水平方向。此后，每次外候选使 R\leftarrow L+R，其纵差严格增加，且外点条件要求该纵差小于 Y，故这种更新只能发生有限次。固定非水平 R 时，连续的内候选使横差严格增加，最终满足式 (11.17)。所以每次方向搜索有限。搜索结束后的栈顶是外方向，下一轮至少成功前进一步，使正整数 Y 至少减一；全部几何过程因而也有限。

注。 栈保留此前搜索得到的相邻方向，供下一次搜索继续使用；式 (11.17) 决定每次搜索的停止位置。

定理 11.11。 几何过程结束后，补上剩余行：

则

证明。 几何过程与剩余行合起来恰计入 y\le m 的全部内点。交换横纵坐标得到 x\le m 的部分；两部分覆盖 xy\le N，交为 [1,m]^2。两个权重分别在该正方形上求和即得所减项。\square

§5　可执行的整数构造

以下 Python 程序使用任意精度整数实现上述过程。输入 N\ge1，输出精确整数 (D_0(N),D_1(N))。下一节按基本算术操作次数分析其复杂度。

from math import isqrt

def merge(L, R):
    p, q, C, U, V, W = L
    r, s, D, E, F, G = R
    assert q * r - p * s == 1
    A = p * (p - 1) // 2
    B = -s * (2 * q + s + 1) // 2
    return (p + r, q + s,
            C + D + p * s - 1,
            U + E + p * D + s * A - p,
            V + F - q * D + p * B + q,
            W + G + p * F - q * E - p * q * D + A * B + p * q)

def edge_weight(X, Y, R):
    p, q, C, U, V, W = R
    rectangle = X * (X - 1) // 2 * (q * (2 * Y - q - 1) // 2)
    return ((X - 1) * q + C - 2,
            rectangle + W + X * V + Y * U + X * Y * C
            - X * Y - (X + p) * (Y - q))

def hyperbola_moments(N):
    m = isqrt(N)
    X, Y = N // m, m + 1
    stack = [(1, 0, 2, 1, 0, 0), (1, 1, 3, 1, -2, -1)]
    Z0 = Z1 = 0

    def outside(R):
        return (X + R[0]) * (Y - R[1]) > N

    while True:
        L = stack.pop()
        while outside(L):
            E0, E1 = edge_weight(X, Y, L)
            Z0 += E0
            Z1 += E1
            X += L[0]
            Y -= L[1]
        if Y ** 3 <= N:
            break

        R = stack[-1]
        while not outside(R):
            L = R
            stack.pop()
            R = stack[-1]

        while True:
            M = merge(L, R)
            if outside(M):
                R = M
                stack.append(R)
            elif (X + M[0]) ** 2 * R[1] >= N * R[0]:
                break
            else:
                L = M

    for y in range(1, Y):
        q = N // y
        Z0 += q
        Z1 += y * q * (q + 1) // 2

    return 2 * Z0 - m * m, 2 * Z1 - (m * (m + 1) // 2) ** 2

例 11.12。 对 N=10，程序返回 (27,170)。直接按乘积分类也得到

这里同时求两个矩量只把每个方向记录的大小增加为常数，不改变几何分支。

§6　操作数与区域计数

设 B=\lfloor N^{1/3}\rfloor。构造在 Y\le B 时转为直接枚举，最后的行处理需 O(B) 次运算。几何部分由沿边前进、方向搜索和出入栈组成。

记 E(N) 为步骤 1 前进的次数，J(N) 为步骤 4 形成中间方向的次数。程序的总操作数为

每次前进、合并均只做常数次运算。每次非初始压栈对应一次中间方向合并；一个记录最多弹出一次，因此所有出入栈也被 J(N)+2 控制。同一方向上的每次前进分别计入 E(N)。每步使 Y 至少减一，故 E(N)\le m。

下面估计切线区域算法中有理方向细分产生的区域数量。

对有理方向区域细分，可把非平凡的处理归属于相邻分数 a/b>c/d>0，其中 ad-bc=1。对应的两个切点横坐标为 \sqrt{Nb/a} 和 \sqrt{Nd/c}。一个同时在两条整数基方向上具有非零宽度的区域必须满足

这是因为两个单位基方向的横坐标跨度为 b+d。

引理 11.13。 满足式 (11.21) 的正整数四元组数量为

证明。 置 z=bc。由 ad=z+1，

所以 4abc^2(b+d)^2\le N。代入 c=z/b、d=(z+1)/a 后，左侧至少为 8z^3，从而 z=O(N^{1/3})。给定 z，选择 b\mid z、a\mid z+1 就唯一确定 c,d，故所求数量至多

设 d_4(n) 为 n 写成四个正整数有序乘积的方式数。逐质数幂比较

引理 11.13 给出相邻分数区域的数量界。本程序的时间与空间分别由式 (11.20) 和 O(J(N)+1) 描述。维护四个矩量只对同一几何过程增加常数因子。

来源。 朱震霆原论文第 108—109 印刷页给出区域方法的 \log^3 N 计数，第 110 页讨论该方法与折线方法的几何比较。

当答案要求取模时，可在 C,U,V,W,Z_0,Z_1 的加乘中约化；坐标、大小比较、整除判断和式 (11.17) 始终在整数上进行。式中的二分之一均来自整数求和公式，应先在整数中约去，或在 2 可逆的模环中使用其逆元。

§7　接回最终的 Möbius 求和

上一章已经得到

本章使其中的 D_0,D_1 成为明确可执行的子程序。实际实现应缓存相同参数 \lfloor N/k^2\rfloor 的结果；对小参数可直接预处理 n^t\tau(n) 的前缀和，省去几何搜索。

逐项调用本章程序时，总成本为

其中 T 由式 (11.20) 给出；开头的 O(\sqrt N) 包括 Möbius 筛。使用上一章的普通双曲线子程序时，总界为 O(\sqrt N\log(N+1))。

推论 11.14。 若精确双曲线子程序的单次算术操作数为 O(x^{1/3}\log^3(x+1))，则式 (11.24) 的逐项求值总成本为 O(\sqrt N)，包括线性筛 Möbius 函数的时间。

证明。 将 k 按 \sqrt N/2^{j+1}<k\le\sqrt N/2^j 分组。该组至多有 O(\sqrt N/2^j) 个元素，且每个参数至多 2^{2j+2}。假定的子程序界在该组内为 O(2^{2j/3}(j+1)^3)，故各组总代价至多

线性筛 Möbius 函数亦满足该时间界。\square

编校附记

本稿重构对象为苏紫苒（minstdfx）《琪露诺的完美数论教室：夏季学期算法竞赛课程》，文件版本为 20260619-2 ver.，共 128 页。篇章和证明体例参照《数据结构入门补充材料_布尔巴基式重构.md》。原讲义的口语说明、图示及代码，改写为明确输入条件、状态、变换和正确性命题；按数学依赖重新安排顺序。为展开原稿略述的算法，另补入局部递推、素数筛的逆过程、稀疏卷积的区间累加，以及格点求和的几何构造。

内容对应

以下页码均为 PDF 的实际页码。

校订与补充

1. 最大公约数用整数理想的非负生成元定义，以容纳全零族；最小公倍数在所用范围内取最小正公倍数。素数、合数与单位分开，不把 0,\pm1 归为合数。因子配对限定于正整数。
2. 同余式约去公因子时同步改变模数。ExCRT 合并维护完整解类；GCD Table 以最大公约数等式检验候选解，并分别处理“起点有界”与“整段位于表内”两种边界。
3. Euler 的非互素指数约简按完整素数幂分量证明，补上指数达到阈值的条件；指数零、模数一及不可逆分母分别处理。ex-Lucas 明确区分去除全部 p 因子的阶乘与只删去 p 倍数的乘积。
4. exBSGS 保留逐轮乘入 a/g 的累计系数，并先检查当前最小指数。单次剥离、模数变化、最小解与终止性以同一不变量证明。
5. 离散对数所在指数环的模数是群的阶；CF1106F 中为 p-1，其线性同余按一般模数求解。素性测试分别给出单侧正确性、随机误差定理与有限范围底数集。Pollard–Rho 在 1<d<n 时返回因子，失败批次回退或更换参数。
6. 第 65、67 页的常值一函数统一写作 \mathbf1，恒等函数写作 \operatorname{id}。约数和的局部因子改为几何和；\prod_p(\nu_p(n)+1) 是约数个数。\varphi 的计数域补为有限完整剩余系，\mu(1)=1 单独明确。
7. 反演适用于一般算术函数；相应的加权恒等式使用完全积性。倍数求和采用有限截断；生成函数以形式级数建立，实级数的收敛性在使用处证明。
8. 第 74 页两重前缀和按其原定义应满足 F_k(x)=T_k(2x)-2T_k(x)。正文按有序对之和的各个系数证明此式，补上原稿缺少的系数 2；在模 2^{32} 下使用不含除法的预处理。
9. 区间素数判定与区间完全分解采用不同的起点。Haitang 实例中候选条件为 S(m)=s、m>s 与 m\mid n-s；需要分解整个短区间，检查成本由该区间全部约数数目控制。
10. 第 85 页的线性筛分支改为“新乘素数 p 是否整除旧下标 i”，结果更新到 ip。一般系数环上的计算通过互素余因子与局部数据构造，复杂度计入素数幂数据的取得代价。
11. 给定两张商值表，一般卷积的直接双曲线求值需 O(N^{3/4}) 次操作；若较小前缀可线性预处理，截断方法给出 O(N^{2/3})。杜教筛使用首项可逆的卷积因子和可快速求值的前缀。
12. 素数权重筛按素数递增、商值递减更新；恢复完全积性前缀时，通过三角变换的逆过程将两种次序反向。Min_25 的递归代价按节点数计算。
13. PN 例题实际恢复的辅助函数为 h；原稿第 114 页末尾将其误写为 g。完全积性选择下的局部差分用望远镜恒等式证明。对相同贡献作排名区间累加，可在 O(\sqrt N\log N) 时间内恢复全体商值。
14. 带权最大公约数矩阵通过整数单位下三角矩阵分解，适用于合数模数。消元的行倍数、对角线及行列式均由同一分解导出。
15. 最后一题的双重和先化为单一前缀和的平方，再作平方自由反演。第十一章展开内层格点算法，并计算外层 Möbius 求和的总成本。

原讲义题单

原 PDF 第 125 页提供以下题单，可用于继续练习：

• gza 的趣味数论 1
• gza 的趣味数论 2
• Leasier 的数论题单
• CF2081G2。