Computing

Hiding data in permutations

Posted on 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

log2(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 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,

CAFEhex = 1100 1010 1111 1110two

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.921fb54442d18hex to binary by converting 1, 9, 2, etc. each to binary.

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

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

π = 11.001001000011111101101010100010001000010110100011000two

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.03ten is

1.000001111010111000010100011110101110000101000111two

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

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.

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.

Regular expression speed and error rates

Posted on 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

  1. Begins with a capital letter, excluding U
  2. Followed by a digit
  3. Followed by a digit or A or B
  4. Optionally followed by a period
  5. 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 + 362 + 363 + 364) = 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.

Related posts

[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 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 J1(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 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 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 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

mk 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.

Related posts

Solving a chess puzzle with Claude and Prolog

Posted on by John

Prolog is the original logic programming language. The name comes from programming in logic. More specifically, the name comes from programmation en logique because the inventor of the language, Philippe Roussel, is French.

Prolog has its advantages and disadvantages. One of the advantages is that the language represents logical problems directly. One of the disadvantages is that the syntax can be quirky. But if an LLM is writing the code, or at least helping to write the code, the syntax doesn’t matter so much.

I wanted to see how well Claude (Sonnet 4.6, medium effort) could solve a chess puzzle by Martin Gardner that I wrote about a little over a year ago. I chose a relatively obscure problem rather than something like the Eight Queens puzzle because an LLM could simply quote one of countless articles on the puzzle.

The puzzle

As I stated in the post last year, the task is to place two rooks, two bishops, and two knights on a 4 by 4 chessboard so that no piece attacks any other.

There are two basic solutions, twelve if you count rotations and reflections as different solutions.

Prolog results

Claude wrote an SWI-Prolog program that I ran with

swipl -g "run, halt" chess_placement.pl

and it gave the following output.

Total solutions (raw): 12
Solutions up to D4 symmetry: 2
Solutions up to D4 symmetry: 2

Solution 1:
B . . N
N . . B
. R . .
. . R .

Solution 2:
. B N .
. . . R
R . . .
. N B .

Claude made a nice image of the solutions.

Note that you can make eight variations of the first solution by a combination of rotations and flips. There are only four variations of the second solution because two rotations put you back where you started.

Prolog source

Here’s the Prolog code that Claude generated.

% ============================================================
%  Non-attacking placement of 2 Rooks, 2 Bishops, 2 Knights
%  on a 4x4 chessboard
%
%  Run with SWI-Prolog:
%    swipl -g "run, halt" chess_placement.pl
%
%  Or interactively:
%    swipl chess_placement.pl
%    ?- run.
%    ?- count_solutions(N).
%    ?- count_unique(N).
% ============================================================

:- use_module(library(lists)).

% ---- Board squares: (Row, Col) with Row, Col in 0..3 ----------------

square(R, C) :- member(R, [0,1,2,3]), member(C, [0,1,2,3]).

% ---- Attack predicates ----------------------------------------------

rook_attacks(R,_,R,_).          % same row
rook_attacks(_,C,_,C).          % same column

bishop_attacks(R1,C1,R2,C2) :-
    DR is abs(R1-R2), DC is abs(C1-C2),
    DR =:= DC, DR > 0.

knight_attacks(R1,C1,R2,C2) :-
    DR is abs(R1-R2), DC is abs(C1-C2),
    ( DR =:= 1, DC =:= 2 ; DR =:= 2, DC =:= 1 ).

% attacks(+Type1, +R1, +C1, +Type2, +R2, +C2)
% True if either piece threatens the other's square.
attacks(rook,   R1,C1, _,     R2,C2) :- rook_attacks(R1,C1,R2,C2).
attacks(_,      R1,C1, rook,  R2,C2) :- rook_attacks(R1,C1,R2,C2).
attacks(bishop, R1,C1, _,     R2,C2) :- bishop_attacks(R1,C1,R2,C2).
attacks(_,      R1,C1, bishop,R2,C2) :- bishop_attacks(R1,C1,R2,C2).
attacks(knight, R1,C1, _,     R2,C2) :- knight_attacks(R1,C1,R2,C2).
attacks(_,      R1,C1, knight,R2,C2) :- knight_attacks(R1,C1,R2,C2).

% ---- All pairs safe -------------------------------------------------

% Check all unordered pairs in a list.
no_attack_pair(_, []).
no_attack_pair(T1-(R1,C1), [T2-(R2,C2)|Rest]) :-
    \+ attacks(T1,R1,C1,T2,R2,C2),
    no_attack_pair(T1-(R1,C1), Rest).

all_pairs_safe([]).
all_pairs_safe([P|Rest]) :-
    no_attack_pair(P, Rest),
    all_pairs_safe(Rest).

% ---- Generate a placement -------------------------------------------
%
%  Placement = [rook-(R1,C1), rook-(R2,C2),
%               bishop-(R3,C3), bishop-(R4,C4),
%               knight-(R5,C5), knight-(R6,C6)]
%
%  We enforce an ordering within each pair to avoid counting
%  same-type duplicates (rook1 vs rook2 swapped).
%  Squares are represented canonically: (R,C) with R*4+C as index.

sq_index(R, C, I) :- I is R*4 + C.

placement(Placement) :-
    % --- Rooks ---
    square(R1r, C1r), square(R2r, C2r),
    sq_index(R1r,C1r,I1r), sq_index(R2r,C2r,I2r),
    I1r < I2r,          % canonical order for the two rooks

    % --- Bishops ---
    square(R1b, C1b), square(R2b, C2b),
    sq_index(R1b,C1b,I1b), sq_index(R2b,C2b,I2b),
    I1b < I2b,          % canonical order for the two bishops

    % --- Knights ---
    square(R1n, C1n), square(R2n, C2n),
    sq_index(R1n,C1n,I1n), sq_index(R2n,C2n,I2n),
    I1n < I2n, % canonical order for the two knights % --- All six squares distinct --- Squares = [(R1r,C1r),(R2r,C2r),(R1b,C1b),(R2b,C2b),(R1n,C1n),(R2n,C2n)], all_distinct_squares(Squares), % --- Build placement list --- Placement = [ rook-(R1r,C1r), rook-(R2r,C2r), bishop-(R1b,C1b), bishop-(R2b,C2b), knight-(R1n,C1n), knight-(R2n,C2n) ], % --- No piece attacks any other --- all_pairs_safe(Placement). all_distinct_squares([]). all_distinct_squares([S|Rest]) :- \+ member(S, Rest), all_distinct_squares(Rest). % ---- Count all solutions -------------------------------------------- count_solutions(N) :- findall(P, placement(P), Ps), length(Ps, N), format("Total solutions (raw): ~w~n", [N]). % ---- D4 symmetry transformations on a 4x4 board --------------------- % Squares (R,C) in 0..3 transform(identity, R, C, R, C). transform(rot90, R, C, C, Nr) :- Nr is 3-R. transform(rot180, R, C, Nr, Nc) :- Nr is 3-R, Nc is 3-C. transform(rot270, R, C, Nc, R) :- Nc is 3-C. transform(flipH, R, C, R, Nc) :- Nc is 3-C. transform(flipV, R, C, Nr, C) :- Nr is 3-R. transform(flipD1, R, C, C, R). transform(flipD2, R, C, Nr, Nc) :- Nr is 3-C, Nc is 3-R. apply_transform(_, [], []). apply_transform(T, [Type-(R,C)|Rest], [Type-(NR,NC)|TRest]) :- transform(T, R, C, NR, NC), apply_transform(T, Rest, TRest). % Canonical form: sort pieces within same-type pairs, then sort the % whole list to get a unique representative. canonical_placement(Placement, Canonical) :- findall(T, member(T,[identity,rot90,rot180,rot270, flipH,flipV,flipD1,flipD2]), Ts), maplist(transform_and_sort(Placement), Ts, AllForms), msort(AllForms, Sorted), Sorted = [Canonical|_]. transform_and_sort(Placement, T, Sorted) :- apply_transform(T, Placement, TPl), msort(TPl, Sorted). % ---- Count solutions up to D4 symmetry ------------------------------ count_unique(N) :- findall(P, placement(P), Ps), maplist(canonical_placement, Ps, Canonicals), list_to_set(Canonicals, Unique), length(Unique, N), format("Solutions up to D4 symmetry: ~w~n", [N]). % ---- Pretty-print a board ------------------------------------------- print_board(Placement) :- forall(member(R, [0,1,2,3]), print_row(R, Placement)), nl. print_row(R, Placement) :- forall(member(C, [0,1,2,3]), print_cell(R, C, Placement)), nl. print_cell(R, C, Placement) :- ( member(rook-(R,C), Placement) -> write('R ')
    ;   member(bishop-(R,C), Placement) -> write('B ')
    ;   member(knight-(R,C), Placement) -> write('N ')
    ;   write('. ')
    ).

% ---- Print all unique solutions -------------------------------------

print_unique_solutions :-
    findall(P, placement(P), Ps),
    maplist(canonical_placement, Ps, Canonicals),
    list_to_set(Canonicals, Unique),
    length(Unique, N),
    format("~nSolutions up to D4 symmetry: ~w~n~n", [N]),
    forall(nth1(I, Unique, Sol),
           ( format("Solution ~w:~n", [I]),
             print_board(Sol) )).

% ---- Top-level entry point ------------------------------------------

run :-
    count_solutions(Raw),
    count_unique(Sym),
    format("~n"),
    print_unique_solutions,
    format("Summary: ~w raw solutions, ~w up to D4 symmetry.~n",
           [Raw, Sym]).