∑MathIntermediate

Multiplicative Functions

Key Points

•
A multiplicative function is an arithmetic function f with f(mn) = f(m)f(n) whenever gcd(m, n) = 1.
•
Classic examples are Euler’s totient $φ,$ Möbius $μ,$ divisor count d (also τ), and divisor sum $σ,$ all computable from prime powers.
•
Dirichlet convolution (f * g)(n) = $\sum_{d ∣ n}$ f(d) g $\frac{n}{d}$ turns arithmetic functions into an algebra with identity ε where $ε (1)$ = 1 and $ε (n > 1)$ = 0.
•
Key identity: μ * 1 = $ε,$ which underlies Möbius inversion and lets you recover f from F(n) = $\sum_{d ∣ n}$ f(d).
•
For multiplicative f, knowing f( $p^{a}$ ) for prime powers determines f(n) for all n via factorization.
•
Efficient computation up to N uses sieves: SPF (smallest prime factor) sieve for O(N) factorization or a linear sieve to compute μ and φ in O(N).
•
Dirichlet convolution over all $n \leq N$ can be computed in O(N log N) by iterating over divisors with nested loops.
•
In competitive programming, store values in 64-bit or modulo arithmetic, handle n = 1 carefully, and reuse precomputed arrays.

Prerequisites

→Prime factorization and the Fundamental Theorem of Arithmetic — Multiplicative functions decompose across distinct primes; understanding factorization is essential.
→Greatest common divisor (gcd) and coprimality — The definition of multiplicativity requires gcd(m, n) = 1.
→Divisors and divisor enumeration — Dirichlet convolution sums over divisors; efficient loops rely on divisor properties.
→Big-O complexity and sieve techniques — Efficient precomputation (SPF, linear sieve) is key to solving large-constraint problems.
→Modular arithmetic and 64-bit integers — Many problems require modulo operations or careful handling to avoid overflow.
→Discrete convolution patterns — Dirichlet convolution resembles other convolution patterns and benefits from similar loop optimizations.

Detailed Explanation

Tap terms for definitions

01Overview

Multiplicative functions are a central class of arithmetic functions in number theory. Informally, they behave nicely with respect to multiplying numbers that share no common prime factors. If two numbers m and n are coprime, a multiplicative function f satisfies f(mn) = f(m)f(n). This property means that once you know f on prime powers p^a, you can determine f(n) for any n by factoring n and multiplying contributions. Famous examples include Euler’s totient φ(n) (the count of integers up to n that are coprime to n), the Möbius function μ(n) (which detects square-free structure and parity of prime factors), the divisor count d(n), and the divisor sum σ(n). These functions exhibit elegant product formulas on prime powers and tie deeply into Dirichlet convolution, an operation that combines arithmetic functions by summing over divisors. The convolution turns the set of arithmetic functions into a commutative ring with identity ε, enabling algebraic manipulation. A cornerstone identity is μ * 1 = ε, which powers Möbius inversion: given F(n) = \sum_{d\mid n} f(d), you can recover f(n) = \sum_{d\mid n} μ(d) F(n/d). Computationally, multiplicative functions can be precomputed up to a bound N with sieve techniques, enabling fast answers to many divisor-sum and coprime-counting problems common in programming contests.

02Intuition & Analogies

Think of building a LEGO model (your number n) from prime bricks. If two LEGO sections (numbers m and n) don’t physically connect (they are coprime), then a multiplicative function f assesses the whole model by independently assessing each section and multiplying the scores: f(mn) = f(m) × f(n). This independence mirrors how prime factorizations split a number into non-overlapping parts. For instance, counting divisors d(n) is like counting ways to pick heights for each stack of identical bricks: for p^a you have a+1 choices (0 through a bricks). The total choices multiply across independent stacks, giving d(n) = ∏(a_i + 1). The Möbius function μ(n) is a yes/no parity checker that says 0 if any stack has height ≥ 2 (the model is not square-free), and otherwise flips sign based on how many stacks you used. Dirichlet convolution is like assembling a project using subassemblies: to compute (f * g)(n), you consider every way to split n into a divisor d and a complementary part n/d; you score the split as f(d)g(n/d), then sum over all splits. The identity ε plays the role of a neutral subassembly: only the trivial split at n = 1 contributes. The identity μ * 1 = ε states that μ ‘undoes’ the effect of summing over divisors; it’s like having an instruction that reverses the “add all subassemblies” operation. Because multiplicative functions reduce to prime powers, coding them efficiently often means: precompute smallest prime factors, factor n quickly, and apply the prime-power formula—much like sorting your LEGO bricks by type for rapid builds.

