An almost periodic function

Posted on by John

This post takes a more abstract view of the previous post. That post looked at the concrete question of whether a number ever has the same sine in radians as in degrees. The relation between radians and degrees is irrelevant except that π/180 is an irrational number.

Suppose α and β are two positive numbers such that α/β is irrational. In the previous post, α = 1 and β = π/180. Then the function

f(x) = sin(αx) − sin(βx)

is almost periodic: it is not periodic, but it comes close to being periodic, as close as you’d like provided you’re willing to look over a sufficiently long range of x‘s.

The identity

sin(αx) − sin(βx) = 2 cos((α + β)x/2) sin((α − β)x/2)

shows that f(x) is the product of two periodic functions but is not periodic itself. The periods of the cosine and sine above never coincide because the ratio of their frequencies is irrational.

The zeros of f are not periodic, though they can be divided into two subsequences that are periodic.

When sine of x degrees equals sine of x radians

Posted on by John

Ordinarily the sine of x radians and the sine of x degrees are very different numbers. Having your calculator in radian mode when it should be in degree mode, or vice versa, results in a major error.

But sometimes it doesn’t matter. A trivial example is when x = 0. A more interesting example is

x = 180π/(180 + π) = 3.08770208….

For that value of x,

sin(x) = sin(x°).

In this article I’ll use the common convention of using radians by default and denoting degrees with ° as above.

Note that

x = πx°/180

and so we are interested in solutions to the equation

sin(x) = sin(πx/180)

Now two angles A and B have the same sine if they differ by a multiple of 2π, or if they’re supplementary (i.e. A = π − B), or both. To put it another way, if A and B have the same sine, they are either equal mod 2π or supplementary mod 2π. This means that

sin(x) = sin(πx/180)

if and only if

x = πx/180 + 2πk

or

x = π − πx/180 + 2πk

for some integer k.

Therefore all solutions have the form

x = 360πk/(180 − π)

or

x = 180π(2k + 1)/(180 + π).

Alternative solution

The derivation above is correct, but it occurred to me later that a simpler argument would be to use the identity

sin(A) − sin(B) = 2 cos((A + B)/2) sin((A − B)/2).

Thus A and B have the same sine if

cos((A + B)/2) = 0

or if

sin((A − B)/2) = 0.

These two possibilities correspond to the two families of solutions above.

Density

When reduced modulo 2π, both families are dense in [0, 2π]. This means that for every y in [−1, 1], there is a number x such that

sin(x) = sin(x°) ≈ y

and we can make the approximation as good as we’d like.

Example 1

For example, today is July 22, so let’s set y = 0.722. We’d like to find a value of x such that the sine of x radians and the sine of x degrees both approximately equal 0.722. And let’s say our approximation tolerance is ε = 0.0001.

We can search for a value of x in the first family of solutions by looking for a value of k with

| sin(360πk/(180 − π)) − 0.722 | < 0.0001

and the smallest such k is 96343 and so

x = 360×96343 π/(180 − π) = 616093.78713621…

will do, and sin(x) = 0.72191…

Example 2

Now let’s set y = 0.2026 and look for a solution in the other family of solutions, and this time let’s set ε = 10−6. The smallest value of k such that

| sin(180π(2k + 1)/(180 + π)) − 0.2026 | < 10−6

is k = 741141. Then

sin( 4576848.310950611 ) = sin( 4576848.310950611° ) = 0.202600139…

Forensic accounting in Python

Posted on by John

I recently had a project in which I had to reverse engineer a data analysis. There was some ambiguity regarding which of several possibilities someone chose for several of the variables, something analogous to the following example.

Suppose you have three numbers with uncertain values with a known, or at least purported, sum. The first number could be 31, 41, or 59; the second could be either 26 or 53; the last could be 58, 97, 93, or 23.

The following code enumerates all 3 × 2 × 4 = 24 possibilities and prints their sums.

from itertools import product

# Example input
possibilities = [(31, 41, 59), (26, 53), (58, 97, 93, 23)]

for combo in product(*possibilities):
    total = sum(combo) 
    print(f"Combination {combo} sums to: {total}")

In this example all the sums are unique, though of course that might not happen in practice. If, for example, you know the sum is 187, you know the three numbers were 41, 53, and 93. If the reported sum is 200, you know some assumption has been violated because none of the possible choices add up to 200.

More forensics posts

Locally everywhere does not imply everywhere

Posted on by John

A couple days ago, Levent Alpöge, a mathematician working at Anthropic, discovered a counterexample to the Jacobian conjecture using Claude Fable 5.

I was curious whether most mathematicians were trying to prove or disprove the conjecture, so I asked Claude.

