Monthly Archives: December 2011

Generating Functions as Cardinality of Set Maps

There is a class of all cardinalities \mathbf{Card}, and it
has elements 0, 1 and operations +, \cdot, and so forth defined on it. Furthermore, there is a map
\mathrm{card}\colon\mathbf{Set}\to\mathbf{Card} which takes
sets to cardinalities such that \mathrm{card}(A\times  B)=\mathrm{card}(A)\cdot\mathrm{card}(B) (and so on).

Ordinary generating functions can be thought of entirely analogously
with set maps \mathbf{Set}\to\mathbf{Set} replacing sets:
There is a class \mathbf{GenFunc} with elements 0,
1, and operations +, \cdot. Furthermore,
there is a (partial) map \mathrm{genFunc}\colon (\mathbf{Set}\to\mathbf{Set})\to\mathbf{GenFunc} such that \mathrm{genFunc}(F\times G)=\mathrm{genFunc}(F)\cdot\mathrm{genFunc}(G) (and so on). Here, F\times G is defined by (F\times G)(S)=F(S)\times G(S). Other operations on set maps (like disjoint union) are similarly defined pointwise.

(This is probably obvious and trivial to anyone who actually works
with generating functions, but it only occurred to me recently, so I
thought I’d write a blog post about it.)

The class \mathbf{GenFunc} is in fact a set, and is just the set of formal power series \{\sum_{i\geq 0} a_i z^i\mid a_i\in\mathbb{N}\}. The partial map \mathrm{genFunc} takes F to \sum_{i\geq 0} a_i z^i just in case F is “canonically isomorphic” (a notion I’ll leave slippery and undefined but that can be made precise) to the map Z\mapsto \coprod_{i\geq 0} \{1,2,\ldots,a_i\}\times Z^i, where \coprod indicates disjoint union.

That provides a semantics for ordinary generating functions. Furthermore, this semantics has a number of features beyond those of cardinality. For example, in addition to respecting \coprod and \times, \mathrm{genFunc} represents composition.

A similar semantics can be provided for exponential generating functions, but it takes a little more work. In particular, we have to single out \lbrack0,1\rbrack as a distinguished set. Let \mathbf{Meas} be the smallest set containing all measurable subsets of \lbrack0,1\rbrack^n for any finite n and which is closed under finite products, countable disjoint unions, and products with sets \{1,\ldots,n\} for finite n\geq 1.

We can define the measure \mu of all sets in \mathbf{Meas} by extending Lebesgue measure in the obvious way (taking the product of a set with \{1,\ldots,n\} will multiply the measure by n). Furthermore, notice that, by construction, every element of every set M in \mathbf{Meas} is a tuple which (after flattening) has all of its elements either natural numbers or elements of \lbrack 0,1\rbrack and has at least one element of \lbrack 0,1\rbrack. Therefore, we can define a pre-ordering on M by comparing the corresponding first elements that are in \lbrack 0,1\rbrack.

The point of all that is that, for M\in\mathbf{Meas}, we can form the set M^n_{<}=\{\langle m_1,\ldots,m_n\rangle\mid m_1 < \cdots < m_n\} which will again be in \mathbf{Meas} and its measure will be \mu(M)^n/n!. The corresponding statement with cardinality is not true since you have to worry about the case when elements in the tuple are equal (\mathrm{card}(X^n_{<}) = \mathrm{card}(X)(\mathrm{card}(X)-1)\cdots(\mathrm{card}(X)-n+1)/n!) but the set of tuples that have duplicates has measure 0, so by working with measure, we can get the equality we want.

Finally, let \mathbf{ExpGenFunc} be the set of formal power series \{\sum_{i\geq 0} a_i/i! z^i\mid a_i\in\mathbb{N}\}. The partial map \mathrm{expGenFunc} takes F to \sum_{i\geq 0} a_i/i! z_i just in case F is “canonically isomorphic” to the map Z\mapsto \coprod_{i\geq 0} \{1,2,\ldots,a_i\}\times Z^i_{<} for all Z in \mathbf{Meas}. Just as before, this map respects +, \cdot, composition, etc.

Note that the exponential generating functions are usually explained via labeled objects and some sort of relabeling operation. This approach weasels out of that by observing that the event that there was a label collision has probability 0, so you can just ignore it.

1 Comment

Filed under Uncategorized

