Go deh!

From all truths to (ir)relevancies

2025-09-14T12:51:00.003+01:00

Following up on my previous post about truth tables, I now ask a subtler question: which inputs actually matter? Some variables, though present, leave no trace on the output. In this post, I uncover those quiet bits — the irrelevant inputs — and learn how to spot them with precision.

Building on the previous

The previous post showed that for a given number of inputs, there is a finite, but rapidly growing, number of possible truth tables. For two inputs, there are 16. For three inputs, there are 256. This leads to a powerful idea: we can create a standardized format for all truth tables and uniquely identify the parts in each one.

The Standardized Truth Table (STT) Format

i[1] i[0] : r ================= 0 0 : 0 0 1 : 0 1 0 : 0
1 1 : 1

Colour Key

Inputs
Input count
Result, r
Result vector: 0001

We have a standardized truth table with $i$ inputs. Each row in the inputs section of the truth table is an increasing binary count from $0$ to $2**^{i} - 1$ . The result column for the truth table has $2**^{i}$ rows, leading to $2^{**(2**^{i})}$ different possible truth table result vectors for $i$ inputs.

Any one possible result, $r$ , of the standardized truth table is read as the binary bits of the output going from the topmost row (where the inputs are all zero) down to the bottom-most row, in order. This single binary number, the "result number" or vector, uniquely identifies the truth table and the boolean function it represents.

Irrelevant Variables in Boolean Expressions

A variable in a boolean expression is considered irrelevant if its value has no effect on the final output. In a truth table, this means that changing the value of that input variable while all other inputs remain the same does not change the output.

While this is easy to understand, it can be difficult to spot for complex functions. This is where the standardized truth table format, STT, comes in handy.

Irrelevancy Calculation

For an $i$ -input truth table with a $2**^{2**^{i}}$ bit result, $r$ , we can efficiently check for irrelevant variables. The key insight is that for any single input variable, say $i_{[n]}$ , the other input bits change in the exact same order when $i_{[n]} = 0$ as they do when $i[_{n]} = 1$ , in the input count region of the STT.

Therefore, if the result bits for the rows where $i_{[n]} = 0$ are the same as the result bits for the rows where $i_{[n]} = 1$ , then the input variable $i[_{n]}$ has no effect on the output and is therefore irrelevant.

Let's illustrate with an example for a 3-input truth table.

Example

i[2]	i[1]	i[0]	:	Output
0	0	0	:	r0
0	0	1	:	r1
0	1	0	:	r2
0	1	1	:	r3
1	0	0	:	r4
1	0	1	:	r5
1	1	0	:	r6
1	1	1	:	r7

Consider checking if $i [1]$ is irrelevant. The values of the other inputs ( $i [2]$ and $i [0]$ ) follow the sequence 00, 01, 10, 11 both when $i [1]$ is 0 and when $i [1]$ is 1.

When $i [1] = 0$ , the corresponding output bits are $r_{0}, r_{1}, r_{4}, r_{5}$ .
When $i [1] = 1$ , the corresponding output bits are $r_{2}, r_{3}, r_{6}, r_{7}$ .

If the sequence of bits $(r_{0}, r_{1}, r_{4}, r_{5})$ is identical to the sequence of bits $(r_{2}, r_{3}, r_{6}, r_{7})$ , then the input variable $i [1]$ has no effect on the output and is therefore irrelevant.

This method allows for a very efficient, algorithmic approach to simplifying boolean expressions.

from collections import defaultdict
import pprint


def is_irrelevant(considered: int,
                  input_count:int,
                  result_vector: int) -> bool:
    """
    Determine whether a specific input variable is irrelevant to the output of the truth table.
    
    Args:
        considered: Index of the input variable to test.
        input_count: Total number of input variables.
        result_vector: Integer representing the output bits of the truth table.
    
    Returns:
        True if the input is irrelevant, False otherwise.
    """
    considered_to_resultbits = {0: [], 1: []}

    for in_count in range(2**input_count):
        considered_bit = (in_count >> considered) & 1
        resultbit = (result_vector >> in_count) & 1
        considered_to_resultbits[considered_bit].append(resultbit)
    
    return considered_to_resultbits[0] == considered_to_resultbits[1]


def find_irrelevant_ttable_inputs(input_count:int,
                                  result_vector: int) -> list[int]:
    """
    Identify which input variables are irrelevant in a standardized truth table.
    
    Args:
        input_count: Number of input variables.
        result_vector: Integer representing the output bits of the truth table.
    
    Returns:
        List of input indices that are irrelevant.
    """
    irrelevant = [i for i in range(input_count)
                  if is_irrelevant(i, input_count, result_vector)]
    
    return irrelevant

Relevant and irrelevant result vectors for STT's

I am interested in the truthtables/boolean expressions of an increasing number of inputs. Previously I was taking all the zero input STT, then all the 1-input, all the 2-input, ...
That had repetitions and irrelevancies. I can now take just the relevant result vectors for each case, or, take the maximum inputs I can handle and sort the result vectors so that those with the most irrelevancies for that maximum number of inputs, come first.

Here's the code I used to investigate these properties of irrelevances:

def print_STT_results_sensitive_to_inputs(max_input_count: int) -> None:
    """
    Print a summary of irrelevance patterns across all standardized truth tables (STTs)
    for input counts from 0 up to `max_input_count - 1`.

    For each input count `i`, the function:
      - Computes all possible result vectors (2**(2**i))
      - Identifies which inputs are irrelevant for each result
      - Summarizes how many result vectors have at least one irrelevant input
      - Prints spacing patterns (diffs) between result indices with irrelevance
      - Checks whether irrelevance indices for input count `i` begin with those from `i-1`


    Args:
        max_input_count: The exclusive upper bound on input counts to analyze.
    """

    imax2irrelevances = {}
    for imax in range(max_input_count):
        print(f"\nINPUTS = {imax}\n==========")
        result_count = 2**2**imax
        print(f"  total_possible_tt_results =", result_count)

        result2irrelevances = {r:[i for i in range(imax)
                                  if is_irrelevant(i, imax, r)]
                               for r in range(result_count)}
        txt = pprint.pformat(result2irrelevances, compact=True, indent=2).replace('\n', '\n  ')
        txt = txt.replace('{', '{\n   ')
        #print("  result2irrelevances = ", txt)

        irrelevances = [k for k, v in result2irrelevances.items()
                        if v]
        
        imax2irrelevances[imax] = irrelevances

        relevances = result_count - len(irrelevances)
        print(f"  {irrelevances = }")
        print(f"    {len(irrelevances) = }, {relevances = }")

        table = chop_line(format_irrelevance_table_full(result2irrelevances, imax)[0], 220)
        print("\nSTT Result Irrelevance vs Input Table"
              f"\n{table}")

        # First-order differences between irrelevance indices
        diff0 = [irrelevances[j] - irrelevances[j-1] 
                 for j in range(1, len(irrelevances))]
        #print(f"    {diff0 = }")

        # Second-order differences between irrelevance indices
        diff1 = [diff0[j] - diff0[j-1] 
                 for j in range(1, len(diff0))]
        #print(f"    {diff1 = }")
    
    print("\n\nIrrelevance indices reflected about the center of the r count?"
          f"\n  {is_irrelevance_reflected(imax2irrelevances)}")  # True so far

    print('Irrelevances for `i` inputs begin with those from `i-1`?')
    previous_prefixed = all((i0 := imax2irrelevances[i-1]) == imax2irrelevances[i][:len(i0)]
                            for i in range(1, max_input_count))
    print(f"  {previous_prefixed}")  # True so far


def is_irrelevance_reflected(ins2irrel: dict[int, list[int]]) -> bool:
    """
    Check whether irrelevance indices are symmetric about the center of the result vector space.

    For each input count `i`, the result vector space has size 2**(2**i).
    This function checks whether for every irrelevance index `r < half_range`,
    its mirror index `full_range - 1 - r` is also marked as irrelevant.

    Args:
        ins2irrel: Mapping from input count to list of result indices with irrelevance.

    Returns:
        True if all irrelevance sets are symmetric about the center, False otherwise.
    """
    for imax, irrelevances in ins2irrel.items():
        irrel_set = set(irrelevances)
        full_range = 2**2**imax
        half_range = full_range / 2

        if not all((full_range - 1 - i) in irrel_set for i in irrel_set if i < half_range):
            return False
    return True


def format_irrelevance_table(result2irrelevances: dict[int, list[int]], input_count: int) -> tuple[str, str]:
    """
    Create a compact table showing which inputs are irrelevant for each result vector `r`.

    Each column is sized to match the width of its corresponding `r` value.
    Only includes columns for `r` values that have at least one irrelevant input.
    Each row corresponds to an input index, with '@' if that input is irrelevant for that `r`, else blank.

    Returns:
        A tuple of (plain text table string, markdown table string)
    """
    if input_count == 0:
        msg = "No inputs to analyze (input_count = 0)."
        md = "| Input | (none) |\n|-------|--------|\n| (none) |        |"
        return msg, md

    filtered_r = [r for r, irrels in result2irrelevances.items() if irrels]
    if not filtered_r:
        msg = "No irrelevant inputs found for any result vector."
        md = "| Input | (none) |\n|-------|--------|\n" + "\n".join(
            f"| i[{i}] |        |" for i in range(input_count)
        )
        return msg, md

    # Determine individual column widths based on r string length
    r_labels = [str(r) for r in filtered_r]
    col_widths = [len(label) for label in r_labels]

    # Header row
    header_plain = " " * 8 + " ".join(label.rjust(w) for label, w in zip(r_labels, col_widths))
    header_md = "| Input | " + " | ".join(label.rjust(w) for label, w in zip(r_labels, col_widths)) + " |"

    # Markdown separator row
    separator_md = "|-------|" + "|".join("-" * (w + 2) for w in col_widths) + "|"

    # Rows
    rows_plain = []
    rows_md = []
    for i in range(input_count):
        label = f"i[{i}]".ljust(8)
        row_plain = label + " ".join("@" .rjust(w) if i in result2irrelevances[r] else " " * w
                                     for r, w in zip(filtered_r, col_widths))
        row_md = "| " + f"i[{i}]".ljust(5) + " | " + " | ".join("@" .rjust(w) if i in result2irrelevances[r] else " " * w
                                                               for r, w in zip(filtered_r, col_widths)) + " |"
        rows_plain.append(row_plain)
        rows_md.append(row_md)

    plain_table = "\n".join([header_plain] + rows_plain)
    markdown_table = "\n".join([header_md, separator_md] + rows_md)

    return plain_table, markdown_table

def format_irrelevance_table_full(result2irrelevances: dict[int, list[int]], input_count: int) -> tuple[str, str]:
    """
    Create a full table showing which inputs are irrelevant for each result vector `r`.

    Includes columns for all result vector indices (0 to max(r)).
    Each row corresponds to an input index, with '@' if that input is irrelevant for that `r`, else blank.

    Returns:
        A tuple of (plain text table string, markdown table string)
    """
    if input_count == 0:
        msg = "No inputs to analyze (input_count = 0)."
        md = "| Input | (none) |\n|-------|--------|\n| (none) |        |"
        return msg, md

    all_r = sorted(result2irrelevances.keys())
    r_labels = [str(r) for r in all_r]
    col_widths = [len(label) for label in r_labels]

    # Header row
    header_plain = " " * 8 + " ".join(label.rjust(w) for label, w in zip(r_labels, col_widths))
    header_md = "| Input | " + " | ".join(label.rjust(w) for label, w in zip(r_labels, col_widths)) + " |"

    # Markdown separator row
    separator_md = "|-------|" + "|".join("-" * (w + 2) for w in col_widths) + "|"

    # Rows
    rows_plain = []
    rows_md = []
    for i in range(input_count):
        label = f"i[{i}]".ljust(8)
        row_plain = label + " ".join("@" .rjust(w) if i in result2irrelevances.get(r, []) else " " * w
                                     for r, w in zip(all_r, col_widths))
        row_md = "| " + f"i[{i}]".ljust(5) + " | " + " | ".join("@" .rjust(w) if i in result2irrelevances.get(r, []) else " " * w
                                                               for r, w in zip(all_r, col_widths)) + " |"
        rows_plain.append(row_plain)
        rows_md.append(row_md)

    plain_table = "\n".join([header_plain] + rows_plain)
    markdown_table = "\n".join([header_md, separator_md] + rows_md)

    return plain_table, markdown_table

def chop_line(table_text: str, max_length: int) -> str:
    """
    Chop each line in a multi-line string to a maximum length.
    If a line exceeds `max_length`, it is truncated to `max_length - 3` and '...' is appended.

    Args:
        table_text: The full multi-line string to process.
        max_length: The maximum allowed line length.

    Returns:
        A new multi-line string with long lines chopped and marked with '...'.
    """
    chopped_lines = []
    for line in table_text.splitlines():
        if len(line) > max_length:
            chopped_lines.append(line[:max_length - 3] + '...')
        else:
            chopped_lines.append(line)
    return '\n'.join(chopped_lines)


print_STT_results_sensitive_to_inputs(5)

Irrelevances: Output

INPUTS = 0
==========
  total_possible_tt_results = 2
  irrelevances = []
    len(irrelevances) = 0, relevances = 2

STT Result Irrelevance vs Input Table
No inputs to analyze (input_count = 0).

INPUTS = 1
==========
  total_possible_tt_results = 4
  irrelevances = [0, 3]
    len(irrelevances) = 2, relevances = 2

STT Result Irrelevance vs Input Table
        0 1 2 3
i[0]    @     @

INPUTS = 2
==========
  total_possible_tt_results = 16
  irrelevances = [0, 3, 5, 10, 12, 15]
    len(irrelevances) = 6, relevances = 10

STT Result Irrelevance vs Input Table
        0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
i[0]    @     @                    @        @
i[1]    @         @          @              @

INPUTS = 3
==========
  total_possible_tt_results = 256
  irrelevances = [0, 3, 5, 10, 12, 15, 17, 34, 48, 51, 60, 63, 68, 80, 85, 90, 95, 102, 119, 136, 153, 160, 165, 170, 175, 187, 192, 195, 204, 207, 221, 238, 240, 243, 245, 250, 252, 255]
    len(irrelevances) = 38, relevances = 218

STT Result Irrelevance vs Input Table
        0 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 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 ...
i[0]    @     @                    @        @                                                                                                  @        @                          @        @                            ...
i[1]    @         @          @              @                                                                                                                                                                            ...
i[2]    @                                         @                                                  @                                                  @                                                  @             ...

INPUTS = 4
==========
  total_possible_tt_results = 65536
  irrelevances = [0, 3, 5, 10, 12, 15, 17, 34, 48, 51, 60, 63, 68, 80, 85, 90, 95, 102, 119, 136, 153, 160, 165, 170, 175, 187, 192, 195, 204, 207, 221, 238, 240, 243, 245, 250, 252, 255, 257, 514, 768, 771, 780, 783, 816, 819, 828, 831, 960, 963, 972, 975, 1008, 1011, 1020, 1023, 1028, 1280, 1285, 1290, 1295, 1360, 1365, 1370, 1375, 1440, 1445, 1450, 1455, 1520, 1525, 1530, 1535, 1542, 1799, 2056, 2313, 2560, 2565, 2570, 2575, 2640, 2645, 2650, 2655, 2720, 2725, 2730, 2735, 2800, 2805, 2810, 2815, 2827, 3072, 3075, 3084, 3087, 3120, 3123, 3132, 3135, 3264, 3267, 3276, 3279, 3312, 3315, 3324, 3327, 3341, 3598, 3840, 3843, 3845, 3850, 3852, 3855, 3888, 3891, 3900, 3903, 3920, 3925, 3930, 3935, 4000, 4005, 4010, 4015, 4032, 4035, 4044, 4047, 4080, 4083, 4085, 4090, 4092, 4095, 4112, 4352, 4369, 4386, 4403, 4420, 4437, 4454, 4471, 4488, 4505, 4522, 4539, 4556, 4573, 4590, 4607, 4626, 4883, 5140, 5397, 5654, 5911, 6168, 6425, 6682, 6939, 7196, 7453, 7710, 7967, 8224, 8481, 8704, 8721, 8738, 8755, 8772, 8789, 8806, 8823, 8840, 8857, 8874, 8891, 8908, 8925, 8942, 8959, 8995, 9252, 9509, 9766, 10023, 10280, 10537, 10794, 11051, 11308, 11565, 11822, 12079, 12288, 12291, 12300, 12303, 12336, 12339, 12348, 12351, 12480, 12483, 12492, 12495, 12528, 12531, 12540, 12543, 12593, 12850, 13056, 13059, 13068, 13071, 13073, 13090, 13104, 13107, 13116, 13119, 13124, 13141, 13158, 13175, 13192, 13209, 13226, 13243, 13248, 13251, 13260, 13263, 13277, 13294, 13296, 13299, 13308, 13311, 13364, 13621, 13878, 14135, 14392, 14649, 14906, 15163, 15360, 15363, 15372, 15375, 15408, 15411, 15420, 15423, 15552, 15555, 15564, 15567, 15600, 15603, 15612, 15615, 15677, 15934, 16128, 16131, 16140, 16143, 16176, 16179, 16188, 16191, 16320, 16323, 16332, 16335, 16368, 16371, 16380, 16383, 16448, 16705, 16962, 17219, 17408, 17425, 17442, 17459, 17476, 17493, 17510, 17527, 17544, 17561, 17578, 17595, 17612, 17629, 17646, 17663, 17733, 17990, 18247, 18504, 18761, 19018, 19275, 19532, 19789, 20046, 20303, 20480, 20485, 20490, 20495, 20560, 20565, 20570, 20575, 20640, 20645, 20650, 20655, 20720, 20725, 20730, 20735, 20817, 21074, 21331, 21588, 21760, 21765, 21770, 21775, 21777, 21794, 21811, 21828, 21840, 21845, 21850, 21855, 21862, 21879, 21896, 21913, 21920, 21925, 21930, 21935, 21947, 21964, 21981, 21998, 22000, 22005, 22010, 22015, 22102, 22359, 22616, 22873, 23040, 23045, 23050, 23055, 23120, 23125, 23130, 23135, 23200, 23205, 23210, 23215, 23280, 23285, 23290, 23295, 23387, 23644, 23901, 24158, 24320, 24325, 24330, 24335, 24400, 24405, 24410, 24415, 24480, 24485, 24490, 24495, 24560, 24565, 24570, 24575, 24672, 24929, 25186, 25443, 25700, 25957, 26112, 26129, 26146, 26163, 26180, 26197, 26214, 26231, 26248, 26265, 26282, 26299, 26316, 26333, 26350, 26367, 26471, 26728, 26985, 27242, 27499, 27756, 28013, 28270, 28527, 28784, 29041, 29298, 29555, 29812, 30069, 30326, 30464, 30481, 30498, 30515, 30532, 30549, 30566, 30583, 30600, 30617, 30634, 30651, 30668, 30685, 30702, 30719, 30840, 31097, 31354, 31611, 31868, 32125, 32382, 32639, 32896, 33153, 33410, 33667, 33924, 34181, 34438, 34695, 34816, 34833, 34850, 34867, 34884, 34901, 34918, 34935, 34952, 34969, 34986, 35003, 35020, 35037, 35054, 35071, 35209, 35466, 35723, 35980, 36237, 36494, 36751, 37008, 37265, 37522, 37779, 38036, 38293, 38550, 38807, 39064, 39168, 39185, 39202, 39219, 39236, 39253, 39270, 39287, 39304, 39321, 39338, 39355, 39372, 39389, 39406, 39423, 39578, 39835, 40092, 40349, 40606, 40863, 40960, 40965, 40970, 40975, 41040, 41045, 41050, 41055, 41120, 41125, 41130, 41135, 41200, 41205, 41210, 41215, 41377, 41634, 41891, 42148, 42240, 42245, 42250, 42255, 42320, 42325, 42330, 42335, 42400, 42405, 42410, 42415, 42480, 42485, 42490, 42495, 42662, 42919, 43176, 43433, 43520, 43525, 43530, 43535, 43537, 43554, 43571, 43588, 43600, 43605, 43610, 43615, 43622, 43639, 43656, 43673, 43680, 43685, 43690, 43695, 43707, 43724, 43741, 43758, 43760, 43765, 43770, 43775, 43947, 44204, 44461, 44718, 44800, 44805, 44810, 44815, 44880, 44885, 44890, 44895, 44960, 44965, 44970, 44975, 45040, 45045, 45050, 45055, 45232, 45489, 45746, 46003, 46260, 46517, 46774, 47031, 47288, 47545, 47802, 47872, 47889, 47906, 47923, 47940, 47957, 47974, 47991, 48008, 48025, 48042, 48059, 48076, 48093, 48110, 48127, 48316, 48573, 48830, 49087, 49152, 49155, 49164, 49167, 49200, 49203, 49212, 49215, 49344, 49347, 49356, 49359, 49392, 49395, 49404, 49407, 49601, 49858, 49920, 49923, 49932, 49935, 49968, 49971, 49980, 49983, 50112, 50115, 50124, 50127, 50160, 50163, 50172, 50175, 50372, 50629, 50886, 51143, 51400, 51657, 51914, 52171, 52224, 52227, 52236, 52239, 52241, 52258, 52272, 52275, 52284, 52287, 52292, 52309, 52326, 52343, 52360, 52377, 52394, 52411, 52416, 52419, 52428, 52431, 52445, 52462, 52464, 52467, 52476, 52479, 52685, 52942, 52992, 52995, 53004, 53007, 53040, 53043, 53052, 53055, 53184, 53187, 53196, 53199, 53232, 53235, 53244, 53247, 53456, 53713, 53970, 54227, 54484, 54741, 54998, 55255, 55512, 55769, 56026, 56283, 56540, 56576, 56593, 56610, 56627, 56644, 56661, 56678, 56695, 56712, 56729, 56746, 56763, 56780, 56797, 56814, 56831, 57054, 57311, 57568, 57825, 58082, 58339, 58596, 58853, 59110, 59367, 59624, 59881, 60138, 60395, 60652, 60909, 60928, 60945, 60962, 60979, 60996, 61013, 61030, 61047, 61064, 61081, 61098, 61115, 61132, 61149, 61166, 61183, 61423, 61440, 61443, 61445, 61450, 61452, 61455, 61488, 61491, 61500, 61503, 61520, 61525, 61530, 61535, 61600, 61605, 61610, 61615, 61632, 61635, 61644, 61647, 61680, 61683, 61685, 61690, 61692, 61695, 61937, 62194, 62208, 62211, 62220, 62223, 62256, 62259, 62268, 62271, 62400, 62403, 62412, 62415, 62448, 62451, 62460, 62463, 62708, 62720, 62725, 62730, 62735, 62800, 62805, 62810, 62815, 62880, 62885, 62890, 62895, 62960, 62965, 62970, 62975, 63222, 63479, 63736, 63993, 64000, 64005, 64010, 64015, 64080, 64085, 64090, 64095, 64160, 64165, 64170, 64175, 64240, 64245, 64250, 64255, 64507, 64512, 64515, 64524, 64527, 64560, 64563, 64572, 64575, 64704, 64707, 64716, 64719, 64752, 64755, 64764, 64767, 65021, 65278, 65280, 65283, 65285, 65290, 65292, 65295, 65297, 65314, 65328, 65331, 65340, 65343, 65348, 65360, 65365, 65370, 65375, 65382, 65399, 65416, 65433, 65440, 65445, 65450, 65455, 65467, 65472, 65475, 65484, 65487, 65501, 65518, 65520, 65523, 65525, 65530, 65532, 65535]
    len(irrelevances) = 942, relevances = 64594

STT Result Irrelevance vs Input Table
        0 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 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 ...
i[0]    @     @                    @        @                                                                                                  @        @                          @        @                            ...
i[1]    @         @          @              @                                                                                                                                                                            ...
i[2]    @                                         @                                                  @                                                  @                                                  @             ...
i[3]    @                                                                                                                                                                                                                ...


Irrelevance indices reflected about the center of the r count?
  True
Irrelevances for `i` inputs begin with those from `i-1`?
  True

OEIS

The sequence of relevances: 2, 2, 10, 218, 64594 is already present as A000371 on OEIS.

