cos(200!)

Posted on 7 August 2026 by John

In a footnote to the previous post, I said that Python’s math library can calculate the logarithm of extremely large numbers but not the cosine. This post will expand on that comment.

In this post I’ll use n = 200! as my example rather than 1000! because this value of N is larger than the largest representable floating point number but small enough to be more convenient to work with.

Suppose someone calculates 200! for you:

78865786736479050355236321393218506229513597768717326329474253324435\
94499634033429203042840119846239041772121389196388302576427902426371\
05061926624952829931113462857270763317237396988943922445621451664240\
25403329186413122742829485327752424240757390324032125740557956866022\
60319041703240623517008587961789222227896237038973747200000000000000\
00000000000000000000000000000000000

You could now calculate log(n) using

n = 7.886578673647905 × 10³⁷⁴

and so

log(n) = log(7.886578673647905 × 10³⁷⁴)
= log(7.886578673647905) + 374 log(10) = 863.2319871924055.

The key thing that makes this possible is that the least significant digits of n only affect the least significant digits of log(n). In the calculation above I kept the first 16 digits of n. Python couldn’t make use of any more digits, and had no need of any more digits, in order to produce the logarithm to machine precision.

Cosine doesn’t work that way. The cosine of n depends on the remainder when n is divided by 2π, and that remainder depends on every single digit of n. I’ll illustrate that below.

Using bc -l and setting the scale to 400, I can calculated n then calculate

cos(n + 10ⁱ)

for i running from 0 to 374, tweaking each digit one at a time. (Except when a digit is a 9 and the addition results in a carry.)

    n = 1
    for (i = 1; i <= 200; i++) n *= i
    scale = 400
    for (i = 1; i <= 374; i++) {
        x = c(n+10^i)
        scale = 16
        print x/1, "\n"
        scale = 400
    }

Here’s what a plot of the results look like.

The value of cos(n) is about −0.985, but the values above are all over the map. We can look at the range by projecting all the points over to the left edge then rotating a quarter turn:

The remarkable thing about this image is that there are a few gaps, i.e. a few values the cosine does not take on.

Here’s a more sophisticated way to look at it. The sequence 10ⁱ mod 2π is dense in [0, 2π], and so by going far enough out in the sequence, we can find a value that shifts the phase of n by any desired amount within any given tolerance.

Every digit in n matters, and changing any digit can change the value of cosine to be essentially any value. You cannot calculate the cosine of an enormous number without using some kind of extended precision arithmetic. There are clever range reduction algorithms that minimize the amount of extended arithmetic necessary, but extended arithmetic cannot be completely eliminated.

Hiding data in permutations

Posted on 27 July 2026 by John

The latest issue of Paged Out! has an article by Stephen Hewitt “An off-line backup of your cryptographic key using playing cards.” The idea is to use a deck of 52 to store a 128-bit cryptographic key. To erase the key, shuffle the deck. Hewitt gives his algorithm for embedding a key, one that can be carried out manually but isn’t maximally efficient.

You could store a 225-bit key as a permutation of 52 cards because

log₂(52!) = 225.581.

But then how would you number permutations so you could go from a number to a particular permutation and later decode the permutation to a number? Is this even practical? For a small number n, you could encode a number k < n by enumerating the first k permutations of a set of n items, and you could decode by enumerating permutations until you find the one you have. But this is completely impractical for large n, such as n = 52.

The process of mapping permutation to an integer is called ranking, and the mapping from an integer to a permutation is called unranking. How efficiently can rankings and unrankings be calculated?

Let n be the number of symbols being permuted. Then there are simple algorithms for ranking and unranking with respect to lexicographical order that have complexity O(n²) and more sophisticated algorithms that have complexity O(n log n). There are also O(n) algorithms that do not preserve lexicographical order.

The Permutations class in SymPy has methods unrank_lex and rank to unrank and rank permutations according to lexicographical order.

The notation the Permutations class uses requires a little explanation. For example, suppose we unrank 2026.

>>> from sympy.combinatorics import Permutation
>>> Permutation.unrank_lex(52, 2026)
Permutation(45, 47, 51, 48, 46, 50)

The output is not a full list of 52 numbers in permuted order; it is only a cycle. The notation refers to the permutation that sends 45 to 47, 47 to 51, …, 50 to 45 and leaves everything else fixed.

If we rank the permutation given above, we get 2026 back.

>>> Permutation.rank(Permutation(45, 47, 51, 48, 46, 50))
2026

Note that we didn’t say how many elements (45, 47, 51, 48, 46, 50) is a permutation of. Because of lexicographical order, the rank would be the same whether we viewed this as a permutation of 52 objects or of more objects.

Now let’s do something larger. Let’s generate a 220-bit number and encode it as a permutation.