Mathematica and Quantifier Elimination

In 1931, Alfred Tarski proved that the real ordered field (\mathbb{R}, 0, 1, +,\times, <) allows quantifier elimination: i.e., every first-order formula is equivalent to one with no quantifiers.  This is implemented in Mathematica’s “Resolve” function.

The Resolve function is called like Resolve[formula,domain] where domain gives the domain for the quantifiers in formula.  Since we’ll always be working over \mathbb{R} in this blog post, let’s set that to be the default at the start.

In[1]:= Unprotect[Resolve]; Resolve[expr_] := Resolve[expr, Reals]; Protect[Resolve];

Now let’s see what quantifier elimination lets you do!

(A couple of caveats first though: First, many of these algorithms are extremely inefficient. Second, I had some trouble exporting the Mathematica notebook, so I basically just copy-and-pasted the text. Apologies if it’s unreadable.)

How many solutions?

Let’s start with just existential formulas.  By eliminating quantifiers from \exists x\, \phi(x,a), we can tell what the conditions are on a such that there’s at least one solution x.  For example:

In[2]:= Resolve[Exists[x, x^2 + b x + c == 0]]
Out[2]= -b^2 + 4 c <= 0

This just tells you that there’s a solution to the quadratic if the discriminant is non-negative.  Let’s turn this into a function:

In[3]:= atLeastOneSolution[formula_, variable_] := Resolve[Exists[variable, formula]]

Now we can verify that cubics always have solutions:

In[4]:= atLeastOneSolution[x^3 + b x^2 + c x + d == 0, x]
Out[4]= True

Now suppose we wanted to find when something has at least two solutions.  Just like resolving \exists x\,\phi(x) told us when there was at least one, \exists x_1,x_2\,x_1\ne x_2\wedge \phi(x_1)\wedge\phi(x_2) will be true exactly when there are at least two.

This is just as easy to program as atLeastOneSolution was, except that when we create the variables x_1 and x_2 we have to be careful to avoid capture (what if one of those two already appeared in \phi?).  Mathematica provides a function called Unique where if you call Unique[], you’re guaranteed to get back a variable that’s never been used before.  With that we can define atLeastTwoSolutions correctly (edit: actually, this isn’t right if the passed-in variable is also bound in the passed-in formula):

In[5]:= atLeastTwoSolutions[formula_, v_] :=
With[{s1 = Unique[], s2 = Unique[]},
Resolve[
Exists[{s1, s2},
s1 != s2 && (formula /. v -> s1) && (formula /. v -> s2)]]]

We can check this by verifying that quadratics have two solutions when the discriminant is strictly positive:

In[6]:= atLeastTwoSolutions[x^2 + b x + c == 0, x]
Out[6]= -b^2 + 4 c < 0

Here’s the condition for the cubic to have at least two solutions:

In[7]:= atLeastTwoSolutions[x^3 + b x^2 + c x + d == 0, x]
Out[7]= c < b^2/3 &&
1/27 (-2 b^3 + 9 b c) - 2/27 Sqrt[b^6 - 9 b^4 c + 27 b^2 c^2 - 27 c^3] <=
d <= 1/27 (-2 b^3 + 9 b c) + 2/27 Sqrt[b^6 - 9 b^4 c + 27 b^2 c^2 - 27 c^3]

Note that (and I believe Resolve always does this) the c<b^2/3 condition given first is sufficient that the later square root is well-defined:
In[8]:= Resolve[ForAll[{b, c}, c < b^2/3 ⇒ b^6 - 9 b^4 c + 27 b^2 c^2 - 27 c^3 > 0]]
Out[8]= True

It’s clear that we can determine when there at least n solutions by a very similar trick: just resolve \exists x_1,\ldots,x_n (x_i\ne x_j,i\ne j)\wedge (\phi(x_i),\forall i).

We’ll first write a helper function to produce the conjunction of inequalities we’ll need:

In[9]:= noneEqual[vars_] :=
And @@ Flatten[Table[If[s1 === s2, True, s1 != s2], {s1, vars}, {s2, vars}]]
In[10]:= noneEqual[{x, y, z}]
Out[10]= x != y && x != z && y != x && y != z && z != x && z != y

And now we’ll write atLeastNSolutions:

In[11]:= atLeastNSolutions[formula_, v_, n_] := With[{sList = Array[Unique[] &, n]},
Resolve[
Exists[sList,
noneEqual[sList] && (And @@ Table[formula /. v -> s, {s, sList}])]]]

Given atLeastNSolutions, we can easily write exactlyNSolutions:

In[12]:= exactlyNSolutions[formula_, v_, n_] :=
BooleanConvert[
atLeastNSolutions[formula, v, n] && ! atLeastNSolutions[formula, v, n + 1]]

I used BooleanConvert instead of Resolve since there won’t be any quantifiers left in the formula, so we just have to do Boolean simplifications.

In[13]:= exactlyNSolutions[x^3 + b x^2 + c x + d == 0, x, 2]
Out[13]= ! 1/27 (-2 b^3 + 9 b c) - 2/27 Sqrt[b^6 - 9 b^4 c + 27 b^2 c^2 - 27 c^3] < d <
1/27 (-2 b^3 + 9 b c) + 2/27 Sqrt[b^6 - 9 b^4 c + 27 b^2 c^2 - 27 c^3] &&
1/27 (-2 b^3 + 9 b c) - 2/27 Sqrt[b^6 - 9 b^4 c + 27 b^2 c^2 - 27 c^3] <=
d <= 1/27 (-2 b^3 + 9 b c) +
2/27 Sqrt[b^6 - 9 b^4 c + 27 b^2 c^2 - 27 c^3] && c < b^2/3
In[14]:= exactlyNSolutions[x^2 + b x + c == 0, x, 1]
Out[14]= -b^2 + 4 c <= 0 && -b^2 + 4 c >= 0

This last calculation shows that a quadratic has exactly one solution exactly when the discriminant is both nonnegative and nonpositive (as you can see, there is no guarantee that the formula will be in it’s simplest form).
We now have a way to test whether a formula with one free variable has n solutions for specific values of n, since exactlyNSolutions will return either True or False if you quantify out the only variable.  For example:

In[15]:= p = x^4 - 3 x^3 + 1
Out[15]= 1 - 3 x^3 + x^4
In[16]:= Plot[Evaluate[p], {x, -3, 3}]

In[17]:= exactlyNSolutions[p == 0, x, 2]
Out[17]= True

It would be nice, however, to have a function which will just tell you how many solutions such a formula has.

In the single-variable polynomial case, we could just try exactlyNSolutions for n=0,1,2,\ldots until we find the right n.  However, there might not be finitely many solutions if the formula involves inequalities or higher-dimension polynomials (e.g., x^2 + y^2 = 1 has infinitely many solutions).

