MoreRSS

site iconJohn D. CookModify

I have decades of consulting experience helping companies solve complex problems involving applied math, statistics, and data privacy.
Please copy the RSS to your reader, or quickly subscribe to:

Inoreader Feedly Follow Feedbin Local Reader

Rss preview of Blog of John D. Cook

Lessons Learned With the Z3 SAT/SMT Solver

2025-03-18 01:50:36

Community best practices are useful for helping use a software product more effectively. I’ve just completed a small project using the Z3 solver. Here are some things I’ve learned:

  • My project involves an optimization problem: for a subset of Boolean variables, maximize the count of how many are true. My specific problem is solved much faster with Z3 by converting to a decision problem: set up a base problem to solve for the count being at least a certain fixed number, and iterate using bisection search to find the highest number satisfied. Bisection has been used for this problem before. Also, certain methods may possibly reduce the number of bisection steps.
  • Using Z3  “tactics” can greatly speed up the solve process. I found a good combination of tactics by trial and error, guided in part by the descriptions of the tactics. ChatGPT was of some help in finding good selections to try. An interesting paper discusses use of Monte Carlo tree search to define a good chain of tactics. The branching factor here is high, perhaps around 1000, though there are some redundancies in this number. Training multi-step MCTS might be expensive, but doing this once to get a good static chain of tactics might be worthwhile.
  • The strength of Z3 is in its extremely broad functionality, more so than its raw compute performance. It would be a daunting task for the Z3 team to fully optimize every possible solve option. I examined some of the SMT solver competitions to find faster codes. CVC5 on one case I tried was about twice as fast as Z3; I’ve seen similar reports in the literature. Presently I don’t find it worth the switching costs to use CVC5. One approach might be to use the very capable tactics engine of Z3 and pass the resulting modified problem to CVC5.
  • The specific formulation of the problem can make a big difference in solver performance. I’ve already seen this in the area of iterative linear solvers, for example diagonal matrix scaling can dramatically help (conjugate gradients) or hurt (multigrid) solver performance. Same thing here. Hence the huge importance in good “preprocessing“ for SAT/SMT solvers. One could wish the solver could handle all this automatically without user intervention. However, these powerful tools must be brandished very carefully for maximum effect.
  • Generally, one should move as much of the problem outside of the solver as possible, since the solver is the long pole in the tent in terms of scalability. For example if there is a Z3 integer that must be limited to a certain range and additionally some values in the interval must be blacklisted, it’s better, if possible, to compress all of the valid values into a single interval, to make testing for validity simpler in the Z3 code.
  • Along these lines: the Z3 tactics for propagating constants are not perfect; thus it can help to manually propagate constants (though this unfortunately makes the code more messy). This propagation can also sometimes allow for removal of unneeded constraints, further speeding up performance. Relatedly, some intriguing work by Benjamin Mikek shows how one can use the LLVM code optimizer to optimize the SMT problem in a way that is complementary to Z3 tactics, achieving significant speedup (for more info see here, here and here). I haven’t tried this but it seems promising.
  • Because of the scalability challenge of SMT solvers, various simplifying heuristics to modify the problem can be helpful. For example: solving a subproblem of the main problem and holding the resulting variables fixed in order to solve the rest of the problem. Or solving a simpler, smaller problem first to determine variable presets for the full problem. With these heuristics, one does not in general find the true global optimum; but the result may be adequate.
  • CPU threading does not work for my case (Z3 Python, macOS). Perfect parallelization of SAT and SMP is an unsolved (and perhaps in some sense not fully solvable) problem. One can naïvely parallelize bisection search by converting to trisection, etc., but this does not give perfect speedup (specif., log(P) speedup on P threads). Improvements to parallel bisection in some cases may be possible. Recent work by Armin Biere and colleagues looks promising; as I read it, near perfect speedup up to eight threads (at least for some problems).
  • Some of the main developers of Z3 are on Stack Overflow and have been active in the past answering questions. This seems very useful.