The sequences of irelevances: 0, 2, 6, 38, 942 had no exact match although A005530 came close

END.

Go Deh celebrations!

2025-08-10T17:22:00.001+01:00

This year it's:

Fifty years since I first learned to program!
Thirty years of programming Python!
And In July, this blog topped one million views since its inception!

I first learned to program by crashing an after-school class I saw the private school geeks disappearing to from Central library in Nottingham. One of the Trent Polytechnic lecturers I now see sometimes on Computerphile - he mentioned how the council payed something towards the cost of the computer on the understanding that they have a class for students, but somehow only the guys from the fee paying Nottingham High School for Boys seemed to get the memo, hmm.
That was s long ago, Paper tape and teletypes.

Forward to 1995: A new job, amongst engineers where the few of them that could program, programmed in Perl. Perl was big within the chip design industry back then, but Aargh! I thought Perl was awful - for me it did not gel! I went looking for an alternative and found this little known language called Python.
The community on comp.lang.python seemed so nice, they valued readability and "One obvious way to do it" over Perls more than one way. I was hooked, and did my little bit of evangelizing. (OK, I was known as that Python guy - but they couldn't fault the code).

Was it eleven years ago that I started the blog? It has been one of two main places where I have shown my work - the other is on site rosettacode.org which, no doubt, has been scraped to death for AI's - it is one of the few places online showing different programming languages solving the same task. I wrote 199 of those RosettaCode tasks over the years, and many more of the examples solving those tasks; (mainly Python and Awk language examples). Most, if not all of those task descriptions were my own work - not copied off other sites - if there are similarities in text then they probably copied from me or the task is so simple that there is going to be convergence).

Lately, as you can tell from the pictures above, I have been dabbling in AI. AI is the future, but who knows how AI will affect that future? I am a programmer who can describe coding tasks to other people for them to complete the task. AI is an ever changing field, I think that I am learning to use it better, as it gets better too - I am both more and less precise in my prompting and finding out what the AI remembers and does best. I experiment with the amount of description in the code to cut and paste into new sessions to cut the mistakes and continue from a clean base. I experiment with coding from pseudo-code. I have found that some abilities are universal: when exploring the use of SQL - which I am a hunt-n-peck coder in, I can still spot redundancies "You recompute the first query as part of the second and the second as part of the third, can't you reuse prior results", lead to the AI introducing common table expressions I think it was, and telling me all about them.

I'm still scripting 😊

All Truth in Truthtables!

2025-07-20T10:28:00.000+01:00

(Best viewed on a larger than phone screen)

To crib from my RosettaCode tasks description and examples:

A truth table is a display of the inputs to, and the output of a Boolean equation organised as a table where each row gives one combination of input values and the corresponding value of the equation.

And as examples:

Boolean expression: A ^ B

A B : A ^ B
0 0 : 0
0 1 : 1
1 0 : 1
1 1 : 0

Boolean expression: S | ( T ^ U )

S T U : S | ( T ^ U )
0 0 0 : 0
0 0 1 : 1
0 1 0 : 1
0 1 1 : 0
1 0 0 : 1
1 0 1 : 1
1 1 0 : 1
1 1 1 : 1

Format

A truth table has a header row of columns showing first the names of inputs assigned to each column; a visual separator - e.g. ':'; then the column name for the output result.

The body of the table, under the inputs section, contains rows of all binary combinations of the inputs. It is usually arranged as each row of the input section being a binary count from zero to 2**input_count - 1

The body of the table, under the result section, contains rows showing the binary output produced from the input configuration in the same row, to the left.

Format used

I am interested in the number of inputs rather than their names so will show vector i with the most significant indices to the left, (so the binary count in the input sections body looks right).

Similarly I am interested in the bits in the result column rather than a name so will just call the result column r.

From one result to many

Here's the invocation, and truth tables produced for some simple boolean operators:

print("\nOR\n")
pp_ttable(2, 0b1110)

OR i[1] i[0] : r ================= 0 0 : 0 0 1 : 1 1 0 : 1 1 1 : 1

print("\nXOR\n")
pp_ttable(2, 0b110)

XOR i[1] i[0] : r ================= 0 0 : 0 0 1 : 1 1 0 : 1 1 1 : 0

print("\nAND\n")
pp_ttable(2, 0b1000)

AND i[1] i[0] : r ================= 0 0 : 0 0 1 : 0 1 0 : 0 1 1 : 1

For those three inputs, we can extend the table to show result columns for OR, XOR and then AND, like this:

print("\nOR, XOR, then AND result *columns*\n")
pp_ttable(2, [14, 6, 8])  # A list of results, (in decimal this time)

OR, XOR, then AND result *columns* i[1] i[0] : r[0] r[1] r[2] =========================== 0 0 : 0 0 0 0 1 : 1 1 0 1 0 : 1 1 0 1 1 : 1 0 1

All Truth

Just how many results are possible?

Well, i = 2 inputs gives 2**i = 4 possible input boolean combinations; so a result column has 2**i = 4 bits.
The number of different result columns is therefore 2**(2**i) = 2**4 = 16

We can show all possible results by successive results being a binary count, but this time by column in the results section, (with the LSB being closest to the header row)

The pp_table function automatically generates all possible results if a second parameter of None is used

print("\nAll Truths of two inputs!\n")
pp_ttable(2, None)

All Truths of two inputs! i[1] i[0] : r[0] r[1] r[2] r[3] r[4] r[5] r[6] r[7] r[8] r[9] r[10] r[11] r[12] r[13] r[14] r[15] ============================================================================================================== 0 0 : 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 : 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 1 0 : 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 1 1 : 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1

We might say that it shows all possible truths for up-to-and-including two inputs. That is because results include outputs not dependent on any or all of those two inputs. For example r[0], and r[15] do not depend on any input as they give constant outputs of 0 and 1, respectively. r[5] is simply ~i[0] and does not depend on i[1]

It's Big, Oh!

The results grow as 2**(2**i), sometimes called double exponential growth! it gets large, quickly!!

i | 2**(2**i) --|---------- 0 | 2 1 | 4 2 | 16 3 | 256 4 | 65_536 5 | 4_294_967_296

The code

Contemplating using AI and needing to get it to understand what I wanted, as well as endless prompting to get it to do what I want, the way I wanted it; I decided on writing it all by myself - I knew I wanted it just so, and the coding would be nothing new to me, just nailing what I wanted to show.

"""
Truthtable prettyprinter

Author: Paddy3118  2025/07/19


# Regions of a truthtable

```text
+----------------- --+---------+
|      Inputs        | Result  |
+------------------- +---------+
| I[2] | I[1] | I[0] |   R     |
|------|------|------|---------|
|  0   |  0   |  0   |   0     |
|  0   |  0   |  1   |   1     |
|  0   |  1   |  0   |   0     |
|  0   |  1   |  1   |   1     |
|  1   |  0   |  0   |   1     |
|  1   |  0   |  1   |   0     |
|  1   |  1   |  0   |   1     |
|  1   |  1   |  1   |   0     |
+--------------------+---------+
```

"""

# %%
def pp_ttable(input_count:int,
              result_vector: int | list[int] | None,
              input_vector_name: str='i',
              result_vector_name: str='r'
              ) -> None:
    r_single = r_mult = False
    match result_vector:
        case int():     # single result
            r_single = True
            r_width = len(f" {result_vector_name}")
            result_vector = [result_vector]
        case list():    # multiple results
            r_mult = True
            r_width = len(f" {result_vector_name}[{len(result_vector) - 1}]")
        case None:      # All possible results
            assert input_count < 4, f"Refuse to write {2**(2**input_count)} result columns."
            r_mult = True
            result_vector = list(range(2**(2**input_count)))
            r_width = len(f" {result_vector_name}[{len(result_vector) - 1}]")
        case _:
            raise Typeerror("result must be an int, a list of ints, or None.")

    col_width = max(len(f" {input_vector_name}[{input_count - 1}]"), r_width)

    rows = []
    header = []
    for i in range(input_count-1, -1, -1):
        col = f"{input_vector_name}[{i:}]"
        header.append(f"{col:{col_width}}")
    header.append(': ')
    if r_single:
        header.append(f"{result_vector_name:^{col_width}}")
    else:
        for i in range(len(result_vector)):
            col = f"{result_vector_name}[{i:}]"
            header.append(f"{col:^{col_width}}")
    rows.append(''.join(header))
    rows.append('=' * len(rows[-1]))    # header row spacer

    result_bits = [f"{res:0{2**input_count}b}"[::-1]
                    for res in result_vector]
    
    for r in range(2**input_count):
        row = []
        col_bits = f"{r:0{input_count}b}"
        for bit in col_bits:
            row.append(f"{bit:^{col_width}}")
        
        row.append(': ')

        for res_bits in result_bits:
            row.append(f"{res_bits[r]:^{col_width}}")
                
        rows.append(''.join(row))


    print('\n'.join(rows))

# %%
print("\nOR\n")
pp_ttable(2, 0b1110)
print("\nXOR\n")
pp_ttable(2, 0b110)
print("\nAND\n")
pp_ttable(2, 0b1000)
    

# %%

print("\nOR, XOR, then AND result *columns*\n")
pp_ttable(2, [14, 6, 8])  # A list of results, (in decimal this time)

# %%
print("\nAll Truths of two inputs!\n")
pp_ttable(2, None)


# %%

END.

Incremental combinations without caching

2025-03-23T08:54:00.001+00:00

Irie server room

Someone had a problem where they received initial data d1, worked on all r combinations of the data initially received, but by the time they had finished that, they checked and found there was now extra data d2, and they need to, in total, process the r combinations of all data d1+d2.

They don't want to process combinations twice, and the solutions given seemed to generate and store the combinations of d1, then generate the combinations of d1+d2 but test and reject any combination that was found previously.

Seems like a straight-forward answer that is easy to follow, but I thought: Is there a way to create just the extra combinations but without storing all the combinations from before?

My methods

It's out with Vscode as my IDE. I'll be playing with a lot of code that will end up deleted, modified, rerun. I could use a Jupyter notebook, but I can't publish them to my blog satisfactorily. I'll develop a .py file but with cells: a line comment of # %% visually splits the file into cells in the IDE adding buttons and editor commands to execute cells and selected code, in any order, on a restartable kernel running in an interactive window that also runs Ipython.

When doodling like this, I often create long lines, spaced to highlight comparisons between other lines. I refactor names to be concise at the time, as what you can take in at a glance helps find patterns. Because of that, this blog post is not written to be read on the small screens of phones.

So, its combinations; binomials, nCr.

A start

First thing was to try and see patterns in the combinations of d1+d2 minus those of just d1.
Dealing with sets of values; sets are unordered, so will at some time need a function to print them in order to aid pattern finding.
Combinations can be large - use combinations of ints then later work with any type.
Initial combinations of 0 to n-1 ints I later found to be more awkward to reason about so changed to work with combinations of 1 to n ints. extending to n+x added ints n+1 to n+x to combinations.

#!/bin/env python3
"""
Check patterns in extending combinations nCr
"""

# %%
from functools import reduce
from itertools import combinations
from pprint import pformat
from typing import Any, Iterable

# %%
def nCr(n, r):
    return combinations(range(1, n+1), r=r)

def pf_set(s: set) -> str:
    "Format set with sorted elements for printing"
    return  f"{{{pformat(sorted(s), width=120, compact=True)[1:-1]}}}"

Diffs

In the following cell I create combinations for r = 3, and n in some range n_range in c and print successive combinations, and differences between successive combinations.

# %%

r = 3
n_range = range(3, 7)
print(f"\n# Investigate combinations nCr for {r=}, and {n_range=}\n")

# generate some combinations
c = {n: set(nCr(n, r)) for n in n_range}

print("c = {")
for key, val in c.items():
    print(f" {key:2}: {pf_set(val)},")
print("    }")


def pp_diffsby(c, n_range, r, delta_n):
    "Print nCr in c diffs by delta_n"
    print(f"\nDiffs by {delta_n}")
    all_n = list(n_range)
    n1 = all_n[0]
    print(f"      {n1}C{r} = {pf_set(c[n1])}")
    for n1, n2 in zip(all_n, all_n[delta_n:]):
        print(f"{n2}C{r} - {n1}C{r} = {pf_set(c[n2] - c[n1])}")


pp_diffsby(c, n_range, r, 1)
pp_diffsby(c, n_range, r, 2)

Cell output:

# Investigate combinations nCr for r=3, and n_range=range(3, 7)

c = {
  3: {(1, 2, 3)},
  4: {(1, 2, 3), (1, 2, 4), (1, 3, 4), (2, 3, 4)},
  5: {(1, 2, 3), (1, 2, 4), (1, 2, 5), (1, 3, 4), (1, 3, 5), (1, 4, 5), (2, 3, 4), (2, 3, 5), (2, 4, 5), (3, 4, 5)},
  6: {(1, 2, 3), (1, 2, 4), (1, 2, 5), (1, 2, 6), (1, 3, 4), (1, 3, 5), (1, 3, 6), (1, 4, 5), (1, 4, 6), (1, 5, 6),
 (2, 3, 4), (2, 3, 5), (2, 3, 6), (2, 4, 5), (2, 4, 6), (2, 5, 6), (3, 4, 5), (3, 4, 6), (3, 5, 6), (4, 5, 6)},
    }

Diffs by 1
      3C3 = {(1, 2, 3)}
4C3 - 3C3 = {(1, 2, 4), (1, 3, 4), (2, 3, 4)}
5C3 - 4C3 = {(1, 2, 5), (1, 3, 5), (1, 4, 5), (2, 3, 5), (2, 4, 5), (3, 4, 5)}
6C3 - 5C3 = {(1, 2, 6), (1, 3, 6), (1, 4, 6), (1, 5, 6), (2, 3, 6), (2, 4, 6), (2, 5, 6), (3, 4, 6), (3, 5, 6), (4, 5, 6)}

Diffs by 2
      3C3 = {(1, 2, 3)}
5C3 - 3C3 = {(1, 2, 4), (1, 2, 5), (1, 3, 4), (1, 3, 5), (1, 4, 5), (2, 3, 4), (2, 3, 5), (2, 4, 5), (3, 4, 5)}
6C3 - 4C3 = {(1, 2, 5), (1, 2, 6), (1, 3, 5), (1, 3, 6), (1, 4, 5), (1, 4, 6), (1, 5, 6), (2, 3, 5), (2, 3, 6), (2, 4, 5),
 (2, 4, 6), (2, 5, 6), (3, 4, 5), (3, 4, 6), (3, 5, 6), (4, 5, 6)}

Patterns

Looking at diffs 4C3 - 3C3 each tuple is like they took 3C2 = {(1,2), (1,3), (2,3)} and tagged the extra 4 on to every inner tuple.
Lets call this modification extending,

5C3 - 4C3 seems to follow the same pattern.

Function extend

# %%

def extend(s: set[tuple[Any]], value: Any) -> set[tuple[Any]]:
    """
    Returns set of tuples of s with each tuple extended by value
    """
    return set((*tpl, value) for tpl in s)

s = {(0, 1), (1, 2), (0, 2)}
print(f"{s = }; {extend(s, 3) = }")
assert extend(s, 3) == {(0, 1, 3), (1, 2, 3), (0, 2, 3)}

Rename nCr to bino and check extend works

nCr was originally working with 0..n-1 and bino was 1..n. Now they both do

# %%

# binomial combinations of ints 1..

def bino(n: int, r: int) -> set[tuple[int]]:
    """
    All combinations of 1..n ints taken r at a time

    bino(4, 3) == {(1, 2, 3), (1, 2, 4), (1, 3, 4), (2, 3, 4)}
    """
    return set(combinations(range(1, n+1), r))


print(f"{(bino(4, 3) == (bino(3, 3) | extend(bino(3, 2), value=4)))  =  }")

Cell output:

(bino(4, 3) == (bino(3, 3) | extend(bino(3, 2), value=4)))  =  True

Pascal

After finding that pattern I went searching for it using Gemini AI. My question was:

show that comb(n+1, r) = comb(n, r) + (n+1)* comb(n, r-1)

The answer said I got my text prompt wrong and mentioned Pascals Rule.

I scanned the page, as I no time for how things were expressed but it seemed reasonable that I had the algorithm right, and that there were relations of some kind for bigger differences in n.

Pascals rule checker

I wrote a function to do the check then exercised it (a few times, now deleted).

# %%
def pascals_rule(n: int, r: int) -> bool:
    "check C(n+1, r) == C(n, r) | extend(C(n, r-1), n + 1)"
    return bino(n + 1, r) == bino(n, r) | extend(bino(n, r - 1), n + 1)

assert pascals_rule(6, 3)

Diff by 1 extra item "done", attempting diff by 2.

Looking back at the diffs by 2 table and looking for patterns I thought I might need different types of extension functions modifying the tuples within sets in different ways - it seemed "mathematical" so...

# %%

# Some functions that may be needed

def extend_mul(s: set[tuple[Any]], e: Iterable) -> set[tuple[Any]]:
    """
    set where each tuple of s is extended in turn, by every item in e

    s = {(1, 2, 3), (1, 2, 4), (1, 3, 4), (2, 3, 4)}
    extend_mul(s, (5, 6)) == {(1, 2, 3, 5), (1, 2, 3, 6),
                              (1, 2, 4, 5), (1, 2, 4, 6),
                              (1, 3, 4, 5), (1, 3, 4, 6),
                              (2, 3, 4, 5), (2, 3, 4, 6)}
    """
    return {(*t, item) for t in s for item in e}

def extend_muli(s: set[tuple[Any]], e: Iterable) -> set[tuple[Any]]:
    """
    set where each tuple of s is extended in turn, by every *item in e

    s = {(1, 2, 3), (1, 2, 4), (1, 3, 4), (2, 3, 4)}
    extend_muli(s, ((5,), (6,))) == {(1, 2, 3, 5), (1, 2, 3, 6),
                              (1, 2, 4, 5), (1, 2, 4, 6),
                              (1, 3, 4, 5), (1, 3, 4, 6),
                              (2, 3, 4, 5), (2, 3, 4, 6)}
    """
    return {(*t, *item) for t in s for item in e}

def extend_add(s: set[tuple[Any]], e: Iterable) -> set[tuple[Any]]:
    """
    set where each tuple of s is extended once, all  items of *e

    s = {(1, 2, 3), (1, 2, 4), (1, 3, 4), (2, 3, 4)}
    extend_add(s, (5, 6)) == {(1, 2, 3, 5, 6), (1, 2, 4, 5, 6),
                              (1, 3, 4, 5, 6), (2, 3, 4, 5, 6)}
    """
    return {(*t, *e) for t in s}

Diff by 2 pattern finding

It was incremental - find a pattern in 3C?, subtract it from 5C3, Find a pattern in 3C? that covers part of the remainder; repeat.
(Where ? <=3).

I also shortened function pf_set to pf so it would take less space when printing formatted expressions in f-strings

# %%

print(f"To simplify formatted printing of sets: {(pf:=pf_set) = }\n")

print("# Looking again at the diffs by 2 i.e. `5C3 - 3C3`")
n, r, x = 3, 3, 2

print(f"5C3 - 3C3 = {pf(bino(5, 3) - bino(3, 3))}")
print(f"\n  There's {pf(bino(3, 3-1)) = } in there with each tuple extended by 4 and by 5 ")
print(f"\n  {pf(extend_mul(bino(3, 3-1), (3+1, 3+2))) = }")
print(f"\n  Whats left: {pf(tmp1 := (bino(5, 3) - bino(3, 3) - extend_mul(bino(3, 3-1), (3+1, 3+2)))) = }")
print(f"\n    Now {pf(bino(3, 3-2)) = }")
print(f"\n    So {pf(tmp2 := (extend_add(bino(3, 3-2), (3+1, 3+2)))) = }")
print(f"\n  Finally: {pf(tmp1 - tmp2) = }")

Cell output:

To simplify formatted printing of sets: (pf:=pf_set) = <function pf_set at 0x7f1bd45365c0>

# Looking again at the diffs by 2 i.e. `5C3 - 3C3`
5C3 - 3C3 = {(1, 2, 4), (1, 2, 5), (1, 3, 4), (1, 3, 5), (1, 4, 5), (2, 3, 4), (2, 3, 5), (2, 4, 5), (3, 4, 5)}

  There's pf(bino(3, 3-1)) = '{(1, 2), (1, 3), (2, 3)}' in there with each tuple extended by 4 and by 5 

  pf(extend_mul(bino(3, 3-1), (3+1, 3+2))) = '{(1, 2, 4), (1, 2, 5), (1, 3, 4), (1, 3, 5), (2, 3, 4), (2, 3, 5)}'

  Whats left: pf(tmp1 := (bino(5, 3) - bino(3, 3) - extend_mul(bino(3, 3-1), (3+1, 3+2)))) = '{(1, 4, 5), (2, 4, 5), (3, 4, 5)}'

    Now pf(bino(3, 3-2)) = '{(1,), (2,), (3,)}'

    So pf(tmp2 := (extend_add(bino(3, 3-2), (3+1, 3+2)))) = '{(1, 4, 5), (2, 4, 5), (3, 4, 5)}'

  Finally: pf(tmp1 - tmp2) = '{}'

Behold Diff by 2

# %%
n, r, x = None, None, None

print(f"\n# lets set some variables and use those for the diffs by 2")
print(f"  {(n:=3), (r:=3) = }")
print(f"  {len(bino(n+2, r)) = }")
print(f"  {bino(n+2, r) == (bino(n, r) | extend_mul(bino(n, r-1), (n+1, n+2)) | extend_add(bino(n, r-2), (n+1, n+2))) = }")

Cell output:

# lets set some variables and use those for the diffs by 2
  (n:=3), (r:=3) = (3, 3)
  len(bino(n+2, r)) = 10
  bino(n+2, r) == (bino(n, r) | extend_mul(bino(n, r-1), (n+1, n+2)) | extend_add(bino(n, r-2), (n+1, n+2))) = True

Pascals rules for increasing diffs

I followed the same method for diffs of three and ended up with these three functions:

# %%

# By similar observation and much checking:
def pascal_rule_1(n, r) -> bool:
    """
    Checks bino(n+1, r) == (bino(n, r)
                            | extend_add(bino(n, r-1), (n+1,)))
    """
    return bino(n+1, r) == (bino(n, r)
                            | extend_add(bino(n, r-1), (n+1,)))


def pascal_rule_2(n, r) -> bool:
    """
    Checks bino(n+2, r) == (bino(n, r)
                            | extend_mul(bino(n, r-1), (n+1, n+2))
                            | extend_add(bino(n, r-2), (n+1, n+2)))
    """
    return bino(n+2, r) == (bino(n, r)
                            | extend_mul(bino(n, r-1), (n+1, n+2))
                            | extend_add(bino(n, r-2), (n+1, n+2)))

def pascal_rule_3(n, r) -> bool:
    """
    Checks bino(n+3, r) == (bino(n, r)
                            | extend_muli(bino(n, r-1), tuple(combinations((n+1, n+2, n+3), 1)))
                            | extend_muli(bino(n, r-2), tuple(combinations((n+1, n+2, n+3), 2)))
                            | extend_muli(bino(n, r-3), tuple(combinations((n+1, n+2, n+3), 3)))
                            )
    """
    extra_n = tuple(range(n+1, n+4))  # n..n+3 inclusive
    return bino(n+3, r) == (bino(n, r)
                            | extend_muli(bino(n, r-1), tuple(combinations(extra_n, 1)))
                            | extend_muli(bino(n, r-2), tuple(combinations(extra_n, 2)))
                            | extend_muli(bino(n, r-3), tuple(combinations(extra_n, 3)))
                            )

# %%
# Simple Checks
assert pascal_rule_1(7, 4)
assert pascal_rule_2(9, 4)
assert pascal_rule_3(11, 4)

Generalised Pascals rule

What can I say, I looked for patterns in the pascals rule functions for discrete diffs and tried to find patterns. I looked deeper into identities between the extend functions.

I finally found the following function that passed my tests, (many not shown).

# %%
# from pascal_rule_3 to pascal_rule_x

def pascal_rule_x(n, r, x) -> bool:
    """
    Checks bino(n+x, r) == union_reduce(<extend_muli selections of bino(n, 0 < r < r)> )

    ie if already used bino(n, r) and along comes x more for n then it shows how to
    calculate bino(n+x, r) without the need for storing bino(n, r)

    """
    extra_n = tuple(range(n+1, n+x+1))  # n..n+x inclusive
    n_r_terms = (extend_muli((bino(n, r-i) if r-i > 0 else {()}),   # extend this
                             tuple(combinations(extra_n, i)))       # by this
                for i in range(min(x, r) + 1))
    reduction = reduce(set.union, n_r_terms, set())
    return bino(n+x, r) == reduction

assert pascal_rule_x(11, 4, 3)  # n, r, x
assert pascal_rule_x(11, 5, 4)
assert pascal_rule_x(3, 2, 3)

I don't like the if r-i > 0 else {()} bit as it doesn't seem elegant. There is probably some identity to be found that would make it disappear but, you know.

Back to the original problem

If comb(d1, r) is processed and then we find an extra d2 items, then we want to process extra_comb(d1, r, d2) where extra_comb does not include or save comb of d1.

We just need to exclude the nCr term in reduction of function pascal_rule_x.

Eventually I arrive at

# %%

first = list(combinations('abcd', 3))
first
# %%
all = list(combinations('abcdef', 3))
# %%
extra = sorted(set(all) - set(first))
extra
# %%

def extra_combs(orig='abcd', r=3, delta='ef'):
    C = combinations  # Less typing
    extra_n = tuple(delta)  # n..n+x inclusive
    n = tuple(orig)
    n_r_terms = (extend_muli((C(n, r-i) if r-i > 0 else {()}),   # extend this
                             tuple(C(extra_n, i)))               # by this
                for i in range(1, min(len(extra_n), r) + 1))     # miss C(n, r)
    reduction = reduce(set.union, n_r_terms, set())
    # set(C(n+extra_n, r)) - reduction == set(C(n, r))
    return reduction

n, r, delta = 'abcd', 3, 'efg'
assert set(combinations(n+delta, r)) \
         == set(combinations(n, r)).union(extra_combs(n, r, delta))

Tests

# %%

# Test comb(n+x, r) == comb(n, r) | extra_combs(n, r, x)

n, r, delta_ = 'abcdefg', 3, 'hijkl'

for r in range(len(n)):
    for delta in (delta_[:i] for i in range(len(delta_))):
        combnx = set(combinations(n+delta, r))
        combn = set(combinations(n, r))
        extra = extra_combs(n, r, delta)
        assert combnx == (combn | extra), f"Whoops! For  {(n, r, delta) = }"
        # checks that extra does not generate any of combn
        assert not extra.intersection(combn), f"Whoops! For  {(n, r, delta) = }"

END.

Intsets by AI

2025-03-05T22:14:00.000+00:00

I will be using sets a lot in a future algorithm and wanted to find out about the speed of operation of different set implementations.

I will be dealing with different sets of strings from a total of around 200_000, (but the strings being added throughout the calculations). Python sets can do that directly, but I can create an integer index of all the strings and then do calculations on sets of integers.

Sets of integers can be translated to ints with bits set to indicate which ints were in the set. Numpy can do this but needs either all the strings up-front, or the maximum integer up front. Python integers are variable length, as the algorithm adds more strings, they can be looked up and either their index returned or an incrementing index assigned and returned if the string is new. if you want to set the seven thousandth and one bit of a Python int then the int will extend as necessary to allow this.

Code by AI

I decided to talk to Gemini AI to create all code and after a couple of hours of asking, refining and my testing I got my code. (Mine as I expended effort in its creation and debugging, even if it was by questioning Gemini).

Result!

I found that intsets can be hundreds of times faster than the same operation with the same values but with sets. (Frozensets in this particular case ase I cache a conversion function).

max_int: 262144, set_count: 5, repeat_count: 1
Average time for intset operations: 0.001064 seconds
Average time for set operations: 0.921396 seconds
Intset operations are approximately 866.31 times faster than set operations.

Code

"""
This module provides functions for generating random sets, performing set operations,
and timing the performance of set operations using both standard Python sets and
integer-based sets (intsets). It includes caching for intset conversions and calculates
the speedup of intset operations over standard set operations.
"""

import random
import time
import itertools
from functools import lru_cache
from typing import FrozenSet, Tuple, List

def generate_random_sets_with_sample(max_int: int, set_count: int) -> List[FrozenSet[int]]:
    """
    Generates a list of random frozensets of integers using random.sample.

    Args:
        max_int: The maximum integer value in the sets.
        set_count: The number of random sets to generate.

    Returns:
        A list of random frozensets of integers.
    """
    all_ints = frozenset(range(max_int + 1))
    random_sets = []

    for _ in range(set_count):
        population = range(max_int + 1)
        sample_size = random.randint(0, max_int + 1)
        current_set = frozenset(random.sample(population, sample_size))
        random_sets.append(current_set)

    random_sets.append(all_ints)  # Append all_ints to random_sets
    return random_sets

def intset_difference(intset1: int, intset2: int) -> int:
    """Calculates the difference of two intsets."""
    return intset1 & ~intset2

def intset_intersection(intset1: int, intset2: int) -> int:
    """Calculates the intersection of two intsets."""
    return intset1 & intset2

def intset_symmetric_difference(intset1: int, intset2: int) -> int:
    """Calculates the symmetric difference of two intsets."""
    return intset1 ^ intset2

def intset_union(intset1: int, intset2: int) -> int:
    """Calculates the union of two intsets."""
    return intset1 | intset2

def set_difference(set1: FrozenSet[int], set2: FrozenSet[int]) -> FrozenSet[int]:
    """Calculates the difference of two frozensets."""
    return set1.difference(set2)

def set_intersection(set1: FrozenSet[int], set2: FrozenSet[int]) -> FrozenSet[int]:
    """Calculates the intersection of two frozensets."""
    return set1.intersection(set2)

def set_symmetric_difference(set1: FrozenSet[int], set2: FrozenSet[int]) -> FrozenSet[int]:
    """Calculates the symmetric difference of two frozensets."""
    return set1.symmetric_difference(set2)

def set_union(set1: FrozenSet[int], set2: FrozenSet[int]) -> FrozenSet[int]:
    """Calculates the union of two frozensets."""
    return set1.union(set2)

@lru_cache(maxsize=None)
def set_to_intset(input_set: FrozenSet[int]) -> int:
    """
    Converts a frozenset of integers to an intset, with caching.

    Args:
        input_set: The frozenset of integers.

    Returns:
        The intset representation of the frozenset.
    """
    result = 0
    for num in input_set:
        result |= (1 << num)
    return result

def time_set_operations(random_sets: List[FrozenSet[int]], repeat_count: int = 1) -> Tuple[float, float]:
    """
    Times set and intset operations multiple times and returns average timings.

    Args:
        random_sets: A list of frozensets.
        repeat_count: The number of times to repeat the timing.

    Returns:
        A tuple containing the average intset operation time and the average set operation time.
    """
    total_intset_time = 0
    total_set_time = 0

    for _ in range(repeat_count):
        intset_time = 0
        set_time = 0

        for set1, set2 in itertools.permutations(random_sets, 2):
            # Time intset operations
            intset1 = set_to_intset(set1)
            intset2 = set_to_intset(set2)

            start_intset = time.time()
            intset_difference(intset1, intset2)
            intset_intersection(intset1, intset2)
            intset_symmetric_difference(intset1, intset2)
            intset_union(intset1, intset2)
            end_intset = time.time()
            intset_time += end_intset - start_intset

            # Time set operations
            start_set = time.time()
            set_difference(set1, set2)
            set_intersection(set1, set2)
            set_symmetric_difference(set1, set2)
            set_union(set1, set2)
            end_set = time.time()
            set_time += end_set - start_set

        total_intset_time += intset_time
        total_set_time += set_time

    avg_intset_time = total_intset_time / repeat_count
    avg_set_time = total_set_time / repeat_count

    return avg_intset_time, avg_set_time

# Example Usage:
max_int = 2**18
set_count = 5
repeat_count = 1

random_sets = generate_random_sets_with_sample(max_int, set_count)

avg_intset_time, avg_set_time = time_set_operations(random_sets, repeat_count)

print(f"max_int: {max_int}, set_count: {set_count}, repeat_count: {repeat_count}")
print(f"Average time for intset operations: {avg_intset_time:.6f} seconds")
print(f"Average time for set operations: {avg_set_time:.6f} seconds")

if avg_intset_time > 0:
    speedup = avg_set_time / avg_intset_time
    print(f"Intset operations are approximately {speedup:.2f} times faster than set operations.")
else:
    print("Set operations time was zero. Unable to calculate speedup.")
    

My Gemini chat

1. python function set_to_intset that given a set containing integers returns an
    int with just those binary bits set to one; and a reverse function called
    intset_to_set that given an int returns a set of ints where the ints are the
    positions of all the bits that are one in the intset argument.
2. methods between two python sets: difference, intersection,
    symmetric_difference and union can be achieved between two intsets using
    logical operators and logical functions. What are they?
3. add a functionwith arguments max_int, set_count, and repeats; generates
    set_count random sets of the ints 0 to max_int and the set of all ints from
    zero to max_int.
4. use random.sample with a randomised sample count
5. population should be the range
6. append all_ints to random_sets then test the stated set functions on all
    permutations of two of the random sets. accumulate the execution time of the
    set functions
7. remove list from variable permutations then do not use an intermediate
    variable permutations, just iterate on its expression.
8. Do similar timings for intset transformations of all random_sets
9. No, keep the method timings as before.
10. keep two time accumulators, one for timing the runs of the set functions; do
    not time the transformations between set and intset; use the other timer for
    timing the intset functions. return two two times
11. do not time the conversions. time the set operations too
12. time_set_operations should take a repeat count and repeat the for loop that
    many times then return the avrage timings for each of sets, intsets
13. do not print the sets. print max_int, set_count, and repeat_count
14. show how much faster intset is
15. speedup calculation is wrong
16. speedup of intset is intset time / set time
17. no it is with respect to set time so divide by set time
18. that is not what I typed!
19. I was wrong about the speedum calculation
20. use it
21. inline the set andintset calculations when timing instead of calling out to
    the set_* and intset_+ functions, (but keep their definitions)
22. use max_int = 3**18 and set_count = 10
23. use max_int = 2**18; set_count = 5; repeat_count = 1
24. use lru_cached versions of function set_to_intset
25. generate random_sets as frozensets directly
26. add function type_hints
27. only have type hinting on function arguments and return values
28. add docstrings to all functions and add a module docstring

End.

Me, Chatgpt, copilot, gemini, and google search classify quadrilaterals

2025-01-01T14:13:00.000+00:00

(Best viewed on a larger than phone screen)

I decided to use AI more in coding, so found some free tools and incorporated them into my coding experience.

My initial goal

Someone on Linkedin posted their classification of types of quadrilaterals with a view to how they would implement each of them using OO principles. I commented on it, but was more interested in creating a program to state what type of quadrilateral it is given.

Starting off

I won't go through all the many questions I put to the AI's but you ask something of the AI, do background searches on results, revise a question, think of implementation, revise the questions to the AI again, ....

Early on I thought the best way to input a quadrilateral would be by entering four ordered points that traced their way around the periphery of the shape so after chatgpt had given a different definition of a quadrilateral I asked:

Is that the same as defining a quadrilateral as "A closed figure constructed from tracing four ordered, unequal, noncollinear points closing by tracing from the last to the first point"

Because their were crucial occasions were lengths were computed and compared as floating point values I got the AI to align all input points onto a grid. After a while, I had to weigh how precise I would need to be to make specific changes to parts of the source in the AI, (and keeping track of those changes), versus just fixing my copy of some function and then just using select parts of AI generated code - with edits if it was from a different AI or something from a search.

Testing

I was using Vscode as my Python editor, I was not writing a Jupyter notebook, but I started using

# %%

comments to split my source into cells that I could easily run/rerun in its interactive and graphical interpreter.

I needed some test data so, with the aid of image searches for types of quadrilaterals, I defined my type of quadrilaterals, and stated what features they have as:

test_cases = [  # ([points], "classification", "description")
    ([(1, 1), (7, 3), (6, 6), (4, 9)], "Convex", "Non crossing, all angles <180"),
    ([(1, 1), (7, 1), (7, 4), (1, 4)], "Rectangle1", "Non-crossing, all ang == 90, opp sides =="),
    ([(1, 1), (1, 7), (4, 7), (4, 1)], "Rectangle2", "Non-crossing, all ang == 90, all sides =="),
    ([(0, 3), (3, 0), (9, 6), (6, 9)], "Rectangle3", "Non-crossing, all ang == 90, opp sides ==, rotated"),
    ([(1, 1), (7, 1), (7, 7), (1, 7)], "Square1", "Non-crossing, all ang == 90, opp sides =="),
    ([(1, 5), (4, 1), (8, 4), (5, 8)], "Square2", "Non-crossing, all ang == 90, all sides ==, rotated"),
    ([(1, 1), (6, 1), (9, 5), (4, 5)], "Rhombus", "Non-crossing, ang !=90, opp sides || and all sides =="),
    ([(1, 1), (6, 1), (9, 10), (4, 10)], "Parallelogram", "Non-crossing, opp sides || and opp sides =="),
    ([(1, 1), (7, 1), (5, 5), (2, 5)], "Trapezium1", "Non-crossing, at least one set of opposite || sides"),
    ([(1, 1), (5, 3), (5, 6), (1, 7)], "Trapezium2", "Non-crossing, at least one set of opposite || sides"),
    ([(1, 4), (5, 1), (7, 4), (5, 7)], "Kite1", "Non crossing, No interior ang >180, two non== sets of == lines that touch"),
    ([(4, 1), (1, 5), (4, 11), (7, 5)], "Kite2", "Non crossing, No interior ang >180, two non== sets of == lines that touch"),
    ([(1, 1), (5, 5), (7, 4), (4, 8)], "Reflex", "Non crossing, one interior ang >180"),
    ([(1, 4), (7, 1), (4, 4), (7, 7)], "Dart1", "Non crossing, one interior ang >180, two non== sets of == lines that touch"),
    ([(4, 1), (1, 7), (4, 4), (7, 7)], "Dart2", "Non crossing, one interior ang >180, two non== sets of == lines that touch"),
    ([(1, 1), (6, 6), (7, 2), (3, 9)], "Crossed", "Crossing"),
    ([(1, 1), (7, 7), (7, 3), (1, 9)], "Truss1", "Crossing, crossing lines ==, none-crossing lines ||"),
    ([(1, 1), (7, 7), (3, 7), (9, 1)], "Truss2", "Crossing, crossing lines ==, none-crossing lines ||"),
    ([(1, 1), (7, 9), (7, 1), (1, 9)], "Bowtie1", "Crossing, crossing lines ==, none-crossing lines || and =="),
    ([(1, 1), (9, 7), (1, 7), (9, 1)], "Bowtie2", "Crossing, crossing lines ==, none-crossing lines || and =="),
    ([(0, 3), (9, 6), (6, 9), (3, 0)], "Bowtie3", "Crossing, crossing lines ==, none-crossing lines || and ==, rotated"),
]

It took some time working out a few of the points necessary as the AI's "can't count" very well so most were done by my hand calculations after giving up on AI help for this. (One gets tired of the confident but wrong replies).

You may not see the truss type of quadrilateral elsewhere, Bowties, to me always have that extra symmetry, but I needed a name for the category so chose truss.

Somewhere, there is a cell defining a function that visualises all the quadrilaterals. I started with chatgpt's help in starting function visualise_test_cases then modified it by hand, as necessary. The picture at the start of the blog is the final matplotlib output. Note: I used several versions of Python in development and some did not have matplotlib so used # %% cell shenanigans, and if 0: blocks so that the matplotlib stuff only runs in the vscode interpreter).

With testing, and debugging came another lot of AI and google searches, and AI describing some of the algorithms used. A google search lead to a StackOverflow q&A's where the second most used answer for finding the angles between two lines used atan2, and with an explanation of its benefits. I decided to use it and asked chatgpt to do this:

in _angle() is it better to use atan2?

That was most of what was needed to update def _angle() to that shown. Because float angles are also being compared, I wrote _round_angle myself and applied it appropriately.

_classify becomes _categorise

I did not like the first, "default", method of classification from chatgpt, which was refined as function _classify which returned only one type of quadrilateral for the given points, and was a bit untidy in how it got there.. I asked several questions to get chatgpt to first calculate properties of the quadrilateral, and refined it to form function _categorise. Putting my tests through it I got all the properties for my named quadrilaterals and pared them and reordered them to form dict quad2features. If the quadrilaterals features are a superset of those shown in quad2features then the quadrilateral is of the type in the key.

OK output

I started to get output that looked correct:

Convex Non crossing, all angles <180 [(1, 1), (7, 3), (6, 6), (4, 9)] QuadrilateralClassifier(points).classify() =['Convex'] Rectangle1 Non-crossing, all ang == 90, opp sides == [(1, 1), (7, 1), (7, 4), (1, 4)] QuadrilateralClassifier(points).classify() =['Convex', 'Parallelogram', 'Rectangle'] Rectangle2 Non-crossing, all ang == 90, all sides == [(1, 1), (1, 7), (4, 7), (4, 1)] QuadrilateralClassifier(points).classify() =['Convex', 'Parallelogram', 'Rectangle'] Rectangle3 Non-crossing, all ang == 90, opp sides ==, rotated [(0, 3), (3, 0), (9, 6), (6, 9)] QuadrilateralClassifier(points).classify() =['Convex', 'Parallelogram', 'Rectangle'] Square1 Non-crossing, all ang == 90, opp sides == [(1, 1), (7, 1), (7, 7), (1, 7)] QuadrilateralClassifier(points).classify() =['Convex', 'Rhombus', 'Square'] Square2 Non-crossing, all ang == 90, all sides ==, rotated [(1, 5), (4, 1), (8, 4), (5, 8)] QuadrilateralClassifier(points).classify() =['Convex', 'Rhombus', 'Square'] Rhombus Non-crossing, ang !=90, opp sides || and all sides == [(1, 1), (6, 1), (9, 5), (4, 5)] QuadrilateralClassifier(points).classify() =['Convex', 'Rhombus'] Parallelogram Non-crossing, opp sides || and opp sides == [(1, 1), (6, 1), (9, 10), (4, 10)] QuadrilateralClassifier(points).classify() =['Convex', 'Parallelogram'] Trapezium1 Non-crossing, at least one set of opposite || sides [(1, 1), (7, 1), (5, 5), (2, 5)] QuadrilateralClassifier(points).classify() =['Convex', 'Trapezium'] Trapezium2 Non-crossing, at least one set of opposite || sides [(1, 1), (5, 3), (5, 6), (1, 7)] QuadrilateralClassifier(points).classify() =['Convex', 'Trapezium'] Kite1 Non crossing, No interior ang >180, two non== sets of == lines that touch [(1, 4), (5, 1), (7, 4), (5, 7)] QuadrilateralClassifier(points).classify() =['Convex', 'Kite'] Kite2 Non crossing, No interior ang >180, two non== sets of == lines that touch [(4, 1), (1, 5), (4, 11), (7, 5)] QuadrilateralClassifier(points).classify() =['Convex', 'Kite'] Reflex Non crossing, one interior ang >180 [(1, 1), (5, 5), (7, 4), (4, 8)] QuadrilateralClassifier(points).classify() =['Reflex'] Dart1 Non crossing, one interior ang >180, two non== sets of == lines that touch [(1, 4), (7, 1), (4, 4), (7, 7)] QuadrilateralClassifier(points).classify() =['Reflex', 'Dart'] Dart2 Non crossing, one interior ang >180, two non== sets of == lines that touch [(4, 1), (1, 7), (4, 4), (7, 7)] QuadrilateralClassifier(points).classify() =['Reflex', 'Dart'] Crossed Crossing [(1, 1), (6, 6), (7, 2), (3, 9)] QuadrilateralClassifier(points).classify() =['Reflex', 'Crossed'] Truss1 Crossing, crossing lines ==, none-crossing lines || [(1, 1), (7, 7), (7, 3), (1, 9)] QuadrilateralClassifier(points).classify() =['Reflex', 'Crossed', 'Truss'] Truss2 Crossing, crossing lines ==, none-crossing lines || [(1, 1), (7, 7), (3, 7), (9, 1)] QuadrilateralClassifier(points).classify() =['Reflex', 'Crossed', 'Truss'] Bowtie1 Crossing, crossing lines ==, none-crossing lines || and == [(1, 1), (7, 9), (7, 1), (1, 9)] QuadrilateralClassifier(points).classify() =['Reflex', 'Crossed', 'Truss', 'Bowtie'] Bowtie2 Crossing, crossing lines ==, none-crossing lines || and == [(1, 1), (9, 7), (1, 7), (9, 1)] QuadrilateralClassifier(points).classify() =['Reflex', 'Crossed', 'Truss', 'Bowtie'] Bowtie3 Crossing, crossing lines ==, none-crossing lines || and ==, rotated [(0, 3), (9, 6), (6, 9), (3, 0)] QuadrilateralClassifier(points).classify() =['Reflex', 'Crossed', 'Truss', 'Bowtie']

The code

Remember, my goal was to learn more about coding with added AI assistance - not to create production ready code.

#!/bin/env python3.11

# %%
"""
Quadrilateral taxonomy

Classification of quadrilaterals defined by a closed figure constructed by
tracing four ordered, distinct, and noncollinear points, connecting consecutive
oints with straight lines, and closing the shape by connecting the last point
to the first.

"""

import math

class QuadrilateralClassifier:

    # Features of quadrilaterals
    quad2features = {
        'Convex':        {'Convex',},
        'Kite':          {'Convex', 'No Parallel Sides', 'Two Touching Equal Sides'},
        'Trapezium':     {'Convex', 'One Pair of Parallel Sides'},
        'Parallelogram': {'Convex', 'Two Pairs of Parallel Sides', 'Two Pairs of Equal Sides'},
        'Rhombus':       {'Convex', 'Two Pairs of Parallel Sides', 'Four Equal Sides'},
        'Rectangle':     {'Convex', 'Two Pairs of Parallel Sides', 'Two Pairs of Equal Sides', '4 Right Angles'},
        'Square':        {'Convex', 'Two Pairs of Parallel Sides', 'Four Equal Sides', '4 Right Angles'},
        'Reflex':        {'Concave'},
        'Dart':          {'Concave', 'No Parallel Sides', 'Two Touching Equal Sides'},
        'Crossed':       {'Crossed'},
        'Truss':         {'Crossed', 'One Pair of Parallel Sides'},
        'Bowtie':        {'Crossed', 'One Pair of Parallel Sides', 'Two Pairs of Equal Sides'},
    }

    def __init__(self, points, grid_size=0.01):
        """
        Initializes the QuadrilateralClassifier with points and grid size.

        Parameters:
            points (list of tuples): List of four ordered points (x, y).
            grid_size (int): Size of the grid to which points will be aligned.
        """
        if len(points) != 4:
            raise ValueError("Exactly four points are required to define a quadrilateral.")

        self.original_points = points
        self.grid_size = grid_size
        self.points = self._align_to_grid(points)

        # Check if points remain four distinct and noncollinear after alignment
        if len(set(self.points)) != 4:
            raise ValueError("Aligned points are not four distinct points.")
        if not self._is_noncollinear():
            raise ValueError("Aligned points are collinear, not forming a quadrilateral.")

    def _align_to_grid(self, points):
        """
        Aligns points to the nearest grid size.

        Parameters:
            points (list of tuples): List of points to align.

        Returns:
            list of tuples: List of aligned points.
        """
        return [(round(x / self.grid_size) * self.grid_size, round(y / self.grid_size) * self.grid_size) for x, y in points]

    def _round_length(self, length: float | int) -> float | int:
        """
        Rounds length to grid/10 size.

        Parameters:
            length: distance to round.

        Returns:
            Rounded length.
        """
        grid = self.grid_size / 10
        return round(length / grid) * grid

    def _round_angle(self, angle: float | int, _unit=0.01) -> float | int:
        """
        Rounds angle to multiples of _unit size.

        Parameters:
            angle: Angle to round.

        Returns:
            Rounded angle.
        """
        return round(angle / _unit) * _unit

    def _is_noncollinear(self):
        """
        Checks if the points are noncollinear.

        Returns:
            bool: True if points are noncollinear, False otherwise.
        """
        def are_points_collinear(p1, p2, p3):
            """
            Checks if three points are collinear.
            """
            x1, y1 = p1
            x2, y2 = p2
            x3, y3 = p3
            return (x2 - x1) * (y3 - y1) == (y2 - y1) * (x3 - x1)

        # Check all combinations of three points
        return not (
            are_points_collinear(self.points[0], self.points[1], self.points[2]) or
            are_points_collinear(self.points[0], self.points[1], self.points[3]) or
            are_points_collinear(self.points[0], self.points[2], self.points[3]) or
            are_points_collinear(self.points[1], self.points[2], self.points[3])
        )

    def _distance(self, p1, p2):
        """
        Calculates the Euclidean distance between two points.

        Parameters:
            p1, p2 (tuple): Two points.

        Returns:
            float: Distance between the points.
        """
        return self._round_length(math.sqrt((p2[0] - p1[0])**2 + (p2[1] - p1[1])**2))

    def _angle(self, p1, p2, p3):
        """
        Calculates the angle (in degrees) at p2 formed by segments (p1, p2) and (p2, p3).

        Parameters:
            p1, p2, p3 (tuple): Three points.

        Returns:
            float: Angle in degrees.
        """
        dx1, dy1 = p1[0] - p2[0], p1[1] - p2[1]
        dx2, dy2 = p3[0] - p2[0], p3[1] - p2[1]
        angle1 = math.atan2(dy1, dx1)
        angle2 = math.atan2(dy2, dx2)
        angle = math.degrees((angle2 - angle1) % (2 * math.pi))
        if angle > 180:
            angle = 360 - angle
        return self._round_angle(angle)
    
    def _is_convex(self):
        """
        Determines if the quadrilateral is convex by analyzing the cross product of vectors formed by consecutive points. 
        A quadrilateral is convex if all cross products have the same sign, indicating consistent rotational direction.

        Returns:
            bool: True if the quadrilateral is convex, False otherwise.
        """
        def cross_product_sign(p1, p2, p3):
            """
            Computes the sign of the cross product of vectors (p1->p2) and (p2->p3).
            """
            return (p2[0] - p1[0]) * (p3[1] - p2[1]) - (p2[1] - p1[1]) * (p3[0] - p2[0])

        signs = []
        for i in range(4):
            p1 = self.points[i]
            p2 = self.points[(i + 1) % 4]
            p3 = self.points[(i + 2) % 4]
            signs.append(cross_product_sign(p1, p2, p3))

        return all(s > 0 for s in signs) or all(s < 0 for s in signs)

    def _is_crossed(self):
        """
        Checks if the quadrilateral is crossed (self-intersecting).

        Returns:
            bool: True if the quadrilateral is crossed, False otherwise.
        """
        def do_segments_intersect(p1, p2, p3, p4):
            """
            Checks if segments (p1, p2) and (p3, p4) intersect.

            This is determined using orientation tests to check if the endpoints of each segment lie on opposite sides of the other segment.
            If the segments intersect, the intersection point lies within the bounds of each segment.

            Parameters:
                p1, p2 (tuple): Endpoints of the first segment.
                p3, p4 (tuple): Endpoints of the second segment.

            Returns:
                bool: True if the segments intersect, False otherwise.
            """
            def orientation(a, b, c):
                val = (b[1] - a[1]) * (c[0] - b[0]) - (b[0] - a[0]) * (c[1] - b[1])
                if math.fabs(val) < 0.000001:
                    return 0
                if val == 0:
                    return 0  # Collinear
                return 1 if val > 0 else -1

            o1 = orientation(p1, p2, p3)
            o2 = orientation(p1, p2, p4)
            o3 = orientation(p3, p4, p1)
            o4 = orientation(p3, p4, p2)

            # General case
            if o1 != o2 and o3 != o4:
                return True

            return False

        def has_crossed_lines(points):
            """Check if a quadrilateral has crossed lines."""
            p1, p2, p3, p4 = points
            return (
                do_segments_intersect(p1, p2, p3, p4) or  # Check if line 1-2 intersects with line 3-4
                do_segments_intersect(p2, p3, p4, p1)     # Check if line 2-3 intersects with line 4-1
            )

        return has_crossed_lines(self.points)

    def _parallel_sides(self):
        """
        Counts the number of pairs of parallel sides.

        Returns:
            int: Number of pairs of parallel sides (0, 1, or 2).
        """
        def is_parallel(p1, p2, p3, p4):
            """
            Checks if lines (p1, p2) and (p3, p4) are parallel.
            """
            return (p2[1] - p1[1]) * (p4[0] - p3[0]) == (p2[0] - p1[0]) * (p4[1] - p3[1])

        parallel_count = 0
        if is_parallel(self.points[0], self.points[1], self.points[2], self.points[3]):
            parallel_count += 1
        if is_parallel(self.points[1], self.points[2], self.points[3], self.points[0]):
            parallel_count += 1
        return parallel_count

    def _classify(self):
        """
        Classifies the type of quadrilateral.

        Returns:
            str: Type of quadrilateral.
        """
        if not self._is_noncollinear():
            return "Not a quadrilateral (points are collinear)"

        # Precompute distances
        d = [
            self._distance(self.points[i], self.points[(i + 1) % 4])
            for i in range(4)
        ]
        diag1 = self._distance(self.points[0], self.points[2])
        diag2 = self._distance(self.points[1], self.points[3])

        # Check if crossed
        if self._is_crossed():
            # Check for bowtie or hourglass
            if d[0] == d[2] and d[1] == d[3]:
                return "Bowtie (Crossed Quadrilateral)"
            if d[0] != d[2] and d[1] != d[3]:
                return "Hourglass (Crossed Quadrilateral)"
            return "Crossed Quadrilateral"

        # Check convexity
        is_convex = self._is_convex()
        if not is_convex:
            # Check for dart
            if d[0] == d[3] and d[1] == d[2]:
                return "Dart (Concave Quadrilateral)"
            return "Concave Quadrilateral"

        # Check properties
        if d[0] == d[2] and d[1] == d[3]:  # Opposite sides equal
            if diag1 == diag2:  # Diagonals equal
                if d[0] == d[1]:
                    return "Square (Convex)"
                return "Rectangle (Convex)"
            if d[0] == d[1]:
                return "Rhombus (Convex)"
            return "Parallelogram (Convex)"
        elif (d[0] == d[2] or d[1] == d[3]) and not (d[0] == d[1] == d[2] == d[3]):
            return "Trapezium (Convex)"
        elif d[0] == d[1] == d[2] == d[3]:
            return "Rhombus (Convex, Degenerate as Square)"

        return "Irregular Quadrilateral (Convex)"


    def _categorise(self):
        """
        Categorises the quadrilaterals properties.

        Returns:
            set of str: Properties of the quadrilateral.
        """
        types = []

        # Check convexity
        is_convex = self._is_convex()
        if is_convex:
            types.append("Convex")
        else:
            types.append("Concave")

        # Check if crossed
        if self._is_crossed():
            types.append("Crossed")

        # Check for angles > 180°
        reflex_angles = [
            self._angle(self.points[i], self.points[(i + 1) % 4], self.points[(i + 2) % 4]) > 180
            for i in range(4)
        ]
        if any(reflex_angles):
            types.append("Reflex")

        # Check number of right angles
        right_angle_count = sum(
            abs(self._angle(self.points[i], self.points[(i + 1) % 4], self.points[(i + 2) % 4]) - 90) < 1e-6
            for i in range(4)
        )
        types.append(f"{right_angle_count} Right Angles")

        # Check number of equal sides
        d = [
            self._distance(self.points[i], self.points[(i + 1) % 4])
            for i in range(4)
        ]
        unique_lengths = set(d)
        if len(unique_lengths) == 1:
            types.append("Four Equal Sides")
        elif len(unique_lengths) == 2:
            if d[0] == d[2] and d[1] == d[3]:
                types.append("Two Pairs of Equal Sides")
            else:
                types.append("Two Touching Equal Sides")

        # Check number of parallel sides
        parallel_sides = self._parallel_sides()
        if parallel_sides == 2:
            types.append("Two Pairs of Parallel Sides")
        elif parallel_sides == 1:
            types.append("One Pair of Parallel Sides")
        elif parallel_sides == 0:
            types.append("No Parallel Sides")

        return set(types)


    def classify(self):
        """
        Aligns points, validates them, and classifies the quadrilateral.

        Returns:
            str: Classification of the quadrilateral.
        """
        features = self._categorise()
        quads = [name
                 for name, needed in self.quad2features.items()
                 if needed.issubset(features)]
        return quads


# %%
# Example usage of visualize_test_cases function
test_cases = [  # ([points], "classification", "description")
    ([(1, 1), (7, 3), (6, 6), (4, 9)], "Convex", "Non crossing, all angles <180"),
    ([(1, 1), (7, 1), (7, 4), (1, 4)], "Rectangle1", "Non-crossing, all ang == 90, opp sides =="),
    ([(1, 1), (1, 7), (4, 7), (4, 1)], "Rectangle2", "Non-crossing, all ang == 90, all sides =="),
    ([(0, 3), (3, 0), (9, 6), (6, 9)], "Rectangle3", "Non-crossing, all ang == 90, opp sides ==, rotated"),
    ([(1, 1), (7, 1), (7, 7), (1, 7)], "Square1", "Non-crossing, all ang == 90, opp sides =="),
    ([(1, 5), (4, 1), (8, 4), (5, 8)], "Square2", "Non-crossing, all ang == 90, all sides ==, rotated"),
    ([(1, 1), (6, 1), (9, 5), (4, 5)], "Rhombus", "Non-crossing, ang !=90, opp sides || and all sides =="),
    ([(1, 1), (6, 1), (9, 10), (4, 10)], "Parallelogram", "Non-crossing, opp sides || and opp sides =="),
    ([(1, 1), (7, 1), (5, 5), (2, 5)], "Trapezium1", "Non-crossing, at least one set of opposite || sides"),
    ([(1, 1), (5, 3), (5, 6), (1, 7)], "Trapezium2", "Non-crossing, at least one set of opposite || sides"),
    ([(1, 4), (5, 1), (7, 4), (5, 7)], "Kite1", "Non crossing, No interior ang >180, two non== sets of == lines that touch"),
    ([(4, 1), (1, 5), (4, 11), (7, 5)], "Kite2", "Non crossing, No interior ang >180, two non== sets of == lines that touch"),
    ([(1, 1), (5, 5), (7, 4), (4, 8)], "Reflex", "Non crossing, one interior ang >180"),
    ([(1, 4), (7, 1), (4, 4), (7, 7)], "Dart1", "Non crossing, one interior ang >180, two non== sets of == lines that touch"),
    ([(4, 1), (1, 7), (4, 4), (7, 7)], "Dart2", "Non crossing, one interior ang >180, two non== sets of == lines that touch"),
    ([(1, 1), (6, 6), (7, 2), (3, 9)], "Crossed", "Crossing"),
    ([(1, 1), (7, 7), (7, 3), (1, 9)], "Truss1", "Crossing, crossing lines ==, none-crossing lines ||"),
    ([(1, 1), (7, 7), (3, 7), (9, 1)], "Truss2", "Crossing, crossing lines ==, none-crossing lines ||"),
    ([(1, 1), (7, 9), (7, 1), (1, 9)], "Bowtie1", "Crossing, crossing lines ==, none-crossing lines || and =="),
    ([(1, 1), (9, 7), (1, 7), (9, 1)], "Bowtie2", "Crossing, crossing lines ==, none-crossing lines || and =="),
    ([(0, 3), (9, 6), (6, 9), (3, 0)], "Bowtie3", "Crossing, crossing lines ==, none-crossing lines || and ==, rotated"),
]

# %%
if 0:
    ...
    
    # %%
    from  matplotlib import pyplot as plt

    def visualize_test_cases(test_cases):
        """
        Visualizes the test cases and labels them with their classification.

        Parameters:
            test_cases (list of tuples): List of tuples, each containing points and classification.
        """
        fig, axs = plt.subplots(3, 7, figsize=(15, 6))
        axs = axs.flatten()

        # Determine the global axis limits
        all_points = [point for points, __classification, _description in test_cases for point in points]
        xs, ys = zip(*all_points)
        x_min, x_max = min(xs), max(xs)
        y_min, y_max = min(ys), max(ys)

        for ax, (points, classification, _description) in zip(axs, test_cases):
            points = points[::]
            points.append(points[0])  # To complete the quadrilateral loop
            xs, ys = zip(*points)
            ax.plot(xs, ys, marker='o')

            ax.set_title(classification)
            ax.set_xticks([])
            ax.set_yticks([])
            ax.set_xlim(x_min - 1, x_max + 1)
            ax.set_ylim(y_min - 1, y_max + 1)

        plt.tight_layout()
        plt.show()

    visualize_test_cases(test_cases)

# %%
head = \
'NAME           CONVEX  REFLEX  CROSSING    >180DEG     =90DEG  SIDES==     PARALLEL'.strip().split()
body = [line.strip().split() for line in '''
Convex          True    False   False           0       -         -           -
Rectangle       True    False   False           0       4       2,pairs     2,pairs
Square          True    False   False           0       4       4           2,pairs
Rhombus         True    False   False           0       0       4           2,pairs
Parallelogram   True    False   False           0       0       2,pairs     2,pairs
Trapezium       True    False   False           -       -       -           1,pair
Kite            True    False   False           0       0       2,Tpairs    0
Reflex          False   True    False           1       0       -           -
Dart            False   True    False           1       0       2,Tpairs    0
Crossed         False   False   True            -       -       -           -
Truss           False   False   True            -       -       -           1,pair
Bowtie          False   False   True            -       -       2,pairs     1,pair

# KEY
#
# CONVEX        not REFLEX and not CROSSED
# REFLEX        One angle greater than 180 degrees
# CROSSING      two lines in figure cross over
# >180DEG       Interior angles greater than 180 degrees
# =90DEG        Interior right-angles
# SIDES==       Number of equal length sides
# PARALLEL      Number of parallel sides
#
# -             Don't care
# Tpairs        Pair of equal length sides Touching at one shared point

'''.strip().splitlines() if line.strip() and not line.startswith('#')]
# %%

#for points, classification, description in test_cases[4:5]:
#for points, classification, description in test_cases[-1:]:
quad2features = {}
for points, classification, description in test_cases:
    quad2features[classification] = QuadrilateralClassifier(points).classify()
    print(f"\n{classification}\n  {description}\n  {points}")
    print(f"  {QuadrilateralClassifier(points).classify() =}")

IS AI Worth it?

I was surprised by just how much AI will do on early prompts, but days are then spent honing and checking, again and again.. Some things, like the has_crossed_lines function just could not be done correctly by chatgpt. It continuousely got the line segments and points of the line segments wrong. In the end I had to cut-n-paste the correct points andintersection checks into my query for it to then "generate" the right function. I had to know what I was doing to have first found the error.

One needs to know more, or find out more, than the AI to spot mistakes. The AI can help with this by asking it to explain parts of its code then searching for corroborating evidence, or for limitations in the algorithms it has chosen.

The AI's can't calculate worth a damn, and this can affect their orderings too. It helps if the AI interface can not only generate code, but also run it. Asking AI to generate points for various types of quadrilateral was error prone, leading me to do it without AI aid.

Free AI accounts have limits! I had to wait till the next day on several occasions after using up my free time.

In summary, I will carry on trying out AI, hopefully both it and I will get better!

There's the easy way...

2024-11-17T23:44:00.002+00:00

Best seen on a larger than landscape phone

Someone blogged about a particular problem:

From: https://theweeklychallenge.org/blog/perl-weekly-challenge-294/#TASK1

Given an unsorted array of integers, `ints`
Write a script to return the length of the longest consecutive elements sequence.
Return -1 if none found. *The algorithm must run in O(n) time.*

The solution they blogged used a sort which meant it could not be O(n) in time, but the problem looked good so I gave it some thought.

Sets! sets are O(1) in Python and are good for looking things up.

What if when looking at the inputted numbers, one at a time, you also looked for other ints in the input that would extend the int you have to form a longer range? Keep tab of the longest range so far and if you remove ints from the pool as they form ranges, when the pool is empty, you should know the longest range.

I added the printout of the longest range too.

My code

def consec_seq(ints) -> tuple[int, int, int]:
    "Extract longest_seq_length, its_min, its_max"
    pool = set(ints)
    longest, longest_mn, longest_mx = 0, 1, 0
    while pool:
        this = start = pool.pop()
        ln = 1
        # check down
        while (this:=(this - 1)) in pool:
            ln += 1
            pool.remove(this)
        mn = this + 1
        # check up
        this = start
        while (this:=(this + 1)) in pool:
            ln += 1
            pool.remove(this)
        mx = this - 1
        # check longest
        if ln > longest:
            longest, longest_mn, longest_mx = ln, mn, mx

    return longest,longest_mn,longest_mx

def _test():
    for ints in[(),
            (69,),
            (-20, 78, 79, 1, 100),
            (10, 4, 20, 1, 3, 2),
            (0, 6, 1, 8, 5, 2, 4, 3, 0, 7),
            (10, 30, 20),
            (2,4,3,1,0, 10,12,11,8,9),  # two runs of five
            (10,12,11,8,9, 2,4,3,1,0),  # two runs of five - reversed
            (2,4,3,1,0,-1, 10,12,11,8,9),  # runs of 6 and 5
            (2,4,3,1,0, 10,12,11,8,9,7),   # runs of 5 and 6
            ]:
        print(f"Input {ints = }")
        longest, longest_mn, longest_mx = consec_seq(ints)

        if longest <2:
            print("  -1")
        else:
            print(f"  The/A longest sequence has {longest} elements {longest_mn}..{longest_mx}")


# %%
if __name__ == '__main__':
    _test()

Sample output

Input ints = ()
  -1
Input ints = (69,)
  -1
Input ints = (-20, 78, 79, 1, 100)
  The/A longest sequence has 2 elements 78..79
Input ints = (10, 4, 20, 1, 3, 2)
  The/A longest sequence has 4 elements 1..4
Input ints = (0, 6, 1, 8, 5, 2, 4, 3, 0, 7)
  The/A longest sequence has 9 elements 0..8
Input ints = (10, 30, 20)
  -1
Input ints = (2, 4, 3, 1, 0, 10, 12, 11, 8, 9)
  The/A longest sequence has 5 elements 0..4
Input ints = (10, 12, 11, 8, 9, 2, 4, 3, 1, 0)
  The/A longest sequence has 5 elements 0..4
Input ints = (2, 4, 3, 1, 0, -1, 10, 12, 11, 8, 9)
  The/A longest sequence has 6 elements -1..4
Input ints = (2, 4, 3, 1, 0, 10, 12, 11, 8, 9, 7)
  The/A longest sequence has 6 elements 7..12

Another Algorithm

What if, you kept and extended ranges untill you amassed all ranges then chose the longest? I need to keep the hash lookup. dict key lookup should also be O(1). What to look up? Look up ints that would extend a range!

If you have an existing (integer) range, say 1..3 inclusive of end points then finding 0 would extend the range to 0..3 or finding one more than the range maximum, 4 would extend the original range to 1..4

So if you have ranges then they could be extended by finding rangemin - 1 or rangemax +1. I call then extends

If you do find that the next int from the input ints is also an extends value then you need to find the range that it extends, (by lookup), so you can modify that range. - use a dict to map extends to their range and checking if an int is in the extends dict keys should also take O(1) time.

I took that sketch of an algorithm and started to code. It took two evenings to finally get something that worked and I had to work out several details that were trying. The main problem was what about coalescing ranges? if you have ranges 1..2 and 4..5 what happens when you see a 3? the resultant is the single range 1..5. It took particular test cases and extensive debugging to work out that the extends2range mapping should map to potentially more than one range and that you need to combine ranges if two of them are present for any extend value being hit.

So for 1..2 the extends being looked for are 0 and 3. For 4..5 the extends being looked for are 3, again, and 6. The extends2ranges data structure for just this should look like:

{0: [[1, 2]],
 3: [[1, 2], [4, 5]],
 6: [[4, 5]]}

The Code #2

from collections import defaultdict


def combine_ranges(min1, max1, min2, max2):
    "Combine two overlapping ranges return the new range as [min, max], and a set of limits unused in the result"
    assert (min1 <= max1 and min2 <= max2          # Well formed
            and ( min1 <= max2 and min2 <= max1 )) # and ranges touch or overlap
    range_limits = set([min1, max1, min2, max2])
    new_mnmx = [min(range_limits), max(range_limits)]
    unused_limits = range_limits - set(new_mnmx)

    return new_mnmx, unused_limits

def consec_seq2(ints) -> tuple[int, int, int]:
    "Extract longest_seq_length, its_min, its_max"
    if not ints:
        return -1, 1, -1
    seen = set()  # numbers seen so far
    extends2ranges = defaultdict(list)  # map extends to its ranges
    for this in ints:
        if this in seen:
            continue
        else:
            seen.add(this)

        if this not in extends2ranges:
            # Start new range
            mnmx = [this, this]    # Range of one int
            # add in the extend points
            extends2ranges[this + 1].append(mnmx)
            extends2ranges[this - 1].append(mnmx)
        else:
            # Extend an existing range
            ranges = extends2ranges[this]  # The range(s) that could be extended by this
            if len(ranges) == 2:
                # this joins the two ranges
                extend_and_join_ranges(extends2ranges, this, ranges)
            else:
                # extend one range, copied
                extend_and_join_ranges(extends2ranges, this, [ranges[0], ranges[0].copy()])

    all_ranges = sum(extends2ranges.values(), start=[])
    longest_mn, longest_mx = max(all_ranges, key=lambda mnmx: mnmx[1] - mnmx[0])

    return (longest_mx - longest_mn + 1), longest_mn, longest_mx

def extend_and_join_ranges(extends2ranges, this, ranges):
    mnmx, mnmx2 = ranges
    mnmx_orig, mnmx2_orig = mnmx.copy(), mnmx2.copy() # keep copy of originals
    mn, mx = mnmx
    mn2, mx2 = mnmx2
    if this == mn - 1:
        mnmx[0] = mn = this  # Extend lower limit of the range
    if this == mn2 - 1:
        mnmx2[0] = mn2 = this  # Extend lower limit of the range
    if this == mx + 1:
        mnmx[1] = mx = this  # Extend upper limit of the range
    if this == mx2 + 1:
        mnmx2[1] = mx2 = this  # Extend lower limit of the range
    new_mnmx, _unused_limits = combine_ranges(mn, mx, mn2, mx2)

    remove_merged_from_extends(extends2ranges, this, mnmx, mnmx2)
    add_combined_range_to_extends(extends2ranges, new_mnmx)


def add_combined_range_to_extends(extends2ranges, new_mnmx):
    "Add in the combined of two ranges's extends"
    new_mn, new_mx = new_mnmx
    for extend in (new_mn - 1, new_mx + 1):
        r = extends2ranges[extend]  # ranges at new limit extension
        if new_mnmx not in r:
            r.append(new_mnmx)

def remove_merged_from_extends(extends2ranges, this, mnmx, mnmx2):
    "Remove original ranges that were merged from extends"
    for lohi in (mnmx, mnmx2):
        lo, hi = lohi
        for extend in (lo - 1, hi + 1):
            if extend in extends2ranges:
                r = extends2ranges[extend]
                for r_old in (mnmx, mnmx2):
                    if r_old in r:
                        r.remove(r_old)
                if not r:
                    del extends2ranges[extend]
    # remove joining extend, this
    del extends2ranges[this]


def _test():
    for ints in[
            (),
            (69,),
            (-20, 78, 79, 1, 100),
            (4, 1, 3, 2),
            (10, 4, 20, 1, 3, 2),
            (0, 6, 1, 8, 5, 2, 4, 3, 0, 7),
            (10, 30, 20),
            (2,4,3,1,0, 10,12,11,8,9),  # two runs of five
            (10,12,11,8,9, 2,4,3,1,0),  # two runs of five - reversed
            (2,4,3,1,0,-1, 10,12,11,8,9),  # runs of 6 and 5
            (2,4,3,1,0, 10,12,11,8,9,7),   # runs of 5 and 6
            ]:
        print(f"Input {ints = }")
        longest, longest_mn, longest_mx = consec_seq2(ints)

        if longest <2:
            print("  -1")
        else:
            print(f"  The/A longest sequence has {longest} elements {longest_mn}..{longest_mx}")


# %%
if __name__ == '__main__':
    _test()

Its Output

 Input ints = ()
  -1
Input ints = (69,)
  -1
Input ints = (-20, 78, 79, 1, 100)
  The/A longest sequence has 2 elements 78..79
Input ints = (4, 1, 3, 2)
  The/A longest sequence has 4 elements 1..4
Input ints = (10, 4, 20, 1, 3, 2)
  The/A longest sequence has 4 elements 1..4
Input ints = (0, 6, 1, 8, 5, 2, 4, 3, 0, 7)
  The/A longest sequence has 9 elements 0..8
Input ints = (10, 30, 20)
  -1
Input ints = (2, 4, 3, 1, 0, 10, 12, 11, 8, 9)
  The/A longest sequence has 5 elements 0..4
Input ints = (10, 12, 11, 8, 9, 2, 4, 3, 1, 0)
  The/A longest sequence has 5 elements 8..12
Input ints = (2, 4, 3, 1, 0, -1, 10, 12, 11, 8, 9)
  The/A longest sequence has 6 elements -1..4
Input ints = (2, 4, 3, 1, 0, 10, 12, 11, 8, 9, 7)
  The/A longest sequence has 6 elements 7..12

This second algorithm gives correct results but is harder to develop and explain. It's a testament to my stubbornness as I thought there was a solution there, and debugging was me flexing my skills to keep them honed.

END.

Recreating the CVM algorithm for estimating distinct elements gives problems

2024-05-27T20:06:00.002+01:00

Someone at work posted a link to this Quanta Magazine article. It describes a novel, and seemingly straight-forward way to estimate the number of distinct elements in a datastream.

Quanta describes the algorithm, and as an example gives "counting the number of distinct words in Hamlet".

Following Quanta

I looked at the description and decided to follow their text. They carefully described each round of the algorithm which I coded up and then looked for the generalizations and implemented a loop over alll items in the stream ....

It did not work! I got silly numbers. I could download Hamlet split it into words, (around 32,000), do len(set(words) to get the exact number of distinct words, (around 7,000), then run it through the algorithm and get a stupid result with tens of digits for the estimated number of distinct words.
I re-checked my implementation of the Quanta-described algorithm and couldn't see any mistake, but I had originally noticed a link to the original paper. I did not follow it at first as original papers can be heavily into maths notation and I prefer reading algorithms described in code/pseudocode.

I decided to take a look at the original.

The CVM Original Paper

I scanned the paper.

I read the paper.

I looked at Algorithm 1 as a probable candidate to decypher into Python, but the description was cryptic. Heres that description taken from the paper:

AI To the rescue!?

I had a brainwave💡lets chuck it at two AI's and see what they do. I had Gemini and I had Copilot to hand and asked them each to express Algorithm 1 as Python. Gemini did something, and Copilot finally did something but I first had to open the page in Microsoft Edge.
There followed hours of me reading and cross-comparing between the algorithm and the AI's. If I did not understand where something came from I would ask the generating AI; If I found an error I would first, (and second and...), try to get the AI to make a fix I suggested.

At this stage I was also trying to get a feel for how the AI's could help me, (now way past what I thought the algorithm should be, just to see what it would take to get those AI's to cross T's and dot I's on a good solution).
Not a good use of time! I now know that asking questions to update one of the 20 to 30 lines of the Python function might fix that line, but unfix another line you had fixed before. Code from the AI does not have line numbers making it difficult to state what needs changing, and where.They can suggest type hints and create the beginnings of docstrings, but, for example, it pulled out the wrong authors for the name of the algorithm.
In line 1 of the algorithm, the initialisation of thresh is clearly shown, I thought, but both AI's had difficulty getting the Python right. eventually I cut-n-pasted the text into each AI, where they confidentially said "OF course...", made a change, and then I had to re-check for any other changes.

My Code

I first created this function:

def F0_Estimator(stream: Collection[Any], epsilon: float, delta: float) -> float:
    """
    ...
    """
    p = 1
    X = set()
    m = len(stream)
    thresh = math.ceil(12 / (epsilon ** 2) * math.log(8 * m / delta))

    for item in stream:
        X.discard(item)
        if random.random() < p:
            X.add(item)
        if len(X) == thresh:
            X = {x_item for x_item in X
                    if random.random() < 0.5}
            p /= 2
    return len(X) / p

I tested it with Hamlet data and it made OK estimates.

Elated, I took a break.

Hacker News

The next evening I decided to do a search to see If anyone else was talking about the algorithm and found a thread on Hacker News that was right up my street. People were discussing those same problems found in the Quanta Article - and getting similar ginormous answers. They had one of the original Authors of the paper making comments! And others had created code from the actual paper and said it was also easier than the Quanta description.

The author mentioned that no less than Donald Knuth had taken an interest in their algorithm and had noted that the expression starting `X = ...` four lines from the end could, thoretically, make no change to X, and the solution was to encase the assignment in a while loop that only exited if len(X) < thresh.

Code update

I decided to add that change:

def F0_Estimator(stream: Collection[Any], epsilon: float, delta: float) -> float:
    """
    Estimates the number of distinct elements in the input stream.

    This function implements the CVM algorithm for the problem of 
    estimating the number of distinct elements in a stream of data.
    
    The stream object must support an initial call to __len__ 

    Parameters:
    stream (Collection[Any]): The input stream as a collection of hashable 
        items.
    epsilon (float): The desired relative error in the estimate. It must be in 
        the range (0, 1).
    delta (float): The desired probability of the estimate being within the 
        relative error. It must be in the range (0, 1).

    Returns:
    float: An estimate of the number of distinct elements in the input stream.
    """
    p = 1
    X = set()
    m = len(stream)
    thresh = math.ceil(12 / (epsilon ** 2) * math.log(8 * m / delta))

    for item in stream:
        X.discard(item)
        if random.random() < p:
            X.add(item)
        if len(X) == thresh:
            while len(X) == thresh:  # Force a change
                X = {x_item for x_item in X
                     if random.random() < 0.5}  # Random, so could do nothing
            p /= 2
    return len(X) / p

thresh

In the code above, the variable thresh, (threshhold), named from Algorithm 1, is used in the Quanta article to describe the maximum storage available to keep items from the stream that have been seen before. You must know the length of the stream - m, epsilon, and delta to calculate thresh.

If you were to have just the stream and thresh as the arguments you could return both the estimate of the number of distinct items in the stream as well as counting the number of total elements in the stream.
Epsilon could be calculated from the numbers we now know.

def F0_Estimator2(stream: Iterable[Any],
                 thresh: int, 
                 ) -> tuple[float, int]:
    """
    Estimates the number of distinct elements in the input stream.

    This function implements the CVM algorithm for the problem of 
    estimating the number of distinct elements in a stream of data.
    
    The stream object does NOT have to support a call to __len__ 

    Parameters:
    stream (Iterable[Any]): The input stream as an iterable of hashable 
        items.
    thresh (int): The max threshhold of stream items used in the estimation.py

    Returns:
    tuple[float, int]: An estimate of the number of distinct elements in the 
        input stream, and the count of the number of items in stream.
    """
    p = 1
    X = set()
    m = 0  # Count of items in stream

    for item in stream:
        m += 1
        X.discard(item)
        if random.random() < p:
            X.add(item)
        if len(X) == thresh:
            while len(X) == thresh:  # Force a change
                X = {x_item for x_item in X
                     if random.random() < 0.5}  # Random, so could do nothing
            p /= 2
            
    return len(X) / p, m

def F0_epsilon(
               thresh: int,
               m: int,
               delta: float=0.05,  #  0.05 is 95%
              ) -> float:
    """
    Calculate the relative error in the estimate from F0_Estimator2(...)

    Parameters:
    thresh (int): The thresh value used in the call TO F0_Estimator2.
    m (int): The count of items in the stream FROM F0_Estimator2.
    delta (float): The desired probability of the estimate being within the 
        relative error. It must be in the range (0, 1) and is usually 0.05
        to 0.01, (95% to 99% certainty).

    Returns:
    float: The calculated relative error in the estimate

    """
    return math.sqrt(12 / thresh * math.log(8 * m / delta))

Testing

def stream_gen(k: int=30_000, r: int=7_000) -> list[int]:
    "Create a randomised list of k ints of up to r different values."
    return random.choices(range(r), k=k)

def stream_stats(s: list[Any]) -> tuple[int, int]:
    length, distinct = len(s), len(set(s))
    return length, distinct

# %% 
print("CVM ALGORITHM ESTIMATION OF NUMBER OF UNIQUE VALUES IN A STREAM")

stream_size = 2**18
reps = 5
target_uniques = 1
while target_uniques < stream_size:
    the_stream = stream_gen(stream_size+1, target_uniques)
    target_uniques *= 4
    size, unique = stream_stats(the_stream)

    print(f"\n  Actual:\n    {size = :_}, {unique = :_}\n  Estimations:")

    delta = 0.05
    threshhold = 2
    print(f"    All runs using {delta = :.2f} and with estimate averaged from {reps} runs:")
    while threshhold < size:
        estimate, esize = F0_Estimator2(the_stream.copy(), threshhold)
        estimate = sum([estimate] +
                    [F0_Estimator2(the_stream.copy(), threshhold)[0]
                        for _ in range(reps - 1)]) / reps
        estimate = int(estimate + 0.5)
        epsilon = F0_epsilon(threshhold, esize, delta)
        print(f"      With {threshhold = :7_} -> "
            f"{estimate = :_}, +/-{epsilon*100:.0f}%"
            + (f" {esize = :_}" if esize != size else ""))
        threshhold *= 8

The algorithm generates an estimate based on random sampling, so I run it multiple times for the same input and report the mean estimate from those runs.

Sample output

CVM ALGORITHM ESTIMATION OF NUMBER OF UNIQUE VALUES IN A STREAM

  Actual:
    size = 262_145, unique = 1
  Estimations:
    All runs using delta = 0.05 and with estimate averaged from 5 runs:
      With threshhold =       2 -> estimate = 1, +/-1026%
      With threshhold =      16 -> estimate = 1, +/-363%
      With threshhold =     128 -> estimate = 1, +/-128%
      With threshhold =   1_024 -> estimate = 1, +/-45%
      With threshhold =   8_192 -> estimate = 1, +/-16%
      With threshhold =  65_536 -> estimate = 1, +/-6%

  Actual:
    ...
 
  Actual:
    size = 262_145, unique = 1_024
  Estimations:
    All runs using delta = 0.05 and with estimate averaged from 5 runs:
      With threshhold =       2 -> estimate = 16_384, +/-1026%
      With threshhold =      16 -> estimate = 768, +/-363%
      With threshhold =     128 -> estimate = 1_101, +/-128%
      With threshhold =   1_024 -> estimate = 1_018, +/-45%
      With threshhold =   8_192 -> estimate = 1_024, +/-16%
      With threshhold =  65_536 -> estimate = 1_024, +/-6%

  Actual:
    size = 262_145, unique = 4_096
  Estimations:
    All runs using delta = 0.05 and with estimate averaged from 5 runs:
      With threshhold =       2 -> estimate = 13_107, +/-1026%
      With threshhold =      16 -> estimate = 3_686, +/-363%
      With threshhold =     128 -> estimate = 3_814, +/-128%
      With threshhold =   1_024 -> estimate = 4_083, +/-45%
      With threshhold =   8_192 -> estimate = 4_096, +/-16%
      With threshhold =  65_536 -> estimate = 4_096, +/-6%

  Actual:
    size = 262_145, unique = 16_384
  Estimations:
    All runs using delta = 0.05 and with estimate averaged from 5 runs:
      With threshhold =       2 -> estimate = 0, +/-1026%
      With threshhold =      16 -> estimate = 15_155, +/-363%
      With threshhold =     128 -> estimate = 16_179, +/-128%
      With threshhold =   1_024 -> estimate = 16_986, +/-45%
      With threshhold =   8_192 -> estimate = 16_211, +/-16%
      With threshhold =  65_536 -> estimate = 16_384, +/-6%

  Actual:
    size = 262_145, unique = 64_347
  Estimations:
    All runs using delta = 0.05 and with estimate averaged from 5 runs:
      With threshhold =       2 -> estimate = 26_214, +/-1026%
      With threshhold =      16 -> estimate = 73_728, +/-363%
      With threshhold =     128 -> estimate = 61_030, +/-128%
      With threshhold =   1_024 -> estimate = 64_422, +/-45%
      With threshhold =   8_192 -> estimate = 64_760, +/-16%
      With threshhold =  65_536 -> estimate = 64_347, +/-6%

Looks good!

Wikipedia

Another day, and I decide to start writing this blog post. I searched again and found the Wikipedia article on what it called the Count-distinct problem.

Looking through it, It had this wrong description of the CVM algorithm:

The, (or a?), problem with the wikipedia entry is that it shows

 $p \leftarrow \frac{p}{2}$

...within the while loop. You need an enclosing if |B| >= s for the while loop and the assignment to p outside the while loop, but inside this new if statement.

STOp PRESS! 02 May 2024: The Wikipedia entry was actually correct.

It's tough!

Both Quanta Magazine, and whoever added the algorithm to Wikipedia got the algorithm wrong.

I've written around two hundred tasks on site Rosettacode.org for over a decade. Others had to read my description and create code in their chosen language to implement those tasks. I have learnt from the feedback I got on talk pages to hone that craft, but details matter. Examples matter. Constructive feedback matters.

END.

Predicting results from small samples.

2024-05-10T20:09:00.000+01:00

I've run simulations, tens of thousands of them at a time, over and over as we developed chips. In one project I noticed that I could predict the final result after only a small number of results were in which allowed me to halt the rest of the simulations, or make advance preparations for the final result.

I looked it up at the time and, indeed, there is an equation where if you want the pass rate of a "large" population to within a given accuracy, it will give you the minimum sample size to use.

To some, that's all gobledygook so I'll try to explain with some code.

Explanation

Lets say you have a large randomised population of pass/fails or ones and zeroes:

from random import sample

population_size = 100_000    # must be large, > 65K?
sample_size = 123            # Actual
p_hat = 0.5                  # Population pass rate, 0.5 == 50%

_ones = int(population_size*p_hat)
_zeroes = population_size - _ones
population = [1] * _ones + [0] * _zeroes

And that we take a sample from it and compute the pass rate of that single, smaller sample

def random_sample() -> list[int]:
    return sample(population, k=sample_size)

pass_rate = (sum(random_sample())  # how many ones
             / sample_size)        # convert to a pass rate

print(pass_rate)  # e.g. 0.59027552674230146

Every time we run that pass_rate expression we get a different value. It is random. We need to run that pass_rate calculation many times to get an idea of what the likely pass_rate would be when using the given sample size:

runs_per_sample = 10_000     # each sample size is run this many times ~2500

pass_rates = [sum(random_sample()) / sample_size
              for _ in range(runs_per_sample)]

I have learnt that the question to ask is not how many times the sample pass rate was the population pass rate, but to define an acceptable margin of error, (say 5%) and say we want to be confident that pass rates be within that margin a certain percentage of the time

epsilon = 0.05               # target margin of error in sample pass rate. 0.05 == +/-5%

p_hat_min, p_hat_max = p_hat * (1 - epsilon), p_hat * (1 + epsilon)
in_range_count = sum(p_hat_min <= pass_rate < p_hat_max
                     for pass_rate in pass_rates)
sample_confidence_level = in_range_count / runs_per_sample
print(f"{sample_confidence_level = }")  # = 0.4054

So for a sample size of 123 we could expect the pass rate of the sample to be within 5% of the actual pass rate of the population, 0.5 or 50%, onlr o.4 or 40% of the time!

We need more!

What is actually done is we state what we think the population pass rate is, p_hat, (choose closer to 50% if you are unsure); the margin of error around p_hat we want, epsilon, usually +/-5% or +/-3%; and the confidence_level in hitting within the pass_rates margin of error.

There are calculators that will then give you n, the size of the sample needed to satisfy those condition.

Doing it myself

I calculated for one specific sample size, above. Obviousely, if I calculated pass_rates over a range of sample_sizes, and increqasing runs_per_sample, I could search out the sample size needed.

That is done in my next program. I have to switch to using the numpy library for its speed and sample_size becomes a range.

When the pass rate confidence levels are calculated I end up with a list of confidence levels for increasing sample sizes that are usually not increasing due to the randomness, e.g.

range_hits = [...0.94, 0.95, 0.93, 0.954, ... 0.95, 0.96, 0.95, 0.96, 0.96, 0.97, ...] # confidence levels

The range of sample_size corresponding to the first occurrence>= the requested confidence level, and the last occorrence of a confidence level < the requested confidence level in then slightly widened and the runs_per_sample increased on another iteration to get a better result.

Here's a sample of the output I get when searching

Sample output

$ time python3 sample_search.py 
2098 <= n <= 2610
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(50, 5000, 512), 250)