Before a counterexample to the Jacobian conjecture was found, did most mathematicians believe it was true or false?

Claude’s response was

The premise of this question isn’t quite right — no counterexample to the Jacobian conjecture has been found. It remains an open problem in mathematics: no one has proven it true, and no one has found a counterexample disproving it. … If you encountered a claim that a counterexample was found, do you have a source for that? I’d be happy to look into it, since that would actually be a major result in algebraic geometry if true.

Of course Claude doesn’t know that it solved the conjecture. It didn’t even solve the conjecture. It was an inanimate tool in the hand of a mathematician, just like a piece of chalk or a dry erase marker.

The middle part of Claude’s response was that mathematicians are (were) divided on whether the conjecture is true. So it was not like the Riemann hypothesis, which most people believe to be true, or the P = NP conjecture, which most people believe to be false.

Now what is the Jacobian conjecture? It says that a polynomial function from ℝn to ℝn with constant, non-zero Jacobian determinant has a polynomial inverse. (The conjecture was stated more generally for fields of characteristic 0, in which the derivatives defining the Jacobian would have to be defined algebraically, not in terms of limits.)

Alpöge came up with a counterexample, a polynomial function from ℝ³ to ℝ³ with constant Jacobian determinant −2. The function is

\begin{align*}F(x,y,z)={}\bigl(~\!& z (1+xy)^3 + y^2 (1+xy) (4+3xy),\\ &y + 3x(1+xy)^2 z + 3xy^2 (4+3xy), \\ &2x - 3x^2 y - x^3 z ~\!\bigr).\end{align*}

It’s a tedious but simple calculus exercise to show that the determinant equals −2 everywhere. The inverse function theorem says that a function is locally invertible at any point where the Jacobian determinant is non-zero, so Alpöge’s function is locally invertible everywhere.

However, the function takes on some values more than once. For example, (0, 0, −1/4) and (1, −3/2, 13/2) both map to (−1/4, 0, 0). Therefore the function is not invertible globally. So not only does the function not have a polynomial inverse, it doesn’t have an inverse even if you allow non-polynomial functions.

Alpöge’s counterexample disproves the Jacobian conjecture for n = 3. It can trivially be extended to all n > 3 by defining the function to be Alpöge’s function for three variables and the identity for the rest. The conjecture remains open for n = 2.

Volume to Area ratio for Regular Solids

Posted on by John

The volume of a sphere of radius r is

V = 4πr³ / 3

and the surface area is

A = 4πr²

and so the ratio of volume to area is

V / A = r / 3.

Surprisingly, the same ratio holds for all regular solids if r is the radius of the largest sphere that can be inscribed inside the regular solid.

For example, if the edge of a cube is a, then r = a/2. The volume is 8r³, the area is 24r², and the ratio is r/3.

The relationship between edge length and radius, and between radius and volume, is more complicated for the four other regular solids (tetrahedron, octahedron, dodecahedron, and icosahedron). However, in each case the ratio of volume to area is r/3.

The proof is surprisingly simple. Pick a face and form a pyramid by connecting each face vertex to the center of the inscribed sphere. The pyramid has height r and volume equal to B/3 where B is the area of the base. If the regular solid has f faces, the volume of the solid is fBr / 3 and the area is fB. So the ratio of volume to area is r/3.

The theorem generalizes to n > 3 dimensions. The formula for the volume of a pyramid in n dimensions is Bh/n where B is the (n − 1)-dimensional volume of the base, and so the ratio of n-dimensional volume of a regular solid to (n − 1)-dimensional volume of its boundary is r/n.

Solving a chess puzzle with Grok 4.5

Posted on by John

I’ve written several posts about using Claude or ChatGPT to generate Prolog or Lean code to solve a chess puzzle. I didn’t think Grok would be up to the task, though I didn’t try it. I’ve heard good things about Grok 4.5, so I gave it a shot. It did great.

Here’s the problem, a variation on the n-queens problem formulated by Martin Gardner.

Can you place five white queens and three black queens on a 5 × 5 board so that no queen of one color attacks a queen of another color?

Here’s the prompt I have Grok 4.5.

Generate SWI Prolog code to solve the following chess puzzle. On a 5 by 5 chessboard, place 5 white queens and 3 black queens so that no queen is attacking a queen of the other color. Have the code enumerate all solutions.

The code ran correctly, giving all eight solutions.

I then asked Grok

Now generate Lean 4 code to solve the same puzzle.

The code did not run on the first or second attempts but ran correctly on the third iteration. That’s better than I’ve seen before. It’s not an apples-to-apples comparison, but it has usually taken many iterations with Claude before Lean code would run.