Resources like Handbook of Satisfiability and the proceedings of various SAT/SMT conferences seem helpful. More information on best practices for non-expert practitioners would be a great help to the community. If you know of any good resources, please share in the comments.

The post Lessons Learned With the Z3 SAT/SMT Solver first appeared on John D. Cook.

Grokking the James Bond chord

2025-03-17 07:43:08

The last chord of the James Bond theme is known as the James Bond chord. If you’ve ever seen a Bond movie, you’ll recognize it. Here’s a quick video showing how to play it on guitar.

I tried asking Grok 3 to generate some LilyPond code for sheet music a few days ago. LilyPond is pretty obscure as programming languages go—it’s not like asking for Python code—so I was skeptical, but Grok did well.

This afternoon I asked Grok to create sheet music for the James Bond chord. Specifically, I asked

  1. What is the “James Bond chord”?
  2. Create sheet music for the James Bond chord using Lilypond.
  3. Move the chord up an octave.

Here’s the output of each step.

What is the chord?

According to Grok,

The “James Bond chord” is a distinctive musical chord famously associated with the James Bond film franchise. It’s not just a single chord but rather a specific voicing and progression that evokes a sense of mystery, danger, and sophistication—perfectly matching the suave spy’s persona. Officially, it’s often identified as an E minor major 9 chord (EmMaj9), though its exact interpretation can vary slightly depending on the context.

The chord was first introduced by composer Monty Norman in the “James Bond Theme” for Dr. No (1962), arranged by John Barry, whose orchestration gave it its iconic swagger. …

First output

Grok produced the following LilyPond code

\version "2.22.0"

\score {
  \new Staff {
    \clef treble
    \key c \major % Default key, no sharps/flats; we'll specify accidentals manually
    \time 4/4
     <e g b dis' fis'>1 % EmMaj9 chord: E G B D# F#
  }
  \layout {}
  \midi {}
}

When compiled this gave the following music.

Grok also produce a MIDI file.

Second output

I wanted to move the chord up an octave for aesthetic reasons, putting the notes inside the staff. Grok complied, changing one line in the code, essentially adding an extra prime mark after each note.

<e' g' b' dis'' fis''>1 % EmMaj9 chord moved up an octave: E' G' B' D#'' F#''

This compiled to the following music.

Problems and prospects

Grok’s not perfect. In another experiment it produced code that wouldn’t compile. But when I told Grok that the code didn’t compile and asked it to try again, it worked.

I tried to remove the time signature, the C symbol. I asked Grok to remove it, and it did not. I asked Grok “How do you get LilyPond to produce music without a time signature?” and it told me two ways, neither of which worked.

I’ve used LilyPond occasionally for years, not to produce full sheets of music but to produce little fragments for blog posts. I’ve always found it a bit mysterious, in part because I jumped in and used it as needed without studying it systematically. There have been times when I thought about including some music notation in a blog post and didn’t want to go to the effort of using LilyPond (or rather the effort of debugging LilyPond if what I tried didn’t work). I may go to the effort more often now that I have a fairly reliable code generator.

Posts using LilyPond

The post Grokking the James Bond chord first appeared on John D. Cook.

Rates of convergence

2025-03-17 03:02:13

The last several posts have looked at counting two kinds of permutations: those that leave no consecutive integers and those that leave no integer fixed. As n grows large, the proportion of permutations of n elements that fall into both classes approaches 1/e. This post will look a little closer and as how fast each proportion approaches 1/e.

(Incidentally, if you’re looking for an problem to explore, take any statement of the form “X converges to Y” and ask how does X converge to Y. Several of the posts on this blog follow that pattern.)

Definitions

To review our definitions, Q(n) is the number of permutations p of 1, 2, 3, …, n + 1 such that there are no consecutive integers, i.e. for no k does p(k + 1) = p(k) + 1. So, for example, 13254 counts, but 15423 does not count because 2 and 3 are consecutive.