03Formal Definition

An arithmetic function is any function f:

N

\to

C

. It is called multiplicative if f(1)

\neq =

0 (often f(1) = 1 by convention) and for all m,n

\in

N

with

g cd

(m,n) = 1, we have f(mn) = f(m)f(n). It is completely multiplicative if the same equality holds for all m,n without the coprimality condition. Every multiplicative function is determined by its values on prime powers: if n =

\prod_{i = 1}^{k}

i^{a_{i}}

, then f(n) =

\prod_{i = 1}^{k}

f(p_

i^{a_{i}}

). Dirichlet convolution of functions f and g is defined by (f * g)(n) =

\sum_{d ∣ n}

f(d)g

\frac{n}{d}

. The function ε (epsilon) is the identity for convolution:

ε (1)

= 1 and

ε (n)

= 0 for

n > 1

. The constant-1 function is denoted 1(n) = 1 for all n. The Möbius function μ is defined by

μ (1)

= 1,

μ (n)

= 0 if n is divisible by a square, and

μ (n)

= (-1)^{k} if n is a product of k distinct primes. Key identities include μ * 1 =

ε,

φ * 1 =

id

, and hence φ = μ *

id

, where

id

(n) = n. For multiplicative f, the Dirichlet series

\sum_{n \geq 1}

f(n)

n^{- s}

factorizes into an Euler product

\prod_{p}

\sum_{a \geq 0}

p^{a}

)

p^{- a s}

, reflecting independence across primes.

04When to Use

Use multiplicative-function techniques when a problem involves sums over divisors, properties depending only on prime-power exponents, or identities expressible via Dirichlet convolution. Common scenarios: computing arrays of φ(n), μ(n), d(n), or σ(n) up to N for many queries; evaluating \sum_{d\mid n} f(d) or inverting such sums via Möbius inversion; counting coprime pairs, reduced fractions, or primitive structures (φ); detecting square-free constraints (μ); and computing divisor-related statistics (d, σ). When constraints are large (N up to 10^6–10^7), precompute with sieves: a linear sieve builds μ and φ in O(N), while an SPF sieve allows O(N) factorization passes giving O(N log N) total for computing general multiplicative f. When asked to compute convolution-like sums for all n ≤ N, use the divisor-iteration pattern with nested loops for O(N log N). If a problem gives F(n) = \sum_{d\mid n} f(d) and asks for f, directly apply Möbius inversion with a precomputed μ array. In analytical or generating-function contexts, use Euler products and Dirichlet series factorization to prove identities or derive average orders.

⚠️Common Mistakes

• Confusing multiplicative with completely multiplicative: f(mn) = f(m)f(n) only when gcd(m, n) = 1 unless stated “completely multiplicative.” For example, φ is multiplicative but not completely multiplicative since φ(p^2) ≠ φ(p)^2. • Forgetting n = 1 base cases: set f(1) correctly (typically 1). The identity ε and inverses (like μ) rely on correct f(1). • Overflow in σ(n) or products: σ(p^a) can exceed 32-bit; use 64-bit (long long) or modular arithmetic per problem statement. • Inefficient per-query factorization: factoring each n with trial division is too slow; precompute SPF for O(1) next-prime steps (overall O(log n) per n) or use a linear sieve when specialized (μ, φ). • Wrong Möbius implementation: μ(n) = 0 if any exponent ≥ 2, not only if n is divisible by p^2 for some specific p you checked; you must factor fully or track exponents. • Misusing convolution loops: (f * g)(n) sums over divisors of n, not over all k ≤ n. The O(N log N) all-n computation is achieved by looping multiples of each d. • Assuming all multiplicative functions share the same prime-power formula: each function has its own p^a formula (e.g., φ(p^a) = p^{a-1}(p-1); σ_k(p^a) = (p^{k(a+1)} - 1)/(p^k - 1)). • Ignoring memory: arrays of size N up to 10^7 need careful memory and I/O management; prefer iterative loops and reserve vectors.

