The Intellectual Wilderness There is nothing more useless than doing efficiently that which should not be done at all.

2019.11.21 17:24

Testing Textually Composed Numbers for Primality

Filed under: Computing,Science & Tech — Tags: , , , , , — zxq9 @ 17:24

Last night on Twitter one of my favorite accounts, @fermatslibrary, put up an interesting post:

Start at 82 and write it down, then 81, write that down, etc. until you reach one, and you’ve written a (huge) prime number. Wow!

This seemed so strange to me, so of course I got curious:

After some doodling around I wrote a script that checks whether numbers constructed by the method above are prime numbers when the number construction is performed in various numeric bases (from 2 to 36, limited to 36 because for now I’m cheating using Erlang’s integer_to_list/2). It prints the results to the screen as the process handling each base finishes and writes a text file “texty_primes-result-[timestamp].eterms” at the end for convenience so you can use file:consult/1 on it later to mess around.

There are a few strange aspects to large primes, one of them being that checking whether or not they are prime can be a computationally intense task (and nobody knows a shortcut to this). To this end I wrote the Miller-Rabin primality test into the script and allow the caller to decide how many rounds of Miller-Rabin to run against the numbers to check them. So far the numbers that have come out have matched what is expected, but once the numbers get extremely large (and they get pretty big in a hurry!) there is only some degree of confidence that they are really prime, so don’t take the output as gospel.

I wrote the program in Erlang as an escript, so if you want to run it yourself just download the script and execute it.
The script can be found here: texty_primes

A results file containing the (very likely) prime constructions in bases 2 through 36 using “count-back from X” where X is 1 to 500 can be found here: texty_primes-result-20191121171755.eterms
Analyzing from 1 to 500 in bases 2 through 36 took about 25 minutes on a mid-grade 8-core system (Ryzen5). There are some loooooooooong numbers in that file… It would be interesting to test the largest of them for primality in more depth.

(Note that while the script runs you will receive unordered “Base X Result” messages printed to stdout. This is because every base is handed off to a separate process for analysis and they finish at very different times somewhat unpredictably. When all processing is complete the text file will contain a sorted list of {Base, ListOfPrimes} that is easier to browse.)

An interesting phenomenon I observed while doing this is that some numeric bases seem simply unsuited to producing primes when numbers are generated in this manner, bases that themselves are prime numbers in particular. Other bases seem to be rather fruitful places to search for this phenomenon.

Another interesting phenomenon is the wide number of numeric bases in which the numbers “21”, “321”, “4321” and “5321” turn out to be prime. “21” and “4321” in particular turn up quite a lot.

Perhaps most strangely of all is that base 10 is not a very good place to look for these kinds of primes! In fact, the “count back from 82” prime is the only one that can be constructed starting between 1 and 500 that I’ve found. It is remarkable that anyone discovered that at all, and also remarkable that it doesn’t happen to start at 14,562 instead of 82 — I’m sure nobody would have noticed this were any number much higher than 82 the magic starting point for constructing a prime this way.

This was fun! If you have any insights, questions, challenges or improvements, please let me know in the comments.

2017.10.21 14:57

Erlang: Naive Matrix Multiplication

Someone asked what was surely a homework question today on StackOverflow about matrix multiplication in Erlang. I set out to answer him in as simple a way as possible, but wound up writing a naive matrix generation and multiplication module.

The code to the module might be of interest to new Erlangers, as it adheres both to the style of zuuid and includes many examples of using a combination of list operations and explicit recursion to cut clutter and make the meaning of otherwise complex operations clear.

Here is the code:

%%% @doc
%%% A naive matrix generation, rotation and multiplication module.
%%% It doesn't concern itself with much checking, so input dimensions must be known
%%% prior to calling any of these functions lest you receive some weird results back,
%%% as most of these functions do not crash on input that go against the rules of
%%% matrix multiplication.
%%%
%%% All functions crash on obviously bad values.
%%% @end 

-module(naive_matrix).
-export([random/2, random/3, rotate/1, multiply/2]).

-type matrix() :: [[number()]].


-spec random(Size, MaxValue) -> Matrix
    when Size     :: pos_integer(),
         MaxValue :: pos_integer(),
         Matrix   :: matrix().
%% @doc
%% Generate a square matrix of dimensions {Size, Size} populated with random
%% integer values inclusive of 1..MaxValue.

random(Size, MaxValue) when Size > 0, MaxValue > 0 ->
    random(Size, Size, MaxValue).


-spec random(X, Y, MaxValue) -> Matrix
    when X        :: pos_integer(),
         Y        :: pos_integer(),
         MaxValue :: pos_integer(),
         Matrix   :: matrix().
%% @doc
%% Generate a matrix of dimensions {X, Y} populated with random integer values
%% inclusive 1..MaxValue.

random(X, Y, MaxValue) when X > 0, Y > 0, MaxValue > 0 ->
    Columns = lists:duplicate(X, []),
    Populate = fun(Col) -> row(Y, MaxValue, Col) end,
    lists:map(Populate, Columns).