Solutions

There are eight solutions, each being a flip or a rotation of one fundamental solution. Here’s a representative.

SWI Prolog code

:- use_module(library(clpfd)).

%% Enumerate all ways to place 5 white queens and 3 black queens
%% on a 5x5 board so that no queen attacks a queen of the opposite colour.
%% Same-colour queens are allowed to attack each other.
%%
%% Positions are numbered 1..25 in row-major order:
%%   1  2  3  4  5
%%   6  7  8  9 10
%%  11 12 13 14 15
%%  16 17 18 19 20
%%  21 22 23 24 25

solve :-
    solution(Whites, Blacks),
    print_solution(Whites, Blacks),
    fail.
solve :-
    writeln('No more solutions.').

solution(Whites, Blacks) :-
    length(Whites, 5),
    length(Blacks, 3),
    append(Whites, Blacks, All),
    All ins 1..25,
    all_distinct(All),
    % canonical order so each set of positions is generated once
    Whites = [W1,W2,W3,W4,W5],
    W1 #< W2, W2 #< W3, W3 #< W4, W4 #< W5,
    Blacks = [B1,B2,B3],
    B1 #< B2, B2 #< B3, % every black is safe from every white maplist(safe_from_all(Whites), Blacks), label(All). safe_from_all(Whites, Black) :- maplist(no_attack(Black), Whites). %% Two queens do *not* attack each other. no_attack(P, Q) :- R1 #= (P-1) div 5 + 1, C1 #= (P-1) mod 5 + 1, R2 #= (Q-1) div 5 + 1, C2 #= (Q-1) mod 5 + 1, R1 #\= R2, % different row C1 #\= C2, % different column abs(R1-R2) #\= abs(C1-C2). % different diagonal print_solution(Whites, Blacks) :- format('Whites: ~w~n', [Whites]), format('Blacks: ~w~n', [Blacks]), forall(between(1, 5, R), ( forall(between(1, 5, C), ( Pos is (R-1)*5 + C, ( memberchk(Pos, Whites) -> write('W ')
                      ; memberchk(Pos, Blacks) -> write('B ')
                      ; write('. ')
                      )
                    )),
             nl )),
    nl.

Lean 4 code

-- Lean 4 solution for the 5×5 black-and-white queens puzzle.
-- 5 white queens + 3 black queens; no opposite-colour attacks.

def absDiff (a b : Nat) : Nat :=
  if a ≥ b then a - b else b - a

def attacks (p q : Nat × Nat) : Bool :=
  p.1 == q.1 || p.2 == q.2 || absDiff p.1 q.1 == absDiff p.2 q.2

/-- All 25 squares (0-based, row-major). -/
def allPositions : List (Nat × Nat) :=
  (List.range 5).flatMap fun r =>
    (List.range 5).map fun c => (r, c)

/-- Combinations of length exactly `k` (order-preserving). -/
def combinations {α : Type} : Nat → List α → List (List α)
  | 0,     _      => [[]]
  | _,     []     => []
  | k + 1, x::ys  =>
      (combinations k ys).map (fun zs => x :: zs) ++
      combinations (k + 1) ys

def isAttackedBy (pos : Nat × Nat) (queens : List (Nat × Nat)) : Bool :=
  queens.any (fun q => attacks q pos)

/-- Every solution as a pair (white squares, black squares). -/
def solve : List (List (Nat × Nat) × List (Nat × Nat)) :=
  (combinations 3 allPositions).foldl (fun acc blacks =>
    let forbidden := allPositions.filter fun p =>
      blacks.contains p || isAttackedBy p blacks
    let available := allPositions.filter fun p =>
      !(forbidden.contains p)
    acc ++ (combinations 5 available).map fun whites => (whites, blacks)
  ) []

/-- Pretty-print one board. -/
def printBoard (whites blacks : List (Nat × Nat)) : IO Unit := do
  for r in List.range 5 do
    let mut row := ""
    for c in List.range 5 do
      let p := (r, c)
      row := row ++
        (if whites.contains p then "W "
         else if blacks.contains p then "B "
         else ". ")
    IO.println row
  IO.println ""

/-- Enumerate and display every solution. -/
def main : IO Unit := do
  let sols := solve
  IO.println s!"Found {sols.length} solutions.\n"
  for p in (List.range sols.length).zip sols do
    let i := p.1
    let (whites, blacks) := p.2
    IO.println s!"Solution {i + 1}:"
    printBoard whites blacks

#eval main

Fitting a regular expression to a list of words

Posted on by John