>>> n = random.getrandbits(225)
>>> a = Permutation.unrank_lex(52, n)
>>> n
40234719030664563684489051530416964877785781669439875437823431388841
>>> a
Permutation(0, 25, 32, 15, 8, 28)(1, 48, 34, 14, 10, 51, 38, 31, 21, 5, 42, 47, 29, 26, 46, 30, 50, 49, 37, 22, 18, 23)(2, 45, 17, 20, 36, 40, 11, 4, 7, 41, 33, 3, 43, 44, 19, 16, 35, 39, 12, 6, 9)
>>> Permutation.rank(a) == n
True

Now just for fun, let’s display the permutation above applied to a standard (French) deck of 52 cards. As explained here, symbols associated with these cards have a range of Unicode values. By printing these values, we can visualize the permuted deck.

Here’s the code that made the image above.

spades = list(range(0x1F0A1, 0x1F0AF))
spades.remove(0x1F0AC) # take out the knight
cards = [s + 16*i for s in spades for i in range(4)]

a = Permutation.unrank_lex(52, n)
p = a(cards)

for i in range(4):
    for j in range(13):
        print(chr(p[13*i + j]), end="")
    print()

The code above is plenty fast, but Permutation has methods rank_nonlex and unrank_nonlex that run in O(n) time, which could be useful for n much larger than 52.

Printing floating point numbers in binary

Posted on 27 July 2026 by John

It’s well known that you can convert the base 16 (hex) representation of an integer to the base 2 (binary) representation by simply converting each digit from hex to binary. For example,

CAFE_hex = 1100 1010 1111 1110_two

I imagine it’s less well known that you can do the same thing with floating point numbers.

I wanted to find the binary representation of a floating point number using Python, and discovered that it has no function to do this. However, there is a method on floats to show a hex representation. For example, here’s the hex representation of π.

>>> import math
>>> (math.pi).hex()
'0x1.921fb54442d18p+1'

Curiously, the p+k part at the end is an exponent of 2, not an exponent of 16. So after we convert 1.921fb54442d18 to binary, we’ll need to multiply by 2, i.e. move the fractional point one space to the right.

So first we convert 1.921fb54442d18_hex to binary by converting 1, 9, 2, etc. each to binary.

1.1001 0010 0001 1111 1011 0101 0100 0100 0100 0010 1101 0001 1000_two

Then after shifting the fraction point to account for the p+1 part we have

π = 11.001001000011111101101010100010001000010110100011000_two

You could use Python’s bin() function to convert the fractional part, interpreted as an integer, to hex, though you may need to pad with 0 bits. For example,

>>> (1.03).hex()
'0x1.07ae147ae147bp+0

>>> bin(0x7ae147ae147)
'0b1111010111000010100011110101110000101000111'

The binary representation of 1.03_ten is

1.000001111010111000010100011110101110000101000111_two

We added a total of five zero bits, four for the 0 after the fractional point and one for converting 7 to 0111_two.

Fitting a regular expression to a list of words

Posted on 19 July 2026 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.

Regular expressions for HCPCS codes

Posted on 17 July 2026 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.

Regular expression speed and error rates

Posted on 17 July 2026 by John

Seven years ago I wrote a post about regular expressions to match diagnosis codes. I wanted to revisit that post looking at speed and error rates.

Regular expressions usually do not exactly match what you’re looking for and nothing else. They have error false positives and false negatives. But they also have advantages, and context determines whether the advantages make the error rates tolerable.

The post mentioned above gave the following regular expression for ICD-10 diagnosis codes.

    [A-TV-Z][0-9][0-9AB]\.?[0-9A-TV-Z]{0,4}

As cryptic as this may look at first glance, it’s straight-forward to interpret. It says that an ICD-10 code

Begins with a capital letter, excluding U
Followed by a digit
Followed by a digit or A or B
Optionally followed by a period
Followed by up to 4 digits or capital letters, excluding U.

Speed

Now suppose you want to scan a text document for ICD-10 codes. One approach would be to use the regex above. Another would be to compare every alphanumeric sequence in the document to a list of ICD-10. Currently this list has 74,719 codes.

I tested both approaches on a 800kb text file. The regex search

    egrep -o '[A-TV-Z][0-9][0-9AB]\.?[0-9A-TV-Z]{0,4}' notes.txt

took 18 milliseconds. Searching against the list of codes

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

took 386 seconds, about six and a half minutes or five orders of magnitude longer.

Error rates

The regex

    [A-TV-Z][0-9][0-9AB]\.?[0-9A-TV-Z]{0,4}