Let D(n) be the number of derangements of 1, 2, 3, …, n, i.e. permutations p such that for no k does p(k) = k. So, for example, 23451 counts, but 21354 does not because 3 remains fixed.

Note that Q(n) counts a certain kind of permutation of n + 1 elements and D(n) counts a certain kind of permutation of n elements. This is a little awkward, but I chose my definition of Q to match OEIS. In order to count the proportion of permutations of n elements with no consecutive elements we should look at  Q(n − 1)/n!. Both Q(n − 1)/n! and D(n)/n! converge to 1/e, but not at the same rate.

D convergence

The convergence of D(n)/n! to 1/e is very fast. In fact,

D(n) = [n! / e]

where [x] is x rounded to the nearest integer. So the difference between D(n)/n! and 1/e is less than 0.5/n!. For n greater than or equal to 18, D(n)/n! = 1/e within the limits of floating point precision. As a rule of thumb, you can see about three significant figures in a plot, so for n greater than about 6 you can’t see the difference between D(n)/n! and 1/e.

Q convergence

The convergence of Q(n − 1)/n! to 1/e is leisurely by comparison.

Here’s why. In this post I showed Q(n − 1) is asymptotic to (n + 1)(n − 1)! / e.

So it is true that Q(n − 1)/n! is asymptotic to 1/e, but it is more accurate to say that it is asymptotic to (n + 1) / ne.

D and Q convergence

If we look at the ratio of (n + 1) Q(n) / n to n! then it converges to 1/e about as fast as the ratio of D(n) to n! does.

Said another way,

(n + 1) Q(n) / n ~ D(n)

The ratio of (n + 1) Q(n) / n to D(n) converges to 1 very quickly, with an error on the order of 10n.

On to something else

I intend to end the series here. My next post will probably be about something else, something less mathematical.

So here’s the complete series.

The post Rates of convergence first appeared on John D. Cook.

Gluons, quarks, letters, and envelopes

2025-03-17 00:17:54

Yesterday I wrote a couple of posts about a combinatorics question that lead to OEIS sequence A000255. That page has this intriguing line:

This comment derives from a family of recurrences found by Malin Sjodahl for a combinatorial problem for certain quark and gluon diagrams (Feb 27 2010)

I love how pulling on a thread can lead you to unexpected places. What in the world does my little permutation problem have to do with particle physics‽

I found Malin Sjödahl’s website, and apparently the citation above refers to this paper from 2009. The author’s site doesn’t list any papers from 2010. Maybe the paper was published electronically in 2009 and in print in 2010.

Counting tensors

The paper is mostly opaque to me since I don’t know particle physics, but at one point Sjödahl says

The problem of finding all such topologies is equivalent to the number of ways of mapping N elements to each other without mapping a single one to itself …

and says that the solution is

N! \sum_{i=0}^N \frac{(-1)^i}{i!} \to \frac{N!}{e}

Sjödahl is not counting physical permutations but the number possible tensors associated with gluons and quark interaction diagrams.

The right-hand side above is essentially the same as the asymptotic estimate for the function Q(n) in the previous post.

I didn’t find the recurrence that the OEIS comment alluded to. Perhaps I’m not looking at the same paper. Perhaps I’m not looking hard enough because I’m skimming a paper whose contents I don’t understand.

The Montmort letter problem

In more mathematical terminology, Sjödahl is counting the number of permutations of N objects with no fixed point, known as the number derangements of a set N objects.

If you divide by the number of possible permutations N! you have the probability that a random permutation moves every object.

This was historically known as the Montmort letter problem, named after Pierre-Remond Montmort who asked the following question. If N letters are randomly assigned to N envelopes, what is the probability that no letter ends up in the correct envelope?

The probability converges to 1/e as N approaches infinity. It approaches 1/e quickly, and so this is a good approximation even for moderate N. More on the rate of convergence in the next post.

Previous posts in this series

The post Gluons, quarks, letters, and envelopes first appeared on John D. Cook.

Example using a generating function to find asymptotic behavior