Suppose you want to search for a list of words. If you’re using grep, you can add the -f flag provide a file of regular expressions, and you can add the -F to tell it that the regular expressions are in fact just words. I did something like this a couple days ago when searching for diagnosis codes.

grep -w -F -o -f icd10codes.txt notes.txt

Now you might want to combine your list of words into a singular regular expression, for efficiency or possibly for some other reason. Apparently ripgrep does this because when I tried replacing grep with ripgrep in the command above I got an error saying “Compiled regex exceeds size limit of 104857600 bytes.”

Beating brute force

Say you wanted to search for the strings “bluecross”, “blueshield”, and “bluey”. You could simply form the brute force regular expression

bluecross|blueshied|bluey

but that doesn’t take advantage of the fact that all three strings begin with “blue.” A smaller regular expression would be

blue(shield|cross|y)

Finding the shortest regular expression that matches a list of words is a hard problem, but finding a regular expression that’s shorter than brute force is not. The Python package trieregex will do this. According to the documentation,

trieregex creates efficient regular expressions (regexes) by storing a list of words in a trie structure, and translating the trie into a more compact pattern.

Let’s try our blue example with trieregex.

import re
from trieregex import TrieRegEx as TRE

words = ['bluecross', 'blueshield', 'bluey']
tre = TRE(*words) 
print(tre.regex())

This produces the same regular expression as above, except it adds ?: to make the parentheses non-capturing.

blue(?:shield|cross|y)

Prefixes versus suffixes

The library builds a trie data structure using common prefixes. That works well in the example above, but the result is disappointing when we have common suffixes rather than common prefixes. The following code

words = ['javascript', 'typescript']
tre = TRE(*words) 
print(tre.regex())

produces the regular expression

(?:javascript|typescript)

which is no better than brute force, whereas we might have hoped for

(?:java|type)script

HCPCS examples

As mentioned at the top of the post, ripgrep failed to search on a list of ICD-10 codes. The list of HCPCS codes is about 10x smaller, and more compressible. Ripgrep was able to fit all HCPCS codes into a single regex and was able to search the test file much faster than grep. The command

grep -w -F -o -f hcpcs.txt notes.txt

took 73.426 seconds to execute, while the command

rg -w -F -o -f hcpsc.txt notes.txt

took 0.078 seconds, three orders of magnitude faster.

The following code will read a list of HCPCS codes from a file and create a regular expression.

tre = TRE()
with open('hcpcs.txt', 'r') as file:
    for line in file:
        tre.add(line.strip())
print(len(tre.regex()))

This shows that the resulting regular expression has 17,198 characters. The file of codes has 8725 five-character codes, so the regex compresses the code characters by roughly a ratio of 5 to 2.

Sum of low squares

Posted on by John

Squares, high and low

Let p be an odd prime number. Then half the numbers from 1 through p − 1 are squares and half are not. That is, for half of numbers 1 ≤ k < p, the equation

x² = k mod p

has a solution. The traditional name for these numbers is “quadratic residues” but we can just say “squares” if the context is clear. So, for example, the numbers 1, 2, and 4 are squares mod 7, and the numbers 3, 5, and 6 are not.

If k is a square mod p we will call is a low square if 0 ≤ k < p/2 and a high square if p/2 < k < p.

Signatures

Now let p > 3 be a prime congruent to 3 mod 4. Add up all the low squares mod p and take the remainder mod p. Call this the signature of p. Here’s Python code to make this explicit.

from sympy import isprime, factorint, is_quad_residue