had a false negative rate of zero at the time it was written. I tested the regex against the current list of codes with the following command.

    egrep -v '[A-TV-Z][0-9][0-9AB]\.?[0-9A-TV-Z]{0,4}' icd10codes.txt

The -v flag reverses the sense of the search, reporting lines that do not match the regular expression. This returned three matches: U070, U071, and U099. So 3 out of 74,719 valid ICD-10 codes were reported as invalid.

Codes beginning with U are designated for provisional/emergency/special purposes, but these three have become essentially permanent. A change in the application of the ICD-10 standard caused an error in the regular expression.

But the change would also have caused an error in code that did an exhaustive search against the list of ICD-10 codes at the time. In fact, every new code not starting with U would also be reported in error. So the regex is actually more future-proof than an exhaustive search. Presumably the simplified regex

    [A-Z][0-9][0-9AB]\.?[0-9A-Z]{0,4}

will remain valid for the foreseeable future.

We’ve looked at false negatives. What about false positives? That depends on context. The false positive rate when searching medical notes is low: a word matching the regex above in a medical record is most likely an ICD-10 code. But the number of conceivable false positives is enormous. If you were searching a file of randomly generated alphanumeric text, the regex matches would overwhelmingly be false positives [1].

The number of strings matching

    [A-Z][0-9][0-9AB]\.?[0-9A-Z]{0,4}

would be

26 × 10 × 12 × (1 + 36 + 36² + 36³ + 36⁴) = 5,390,127,600.

Out of over five billion strings matching the regular expression, only around 75,000 are valid ICD-10 codes. So a naive theoretical calculation would say the false positive rate is 99.9986%, whereas in practice the false positive rate is very low, though there’s no way to say a priori exactly how low.

[1] You could argue that all positives would be false positives in this context because you’re looking at noise. You couldn’t find an ICD code, though you could find a string of characters that coincides with an ICD code. That may sound like a pedantic distinction, but it matters in the context of evaluating deidentification quality: you want to find instances of PHI, not instances of strings that match the character pattern of PHI.

Who you gonna believe: Grok or the docs?

Posted on 29 June 2026 by John

The calculator utility bc has a minimal math library. For example, there’s no tangent function because you’re expected take the ratio of sine and cosine. (The Gnu version of bc does have a function for tangent, but the POSIX version does not.) And yet bc includes support for Bessel functions J(x).

The bc function j takes two arguments. Is the first argument n or x? Grok said the function arguments are j(n,x). I thought I should run man bc just to make sure, and it said

j(x, n) Returns the bessel integer order n (truncated) of x.

So Grok says j(n,x) and the documentation that ships with the software says j(x,n). Which one should you believe? Neither! You should run a little test.

~$ bc -l
>>> j(1, 0)
0
>>> j(0, 1)
.76519768655796655144

Now J₁(0) = 0, so apparently the first argument is the order n. Grok was right and the man page was wrong.

Groucho Marx saysing

As further confirmation, let’s see which argument is truncated.

>>> j(1.2, 3.4)
.17922585168150711099
>>> j(1, 3.4)
.17922585168150711099
>>> j(1.2, 3)
.33905895852593645892

The first argument is truncated to an integer value, so that’s the order n.

Turns out there’s a bug in the man page. The man page text above comes from running man bc on my Macbook. On my Linux box, the documentation is correct. It says

j(n,x) The Bessel function of integer order n of x.

The software produces the same results on both computers. It’s just a documentation bug.

The version running on my Macbook is the version that ships with the OS. It’s not the Gnu version, though the documentation says “This bc is compatible with both the GNU bc and the POSIX bc spec.” It has a function t for tangent, for example, which a POSIX version does not. But if you run bc --standard -l attempting to call t produces an error.

Brace expansion tree

Posted on 27 June 2026 by John

Here’s a crazy bash one-liner I found via an article by Peter Krumins:

echo {w,t,}h{e{n{,ce{,forth}},re{,in,fore,with{,al}}},ither,at}

This prints 30 English words:

when, whence, whenceforth, where, wherein, wherefore, wherewith, wherewithal, whither, what, then, thence, thenceforth, there, therein, therefore, therewith, therewithal, thither, that, hen, hence, henceforth, here, herein, herefore, herewith, herewithal, hither, hat

This post will explain how the one-liner works.

Bash brace expansion iterates through all possibilities listed within curly braces, with possibilities separated by a comma. Note that the comma is a separator and not a terminator. And so, for example, the expression {w,t,} is effectively {w,t,""}.

When bash sees two brace expressions, these expand to the cartesian product of the two expressions. For example,

echo {A,B}{1,2,3}

produces

A1 A2 A3 B1 B2 B3

In the expression above we have

{w,t,}h{e…,ither,at}

So the expansion will enumerate all possibilities of {w,h,} multiplied by all possibilities of {e…,ither,at} where e… is itself a brace expression.