2025-03-16 07:31:46

The previous post looked at how to compute Q(n), the number of permutations of 1, 2, 3, …, n + 1 that contain no consecutive integers. We found a way to numerically compute Q(n) but no analytic expression that would let us compute asymptotics.

The sequence Q(n) is sequence A000255 in OEIS, and OEIS gives the exponential generating function of Q:

\sum_{n=0}^\infty \frac{Q(n)}{n!} = \frac{\exp(-x)}{(1-x)^2}

We can use this function and Theorem 5.2.1 from [1] to get the asymptotic form of Q(n). According to that theorem, the coefficient of xn in our generating function is asymptotically the same as the coefficient of xn in the principle part at the singularity at 1. This principle part is

\frac{1}{e(1-x)^2} + \frac{1}{e(1-x)} = \frac{1}{e}\left(2 + 3x + 4x^2 + 5x^3 + \cdots\right)

and so the coefficient of xn is (n + 2)/e.

So

Q(n) \sim \frac{n + 2}{e} n!

and Q(n) / (n + 1)! approaches 1/e for large n, as conjectured in the previous post.

[1] Herbert S. Wilf. Generatingfunctionology. Second edition.

The post Example using a generating function to find asymptotic behavior first appeared on John D. Cook.

Permutations with no consecutive elements

2025-03-16 00:32:23

I was working on something this week that made me wonder how many permutations break up all consecutive elements. For example, if we’re looking at permutations of abcde then acebd counts, but acdbe does not. I’d like to count the number of such permutations, and estimate for large N the number of permutations of N elements with no consecutive terms.

A brief search lead to OEIS sequence A000255 which gives a recursive formula for the answer. Before looking at that, let’s define some terms and do some brute force calculation to get a feel for the problem.

Brute force calculation

Given a positive integer n, define Q(n) to be the number of partitions p of 1, 2, 3, …, n + 1 such that for no k does p(k + 1) = p(k) + 1.

from itertools import permutations

def no_consec(perm):
    for k in range(len(perm) - 1):
        if perm[k+1] == perm[k] + 1:
            return False
    return True

def Q(n):
    count = 0
    for p in permutations(range(n+1)):
        if no_consec(p):
            count += 1
    return count

We could use this code to compute Q(n) for moderate-sized n, but the time required to compute Q(n) this way is proportional to n! and so that won’t scale well.

Recurrence relation

OEIS sequence A000255 gives the values of what we’re calling Q(n), and the first few values match the output of the code above. But OEIS does better, giving a recursive solution for Q:

Q(n) = n Q(n − 1) + (n − 1) Q(n − 2)

for n > 1 with initial conditions Q(0) = Q(1) = 1.

As with Fibonacci numbers, you could implement this easily as a recursive function

def Q2(n):
    if n in[0, 1]:
        return 1
    return n*Q2(n-1) + (n-1)*Q2(n-2)

but it is much faster to implement iteratively.

def Q3(n):
    q = [1]*(n+2)
    for k in range(2, n+1):
        q[k] = k*q[k-1] + (k-1)*q[k-2]
    return q[n]

You could make the recursive implementation faster with memoization, using lru_cache from functools. However, memoization does not prevent recursion; it just speeds up recursion using a cache, and you can still get an error saying maximum recursion depth exceeded.

Asymptotics

We can calculate Q(n) for large n, but we don’t have an analytic expression [1] that will let us find the asymptotic behavior. I intended to work that out in my next post.

Empirically, it seems Q(n) approaches (n + 1)!/e as n goes to infinity, but I don’t have a proof of that at the time of writing. If this is correct, this means that in the limit the probability of a permutation having no fixed point equals the probability of it having no consecutive values.

Update: Proof

***

[1] There is an analytic expression for Q(n) in terms of a hypergeometric function

Q(n) = (−1)n2F0(2, −n; ; 1)

but this expression isn’t helpful in finding the asymptotic behavior.

The post Permutations with no consecutive elements first appeared on John D. Cook.