def signature(p):
    assert(p > 3)
    assert(isprime(p))
    assert(p % 4 == 3)
    s = 0
    for k in range(1, 1 + p//2):
        if is_quad_residue(k, p):
            s += k
    return s % p

Inverse signatures

Surprisingly, the signature of each p is unique. Given the signature of p, you can uniquely determine p, and in fact you can do so easily. I ran across this in a paper [1] that presented the results in the form of a parlor trick: have someone pick a prime p such that p = 3 mod 4 and ask them to compute its signature, the sum of the low squares mod p. Then you can quickly tell them what their choice of p was.

Given a signature s, the corresponding prime p is the largest prime factor of 16s + 1.

Not only that,

p = (16s + 1)/m

where m is the smallest of the numbers {3, 7, 11, 15} such that the fraction above is a prime number. In term of Python code, both the following functions should invert the signature of p.

def inverse_signature1(s):
    n = 16*s + 1
    return max(factorint(n).keys())

def inverse_signature2(s):
    n = 16*s + 1
    for m in [3, 7, 11, 15]:
        if n % m == 0 and isprime(n // m):
            return n // m

The following code demonstrates that this is the case for numbers less than 1,000.

for n in range(7, 1000, 4):
    if isprime(n):
        s = signature(n)
        assert(n == inverse_signature1(s))
        assert(n == inverse_signature2(s))

[1] David M. Bloom. A Quadratic Residues Parlor Trick. Mathematics Magazine, Vol. 71, No. 3 (Jun., 1998), pp. 201–203.

Visualizing Medical Code Hierarchy

Posted on by John

Quick follow up to the previous two posts on ICD-10 codes and HCPCS codes. This post uses Python’s squarify library to create treemaps visualizing how many codes begin with each letter.

Here’s the treemap for HCPCS codes.

HCPCS code treemap

And here’s the treemap for ICD-10 codes.

HCPCS code treemap

The sizes of the squares are proportional to the number of codes beginning with that letter. Note that they are not necessarily proportional to how often codes are used.

The HCPCS map omits R and U because these are tiny relative to the rest. The ICD-10 map omits U for the same reason.

Here’s the code that was used to create the HCPCS map.

import matplotlib.pyplot as plt
import squarify

# HCPCS
data = {
    "G": 2010,
    "J": 1232,
    "L": 940,
    "A": 862,
    "E": 671,
    "Q": 639,
    "C": 619,
    "S": 533,
    "M": 506,
    "V": 212,
    "K": 175,
    "T": 114,
    "H": 94,
    "P": 59,
    "B": 51,
  # "U": 5,
  # "R": 3,
}

labels = list(data.keys())
sizes = list(data.values())

# Labels are just the letters (no counts)
display_labels = labels

# Color map — one distinct color per box
colors = plt.cm.tab20.colors[: len(labels)]

fig, ax = plt.subplots(figsize=(12, 8))
squarify.plot(
    sizes=sizes,
    label=display_labels,
    color=colors,
    alpha=0.85,
    ax=ax,
    text_kwargs={"fontsize": 30, "weight": "bold"},
    pad=True,
)
ax.axis("off")

plt.tight_layout()
plt.savefig("treemap.png", dpi=72)
plt.show()

The code to create the ICD-10 map differs only in its data.

# ICD-10
data = {
    "S": 31052, 
    "T": 10090, 
    "M":  6665, 
    "V":  4086, 
    "H":  3330, 
    "O":  2437, 
    "Y":  1590, 
    "I":  1427, 
    "Z":  1411, 
    "W":  1290, 
    "C":  1226, 
    "L":  1000, 
    "E":   971, 
    "Q":   894, 
    "F":   871, 
    "K":   857, 
    "N":   836, 
    "D":   824, 
    "R":   773, 
    "G":   700, 
    "A":   573, 
    "X":   495, 
    "B":   495, 
    "P":   463, 
    "J":   360, 
  # "U":     3,
}

Regular expressions for HCPCS codes

Posted on by John

Since I revisited my old post on ICD code matching, I thought I’d revisit by post on HCPCS codes too.

HCPCS stands for Healthcare Common Procedure Coding System, and is pronounced “hick picks.” When most people say HCPCS, they technically mean HCPCS Level II, and that’s what I mean here.

The format of a HCPCS code is simple: one letter and four digits. In regex terms,

    [A-Z]\d{4}

Not all letters are used, so you can get more specific and say

    [A-CEGHJ-MP-V][0-9]\d{4}

Some sources say no codes begin with U, but there are currently five codes that begin with U.

When I was doing some research on HCPCS codes recently using AI, I was told there is a D code for dentistry, but that was a hallucination.

HCPCS codes can also have modifiers. These consist of a letter and either a letter or digit:

    [A-Z][A-Z0-9]

Not all letters actually appear in modifiers—I, O, W, and Y are missing—so you could be more specific. At the time of writing there are 384 official modifiers.

Modifiers are often stored in a separate column in a database, but in text you’ll see a HCPCS code optionally followed by a dash and a modifier. So a regex to match HCPCS codes with possible modifiers would be

    [A-CEGHJ-MP-V][0-9]\d{4}(-[A-Z][A-Z0-9])?

This regex will have some false positives, but it should not have false negatives: every real HCPCS code should match.

Of course you could search against a complete list of HCPCS codes. This would be more accurate and slower. I did a test similar to the one in the previous post and found a search with the regex above took 20 milliseconds, while a search against the list of HCPCS codes took 46 seconds.

However, the regex searched for possible modifiers and the exhaustive search only looked for unmodified HCPCS codes. A complete list of HCPCS codes with possible modifiers would be tedious to create because some combinations of codes and modifiers make no sense. And I imagine that some combinations that would make sense are not used in practice.