A diagram will help a lot.

The brace expansion does a depth-first traversal of this tree.

Regular expressions that work “everywhere”

Posted on 23 June 2026 by John

The most frustrating aspect of regular expressions is that implementations vary. Features supported in one tool may not be supported at all in another tool, or they may be supported with slightly different syntax.

I learned regular expressions in the context Perl, a maximalist regex environment. This led to frustration when features I expect to work are missing [1]. One way around this is to use Perl analogs of other tools, but this is very non-standard. I want to be able to send colleagues and clients code that works out of the box.

As I mentioned in my post on computational survivalism, I occasionally need to work on computers that I cannot install software on. So a better approach is to identify a subset of regex features that work everywhere. The stricter your definition of “everywhere” the less this includes. The strictest subset would be

literals
character classes […]
the special characters . * ^ $

A more relaxed definition of “everywhere” would be the tools you most care about. Currently the tools I most want to use with regular expressions are sed, awk, grep, and Emacs.

Awk as lowest common denominator

If you use the Gnu versions of sed, awk, and grep, and use the -E option with sed and grep, then the list of common features is bigger. The regular expression features of the three tools are similar, and awk’s features are supported in the other tools, with one exception: word boundaries in awk are \< and \> rather than \b and \B.

I wrote about Awk’s regex features here.

Emacs as the oddball

Emacs supports analogs of most of awk’s regex features. However, the characters

    + ? ( ) { } |

all require a backslash in front in order to act like the awk counterparts. Also, the analog of \s and \S in awk is \s- and \S- in Emacs.

Instead of meaning space or nonspace, \s and \S in Emacs begin a (negated) character class, and one of those classes is - for space. But there are many others. For example, \s. stands for a punctuation character and \S. stands for a non-punctuation character.

What works everywhere

So for my definition of “everywhere,” with the caveats mentioned above, the following features work everywhere. YMMV.

    .
    ^, $
    […], [^…]
    *
    \w, \W, \s, \S
    \1 - \9 backreferences
    \b \B
    ? + 
    | alternation
    {n,m} for counting matches
    (...) capturing

One footnote is that gawk supports backreferences in replacement strings but not in regular expressions per se.

[1] To some extent, basic Perl features work elsewhere and advanced features do not, depending on your idea of what is basic or advanced. I think of look-around features as advanced, and that tracks. But I think of \d for digits as basic, but that’s not supported in many regex flavors.

RSA munitions T-shirt

Posted on 13 June 2026 by John

Back when the US government classified strong encryption as “munitions,” RSA public key cryptography was illegal to export. In 1995, Adam Back protested this by creating a terse, obfuscated implementation of RSA in Perl code and used it as an email signature.

The code was also printed on T-shirts. The shirt was classified as munitions because it contained source code for strong encryption. More on the shirt here.

Adam Back's munitions T-shirt

This was the code:

#!/bin/perl -s-- -export-a-crypto-system-sig -RSA-3-lines-PERL
$m=unpack(H.$w,$m."\0"x$w),$_=`echo "16do$w 2+4Oi0$d*-^1[d2%Sa
2/d0<X+d*La1=z\U$n%0]SX$k"[$m*]\EszlXx++p|dc`,s/^.|\W//g,print
pack('H*',$_)while read(STDIN,$m,($w=2*$d-1+length$n&~1)/2)

My initial intention was to unpack the code, explaining each piece in detail. I don’t have the time or patience for that, and I imagine many readers don’t either. For more of a blow-by-blow explanation, see this commentary from 1995.

dc

In the middle of the code is

    echo ... | dc

This is the most dense and most important part of the code. Perl calls the dc calculator to do the arbitrary precision arithmetic that RSA encryption requires.

I’ve written about bc several times. bc (“basic calculator”) was a originally a more user-friendly wrapper around the reverse-Polish dc (“desktop calculator”). dc is still part of every Unix and Unix-like system, but I imagine bc is far more popular.

The important feature of dc for this post is that it is stack-based, meaning that users would push data and commands on to the stack and pop results off the stack. A sequence of commands that might be understandable when interactively using dc would look cryptic in a transcript. This is part of what makes the code so cryptic.

I’ll parse just a tiny bit of the dc code to give a flavor of what it does. The first four characters 16do instructs dc to push 16 on to the stack, duplicate it, and set the output radix to 16, i.e. these four characters tell dc to work in hexadecimal.

Believe it or not, the dc code is computing

m^k mod n

using fast exponentiation, which is the key step in the RSA algorithm.

Textbook RSA

Note that Adam Back’s code is computing what we would now call textbook RSA, not RSA as it has been refined over the years and is currently implemented.

Computing