-spec row(Size, MaxValue, Acc) -> NewAcc
    when Size     :: non_neg_integer(),
         MaxValue :: pos_integer(),
         Acc      :: [pos_integer()],
         NewAcc   :: [pos_integer()].
%% @private
%% Generate a single row of random integers.

row(0, _, Acc) ->
    Acc;
row(Size, MaxValue, Acc) ->
    row(Size - 1, MaxValue, [rand:uniform(MaxValue) | Acc]).


-spec rotate(matrix()) -> matrix().
%% @doc
%% Takes a matrix of {X, Y} size and rotates it left, returning a matrix of {Y, X} size.

rotate(Matrix) ->
    rotate(Matrix, [], [], []).


-spec rotate(Matrix, Rem, Current, Acc) -> Rotated
    when Matrix  :: matrix(),
         Rem     :: [[number()]],
         Current :: [number()],
         Acc     :: matrix(),
         Rotated :: matrix().
%% @private
%% Iterates doubly over a matrix, packing the diminished remainder into Rem and
%% packing the current row into Current. This is naive, in that it assumes an
%% even matrix of dimentions {X, Y}, and will return one of dimentions {Y, X}
%% based on the length of the first row, regardless whether the input was actually
%% even.

rotate([[] | _], [], [], Acc) ->
    Acc;
rotate([], Rem, Current, Acc) ->
    NewRem = lists:reverse(Rem),
    NewCurrent = lists:reverse(Current),
    rotate(NewRem, [], [], [NewCurrent | Acc]);
rotate([[V | Vs] | Rows], Rem, Current, Acc) ->
    rotate(Rows, [Vs | Rem], [V | Current], Acc).


-spec multiply(ValueA, ValueB) -> Product
    when ValueA  :: number() | matrix(),
         ValueB  :: number() | matrix(),
         Product :: number() | matrix().
%% @doc
%% Accept any legal combination of scalar and matrix values to be multiplied.
%% The correct operation will be chosen based on input values.

multiply(A, B) when is_number(A), is_number(B) ->
    A * B;
multiply(A, B) when is_number(A), is_list(B) ->
    multiply_scalar(A, B);
multiply(A, B) when is_list(A), is_list(B) ->
    multiply_matrix(A, B).


-spec multiply_scalar(A, B) -> Product
    when A       :: number(),
         B       :: matrix(),
         Product :: matrix().
%% @private
%% Simple scalar multiplication of a matrix.

multiply_scalar(A, B) ->
    multiply_scalar(A, B, []).


-spec multiply_scalar(A, B, Acc) -> Product
    when A       :: number(),
         B       :: matrix(),
         Acc     :: matrix(),
         Product :: matrix().
%% @private
%% Scalar multiplication is implemented here as an explicit recursion over
%% a list of lists, each element of which is subjected to a map operation.

multiply_scalar(A, [B | Bs], Acc) ->
    Row = lists:map(fun(N) -> A * N end, B),
    multiply_scalar(A, Bs, [Row | Acc]);
multiply_scalar(_, [], Acc) ->
    lists:reverse(Acc).


-spec multiply_matrix(A, B) -> Product
    when A       :: matrix(),
         B       :: matrix(),
         Product :: matrix().
%% @doc
%% Multiply two matrices together according to the matrix multiplication rules.
%% This function does not check that the inputs are actually proper (regular)
%% matrices, but does check that the input row/column lengths are compatible.

multiply_matrix(A = [R | _], B) when length(R) == length(B) ->
    multiply_matrix(A, rotate(B), []).


-spec multiply_matrix(A, B, Acc) -> Product
    when A       :: matrix(),
         B       :: matrix(),
         Acc     :: matrix(),
         Product :: matrix().
%% @private
%% Iterate a row multiplication operation of each row of A over matrix B until
%% A is exhausted.

multiply_matrix([A | As], B, Acc) ->
    Prod = multiply_row(A, B, []),
    multiply_matrix(As, B, [Prod | Acc]);
multiply_matrix([], _, Acc) ->
    lists:reverse(Acc).


-spec multiply_row(Row, B, Acc) -> Product
    when Row     :: [number()],
         B       :: matrix(),
         Acc     :: [number()],
         Product :: [number()].
%% @private
%% Multiply each row of matrix B by the input Row, returning the list of resulting sums.

multiply_row(Row, [B | Bs], Acc) ->
    ZipProd = lists:zipwith(fun(X, Y) -> X * Y end, Row, B),
    Sum = lists:sum(ZipProd),
    multiply_row(Row, Bs, [Sum | Acc]);
multiply_row(_, [], Acc) ->
    Acc.

Hopefully reading that on a blog won’t drive anyone too nuts. I’ll probably include an expanded version of that (or something related) in a convenience library eventually. Unless I forget. Meh.

Powered by WordPress