2013 <= n <= 2525
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(2013, 2694, 256), 500)

1501 <= n <= 2781
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1501, 3037, 256), 500)

1757 <= n <= 2013
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(221, 4061, 256), 500)

1714 <= n <= 1970
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1714, 2055, 128), 1000)

1458 <= n <= 2098
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1458, 2226, 128), 1000)

1586 <= n <= 1714
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(818, 2738, 128), 1000)

1564 <= n <= 1692
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1564, 1735, 64), 2000)

1500 <= n <= 1564
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1436, 1820, 64), 2000)

1553 <= n <= 1585
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1489, 1575, 32), 4000)

1547 <= n <= 1579
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1547, 1590, 16), 8000)

1547 <= n <= 1579
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1515, 1611, 16), 8000)

1541 <= n <= 1581
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1541, 1584, 8), 16000)

1501 <= n <= 1533
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1501, 1621, 8), 16000)

1501 <= n <= 1533
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1389, 1538, 8), 16000)

1503 <= n <= 1575
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1495, 1677, 8), 16000)

1503 <= n <= 1535
  For p_hat, epsilon, confidence_level =(0.5, 0.05, 0.95)
  Using population_size, sample_size, runs_per_sample =(100000, range(1491, 1587, 4), 32000)

 Use a sample n = 1535 to predict population pass rate of 50.0% +/-5% with a confidence level of 95%.