Key Formulas

Multiplicativity

f (mn) = f (m) f (n) for g cd (m, n) = 1

Explanation: The value on coprime products factors as a product of values. This lets you build f(n) from its values on prime powers.

Dirichlet Convolution

(f * g) (n) = d ∣ n \sum f (d) g (\frac{n}{d})

Explanation: Defines how to combine two arithmetic functions by summing over divisor pairs. It is associative, commutative, and has identity $ε .$

Convolution Identity

ε (n) = [n = 1] = {1, 0, n = 1 n > 1

Explanation: This function acts as the identity for Dirichlet convolution. Convolving any function with ε yields the original function.

Möbius Identity

μ * 1 = ε

Explanation: The Möbius function is the Dirichlet inverse of the constant-1 function. It underlies Möbius inversion.

Möbius Inversion

F (n) = d ∣ n \sum f (d) \Rightarrow f (n) = d ∣ n \sum μ (d) F (\frac{n}{d})

Explanation: If F is the divisor-sum transform of f, then f can be recovered using $μ .$ This is essential for inverting many divisor-sum relations.

Prime-Power Formulas

ϕ (p^{a}) = p^{a - 1} (p - 1), σ (p^{a}) = \frac{p ^{a + 1} - 1}{p - 1}, d (p^{a}) = a + 1

Explanation: Closed forms on prime powers for classic multiplicative functions. Combined multiplicatively they yield values for all n.

Totient Identity

d ∣ n \sum ϕ (d) = n \Leftrightarrow (ϕ * 1) (n) = id (n)

Explanation: The sum of φ over all divisors of n equals n. Equivalently, φ convolved with 1 equals the identity-by-value function.

Totient via Möbius

ϕ = μ * id

Explanation: Totient can be recovered as the convolution of μ with the identity function id(n) = n. Useful for proofs and computations.

Euler Product for Multiplicative f

n \geq 1 \sum \frac{f ( n )}{n ^{s}} = p \prod (a \geq 0 \sum \frac{f ( p ^{a} )}{p ^{a s}})

Explanation: The Dirichlet series of a multiplicative function factorizes over primes. This encodes independence across prime powers.

Riemann Zeta Euler Product

ζ (s) = n \geq 1 \sum \frac{1}{n ^{s}} = p \prod \frac{1}{1 - p ^{- s}}

Explanation: A special case where f(n)=1. Useful as a reference point for understanding Euler products and divisor sums.

Complexity Analysis

Precomputation dominates most uses. An SPF (smallest prime factor) linear sieve runs in O(N) time and O(N) space, producing both the prime list and spf[] so that factoring any

n \leq N

takes O(number of prime factors) = O(log n) divisions. Computing a generic multiplicative function f for all

n \leq N

by factoring each n via spf[] costs O(N log N) overall in the worst case (harmonic sum of logs). Specialized linear sieves can compute certain functions like φ and μ in true O(N) by using recurrence rules when building composites (no per-n factorization required). For Dirichlet convolution over all 1 ≤

n \leq N

, the standard “divisor iteration” pattern loops over divisors d and their multiples

m = k \cdot d

, updating h[m] += f[d]·g[k]. The double loop runs in Θ(∑_{

d=1}^{N}

⌊N/d⌋) = Θ(N log N) time and O(N) space. The divisor zeta transform F(n) =

\sum_{d ∣ n}

f(d) and its inverse via μ follow the same complexity. Memory usage is O(N) for arrays like

φ,

μ,

σ,

spf, and primes; at

N = 1 0^{7}

this means several hundred MB if many arrays are kept simultaneously, so reuse and free arrays when possible. In modular arithmetic settings, multiplications inside convolution remain O(1), but be careful to reduce modulo often to avoid overflow. In 64-bit sums (e.g.,

σ),

sums may exceed 32-bit even for moderate N, so prefer long long.

Code Examples

Compute φ, μ, d, σ for all n ≤ N using a linear SPF sieve and prime-power formulas

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 // Linear sieve to compute SPF (smallest prime factor) for 1..N in O(N)
5 struct SPF {
6     int N;
7     vector<int> spf, primes;
8     SPF(int n): N(n), spf(n+1, 0) {
9         if (N >= 1) spf[1] = 1;
10         for (int i = 2; i <= N; ++i) {
11             if (spf[i] == 0) { spf[i] = i; primes.push_back(i); }
12             for (int p : primes) {
13                 long long v = 1LL * p * i;
14                 if (v > N) break;
15                 spf[(int)v] = p; // first time set gives smallest prime factor
16                 if (p == spf[i]) break; // ensure each composite is written once per minimal prime
17             }
18         }
19     }
20 
21     // Factorize n into vector of (prime, exponent) using spf[]; O(#prime factors)
22     vector<pair<int,int>> factorize(int n) const {
23         vector<pair<int,int>> fac;
24         if (n == 1) return fac;
25         while (n > 1) {
26             int p = spf[n];
27             int a = 0;
28             while (n % p == 0) { n /= p; ++a; }
29             fac.push_back({p, a});
30         }
31         return fac;
32     }
33 };
34 
35 int main() {
36     ios::sync_with_stdio(false);
37     cin.tie(nullptr);
38 
39     int N = 1000000; // adjust per problem constraints
40     SPF sieve(N);
41 
42     vector<long long> phi(N+1), mu(N+1), d(N+1), sigma(N+1);
43     phi[1] = 1; mu[1] = 1; d[1] = 1; sigma[1] = 1;
44 
45     for (int n = 2; n <= N; ++n) {
46         auto fac = sieve.factorize(n);
47         // φ(n) = n * Π_{p|n} (1 - 1/p)
48         long long ph = n;
49         // μ(n) = 0 if any exponent >= 2; else (-1)^{k}
50         long long mval = 1;
51         // d(n) = Π (a_i + 1)
52         long long divcnt = 1;
53         // σ(n) = Π (p^{a+1}-1)/(p-1)
54         long long divsum = 1;
55         for (auto [p, a] : fac) {
56             ph = ph / p * (p - 1);
57             if (a >= 2) mval = 0; // square factor detected
58             // flip sign per distinct prime if still square-free
59             // done after loop if mval != 0
60             divcnt *= (a + 1);
61             // compute p^{a+1}
62             long long pa1 = 1;
63             for (int i = 0; i < a + 1; ++i) pa1 *= p;
64             divsum *= (pa1 - 1) / (p - 1);
65         }
66         // set sign for μ if square-free
67         if (mval != 0) {
68             int k = (int)fac.size();
69             mval = (k % 2 == 0) ? 1 : -1;
70         }
71         phi[n] = ph;
72         mu[n] = mval;
73         d[n] = divcnt;
74         sigma[n] = divsum;
75     }
76 
77     // quick correctness spot-checks
78     cout << "phi(12) = " << phi[12] << " (expect 4)\n";
79     cout << "mu(12) = " << mu[12] << " (expect 0)\n";
80     cout << "d(12) = " << d[12] << " (expect 6)\n";
81     cout << "sigma(12) = " << sigma[12] << " (expect 28)\n";
82 
83     return 0;
84 }
85

We build an O(N) linear sieve that stores the smallest prime factor for every n. Each n is then factorized in O(number of prime factors). Using prime-power closed forms, we compute φ, μ, d, and σ for all n ≤ N. This approach is flexible: any multiplicative function can be added by specifying its prime-power rule.

Time: O(N) to build SPF, plus O(N log N) to compute functions via factorization across all nSpace: O(N)

Compute Dirichlet convolution h = f * g for all n ≤ N, and verify μ * 1 = ε

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 int main() {
5     ios::sync_with_stdio(false);
6     cin.tie(nullptr);
7 
8     int N = 100000; // tune per constraints
9 
10     // Example: build mu via a linear sieve specialized for mu (optional fast path)
11     vector<int> mu(N+1, 0), primes;
12     vector<int> is_comp(N+1, 0);
13     mu[1] = 1;
14     for (int i = 2; i <= N; ++i) {
15         if (!is_comp[i]) { primes.push_back(i); mu[i] = -1; }
16         for (int p : primes) {
17             long long v = 1LL * i * p; if (v > N) break;
18             is_comp[(int)v] = 1;
19             if (i % p == 0) { mu[(int)v] = 0; break; }
20             else mu[(int)v] = -mu[i];
21         }
22     }
23 
24     vector<long long> f(N+1, 1), g(N+1, 0), h(N+1, 0);
25     // f = 1 (constant-one function), g = mu
26     for (int n = 1; n <= N; ++n) g[n] = mu[n];
27 
28     // Compute h = f * g via divisor-iteration: O(N log N)
29     for (int d = 1; d <= N; ++d) {
30         for (int m = d; m <= N; m += d) {
31             h[m] += f[d] * g[m / d];
32         }
33     }
34 
35     // Verify epsilon: h[1] = 1, and h[n] = 0 for n > 1
36     cout << "(mu * 1)(1) = " << h[1] << " (expect 1)\n";
37     bool ok = true;
38     for (int n = 2; n <= N; ++n) if (h[n] != 0) { ok = false; break; }
39     cout << "All n>1 zero? " << (ok ? "YES" : "NO") << "\n";
40     return 0;
41 }
42

This demonstrates the all-n Dirichlet convolution using the classic nested loop over divisors. With f ≡ 1 and g = μ, the result must be the identity ε: only n = 1 yields 1; all others are 0.

Time: O(N log N)Space: O(N)

Möbius inversion: recover f from F(n) = ∑_{d|n} f(d) using μ

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 int main() {
5     ios::sync_with_stdio(false);
6     cin.tie(nullptr);
7 
8     int N = 200000; // adjust per constraints
9 
10     // Build mu by linear sieve
11     vector<int> mu(N+1, 0), primes; vector<int> is_comp(N+1, 0);
12     mu[1] = 1;
13     for (int i = 2; i <= N; ++i) {
14         if (!is_comp[i]) { primes.push_back(i); mu[i] = -1; }
15         for (int p : primes) {
16             long long v = 1LL * i * p; if (v > N) break;
17             is_comp[(int)v] = 1;
18             if (i % p == 0) { mu[(int)v] = 0; break; }
19             else mu[(int)v] = -mu[i];
20         }
21     }
22 
23     // Suppose the hidden f is: f(n) = 1 if n is odd, else -1 (arbitrary example)
24     vector<long long> f_true(N+1, 0), F(N+1, 0), f_rec(N+1, 0);
25     for (int n = 1; n <= N; ++n) f_true[n] = (n % 2 ? 1 : -1);
26 
27     // Forward divisor zeta transform: F(n) = sum_{d|n} f(d)
28     for (int d = 1; d <= N; ++d) {
29         for (int m = d; m <= N; m += d) {
30             F[m] += f_true[d];
31         }
32     }
33 
34     // Inversion: f = mu * F
35     for (int d = 1; d <= N; ++d) {
36         if (mu[d] == 0) continue;
37         for (int m = d; m <= N; m += d) {
38             f_rec[m] += 1LL * mu[d] * F[m / d];
39         }
40     }
41 
42     // Check accuracy on a few values
43     bool ok = true;
44     for (int n = 1; n <= N; ++n) {
45         if (f_rec[n] != f_true[n]) { ok = false; break; }
46     }
47     cout << "Recovered f equals original? " << (ok ? "YES" : "NO") << "\n";
48     return 0;
49 }
50

We create an arbitrary function f, compute its divisor-sum transform F(n) = ∑_{d|n} f(d), and then recover f using Möbius inversion: f = μ * F. The divisor-iteration pattern computes both forward and inverse transforms in O(N log N).

Time: O(N log N)Space: O(N)

True O(N) linear sieve for φ and μ with a sanity identity check ∑_{d|n} φ(d) = n

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 int main() {
5     ios::sync_with_stdio(false);
6     cin.tie(nullptr);
7 
8     int N = 1000000; // choose as needed; O(N) time and space
9     vector<int> primes; primes.reserve(N/10);
10     vector<int> is_comp(N+1, 0);
11     vector<int> phi(N+1, 0), mu(N+1, 0);
12 
13     phi[1] = 1; mu[1] = 1;
14     for (int i = 2; i <= N; ++i) {
15         if (!is_comp[i]) { primes.push_back(i); phi[i] = i - 1; mu[i] = -1; }
16         for (int p : primes) {
17             long long v = 1LL * i * p; if (v > N) break;
18             is_comp[(int)v] = 1;
19             if (i % p == 0) {
20                 // p divides i: update using recurrence
21                 phi[(int)v] = phi[i] * p;
22                 mu[(int)v] = 0;
23                 break; // ensure linear-time by stopping at first prime dividing i
24             } else {
25                 // p does not divide i
26                 phi[(int)v] = phi[i] * (p - 1);
27                 mu[(int)v] = -mu[i];
28             }
29         }
30     }
31 
32     // Spot-check: sum_{d|n} phi(d) == n for a random sample
33     std::mt19937 rng(712367);
34     for (int t = 0; t < 5; ++t) {
35         int n = (rng() % (N-1)) + 1;
36         long long s = 0;
37         for (int d = 1; d * d <= n; ++d) if (n % d == 0) {
38             s += phi[d];
39             if (d * 1LL * d != n) s += phi[n/d];
40         }
41         cout << "Check n=" << n << ": sum phi(d)=" << s << ", n=" << n
42              << (s == n ? " OK" : " FAIL") << "\n";
43     }
44 
45     return 0;
46 }
47

This is the classic linear sieve to compute φ and μ simultaneously in O(N). The recurrence distinguishes whether the next prime multiplies an i already divisible by that prime. A quick check verifies the identity (φ * 1)(n) = id(n).

Time: O(N)Space: O(N)

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	// Linear sieve to compute SPF (smallest prime factor) for 1..N in O(N)
5	struct SPF {
6	int N;
7	vector<int> spf, primes;
8	SPF(int n): N(n), spf(n+1, 0) {
9	if (N >= 1) spf[1] = 1;
10	for (int i = 2; i <= N; ++i) {
11	if (spf[i] == 0) { spf[i] = i; primes.push_back(i); }
12	for (int p : primes) {
13	long long v = 1LL * p * i;
14	if (v > N) break;
15	spf[(int)v] = p; // first time set gives smallest prime factor
16	if (p == spf[i]) break; // ensure each composite is written once per minimal prime
17	}
18	}
19	}
20
21	// Factorize n into vector of (prime, exponent) using spf[]; O(#prime factors)
22	vector<pair<int,int>> factorize(int n) const {
23	vector<pair<int,int>> fac;
24	if (n == 1) return fac;
25	while (n > 1) {
26	int p = spf[n];
27	int a = 0;
28	while (n % p == 0) { n /= p; ++a; }
29	fac.push_back({p, a});
30	}
31	return fac;
32	}
33	};
34
35	int main() {
36	ios::sync_with_stdio(false);
37	cin.tie(nullptr);
38
39	int N = 1000000; // adjust per problem constraints
40	SPF sieve(N);
41
42	vector<long long> phi(N+1), mu(N+1), d(N+1), sigma(N+1);
43	phi[1] = 1; mu[1] = 1; d[1] = 1; sigma[1] = 1;
44
45	for (int n = 2; n <= N; ++n) {
46	auto fac = sieve.factorize(n);
47	// φ(n) = n * Π_{p\|n} (1 - 1/p)
48	long long ph = n;
49	// μ(n) = 0 if any exponent >= 2; else (-1)^{k}
50	long long mval = 1;
51	// d(n) = Π (a_i + 1)
52	long long divcnt = 1;
53	// σ(n) = Π (p^{a+1}-1)/(p-1)
54	long long divsum = 1;
55	for (auto [p, a] : fac) {
56	ph = ph / p * (p - 1);
57	if (a >= 2) mval = 0; // square factor detected
58	// flip sign per distinct prime if still square-free
59	// done after loop if mval != 0
60	divcnt *= (a + 1);
61	// compute p^{a+1}
62	long long pa1 = 1;
63	for (int i = 0; i < a + 1; ++i) pa1 *= p;
64	divsum *= (pa1 - 1) / (p - 1);
65	}
66	// set sign for μ if square-free
67	if (mval != 0) {
68	int k = (int)fac.size();
69	mval = (k % 2 == 0) ? 1 : -1;
70	}
71	phi[n] = ph;
72	mu[n] = mval;
73	d[n] = divcnt;
74	sigma[n] = divsum;
75	}
76
77	// quick correctness spot-checks
78	cout << "phi(12) = " << phi[12] << " (expect 4)\n";
79	cout << "mu(12) = " << mu[12] << " (expect 0)\n";
80	cout << "d(12) = " << d[12] << " (expect 6)\n";
81	cout << "sigma(12) = " << sigma[12] << " (expect 28)\n";
82
83	return 0;
84	}
85