How can we tell if a formula has infinitely many solutions?  Well, the fact that \mathbb{R} has quantifier elimination implies that \{x\in\mathbb{R}\mid \phi(x)\} for \phi with just x free must be a finite union of points and open intervals (since the only quantifier free terms are t_1=t_2 and t_1<t_2. Therefore \{x\in\mathbb{R}\mid\phi(x)\} is infinite iff it contains a non-empty open interval, i.e., iff \exists a,b\,a<b\wedge\forall x\,(a<x<b\implies\phi(x)).

In[18]:= infinitelyManySolutions[formula_, v_] := With[{a = Unique[], b = Unique[]},
Resolve[Exists[{a, b}, a < b && ForAll[v, a < v < b ⇒ formula]]]]

To test:

In[19]:= infinitelyManySolutions[Exists[y, x^2 + y^2 == 1], x]
Out[19]= True

Now we can write numberOfSolutions and be assured that it will always (theoretically) terminate for any formula with a single free variable:

In[20]:= numberOfSolutions[formula_, v_] :=
If[infinitelyManySolutions[formula, v], Infinity,
Block[{n = 0},
While[! exactlyNSolutions[formula, v, n], n++];
n]]

A few examples:

In[21]:= numberOfSolutions[p == 0, x]
Out[21]= 2
In[22]:= numberOfSolutions[p > x^2, x]
Out[22]= ∞
In[23]:= numberOfSolutions[p > x^6 + 5, x]
Out[23]= ∞
In[24]:= numberOfSolutions[p > x^6 + 6, x]
Out[24]= 0
In[26]:= Plot[{p, x^6 + 5, x^6 + 6}, {x, -1.6, -1},
PlotLegend -> {HoldForm[p], x^6 + 5, x^6 + 6}, LegendPosition -> {1, 0},
ImageSize -> Large]

Up to now, all our functions have taken single variables, but we can accomodate tuples of variables as well.  First, we’ll define the analogue of noneEqual to produce the formula asserting that none of the given tuples are equal (recall that two tuples are unequal iff a pair of corresponding components is unequal):

In[27]:= noTuplesEqual[tuples_] := And @@ Flatten[Table[If[t1 === t2, True,
Or @@ MapThread[#1 != #2 &, {t1, t2}]], {t1, tuples}, {t2, tuples}]]
In[28]:= noTuplesEqual[{{x[1], y[1]}, {x[2], y[2]}}]
Out[28]= (x[1] != x[2] || y[1] != y[2]) && (x[2] != x[1] || y[2] != y[1])

Now we can add rules to our old function to deal with tuples of variables as well:

In[29]:= atLeastNSolutions[formula_, variables_List, n_] := With[
{sList = Array[Unique[] &, {n, Length[variables]}]},
Resolve[
Exists[Evaluate[Flatten[sList]],
noTuplesEqual[sList] &&

And @@
Table[
formula /. MapThread[Rule, {variables, tuple}], {tuple, sList}]]]];

We can extend infinitelyManySolutions by observing that a formula \phi(x_1,\ldots,x_n) has infinitely many solutions iff some projection \exists x_1,\ldots,x_{i-1},x_{i+1},\ldots,x_n \phi(x_1,\ldots,x_n) does.

In[30]:= infinitelyManySolutions[formula_, variables_List] := Or @@ Table[
infinitelyManySolutions[Exists[Select[variables, ! (# === v) &], formula],
v], {v, variables}]
In[33]:= ContourPlot[{x^2 + y^3 - 2, x^2 + y^2/4 - 2}, {x, -3, 3}, {y, -3, 3}]

In[34]:= exactlyNSolutions[x^2 + y^3 == 2 && x^2 + y^2/4 == 2, {x, y}, 2]
Out[34]= False

(There are actually four solutions. This example of a set equations for which it’s difficult to tell how many solutions there are by graphing is from Stan Wagon)

Solving Polynomial Equations

In the last section, we saw how to use quantifier elimination to find out how many roots there are.  But how can you actually find the roots?

In a certain sense, you’ve already found them just when you identified how many there are!  To “find” a root in this sense, you just introduce a new symbol for it, and have some means for answering questions about its properties.  Given some property \phi(x), if you want to determine if it holds of the 6th root of some polynomial p with 17 roots, then you just have to decide \exists x_1,\ldots,x_{17}\,(x_i<x_j,i<j)\wedge(p(x_i),\forall i)\wedge\phi(x_6).

We can implement this by a function withSpecificRoot, that takes a variable, the formula it’s supposed to be a solution to, which of the roots it’s a solution to, and a formula in which you want to use this root:

In[35]:= withSpecificRoot[variable_, rootFormula_, whichRoot_, totalRoots_, formula_] :=

With[{roots = Array[Unique[] &, totalRoots]},
Resolve[
Exists[Evaluate[roots~Join~{variable}],
Less[Sequence @@ roots] &&
variable ==
roots[[whichRoot]] &&
(And @@
Table[(rootFormula /. variable -> root), {root, roots}]) && formula]]]

We can tell where various roots are with respect to already-known real numbers:
In[36]:= withSpecificRoot[x, x^2 - 3 == 0, 1, 2, x < 3]
Out[36]= True
In[37]:= withSpecificRoot[x, p == 0, 1, 2, x < 1]
Out[37]= True
In[38]:= withSpecificRoot[x, p == 0, 2, 2, x < 1]
Out[38]= False

We can also compute relationships between roots like \sqrt{2}+\sqrt{3}=\sqrt{5+2\sqrt{6}}:

In[39]:= withSpecificRoot[sqrt6, sqrt6^2 == 6, 2, 2,
withSpecificRoot[lhs, lhs^2 == 5 + 2 sqrt6, 2, 2,
withSpecificRoot[sqrt3, sqrt3^2 == 3, 2, 2,
withSpecificRoot[sqrt2, sqrt2^2 == 2, 2, 2,
lhs == sqrt3 + sqrt2
]]]]
Out[39]= True

That’s all I have time for now, but I hope to write another blog post on the subject soon!

Leave a comment

Filed under Uncategorized