real    3m49.023s
user    3m49.022s
sys     0m1.361s

My Code

# -*- coding: utf-8 -*-
"""
Created on Wed May  8 14:04:17 2024

@author: paddy
"""
# %%
from random import sample
import numpy as np


def sample_search(population_size, sample_size, p_hat, epsilon, confidence_level, runs_per_sample) -> range:
    """
    Arguments with example values:
        population_size = 100_000           # must be large, > 65K?
        sample_size = range(1400, 1750, 16) # (+min, +max, +step)
        p_hat = 0.5                         # Population pass rate, 0.5 == 50%
        epsilon = 0.05                      # target margin of error in sample pass rate. 0.05 == +/-5%
        confidence_level = 0.95             # sample to be within p_hat +/- epsilon, 0.95 == 95% of the time.
        runs_per_sample = 10_000            # each sample size is run this many times ~2500
    
    Return:
        min,max range for the sample size, n, satisfying inputs.
    
    """
    
    def create_1_0_array(population_size=100_000, p_hat=0.5) -> np.ndarray:
        "Create numpy array of ones and zeroes with p_hat% as ones"
        ones = int(population_size*p_hat + 0.5)
        array10 = np.zeros(population_size, dtype=np.uint8)
        array10[:ones] = 1
        
        return array10


    def rates_of_samples(population: np.ndarray, sample_size_range: range, runs_per_sample: int
                        ) -> list[list[float]]:
        "Pass rates for range of sample sizes repeated runs_per_sample times."
        # many_samples_many_rates = [( np.random.shuffle(population),         # shuffle every *run*
        #                             [population[:s_count].sum() / s_count  # The pass rate for samples
        #                             for s_count in sample_size_range]   
        #                         )[1]                                     # drop the shuffle
        #                         for _ in range(runs_per_sample)]         # Every run

        many_samples_many_rates = [[ np.random.shuffle(population),         # shuffle every *sample*
                                     [population[:s_count].sum() / s_count  # The pass rate for samples
                                     for s_count in sample_size_range]   
                                   ][1]                                     # drop the shuffle
                                   for _ in range(runs_per_sample)]         # Every run

        return list(zip(*many_samples_many_rates))      # Transpose to by_sample_size_then_runs


    population = create_1_0_array(population_size, p_hat)
    by_sample_size_then_runs = rates_of_samples(population, sample_size, runs_per_sample)

    # Pass rates within target

    target_pass_range = tmin, tmax = p_hat * (1 - epsilon), p_hat * (1 + epsilon)  # Looking for rates within the range

    range_hits = [sum(tmin <= sample_pass_rate < tmax for sample_pass_rate in single_sample_size)
                for single_sample_size in by_sample_size_then_runs]
    runs_for_confidence_level = confidence_level * runs_per_sample

    for n_min, conf in zip(sample_size, range_hits):
        if conf >= runs_for_confidence_level:
            break
    else:
        n_min = sample_size.start

    for n_max, conf in list(zip(sample_size, range_hits))[::-1]:
        if conf <= runs_for_confidence_level:
            n_max += sample_size.step  # back a step
            break
    else:
        n_max = sample_size.stop
    
    if (n_min + sample_size.step) >= n_max and sample_size.step > 1:
        # Widen
        n_max = n_max + sample_size.step + 1

    return range(n_min, n_max, sample_size.step)


def min_max_mid_step(from_range: range) -> tuple[int, int, float, int]:
    "Extract from **increasing** from_range the min, max, middle, step"
    mn, st = from_range.start, from_range.step
    
    # Handle range where start == stop
    mx = from_range.stop
    for mx in from_range:
        pass
    
    md = (mn + mx) / 2
    
    return mn, mx, md, st
    
def next_sample_size(new_samples, last_samples,
                     runs_per_sample,
                     widener=1.33  # Widen range by
                    ):
    n_min, n_max, n_mid, n_step = min_max_mid_step(new_samples)
    l_min, l_max, l_mid, l_step = min_max_mid_step(last_samples)
    
    # Next range of samples computes in names with prefix s_
    
    increase_runs = True
    if n_max == l_max:
        # Major expand of high end
        s_max = l_max + (l_max - l_min)
        increase_runs = False
    else:
        s_max = (n_mid +( n_max - n_mid)* widener)

    if n_min == l_min:
        # Major expand of low end
        s_min = max(1, l_min + (l_min - l_max))
        increase_runs = False
    else:
        s_min = (n_mid +( n_min - n_mid)* widener)
    
    s_min, s_max = (max(1, int(s_min)), int(s_max + 0.5))
    s_step = n_step
    if s_min == s_max:
        if s_min > 2:
            s_min -= 1
        s_max += 1
        
    if increase_runs or n_max == n_min:
        runs_per_sample *= 2
        if n_max == n_min:
            s_step = 1
        else:
            s_step = max(1, (s_step + 1) // 2)  # Go finer
    
    next_sample_range = range(max(1, int(s_min)), int(s_max + 0.5), s_step)
    
    return next_sample_range, runs_per_sample

# %%

if __name__ == '__main__':

    population_size = 100_000           # must be large, > 65K?
    sample_size = range(50, 5_000, 512) # Increasing!
    p_hat = 0.50                        # Population pass rate, 0.5 == 50%
    epsilon = 0.05                      # target margin of error in sample pass rate. 0.05 == +/-5%
    confidence_level = 0.95             # sample to be within p_hat +/- epsilon, 0.95 == 95% of the time.
    runs_per_sample = 250               # each sample size is run this many time at start, ~250
    max_runs_per_sample = 35_000

    while runs_per_sample < max_runs_per_sample:
        new_range = sample_search(population_size, sample_size, p_hat, epsilon, confidence_level, runs_per_sample)
        n_min, n_max, n_mid, n_step = min_max_mid_step(new_range)
        print(f"{n_min} <= n <= {n_max}")
        print(f"  For {p_hat, epsilon, confidence_level =}\n"
            f"  Using {population_size, sample_size, runs_per_sample =}\n")
            
        sample_size, runs_per_sample = next_sample_size(new_range, sample_size, runs_per_sample)

print(f" Use a sample n = {n_max} to predict population pass rate of {p_hat*100.:.1f}% +/-{epsilon*100.:.0f}% "
      f"with a confidence level of {confidence_level*100.:.0f}%.")

END.

Searching OEIS tables

2024-04-21T15:30:00.002+01:00

A few months ago I submitted a series to OEIS* that was accepted; yes, but OEIS does not seem to leave my series searchable!

*OEIS is the Online Encyclopedia of Integer Series. I guess table is not in the name, but...

(best viewed on larger than a portrait phone)

Let me explain.

The documentation for OEIS, explains that if you have a 2D triangle or table of values rather than a one dimensional strict series, then one should antidiagonalise the data and submit the series produced.

They give as an example A003987 . This gives this table:

Table begins
   0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, ...
   1,  0,  3,  2,  5,  4,  7,  6,  9,  8, 11, 10, ...
   2,  3,  0,  1,  6,  7,  4,  5, 10, 11,  8, ...
   3,  2,  1,  0,  7,  6,  5,  4, 11, 10, ...
   4,  5,  6,  7,  0,  1,  2,  3, 12, ...
   5,  4,  7,  6,  1,  0,  3,  2, ...
   6,  7,  4,  5,  2,  3,  0, ...
   7,  6,  5,  4,  3,  2, ...
   8,  9, 10, 11, 12, ...
   9,  8, 11, 10, ...
  10, 11,  8, ...
  11, 10, ...
  12, ...
  ...

The above gets antidiagonalised to the series beginning:

  0, 1, 1, 2, 0, 2, 3, 3, 3, 3, 4, 2, 0, 2, 4, 5, 5, 1, 1, 5, 5, 6, 4,
  6, 0, 6, 4, 6, 7, 7, 7, 7, 7, 7, 7, 7, 8, 6, 4, 6, 0, 6, 4, 6, 8, 9,
  9, 5, 5, 1, 1, 5, 5, 9, 9, 10, 8, 10, 4, 2, 0, 2, 4, 10, 8, 10, 11,
  11, 11, 11, 3, 3, 3, 3, 11, 11, 11, 11, 12, 10, 8, 10, 12, 2, 0, 2,
  12, 10, 8, 10, 12, 13, 13, 9, 9, 13, 13, 1, 1, 13, 13, 9, 9, 13, 13

Searching...

If we search for a sub-sequence of the antidiagonalised table, we can find the correct entry.

If, however, we search for a row of the values from the table, A003987 is not found!
The values chosen: to search for: 5,4,7,6,1,0,3 appear near the end of the table which shows that that row of numbers should be followed by a 2.
The table shows 13*13 / 2 ~ 85 values. OEIS has a list of 104 values, so it has the data to search through.

No intuitive search of OEIS tables

It seems to me that the most intuitive way to search a table of values is by row, left to right. There are other ways to search a table, (assuming an origin at top left and the table extends to the right and down):

By row, L2R. →, R2L ←
By Column Top2Bottom, ↓, B2T ↑
By 45 degree diagonals, ↘, ↖, ↙, ↗

OEIS doesn't seem to do these searches on tabular data.

Regenerating a 2D table from antidiagonalised data.

I did play around and created some code to recreate a table as a list of row-lists, in Python, given an OEIS B-file. The options handling is a work in progress, but the main part was being able to generate the table.

# -*- coding: utf-8 -*-
# %%
"""
adia_to_table.py file

Generate 2D table from OEIS type anti-diagonalised sequence held in file


Created on Thu Apr  4 18:16:07 2024

@author: paddy3118
"""

from itertools import zip_longest
import math
# from pprint import pp
from antidiagonals import antidiag_triangle_indices


def read_bfile(bfname: str) -> tuple[int | None,  # first index
                                     list[int]]:   # anti-diag values
    """
Read B-file bfname

## B-File format:
* Ref: https://oeis.org/SubmitB.html
* Blank lines ignored
* Lines beginning  with # ignored
* Lines of two, space-separated integers <index> <indexed-value>

It is assumed the index increments by one on subsequent lines.
    """
    first_index, values = None, []
    with open(bfname) as bfile:
        for line in bfile:
            ln = line.strip()
            if not ln or ln.startswith('#'):
                continue
            index, value = [int(field) for field in ln.split()[:2]]
            if first_index is None:
                first_index = index
            values.append(value)

    return first_index, values


def antidiag_to_table(sequence: list[int]) -> list[list[int]]:
    """
    Convert anti-diagonalised sequence back to infinite 2D table.

    Parameters
    ----------
    sequence : list[int]
        Anti-diagonalised values from table.

    Returns
    -------
    list[list[int]]
        Table of rows of ints.

    Table rows will fill in from successive sequence values like this:

     1  2  4  7 11 16 ...
     3  5  8 12 17
     6  9 13 18
    10 14 19
    15 20
    21
    .
    .
    .
    """

    # 1, 3, 6, 10, 15, 21, ... rows*(rows+1) / 2

    # min columns in triangular table generation. ~= min rows
    size = len(sequence)                                  # = rows*(rows+1)/2
    rows = math.ceil((-1 + math.sqrt(4 * 2 * size)) / 2)  # solve for rows
    # cols = rows  # last row may be deleted if last anti-diag is part filled.
    # print(f"{(size, cols) = }")

    # Empty (triangular) table of None's
    table = [[None] * (rows - i)
             for i in range(rows)]

    indices = antidiag_triangle_indices()
    col = 0  # for if sequence is empty
    for val, (row, col) in zip(sequence, indices):
        table[row][col] = val

    # Remove unfilled part of last anti-diag of table
    while col > 0:
        row, col = next(indices)
        table[row].pop(-1)
    # remove last row if present and empty
    if table and not table[-1]:
        table.pop(-1)

    return table


def pp_table(table: list[list[int]]) -> None:
    "Pretty-print table of differring row lengths"
    if not table:
        return
    col_width = max(max(len(str(val)) for val in row) for row in table)
    for row in table:
        print(''.join(f"{val:{col_width}}" for val in row))


def transpose(table: list[list[int]]) -> list[list[int]]:
    "Table of rows to x<>y transposed table of new rows"
    fv = math.nan
    tr = [list(new_row)
          for new_row in zip_longest(*table, fillvalue=fv)]
    # remove fillvalues in triangular transposition
    for row in tr:
        try:
            row[row.index(fv):] = []
        except ValueError:
            continue

    return tr


if __name__ == "__main__":
    print("# TEST FILL BY ANTI-DIAGONAL\n")
    for n in range(0, 8):
        print(f"{n = }:\n")
        ans = antidiag_to_table(list(range(1, n+1)))
        pp_table(ans)
        print()

    fname, m = 'b365096.txt', 505
    print(f"\n\n# Data from {fname}, first {m} values:\n")
    ad = read_bfile(fname)
    ans = antidiag_to_table(ad[1][:m])
    pp_table(ans)
    print("\n## Transposed:\n")
    pp_table((tr:=transpose(ans)))

And antidiagonals.py is this:

# -*- coding: utf-8 -*-
"""
Anti-diagonals:

0,0 0,1 0,2 0,3
1,0 1,1 1,2 1,3
2,0 2,1 2,2 2,3
3,0 3,1 3,2 3,3

Of Square:
0,0  0,1 1,0  0,2 1,1 2,0  0,3 1,2 2,1 3,0   1,3 2,2 3,1  2,3 3,2  3,3

of Infinite table:
0,0  0,1 1,0  0,2 1,1 2,0  0,3 1,2 2,1 3,0   0,4 1,3 2,2 3,1 4,0 ...


Created on Mon Aug 21 13:36:31 2023

@author: paddy

"""
# %% Triangles

from itertools import islice


def antidiag_triangle_indices() -> tuple[int, int]:
    x = y = 0
    while True:
        yield (x, y)
        x, y = (x+1, y-1) if y else (0, x+1)


list(islice(antidiag_triangle_indices(), 15))

# %% Rectangles

from itertools import islice


def antidiag_rectangle_indices(sizex: int=4, sizey: int=4) -> tuple[int, int]:
    x = y = 0
    while True:
        yield (x, y)
        if (x, y) == (sizex - 1, sizey - 1):
            break
        x, y = (x+1, y-1)
        if x == sizex or y < 0:
            u = x + y + 1
            x, y = (0, u) if u < sizey else (u - sizey + 1, sizey - 1)

list(antidiag_rectangle_indices(3, 4))

END.

Finding a sub-list within a list, in Python

2024-03-31T10:03:00.005+01:00

Existing?

As part of a larger project, I thought I might need to search for a sub-list within a given list, and because I am lazy i did a quick google and did not like the answers I found.I started with the thought that the best algorithm for me would be to start searching from the index of the first item in the sublist and so on, but none of the googled answers used list.index.

I decided then to create my own

My version

Well I want to use list.index. If the item is not in the list then it raises an error, so I'll need a try-except block too.

I look for successive first item from the sub-list in the list and if found, accumulate the index in the answer and move on to search for the next match.

It seemed easy to add flags to:

Stop after finding a first index of the sub-list in the list.
Allow overlapping matches or not. [1,0,1] is found twice in [1,0,1,0,1] at indices 0 and 2, but only once if overlapping is not allowed

#!/bin/env python3
#%%
from typing import Any


"""
Find instance of a sub-list in a list
"""

def index_sublist(lst: list[Any],
                  sublst: list[Any],
                  only_first: bool=False,
                  non_overlapping=False,
                  ) -> list[int]:
    "Find instance of a (non-empty), sub-list in a list"
    if not sublst:
        raise ValueError("Empty sub-list")
    if not lst:
        return []
    
    first, ln = sublst[0], len(sublst)
    
    if len(lst) < ln:  # Reddit: check sizes
        return []
    
    ans, i = [], 0
    while True:
        try:
            i = lst.index(first, i)
        except ValueError:
            break
        if lst[i: i+ln] == sublst:
            ans.append(i)
            if only_first:  #Reddit: Only_first calc now fixed 
                break
        i += ln if non_overlapping else 1
    
    return ans

#%%
def test():
    assert index_sublist([], [1], only_first=False) == []
    assert index_sublist([1], [1], only_first=False) == [0]
    assert index_sublist([1,0,1], [1], only_first=False) == [0, 2]
    assert index_sublist([2,1,0,1], [1], only_first=True) == [1]
    assert index_sublist([2,1,0,1], [1, 3], only_first=False) == []
    
    assert index_sublist([1,0,1,0,1], [1,0,1],
                         only_first=False,
                         non_overlapping=False) == [0, 2]
    assert index_sublist([1,0,1,0,1], [1,0,1],
                         only_first=False,
                         non_overlapping=True) == [0]

    # Reddit: Extra checks
    assert index_sublist([2,1,0,1,2,1,2], [1, 2],
                         only_first=True) == [3]
    assert index_sublist([2,1,0,1,2,1,2], [1, 2],
                         only_first=False) == [3,5]


#%%
if __name__ == '__main__':
    test()

# %%

STOP PRESS!

A Reddit commenter spotted a bug that I fixed and added a couple of extra tests for.

End.

How not to check for a key in a dictionary.

2023-12-28T19:03:00.000+00:00

I skim the Linkedin Python group and sometimes comment.

A few days there was a poll asking for the way to check if a key is in a Python dictionary and as I write, more than half of the 900+ respondents have chosen dict.get rather than key in dict, which is the correct answer!

When pushed, it seems they are relying on dict.get returning None if the key is not in the dict, but fail to see tha if the key being tested is in the dict, but has a value of None then their test fails.

Here's some code:

# %%
mydict1 = {'a': 1, 'b': 2}
print(mydict1.get('c'))  # -> None

# %%

def key_in_dict1(key_, dict_: dict) -> bool:
    "FAULTY IMPLEMENTATION"
    return dict_.get(key_) != None

# This works
mydict1 = {'a': 1, 'b': 2}
print(key_in_dict1('a', mydict1))  # -> True
print(key_in_dict1('x', mydict1))  # -> False
# %%

# But adding a key of x with value None gives
# the wrong result
mydict2 = {'a': 1, 'b': 2, 'x': None}
print(key_in_dict1('x', mydict2))  # -> False

# %%

# The correct way is to use 'in'
def key_in_dict2(key_, dict_: dict) -> bool:
    "Pythonic IMPLEMENTATION"
    return key_ in dict_

# Tests:
print(key_in_dict2('a', mydict1))  # -> True
print(key_in_dict2('x', mydict1))  # -> False

# %%

# And now for keys with None as a value:
print(key_in_dict2('x', mydict2))  # -> True

Ugly, code-golf solution using dict.get()

You can alter the default value returned from dict.get if a key is not found. This allows one to use two calls of dict.get with different defaults to be used to give a correct solution:

def key_in_dict3(key_, dict_: dict) -> bool:
    "CODE-GOLF IMPLEMENTATION USING GETs"
    default1 = None
    default2 = ""
    return default1 != default2 \
        and ( dict_.get(key_, default1) != default1
             or dict_.get(key_, default2) != default2)

# Tests:
print(key_in_dict3('a', mydict1))  # -> True
print(key_in_dict3('x', mydict1))  # -> False
print(key_in_dict3('x', mydict2))  # -> True

Of course, don't use key_in_dict3, use version 2, i.e. key in dict.

END.

Google Bard vs Microsoft Bing Chat: Vanity questions

2023-11-24T17:18:00.001+00:00

I ask two questions about my activities online and compare the outputs from the two AI tools.

Google Bard

create a rosettacode.org task in the style of user paddy3118

Sure, here is a Rosetta Code task in the style of user paddy3118:

Task:

Write a program that generates the Godeh series. The Godeh series is an infinite sequence of integers that is defined by the following recursive formula:

godeh(n) = {
  1, if n = 1
  godeh(n // 2), if n is even
  godeh(n // 2) + godeh(n // 2 + 1), if n is odd
}

Input:

An integer n.

Output:

The first n terms of the Godeh series.

Example:

Input: 10

Output:
1
1
2
1
3
2
3
1
4
3

Solution:

def godeh(n):
  if (n == 1):
    return 1
  elif (n % 2 == 0):
    return godeh(n // 2)
  else:
    return godeh(n // 2) + godeh(n // 2 + 1)

n = 10
for i in range(1, n + 1):
  print(godeh(i))

Use code with caution. Learn more

Paddy3118's style:

Paddy3118's programs are typically well-written and easy to understand. They are also often accompanied by clear and concise explanations of the code. In this task, Paddy3118 would likely provide a brief explanation of the Godeh series and the recursive formula used to define it. He would then provide a Python implementation of the program, along with an example of how to use it.

Here is an example of a task that Paddy3118 might write:

Task:

Write a program that generates the prime numbers up to a specified input.

Input:

An integer n.

Output:

A list of all the prime numbers up to n.

Example:

Input: 100

Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

Solution:

def is_prime(n):
  if (n <= 1):
    return False
  if (n <= 3):
    return True
  if (n % 2 == 0 
or n % 3 == 0):
    return 
False
  i = 5
  
while (i * i <= n):
    if (n % i == 0 
or n % (i + 2) == 0):
      return 
False
    i += 6
  
return 
True

def 
prime_numbers_up_to(n):
  primes = []
  for i in 
range(2, n + 1):
    if (is_prime(i)):
      primes.append(i)
  return primes

n = 100
primes = prime_numbers_up_to(n)
print(primes)

Use code with caution. Learn more

Microsoft Bing Chat

create a rosettacode.org task in the style of user paddy3118

…

I’m sorry, but I am not able to create a Rosetta Code task for you. However, I can provide you with some information about Paddy3118, a prolific contributor to Rosetta Code who has created many tasks in the past. Paddy3118 is a pseudonym for Donald McCarthy, who works with Electronic Design Automation tools creating design flows, filling in the gaps with novel tools, and generally helping to design and verify integrated circuits ¹. He has created 392 new pages on Rosetta Code, of which 148 were not Talk: pages, from approximately 4523 edits ¹. You can find more information about Paddy3118 and his contributions to Rosetta Code on his user page ¹. I hope this information helps!

Learn more

1rosettacode.org 2rosettacode.org 3rosettacode.org 4rosettacode.org 5en.wikipedia.org

8of30

Bard vs Bing

The actual formatting does not copy-paste well in either case.

Bard attempts to write a task in my style, and gives some tips on my style.
Bing cannot create the task, but does extract some info on me.
Clearly a Bard win.
Flattery will get them nowhere! (But keep it up mind you).

Google Bard #2

what does paddy3118 do online?

Paddy3118 is an active participant in the online programming community, particularly in areas related to mathematics and algorithms. They have authored numerous programs and solutions to coding challenges, demonstrating expertise in various programming languages and problem-solving techniques. Paddy3118's contributions can be found across various platforms, including:

GitHub: Paddy3118 has a notable presence on GitHub, where they have authored over 50 programs covering a wide range of topics, including:

Mathematical algorithms: Implementing various mathematical concepts and algorithms, such as the Godeh series generator, Collatz conjecture tester, and prime number generator.
Problem-solving challenges: Tackling coding challenges from platforms like Codewars, showcasing their ability to solve complex problems efficiently.
Data manipulation and analysis: Demonstrating proficiency in handling and analyzing data through programs like the factorial calculator and Fibonacci sequence generator.

Their GitHub repositories often include detailed explanations of the code, making them valuable learning resources for other programmers.

Stack Overflow: Paddy3118 is an active contributor on Stack Overflow, a question-and-answer website for programmers. They have answered over 1,000 questions with an average score of 4.5 out of 5, indicating their deep understanding of programming concepts and their willingness to share their knowledge with others.

Codewars: Paddy3118 is a skilled competitor on Codewars, a coding challenge platform. They have solved over 1,000 coding challenges and currently rank in the top 1% of users. This achievement highlights their ability to think creatively and apply their programming expertise to solve complex problems.

GitHub Discussions: Paddy3118 actively engages in discussions on GitHub Discussions, a forum where users discuss technical topics. They provide helpful advice, share insights, and contribute to the overall knowledge base of the community.

Overall, Paddy3118 is a valuable member of the online programming community, recognized for their contributions to various platforms and their willingness to share their knowledge and expertise with others.

Microsoft Bing Chat #2

what does paddy3118 do online?

…

I’m sorry, but I could not find any information about what Paddy3118 does online. However, I can tell you that Paddy3118 is a pseudonym for Donald McCarthy, who works with Electronic Design Automation tools creating design flows, filling in the gaps with novel tools, and generally helping to design and verify integrated circuits ¹. He has created 392 new pages on Rosetta Code, of which 148 were not Talk: pages, from approximately 4523 edits ¹. You can find more information about Paddy3118 and his contributions to Rosetta Code on his user page ¹. I hope this information helps!

Learn more

1you.dj 2photopea.com 3silver.urih.com

10of30

Bard vs Bing #2

Bard has found me on numerous sites and attempted to summarise my contributions to them and comments that I like to help
Bing states up front that it could not find me but then quotes an extract based on my RosettaCode talk page.
None mention Python!?!
Another win for Bard.

At the moment, just using these easy-access AI tools, Bard seems to give better results.

It's still early days however!

END.

"Staring at assorted perms and Python", or "The Ranking and Unranking of Lexical Permutations"

2023-08-26T09:48:00.002+01:00

I have been playing with encryptions and then series that involved permutations. Submitting my Godeh series to the OEIS made me look closer at the permutations that Python produces, and at maybe reducing those permutations to a single integer.

(p.s. Post is best read on a larger than portrait phone screen)

Here are the start of the permutations used in the Godeh series rows:

(0,)
(0, 1)
(1, 0)
(0, 1, 2)
(0, 2, 1)
(1, 0, 2)
(1, 2, 0)
(2, 0, 1)
(2, 1, 0)
(0, 1, 2, 3)
(0, 1, 3, 2)
(0, 2, 1, 3)
(0, 2, 3, 1)
(0, 3, 1, 2)
(0, 3, 2, 1)
(1, 0, 2, 3)
(1, 0, 3, 2)
(1, 2, 0, 3)
(1, 2, 3, 0)
(1, 3, 0, 2)
(1, 3, 2, 0)
(2, 0, 1, 3)
(2, 0, 3, 1)
(2, 1, 0, 3)
(2, 1, 3, 0)
(2, 3, 0, 1)
(2, 3, 1, 0)
(3, 0, 1, 2)
(3, 0, 2, 1)
(3, 1, 0, 2)
(3, 1, 2, 0)
(3, 2, 0, 1)
(3, 2, 1, 0)

It can be generated by this:

from itertools import permutations, chain
from math import factorial
from typing import Sequence
from pprint import pp

# %% perms in Python generation order
generation_order = []
for width in range(1, 5):
    for perm in permutations(range(width)):
        print(perm)
        generation_order.append(perm)

I truncate at all the perms for four items, but you get the idea - it could extend in this way without limit.

Ordered by...

Successive groups of permutations are for incrementing items permuted. Within each group, of permutations of the same amount of items, the order is lexicographical, i.e. the Python sorted order,

I stared at, and messed around with the ordered perms of four items above, looking for patterns. But lets just confirm the order is sort by width then lexicographical:

wtl = sorted(set(generation_order),
             key=lambda x:(len(x), x))
if generation_order == wtl:
    print(' == width_then_lexicographical ordering')
else:
    print(' != width_then_lexicographical ordering of:')
    pp(wtl)

# Prints:  == width_then_lexicographical ordering

Ranks

I remembered that I had created a Rank of a permutation task some time ago, (2012), and that there was something "off" with the ranking.

Sure enough the issue was that a rank number gets transformed into a unique permutation, and vice-versa, knowing the start configuration and given any permutation of it then yoy can generate its rank number - But the order of permutations for increasing rank number does not have to be that lexical sort used by Python, and indeed the Myrvold & Ruskey or Trottor & Johnson algorithms I converted to Python are not in lexical order.

So now I know I needed some way to rank and un-rank a lexically ordered set of permutations, as generated by Python.

Something to look at and talk about

Lets assign ranks to Python perms and column headings to make it easier to talk about:

from string import ascii_uppercase

n = 4  # large enough to hopefully show patterns, but not too large
print(' # :  ' + '  '.join(char for char in ascii_uppercase[:n]))
for rank, perm in enumerate(permutations(range(n))):
    print(f"{rank:>2} : {perm}")
"""
Outputs:

 # :  A  B  C  D
 0 : (0, 1, 2, 3)
 1 : (0, 1, 3, 2)
 2 : (0, 2, 1, 3)
 3 : (0, 2, 3, 1)
 4 : (0, 3, 1, 2)
 5 : (0, 3, 2, 1)
 6 : (1, 0, 2, 3)
 7 : (1, 0, 3, 2)
 8 : (1, 2, 0, 3)
 9 : (1, 2, 3, 0)
10 : (1, 3, 0, 2)
11 : (1, 3, 2, 0)
12 : (2, 0, 1, 3)
13 : (2, 0, 3, 1)
14 : (2, 1, 0, 3)
15 : (2, 1, 3, 0)
16 : (2, 3, 0, 1)
17 : (2, 3, 1, 0)
18 : (3, 0, 1, 2)
19 : (3, 0, 2, 1)
20 : (3, 1, 0, 2)
21 : (3, 1, 2, 0)
22 : (3, 2, 0, 1)
23 : (3, 2, 1, 0)
"""

Patterns

I can't go through the many false trails I had in coming to my result, (there were several), but I finally noticed that The digits in column A stayed the same for 6 successive rows and six is three factorial. Column B stays the same for 2 successive rows i.e. 2! , I guessed 1! for column C and D is what's left.

A lot of deleted code later I found that, kinda, after you pop the first digit out of the starting perm at rank 0, which would be rank // (4 - 1)!; you could take the second digit out using indices using decreasing factorials, but also modulo the decreasing numbers left to index.

Hard to write a textual description of, but I was very chuffed to work out the following on my own:

def nth_perm(initial: Sequence[int], n: int) -> tuple[int]:
    init = list(initial)
    a = len(initial)
    fac_divs = tuple((n // factorial(a - j)) % (a - j + 1)
                     for j in range(1, a))

    return tuple([init.pop(indx) for indx in fac_divs]
                 + init)

nth_perm(range(4), 11)

And indeed the nth perm of (0, 1, 2 ,3) produced for n = 11 was (1, 3, 2, 0).

Check

I decided to check my perm-from-rank generator against the Python permutations:

start = range(4)
for i, perm in enumerate(permutations(range(4))):
    assert perm == nth_perm(start, i)

Generate ranking from arbitrary perm

This had its own aha moments, and eventually I created the following code generating the rank of the sorted perm:

def s_perm_rank(p: Sequence[int]) -> int:
    """
    Ranking of perm p in the sorted sequence of perms of the ints 0 .. len(p)-1
    p must be a permutation of the integers 0 .. (len(p)-1)
    """
    this = list(p)
    a = len(this)
    init = list(range(a))  # Perm of rank 0
    assert set(this) == set(init), "p must be perm of the ints 0 .. (len(p)-1)"

    n, f = 0, a
    while this:
        f -= 1
        n += init.index(this[0]) * factorial(f)
        init.remove(this.pop(0))

    return n

Roundtrip Check

To check I roundtrip from rank to perm to rank, and also check against the Python permutations ranking

n = 4
init = range(n)
for i, py_perm in enumerate(permutations(range(n))):
    s_perm = nth_perm(init, i)
    print(f"{i:>2}: {py_perm=} {py_perm==s_perm=} {i==s_perm_rank(s_perm)=}")

Output:

 0: py_perm=(0, 1, 2, 3) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
 1: py_perm=(0, 1, 3, 2) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
 2: py_perm=(0, 2, 1, 3) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
 3: py_perm=(0, 2, 3, 1) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
 4: py_perm=(0, 3, 1, 2) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
 5: py_perm=(0, 3, 2, 1) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
 6: py_perm=(1, 0, 2, 3) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
 7: py_perm=(1, 0, 3, 2) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
 8: py_perm=(1, 2, 0, 3) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
 9: py_perm=(1, 2, 3, 0) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
10: py_perm=(1, 3, 0, 2) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
11: py_perm=(1, 3, 2, 0) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
12: py_perm=(2, 0, 1, 3) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
13: py_perm=(2, 0, 3, 1) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
14: py_perm=(2, 1, 0, 3) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
15: py_perm=(2, 1, 3, 0) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
16: py_perm=(2, 3, 0, 1) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
17: py_perm=(2, 3, 1, 0) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
18: py_perm=(3, 0, 1, 2) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
19: py_perm=(3, 0, 2, 1) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
20: py_perm=(3, 1, 0, 2) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
21: py_perm=(3, 1, 2, 0) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
22: py_perm=(3, 2, 0, 1) py_perm==s_perm=True i==s_perm_rank(s_perm)=True
23: py_perm=(3, 2, 1, 0) py_perm==s_perm=True i==s_perm_rank(s_perm)=True

It works!

Tidy

I searched online, read more mathematical articles in and around the subject of permutations and decided to tidy things up with renamings and changes to function arguments. For example, I found that lex is used as a short replacement for lexicographic ordering, and that a perm is usually generated from the integer width of the perm and the rank required, rather than giving the perm of rank 0 and the rank

The final, (?), version is:

#%% Tidy
def lex_perm_unrank(width: int, rank: int) -> tuple[int]:
    """
    Generate the lexicographic-permutation of ints 0 .. width-1 of given rank.
    
    Author: Donald S. McCarthy "Paddy3118"
    Date: August 2023
    """    
    initial = list(range(width))
    indices = [(rank // factorial(width - j)) % (width - j + 1)
               for j in range(1, width+1)]

    return tuple(initial.pop(index) for index in indices)


# %% Check
for w in range(7):
    for py_rank, py_perm in enumerate(permutations(range(w))):
        lex_perm = lex_perm_unrank(w, py_rank)
        assert lex_perm == py_perm
print("lex_perm_unrank consistent with Python permutations order.")


# %% Tidier
def lex_perm_rank(p: Sequence[int]) -> int:
    """
    Rank of perm p in the lexicographic ordered perms of ints 0 .. len(p)-1
    p must be a permutation of the integers 0 .. (len(p)-1)
    
    Author: Donald S. McCarthy "Paddy3118"
    Date: August 2023
    """
    perm = list(p)
    width = len(perm)
    initial = list(range(width))  # Perm of rank 0
    assert set(perm) == set(initial), \
           "p must be a permutation of the integers 0 .. (len(p)-1)"

    rank, f = 0, width
    while perm:
        f -= 1
        rank += initial.index(perm[0]) * factorial(f)
        initial.remove(perm.pop(0))

    return rank

# %% Check both
for w in range(7):
    for py_rank, py_perm in enumerate(permutations(range(w))):
        lex_perm = lex_perm_unrank(w, py_rank)
        lex_rank = lex_perm_rank(lex_perm)
        assert lex_perm == py_perm and lex_rank == py_rank
print("lex_perm_rank/unrank consistent with Python permutations order.")

END.

The Godeh Series, Python, and OEIS

2023-08-21T21:34:00.000+01:00

(Best viewed on screens larger than a portrait phone)

In my previous post I generalised a form of encryption to do this:

Task Encrypt grouping chars by mod M, N times:
    Given an input string S; Some order of the K integers 0..K-1 called
    M, where K <= length_of(S); and an integer N > 0

    First form groups:
        Form group G[0] by concatenating every K'th item of S starting from index 0
        Form group G[1] by concatenating every K'th item of S starting from index 1
        ...
        Form group G[K-1] by concatenating every K'th item of S starting from index K-1

    Group Concatenation:
        Concatenate all the groups G in the order given M to form the partial result P
        P(1) = concatenate_left_to_right(G[i] for i in M)
    Repetition:
        Substituting P for S, repeat the above process N times to form the final
        encrypted string P(N)

    Corner cases:
        If S is empty or N is 0 then return S.
        If M is empty return the empty version of sequence S (empty string).

Examples:
    encrypt6('012345', [1, 0], 0) => '012345'
    encrypt6('012345', [1, 0], 1)          => '135024'
    encrypt6('012345', [1, 0], 2)                   => '304152'
    encrypt6('012345', [1, 0], 3)                            => '012345'
    encrypt6('012345', [1, 0], 4)                                     => '135024'

    encrypt6('01234567', [2, 1, 0], 0) => '01234567'
    encrypt6('01234567', [2, 1, 0], 1)          => '25147036'
    encrypt6('01234567', [2, 1, 0], 2)                   => '10576243'
    encrypt6('01234567', [2, 1, 0], 3)                            => '52063174'
    encrypt6('01234567', [2, 1, 0], 4)                                     => '01234567'

Created from Sun Aug  6 09:10:10 2023

@author: paddy3118

Repetitions

In the example of we see that for n = 3, the output of encrypt6('012345', [1, 0], 3) starts repeating for higher n.

I have found, but have not proved, that

encrypt6(S, M, n) == encrypt6(S, M, (n modulo x) )

for some x.

I.e. outputs will eventually equal the input to form a repettitive cycle of outputs.

The length of the cycle, x before repetition does not seem to be simple.

The Godeh series

This series calculates the length of the cycles of encrypt6 untill the input is produced again.

Godeh series - Original definition by Author
    1. Given:
    * An ordered list of S distinct items, S >= 1.
    * Some permutation of the K integers 0..K-1 called M, where 1 <= K <= S.

    2. Form groups:
    * Form group G[0] by concatenating every K'th item of S starting from index 0
    * Form group G[1] by concatenating every K'th item of S starting from index 1
      ...
    * Form group G[K-1] by concatenating every K'th item of S starting from index K-1

    3. Concatenate groups:
    * Concatenate all the groups G in the order given by M to form the first partial
      result P(1)
        P(1) = concatenate_left_to_right(G[i] for i in M)

    4. Repetition:
        Substituting P for S, repeat the above process until the order of items in P
        equals the initial order, S.

    5. Result
    * Godeh(S, M) = the number of repetitions needed of steps 1-through-4

Created on Thu Aug 17 19:40:32 2023

@author: paddy3118

Python code to generate one term

from collections.abc import Sequence
from collections import defaultdict
from itertools import zip_longest, permutations, chain, islice
from functools import cache, lru_cache
from pprint import pp
from string import ascii_letters, digits
from time import sleep
from typing import Callable, Any, overload
import sys
sys.path.append(r'C:\Users\paddy\Google Drive\Code\oeisify')
from oeis_request import oeis_request
from antidiagonals import antidiag_triangle_indices


def godeh(s: int, m: Sequence[int]) -> int:
    "Godeh function Godeh(s, m)"
    k = len(m)
    assert s >= k
    assert set(range(k)) == set(m), \
           f"Sequence m of length {k} should contain a permutation of all the " \
           f"numbers 0..{k-1} inclusive."

    s_init = list(range(s))
    n, s = 0, None
    while s != s_init:
        if n == 0:
            s = s_init
        n += 1
        s = sum((s[offset::k] for offset in m),
                start=[])

    return n

Table of values

Each row should represent M as "successive permutations" - i.e perms of (0,) then perms of (0,1), perms of (0,1,2) ...

Successive values in each row are then godeh(S,M) where S = len(M), len(M)+1, len(M+2), ...

Code:

GodehTableType = dict[tuple[int], list[int]]

def godeh_table(max_s: int=15, max_permed_ints: int=3) -> GodehTableType:
    """
    Create a table of Godeh series values

    Parameters
    ----------
    max_s : int, optional
        Godeh(s, m) for s=len(m) to len(m) + max_s, for all m. The default is 10.
    max_permed_ints : int, optional
        All perms of (0), (0, 1), (0, 1, 2), ... (0, 1, ..max_permed_ints).
        The default is 3.

    Returns
    -------
    dict[
         tuple[int],    # M
         list[int]]     # Godeh(S, M) for S = len(M), len(M)+1, ...
        Tablulated Godeh values.

    """
    table = defaultdict(list)
    for m in chain(*(permutations(range(i)) for i in range(1, max_permed_ints+1))):
        for s in range((k:=len(m)), max_s + k):
            table[(m)].append(godeh(s, m))
    return dict(table)

godeh_series = godeh_table(20, 3)
# %%
pp(godeh_series, width=222)

Output:

{(0,): [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
 (0, 1): [1, 2, 2, 4, 4, 3, 3, 6, 6, 10, 10, 12, 12, 4, 4, 8, 8, 18, 18, 6],
 (1, 0): [2, 2, 4, 4, 3, 3, 6, 6, 10, 10, 12, 12, 4, 4, 8, 8, 18, 18, 6, 6],
 (0, 1, 2): [1, 3, 4, 4, 6, 2, 2, 6, 5, 5, 11, 6, 6, 15, 16, 16, 52, 4, 4, 38],
 (0, 2, 1): [2, 2, 2, 4, 6, 6, 4, 4, 4, 21, 3, 3, 30, 4, 4, 90, 18, 18, 24, 5],
 (1, 0, 2): [2, 2, 4, 4, 6, 4, 4, 4, 21, 21, 3, 30, 30, 8, 90, 90, 18, 24, 24, 10],
 (1, 2, 0): [3, 3, 4, 6, 6, 4, 6, 6, 5, 11, 11, 6, 15, 15, 16, 52, 52, 4, 38, 38],
 (2, 0, 1): [3, 4, 4, 6, 2, 2, 6, 5, 5, 11, 6, 6, 15, 16, 16, 52, 4, 4, 38, 11],
 (2, 1, 0): [2, 2, 4, 6, 6, 4, 4, 4, 21, 3, 3, 30, 4, 4, 90, 18, 18, 24, 5, 5]}

Duplicate rows?

I suspect that there are duplicate rows of values so generate a larger table and search for row duplication, (I use text searches but it gets the job done).

Code:

def find_reps_in_godeh(g: GodehTableType,
                       _slice_start: int=10,
                       _slice_end: int=18) -> None:
    """
    Prints information on duplicated Godeh table rows.

    Parameters
    ----------
    g : GodehTableType
        Table of Godeh(S,M) values.
    _slice_start : int, optional
        Defines slice of a rows values to search for in other rows of g.
        The default is 10.
    _slice_end : int, optional
        Defines slice of a rows values to search for in other rows of g.
        The default is 18.

    Returns
    -------
    None
    """
    # Uses a string search on stringified row values between commas
    searchable = {m: ',' + str(values)[1:-1].replace(' ','') + ','
                  for m, values in g.items()}
    keys, rows = list(searchable), list(searchable.values())
    # Search for every row in all "higher" rows
    maxrows = len(rows)
    for i in range(maxrows - 1):
        isearch = ',' + str(g[keys[i]][10:17])[1:-1].replace(' ','') + ','
        for k in range(i+1, maxrows):
            if isearch in (rowk:=rows[k]):
                offset = _slice_start - rowk[:rowk.index(isearch)].count(',')
                print(f"Found {keys[i]} in {keys[k]} at {offset = }")


godeh_series = godeh_table(40, 5)  # more values
find_reps_in_godeh(godeh_series, 10, 20)

Output:

Found (0, 1) in (1, 0) at offset = 1
Found (0, 1, 2) in (2, 0, 1) at offset = 1
Found (0, 2, 1) in (2, 1, 0) at offset = 1
Found (0, 1, 2, 3) in (3, 0, 1, 2) at offset = 1
Found (0, 1, 3, 2) in (3, 0, 2, 1) at offset = 1
Found (0, 2, 1, 3) in (3, 1, 0, 2) at offset = 1
Found (0, 2, 3, 1) in (3, 1, 2, 0) at offset = 1
Found (0, 3, 1, 2) in (3, 2, 0, 1) at offset = 1
Found (0, 3, 2, 1) in (3, 2, 1, 0) at offset = 1
Found (0, 1, 2, 3, 4) in (4, 0, 1, 2, 3) at offset = 1
Found (0, 1, 2, 4, 3) in (4, 0, 1, 3, 2) at offset = 1
Found (0, 1, 3, 2, 4) in (4, 0, 2, 1, 3) at offset = 1
Found (0, 1, 3, 4, 2) in (4, 0, 2, 3, 1) at offset = 1
Found (0, 1, 4, 2, 3) in (4, 0, 3, 1, 2) at offset = 1
Found (0, 1, 4, 3, 2) in (4, 0, 3, 2, 1) at offset = 1
Found (0, 2, 1, 3, 4) in (4, 1, 0, 2, 3) at offset = 1
Found (0, 2, 1, 4, 3) in (4, 1, 0, 3, 2) at offset = 1
Found (0, 2, 3, 1, 4) in (4, 1, 2, 0, 3) at offset = 1
Found (0, 2, 3, 4, 1) in (4, 1, 2, 3, 0) at offset = 1
Found (0, 2, 4, 1, 3) in (4, 1, 3, 0, 2) at offset = 1
Found (0, 2, 4, 3, 1) in (4, 1, 3, 2, 0) at offset = 1
Found (0, 3, 1, 2, 4) in (4, 2, 0, 1, 3) at offset = 1
Found (0, 3, 1, 4, 2) in (4, 2, 0, 3, 1) at offset = 1
Found (0, 3, 2, 1, 4) in (4, 2, 1, 0, 3) at offset = 1
Found (0, 3, 2, 4, 1) in (4, 2, 1, 3, 0) at offset = 1
Found (0, 3, 4, 1, 2) in (4, 2, 3, 0, 1) at offset = 1
Found (0, 3, 4, 2, 1) in (4, 2, 3, 1, 0) at offset = 1
Found (0, 4, 1, 2, 3) in (4, 3, 0, 1, 2) at offset = 1
Found (0, 4, 1, 3, 2) in (4, 3, 0, 2, 1) at offset = 1
Found (0, 4, 2, 1, 3) in (4, 3, 1, 0, 2) at offset = 1
Found (0, 4, 2, 3, 1) in (4, 3, 1, 2, 0) at offset = 1
Found (0, 4, 3, 1, 2) in (4, 3, 2, 0, 1) at offset = 1
Found (0, 4, 3, 2, 1) in (4, 3, 2, 1, 0) at offset = 1

Observation, Duals!

The first line is saying that M = (1, 0) values are the same as those for M = (0, 1); just with an offset of being one term over.

There are only row duals found. no rows data is replicated in more than one other row.
All the offsets are one.
Duals are between rows of the same size perms in M; with the Leftmost value of M being max(M) paired with an M where the leftmost value of M is min(M).

OEIS Searching by hand: Abulsme

I put a few selections of the table rows into OEIS by hand and found that, (but not proved), I can generate all the Abulsme function values from my Godeh function!

def a105272_table(max_s: int=10, max_permed_ints: int=40):
    """
    Abulsme function by different means: http://www.abulsme.com/function.html
    == Godeh(S, M) where All M's are ints sorted high to low
    """
    table = defaultdict(list)
    for m in (range(i)[::-1] for i in range(1, max_permed_ints+1)):
        for s in range(k:=len(m), max_s + k):
            table[k].append(godeh(s, m))
    return dict(table)

a105272 = a105272_table(20, 15)
# %%
pp(a105272, width=222)

Table of Abulsme values produced

{1: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
 2: [2, 2, 4, 4, 3, 3, 6, 6, 10, 10, 12, 12, 4, 4, 8, 8, 18, 18, 6, 6],
 3: [2, 2, 4, 6, 6, 4, 4, 4, 21, 3, 3, 30, 4, 4, 90, 18, 18, 24, 5, 5],
 4: [2, 2, 4, 7, 3, 3, 8, 10, 6, 6, 20, 4, 4, 4, 132, 30, 3, 3, 24, 252],
 5: [2, 2, 4, 7, 15, 5, 5, 12, 40, 45, 4, 4, 12, 8, 16, 6, 6, 16, 28, 4],
 6: [2, 2, 4, 14, 6, 10, 12, 12, 7, 15, 56, 90, 9, 9, 51, 21, 105, 23, 5, 5],
 7: [2, 2, 4, 14, 6, 12, 30, 4, 4, 20, 15, 70, 30, 18, 10, 10, 48, 72, 24, 25],
 8: [2, 2, 4, 14, 6, 13, 13, 24, 8, 8, 40, 70, 91, 18, 18, 22, 20, 20, 44, 312],
 9: [2, 2, 4, 14, 6, 13, 15, 15, 63, 9, 9, 20, 18, 85, 30, 24, 60, 44, 3, 3],
 10: [2, 2, 4, 14, 6, 26, 16, 10, 18, 12, 6, 6, 42, 21, 22, 42, 68, 165, 182, 30],
 11: [2, 2, 4, 14, 6, 26, 16, 60, 78, 48, 20, 22, 22, 12, 25, 390, 90, 21, 360, 42],
 12: [2, 2, 4, 14, 6, 26, 16, 65, 120, 84, 21, 36, 20, 20, 44, 132, 120, 24, 100, 108],
 13: [2, 2, 4, 14, 6, 26, 16, 65, 168, 168, 120, 20, 25, 9, 9, 60, 24, 187, 616, 84],
 14: [2, 2, 4, 14, 6, 26, 16, 130, 24, 60, 130, 144, 20, 26, 28, 28, 15, 60, 84, 84],
 15: [2, 2, 4, 14, 6, 26, 16, 130, 24, 204, 24, 144, 170, 170, 40, 10, 10, 30, 510, 34]}

Auto-check OEIS

I hooked in my oeis search routine and searched for each row of the Godeh table of values.

Code:

def oeis_check_values(table: dict[Any, list[int]],
                      _slice_start: int=4,
                      _slice_end: int=12,
                      _secs_between_searches: int=5,
                      _show_not_found=False
                      ) -> None:
    "Search for each row of table in OEIS"
    
    for key, vals in table.items():
        oeis = oeis_request(str(vals[_slice_start:_slice_end])[1:-1])
        found = [k for k in oeis if k[0] == 'A']
        if found:
            print(f"{key}: is found in OEIS {', '.join(found)}")
        elif _show_not_found:
            print(f"{key}: Not found in oeis.")
            
        sleep(_secs_between_searches)

godeh_series = godeh_table(40, 5)  # more values
oeis_check_values(godeh_series, _secs_between_searches=2)

It takes some time to run, but eventually I get these results:

Output:

(0,): is found in OEIS A000012, A055642, A135010, A057427, A138121, A000030, A047999, A071625, A079944, A228351
(0, 1): is found in OEIS A024222
(1, 0): is found in OEIS A024222
(0, 2, 1): is found in OEIS A118960
(2, 1, 0): is found in OEIS A118960
(0, 3, 2, 1): is found in OEIS A120280
(3, 2, 1, 0): is found in OEIS A120280
(0, 4, 3, 2, 1): is found in OEIS A120363
(4, 3, 2, 1, 0): is found in OEIS A120363
(0, 5, 4, 3, 2, 1): is found in OEIS A120654
(5, 4, 3, 2, 1, 0): is found in OEIS A120654

Ignoring (0,) which is just all the ones; All the other matches are selections from the Abulsme function,

a105272, when M is either reverse sorted, or its dual.

OEIS ENTRY?

It seems that the Godeh Series is more general than a105272 so I thought I'd submit it to OEIS. It's submitted and allocated as A365096. Fingers crossed!

Antidiagonalisation

When you submit a table of values the DATA entry must be in antidiagonalisation format. I did a quick module to provide that functionality:

Module antidiagonals.py:

# -*- coding: utf-8 -*-
"""
Anti-diagonals:

0,0 0,1 0,2 0,3
1,0 1,1 1,2 1,3
2,0 2,1 2,2 2,3
3,0 3,1 3,2 3,3

Of Square:
0,0  0,1 1,0  0,2 1,1 2,0  0,3 1,2 2,1 3,0   1,3 2,2 3,1  2,3 3,2  3,3

of Infinite table:
0,0  0,1 1,0  0,2 1,1 2,0  0,3 1,2 2,1 3,0   0,4 1,3 2,2 3,1 4,0 ...


Created on Mon Aug 21 13:36:31 2023

@author: paddy

"""
# %% Triangles

from itertools import islice


def antidiag_triangle_indices() -> tuple[int, int]:
    x = y = 0
    while True:
        yield (x, y)
        x, y = (x+1, y-1) if y else (0, x+1)

list(islice(antidiag_triangle_indices(), 15))


# %% Rectangles

from itertools import islice


def antidiag_rectangle_indices(sizex: int=4, sizey: int=4) -> tuple[int, int]:
    x = y = 0
    while True:
        yield (x, y)
        if (x, y) == (sizex - 1, sizey - 1):
            break
        x, y = (x+1, y-1)
        if x == sizex or y < 0:
            u = x + y + 1
            x, y = (0, u) if u < sizey else (u - sizey + 1, sizey - 1)


list(antidiag_rectangle_indices(3, 4))

Antidiagonalisation test

def antidiagonalise(table: dict[Any, list[int]]) -> int:
    maxx = len(xkeys:=list(table))
    maxy = min(len(row) for row in table.values())
    antigen = antidiag_triangle_indices().__next__
    x, y = antigen()
    while x < maxx and y < maxy:
        yield table[xkeys[x]][y]
        x, y = antigen()

t = {f"{x = }": [(x, y) for y in range(7)] for x in range(3)}

print('# Sample table')
pp(t, width=100, sort_dicts=False)

print("\nAntidiagonalised\n", list(antidiagonalise(t)))

Test output:

# Sample table
{'x = 0': [(0, 0), (0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6)],
 'x = 1': [(1, 0), (1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6)],
 'x = 2': [(2, 0), (2, 1), (2, 2), (2, 3), (2, 4), (2, 5), (2, 6)]}

Antidiagonalised
 [(0, 0), (0, 1), (1, 0), (0, 2), (1, 1), (2, 0), (0, 3), (1, 2), (2, 1)]

Better formatted table for comments

godeh_series2 = godeh_table(20, 3)

for M, v in godeh_series2.items():
    m = f"{M = }"
    print(f"{m:>13s}: {str(v)[1:-1].replace(' ', '')},...")

Output:

     M = (0,): 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...
   M = (0, 1): 1,2,2,4,4,3,3,6,6,10,10,12,12,4,4,8,8,18,18,6,...
   M = (1, 0): 2,2,4,4,3,3,6,6,10,10,12,12,4,4,8,8,18,18,6,6,...
M = (0, 1, 2): 1,3,4,4,6,2,2,6,5,5,11,6,6,15,16,16,52,4,4,38,...
M = (0, 2, 1): 2,2,2,4,6,6,4,4,4,21,3,3,30,4,4,90,18,18,24,5,...
M = (1, 0, 2): 2,2,4,4,6,4,4,4,21,21,3,30,30,8,90,90,18,24,24,10,...
M = (1, 2, 0): 3,3,4,6,6,4,6,6,5,11,11,6,15,15,16,52,52,4,38,38,...
M = (2, 0, 1): 3,4,4,6,2,2,6,5,5,11,6,6,15,16,16,52,4,4,38,11,...
M = (2, 1, 0): 2,2,4,6,6,4,4,4,21,3,3,30,4,4,90,18,18,24,5,5,...

END.

Generalising a simple encryption using Python

2023-08-20T11:32:00.002+01:00

(Best viewed on screens larger than a portrait phone)

Recap

In my last but one post I coded solutions to the task:

Implement a pseudo-encryption algorithm which given a string S and an integer N 
concatenates all the odd-indexed characters of S with all the even-indexed 
characters of S, this process should be repeated N times.

First implementation

My first implementation didn't use index variables to walk through strings, but instead used slicing and zip_longest:

from itertools import zip_longest

def encrypt(s: str, n: int) -> str:
    if not s or n <= 0:
        return s
    for _ in range(n):
        s = s[1::2] + s[0::2]
    return s


def decrypt(s: str, n: int) -> str:
    if not s or n <= 0:
        return s
    l2 = len(s) // 2
    for _ in range(n):
        s = ''.join(''.join(x)
                    for x in
                    zip_longest(s[l2:], s[:l2],
                                fillvalue=''))
    return s

N modulo

In that same post I found that encryption repeats after some value found to be associated with
https://oeis.org/A002326 So if n is ls large, work will be reduced.

from functools import cache

@cache
def a002326(n: int) -> int:
    "Least m > 0 such that 2n+1 divides 2^m-1"
    n21 = 2 * n + 1
    m, m2 = 1, 2
    while (m2 - 1) % n21:
        m, m2 = m+1, m2 * 2
    return m

def encrypt2(s: str, n: int) -> str:
    if not s or n <= 0:
        return s
    l2 = len(s) // 2
    for _ in range(n % a002326(l2)):  # Modulo reduction on large n
        s = s[1::2] + s[0::2]
    return s


def decrypt2(s: str, n: int) -> str:
    if not s or n <= 0:
        return s
    l2 = len(s) // 2
    for _ in range(n % a002326(l2)):  # Modulo reduction on large n
        s = ''.join(''.join(x)
                    for x in
                    zip_longest(s[l2:], s[:l2],
                                fillvalue=''))
    return s

New Stuff

this marks the end of the previous posts precis.

Task generalisation #1: Reverse

Reverse the concatenation of the original.

The new task descripption:

New Task:
     given a string S and an integer N concatenates all the EVEN-indexed
     characters of S with all the ODD-indexed characters of S, this process
     should be repeated N time. (Prevous was ODD-EVEN).

Examples:
    encrypt('012345', 1) => '024135'
    encrypt('012345', 2)             => '043215'
    encrypt('012345', 3)                         => '031425'

    encrypt('0123456', 1) => '0246135'
    encrypt('0123456', 2)              => '0415263'
    encrypt('0123456', 3)                           => '0123456'

Its code

def encrypt3(s: str, n: int) -> str:
    if not s or n <= 0:
        return s
    for _ in range(n):
        s = s[::2] + s[1::2]
    return s


def decrypt3(s: str, n: int) -> str:
    if not s or n <= 0:
        return s
    l2 = (len(s) + 1) // 2  # Odd char now in first half for odd len(s)
    for _ in range(n):
        s = ''.join(''.join(x)
                    for x in
                    zip_longest(s[:l2], s[l2:],
                                fillvalue=''))
    return s

N modulo for Reverse

I looked again for a value for n after which encryptions repeat and found a similar dependency on a002326 but with a subtle change in offset

def encrypt4(s: str, n: int) -> str:
    if not s or n <= 0:
        return s
    l2 = (l1:=len(s)) // 2
    for _ in range(n % a002326((l1 - 1) // 2)):
        s = s[::2] + s[1::2]
    return s


def decrypt4(s: str, n: int) -> str:
    if not s or n <= 0:
        return s
    l2 = ((l1:=len(s)) + 1) // 2  # Odd char now in first half for odd len(s)
    for _ in range(n % a002326((l1 - 1) // 2)):
        s = ''.join(''.join(x)
                    for x in
                    zip_longest(s[:l2], s[l2:],
                                fillvalue=''))
    return s

Three groups

A task description like the first, that gave us encrypt1, but this time for three groups.

I am thinking about later generalisations and realise that after splitting into three, there could be 3! = 6 possible ways of concatenating groups, so call this encrypt_210 as 210 is the order used here:

New Task encrypt_210:
    Given a string S and an integer N select chars at positions div 3
    and concatenate all the remainder-2, then rem-1 then rem-0 ordered
    groups of character. This process should be repeated N times to
    form the encrpted string.
    If N is 0 or the string is null return the input string.

Examples:
    encrypt('012345', 0) => '012345'
    encrypt('012345', 1) 03,14,25 => '251403'
    encrypt('012345', 2)          24,50,13 => '135024'

Decode is encode!!

I realised that if you encode individual items to their output positions then the decode is always just the reverse mapping of encoded to starting state.

The above is easiest done by first encoding a list of integers then doing reverse indexing.

I decided to extend my knowledge of typing too. (Although I want to see if it helps with readability rather than running Mypy).

from typing import overload

@overload
def encrypt5(s: str, n: int) -> str: ...  # normal
@overload
def encrypt5(s: list[int], n: int) -> list[int]: ...  # internal, for decryption
def encrypt5(s, n) -> str:
    if not s or n <= 0:
        return s
    for _ in range(n):
        s = s[2::3] + s[1::3] + s[0::3]
    return s

def decrypt5(s: str, n: int) -> str:
    """
    In which I realise that for this _class_ of problem, decryption is reversing the encrpt
    """
    base = list(range(len(s)))
    enc = encrypt5(base, n)
    if not enc:
        return ''

    return ''.join(s[enc.index(i)] for i in base)

Generalised groupings

I realized that these encodings might be described in terms of cards and dealings:

You have a pack of S cards and dealt them into K piles, left to right, starting again at the leftmost pile when you last dealt at the rightmost until you dealt the last card.
Number the piles, left-to-right from zero to K-1,
You also have a permutation M of the integers 0 to K-1 and you concatenate the dealt piles in the order given in M to form the first reordering of S.
Repeat the process N times to form the final ordering of S

The picture at the beginning is for K == 3 and using Ace as 1, starting from Ace,2,3,...

The following picture is the first concatenation arrangement of M = 2, 1, 0:

My more formal definition of the generalised case:

Task Encrypt grouping chars by mod M, N times:
    Given an input string S; Some order of the K integers 0..K-1 called
    M, where K <= length_of(S); and an integer N > 0
    
    First form groups:
        Form group G[0] by concatenating every K'th item of S starting from index 0
        Form group G[1] by concatenating every K'th item of S starting from index 1
        ...
        Form group G[K-1] by concatenating every K'th item of S starting from index K-1
        
    Group Concatenation:
        Concatenate all the groups G in the order given M to form the partial result P
        P(1) = concatenate_left_to_right(G[i] for i in M)
    Repetition:
        Substituting P for S, repeat the above process N times to form the final
        encrypted string P(N)
    
    Corner cases:
        If S is empty or N is 0 then return S.
        If M is empty return the empty version of sequence S (empty string).

Examples:
    encrypt6('012345', [1, 0], 0) => '012345'
    encrypt6('012345', [1, 0], 1)          => '135024'
    encrypt6('012345', [1, 0], 2)                   => '304152'
    encrypt6('012345', [1, 0], 3)                            => '012345'
    encrypt6('012345', [1, 0], 4)                                     => '135024'

    encrypt6('01234567', [2, 1, 0], 0) => '01234567'
    encrypt6('01234567', [2, 1, 0], 1)          => '25147036'
    encrypt6('01234567', [2, 1, 0], 2)                   => '10576243'
    encrypt6('01234567', [2, 1, 0], 3)                            => '52063174'
    encrypt6('01234567', [2, 1, 0], 4)                                     => '01234567'

Past encryptions in its terms

Note that:
    encrypt(S, N), encrypt2(S, N) do encrypt6(S, [1, 0], N)
    encrypt3(S, N), encrypt4(S, N) do encrypt6(S, [0, 1], N)
    and
    encrypt5(S, N), does encrypt6(S, [2, 1, 0], N)

Generalised Code

from typing import Sequence

def encrypt6(s: Sequence,               # Supporting indexing and __add__ 
            m: Sequence[int],           # perm of ints 0..k-1, k <= len(s)
            n: int
            ) -> Sequence: 
    if not s or n == 0:
        return s
    if not m:
        return type(s)()  # Empty sequence
    assert n >= 0
    assert (k:=len(m)) <= len(s)
    assert set(range(k)) == set(m), \
               f"List {m = } should contain a permutation of all the numbers 0..{k-1} inclusive."
    
    s_type = type(s)  #[0])
    if issubclass(s_type, str):
        for _ in range(n):
            s = ''.join(s[offset::k] for offset in m)
    else:  # non-strings
        empty_val = s_type()
        for _ in range(n):
            s = sum((s[offset::k] for offset in m),
                    start=empty_val)
  
    return s

def decrypt6(s: Sequence,               # Supporting indexing and __add__ 
            m: Sequence[int],           # perm of ints 0..k-1, k <= len(s)
            n: int
            ) -> Sequence:
    "By a reverse of the encoding"
    
    # Where the encoding maps ordered items
    base = list(range(len(s)))
    enc = encrypt6(base, m, n)
    
    if not enc:
        return enc
    
    # Return the revered mapping
    s_type = type(s)
    if issubclass(s_type, str):
        return ''.join(s[enc.index(i)] for i in base)
    else:
        empty_val = s_type()
        return    sum((s[enc.index(i)] for i in base),
                      start=empty_val)

END BIT.

I have looked into the number sequences produced by encrypt6, but this is already a long post so I'll leave that 'till next time

OEISify

2023-08-19T10:04:00.000+01:00

Best read on larger than a phone screen

The title is a rather mean pun on ossify, but it stuck.

I had been creating series of integers from prior work, and had been looking them up, manually, on

The On-Line Encyclopedia of Integer Sequences, OEIS.

Work progressed and I was creating a table of related sequences and looked to automate the OEIS searches.

I found a description of the textual internal format here which seemed steraight-forward and "awkable" - meaning that the format is written in a way that a short awk program could parse. I did find a couple of old python libraries that said they did this, but I wanted to flex my parsing skills and this seemed an enticing.

Mile high looking down

Checking that internal format, the example result seem to be of the form:

%<TAG> A<number> <value>

To show search results in the internal format you need to make a slight change to a url:

https://oeis.org/search?fmt=text&q=3%2C11%2C26%2C50%2C85 (test)

Reading the internal format in more detail

I saw that the %S %T and %U lines are all giving members of the sequence - there seems to be no reason to keep a distinction, I could parse and return all the values in order.

The %p, %t, %o lines when you concatenate all of the same type, give program sources in three different languages. They can be treated similarly.

the %k keywords are best returned as a list of strings.

The other keywords, if they appear more than once in a sequence definition then subsequent values can append as new lines, accumulating a final value.

Extra Info given

When you look at a search result in text format you get a Search: and Showing line I might want to return this data and maybe also derive a Next entry from the details of Showing that is what is needed to generate the next page of results for the given search.

The Regexp

Yep, I started with a regexp, so put a lot of sample text into regex101 as TEST STRING then generated the regexp used in the program.

The Output

A Python dict with all string keys. Each found A-series will be a nested dict of keys corresponding to the %-letter tags.

Keys Search, Showing and Next are about the search and so appear in the top dict.

The Code: oeis_request.py

#!/bin/env python3
# -*- coding: utf-8 -*-
"""

Module to search and return sequences from:
    The On-Line Encyclopedia of Integer Sequences® (OEIS®)
    https://oeis.org/
    https://oeis.org/eishelp1.html

Standalone Usage:
    oeis_request.py <query> <start>
    oeis_request.py --help
    oeis_request.py --test


Created on Sat Aug 12 07:57:53 2023

@author: paddy3118, paddy3118.blogspot.com

"""

from pprint import pp
import re
import ssl
import sys
from typing import Any
import urllib.request
import urllib.parse
import json  # pylint: disable=W0611


URL = 'https://oeis.org/search'

DBG = not True

finditer = re.compile(r"""
    # https://regex101.com/r/iNL0ve/2

    (?: Search: \s+ (?P<SEARCH> .* ))
    | (?: (?P<SHOWING> ^Showing\s.*))
    | (?: ^% (?P<_TAG>  (?P<SEQ>  [STU] ) |  (?P<PROG>  [pto] ) |  (?P<REST> [A-Za-z]) ) \s+
             (?P<ANUMBER> A\d+) \s+ (?P<VALUE> .*) )
    """, re.MULTILINE | re.VERBOSE).finditer

OEISTYPE = dict[str, str | dict[str, str | list[int | str]]]

def oeis_request(query: str, start: int=0) -> OEISTYPE:
    """
    Search OEIS for query returning a page of results from position start as a Python dict


    Parameters
    ----------
    query : str
        The OEIS query in a format consistent with https://oeis.org/hints.html.
    start : int, optional
        page offset when multiple pages of results are available. The default is 0.

    Returns
    -------
    dict[str,                   # key0, top level keys: SEARCH, SHOWING, NEXT
         str                    # string value,
         | dict[                # Or a nested dict for each sequence definition:
                str,                # key1, sequence level string keys,
                str |               # And either a string value,
                list[int            # A list for the int Sequence,
                     | str]]]       # or a list of the str Keywords; as value.

        Returned Python datastructure representing the (page of) sequences found.
        Top level keys are the 'SEARCH' and what results are 'SHOWING', together
        with the sequence numbers of returned sequences like 'A123456' - which
        all have one sub-dict as its value. Top level key 'NEXT', if present,
        has as its value the query and start value pre-computed to retrieve the
        next page of results from oeis_request.

        All possible keys of a sequence sub-dict are the values of global variable
        RECORD_MAP
    """

    # https://stackoverflow.com/a/60671292/10562
    # Security risk as ssl cert checking is skipped for this function
    # pylint: disable=protected-access
    tmp, ssl._create_default_https_context = ssl._create_default_https_context, \
                                             ssl._create_unverified_context
    # pylint: enable=protected-access

    data = {}
    data['fmt'] = 'text'
    data['q'] = query
    if start:
        data['start'] = start
    url_values = urllib.parse.urlencode(data)
    url = URL + '?' + url_values
    if DBG:
        sys.stderr.write(url+'\n')
        sys.stderr.flush()
    with urllib.request.urlopen(url) as req:
        txt =req.read().decode('utf8')

    # Restored:
    ssl._create_default_https_context = tmp    # pylint: disable=protected-access

    return oeis_parser(txt)

def oeis_parser(text: str) -> dict[str, Any]:  # pylint: disable=too-many-branches
    """
    Parse text formatted like https://oeis.org/eishelp1.html

    Parameters
    ----------
    text : str
        text returned from OEIS search.

    Returns
    -------
    dict[str, Any]
        parsed data

    """
    data = {}

    for matchobj in finditer(text):
        group = matchobj.groupdict()
        try:
            anumber = group['ANUMBER']
            if anumber:  # Skip None
                adict = data.setdefault(anumber, {})
        except KeyError:
            anumber = ''
        try:
            value = group['VALUE']
        except KeyError:
            value = ''

        if val:=group[key0:='SEARCH']:
            data[key0] = val
        elif val:=group[key0:='SHOWING']:
            data[key0] = val

        elif val:=group[key1:='SEQ']:
            # Concat all %S $T and %U lines as list[str]
            adict.setdefault(key1, []).extend(int(v) for v in
                                             value.strip().strip(',').split(','))
        elif val:=group['PROG']:
            # Concat all individual %p %t %o lines as list[str]
            adict.setdefault(RECORD_MAP[val], []).append(value.strip())
        elif (val:=group['REST']) and val == 'K':
            # Split and concat Keywords, %K lines as list[str]
            adict.setdefault(RECORD_MAP[val], []).extend(value.strip().split(','))
        elif val:=group['REST']:
            # Concat other %_ lines as list[str]
            adict.setdefault(RECORD_MAP[val], []).append(value.strip())
        else:
            # Should never arrive here from the regexp!
            assert False, f"Got {group = }"

    # fixup
    for key, subd in data.items():
        if key[0] == 'A' and issubclass(type(subd), dict):
            for record, value in subd.items():
                if record != 'Keywords' and issubclass(type(value), list) \
                   and issubclass(type(value[0]), str):
                    subd[record] = '\n'.join(value)

    next_from_showing(data)

    return data

def next_from_showing(parsed: OEISTYPE)-> None:
    """
    Parse "SHOWING" key to generate "NEXT" value that would request the next page

    If already at the end then no next key is in data
    parsed is altered in-place

    Parameters
    ----------
    parsed : OEISTYPE
        Parsed data to insert a NEXT key/value.

    Returns
    -------
    None.
        argument updated in-place

    """
    matchobj = re.match(r"^Showing\s+(\d+)-(\d+)\s+of\s+(\d+)\s*$", parsed.get('SHOWING', ''))
    if matchobj:
        _start, stop, end = [int(x) for x in matchobj.groups()]
        if stop < end:
            search = parsed['SEARCH'].split(':', 1)[-1]
            parsed['NEXT'] = [search, stop]
            return

    if 'NEXT' in parsed:
        del parsed['NEXT']


record_info = """
%I A000001 Identification line (required)
%S A000001 First line of sequence (required)
%T A000001 2nd line of sequence.
%U A000001 3rd line of sequence.
%N A000001 Name (required)
%D A000001 Reference Detailed reference line.
%H A000001 Link to other site.
%F A000001 Formula .
%Y A000001 Cross-references to other sequences.
%A A000001 Author (required)
%O A000001 Offset (required)
%E A000001 Etc Extensions, errors, etc.
%e A000001 Examples examples to illustrate initial terms.
%p A000001 Maple program.
%t A000001 Mathematica program.
%o A000001 OtherProgram in another language.
%K A000001 Keywords (required)
%C A000001 Comments.
""".strip().splitlines()

RECORD_MAP = {words[0][1]: words[2]
              for line in record_info
              for words in [line.split()]}
# pp(RECORD_MAP, sort_dicts=False)

HELP = f"""\
{__doc__}

function oeis_request
=====================

{oeis_request.__doc__}
"""

if __name__ == '__main__':
    if '--help' in sys.argv:
        print(HELP)
        sys.exit(0)

    if '--test' in sys.argv[1:]:
        _req, _start = '1,2,3,4,5,6,6,7,7,8', 0
    elif len(sys.argv) == 3:
        _req, _start = sys.argv[1:]  # pylint: disable=W0632
    else:
        print(HELP)
        sys.exit(1)

    _data = oeis_request(_req, _start)  # pylint: disable=E0601
    pp(_data, width=512, sort_dicts=False)
    #print(json.dumps(_data, indent=2))

Test output

$ ./oeis_request.py 3,11,26,50,85 0 
{'SEARCH': 'seq:3,11,26,50,85',
 'SHOWING': 'Showing 1-1 of 1',
 'A051925': {'Identification': '#84 Jun 26 2022 03:06:23',
             'SEQ': [0, 0, 3, 11, 26, 50, 85, 133, 196, 276, 375, 495, 638, 806, 1001, 1225, 1480, 1768, 2091, 2451, 2850, 3290, 3773, 4301, 4876, 5500, 6175, 6903, 7686, 8526, 9425, 10385, 11408, 12496, 13651, 14875, 16170, 17538, 18981, 20501, 22100, 23780],
             'Name': 'a(n) = n*(2*n+5)*(n-1)/6.',
             'Comments.': 'Related to variance of number of inversions of a random permutation of n letters.\n'
                          'Zero followed by partial sums of A005563. - _Klaus Brockhaus_, Oct 17 2008\n'
                          'a(n)/12 is the variance of the number of inversions of a random permutation of n letters. See evidence in Mathematica code below. - _Geoffrey Critzer_, May 15 2010\n'
                          'The sequence is related to A033487 by A033487(n-1) = n*a(n) - Sum_{i=0..n-1} a(i) = n*(n+1)*(n+2)*(n+3)/4. - _Bruno Berselli_, Apr 04 2012\n'
                          "Deleting the two 0's leaves row 2 of the convolution array A213750. - _Clark Kimberling_, Jun 20 2012\n"
                          'For n>=4, a(n-2) is the number of permutations of 1,2...,n with the distribution of up (1) - down (0) elements 0...0110 (the first n-4 zeros), or, the same, a(n-2) is up-down coefficient {n,6} (see comment in A060351). - _Vladimir Shevelev_, Feb 15 2014',
             'Reference': 'V. N. Sachkov, Probabilistic Methods in Combinatorial Analysis, Cambridge, 1997.',
             'Link': 'Vincenzo Librandi, <a href="/A051925/b051925.txt">Table of n, a(n) for n = 0..1000</a>\nJ. Wang and H. Li, <a href="http://dx.doi.org/10.1016/S0012-365X(01)00301-6">The upper bound of essential chromatic numbers of hypergraphs</a>, Discr. Math. 254 (2002), 555-564.\n<a href="/index/Rec#order_04">Index entries for linear recurrences with constant coefficients</a>, signature (4,-6,4,-1).',
             'Formula': 'a(n) = A000330(n) - n. - _Andrey Kostenko_, Nov 30 2008\nG.f.: x^2*(3-x)/(1-x)^4. - _Colin Barker_, Apr 04 2012\na(n) = 4*a(n-1) - 6*a(n-2) + 4*a(n-3) - a(n-4). - _Vincenzo Librandi_, Apr 27 2012\nE.g.f.: (x^2/6)*(2*x + 9)*exp(x). - _G. C. Greubel_, Jul 19 2017',
             'Mathematica': 'f[{x_, y_}] := 2 y - x^2; Table[f[Coefficient[ Series[Product[Sum[Exp[i t], {i, 0, m}], {m, 1, n - 1}]/n!, {t, 0, 2}], t, {1, 2}]], {n, 0, 41}]*12 (* _Geoffrey Critzer_, May 15 2010 *)\nCoefficientList[Series[x^2*(3-x)/(1-x)^4,{x,0,50}],x] (* _Vincenzo Librandi_, Apr 27 2012 *)',
             'OtherProgram': '(PARI) {print1(a=0, ","); for(n=0, 42, print1(a=a+(n+1)^2-1, ","))} \\\\ _Klaus Brockhaus_, Oct 17 2008\n(Magma) I:=[0, 0, 3, 11]; [n le 4 select I[n] else 4*Self(n-1)-6*Self(n-2)+4*Self(n-3)-Self(n-4): n in [1..50]]; // _Vincenzo Librandi_, Apr 27 2012',
             'Cross-references': 'Cf. A000330, A005563, A033487.',
             'Keywords': ['nonn', 'easy'],
             'Offset': '0,3',
             'Author': '_N. J. A. Sloane_, Dec 19 1999'}}
$

A test with paged results

Highlighting the NEXT value of what is needed to get the next page of results

$ ./oeis_request.py 1,2,3,4,5 0 
{'SEARCH': 'seq:1,2,3,4,5',
 'SHOWING': 'Showing 1-10 of 7513',
 'A000027': {'Identification': 'M0472 N0173 #637 Aug 14 2023 15:10:40',
             ...
 'A007953': {'Identification': '#280 Jun 18 2023 11:41:19',
             ...
 'A000961': 
 ...
 'A002260': {'Identification': '#205 Feb 03 2023 18:43:52',
             ...
 'NEXT': ['1,2,3,4,5', 10]}
$