barak/suffix-primes.hs

## suffix-primes.hs
{-# LANGUAGE UnicodeSyntax #-}
{-# LANGUAGE ParallelListComp #-}
{-# LANGUAGE BangPatterns #-}

-- A suffix prime is a prime number all of whose suffixes are also
-- prime.  No zero digits allowed.  Here we find all, in a given base.

import Prelude.Unicode

import Data.List (intersperse, foldl')
import Data.Char (isDigit)
import Data.Tuple.HT (mapPair)                    -- cabal install utility-ht
import Math.NumberTheory.Primes (isPrime, primes) -- cabal install arithmoi
import System.Environment (getArgs, getProgName)
import System.Exit (exitFailure)

-- Try to raise the abstraction bar a bit here
-- import Control.Monad.SearchTree -- cabal install tree-monad
-- import Control.Parallel.TreeSearch -- cabal install parallel-tree-search

-- Depth-First Search should use less space when seeking the largest.
-- Note that the time should be about the same, since it is dominated
-- by isPrime which must be called on the same set of numbers, albeit
-- in a different order.

-- Breadth-First Search yields numerically sorted list
suffixPrimes ∷ Integer → [Integer]
suffixPrimes b = concat $ takeWhile ((¬) ∘ null) spl
  where spl = [filter isPrime [a+p | a ← map (⋅b^n) [1..b-1], p ← ps]
              | ps ← [0]:spl | n ← [(0∷Int)..]]

-- Depth-First Search yields alphabetically (by reversed base-b
-- digits) sorted list.
suffixPrimes' ∷ Integer → [Integer]
suffixPrimes' b = loop (0∷Int) 0
  where loop len !start
          = concat [p:loop (len+1) p
                   | d←[1..b-1], let p = d ⋅ b^len + start, isPrime p]

-- Depth-First Search for unique (non-extendable) suffix primes.
uniqSuffixPrimes' ∷ Integer → [Integer]
uniqSuffixPrimes' b = loop (0∷Int) 0
  where loop len !start
          = concat [if null ps then [p] else ps
                   | d←[1..b-1],
                     let p = d ⋅ b^len + start,
                     isPrime p,
                     let ps=loop (len+1) p]

-- Extract the length and last element of a long list in a fashion
-- allowing the storage for the beginning of the list to be reclaimed
-- as the computation proceeds.  Unlike
--  (length xs, last xs)
-- which has poor space properties.
lenLast ∷ [a] → (Int, Maybe a)
lenLast = foldl' (\(n,_) x → (n+1, Just x)) (0, Nothing)

lenMax ∷ Ord a ⇒ [a] → (Int, Maybe a)
lenMax = foldl' (flip p1) (0, Nothing)
  where
    p1 x = mapPair (succ, return ∘ maybe x (max x))

-- Breadth-First Search:
suffixPrimesCountLargest ∷ Integer → (Int, Integer)
suffixPrimesCountLargest b = (count, largest)
  where (count, Just largest) = lenLast (suffixPrimes b)

-- Depth-First Search:
suffixPrimesCountLargest' ∷ Integer → (Int, Integer)
suffixPrimesCountLargest' b = (count, largest)
  where (count, Just largest) = lenMax (suffixPrimes' b)

-- Depth-first Search for count, unique count, and largest.
suffixPrimesCountUniqueLargest ∷ Integer → (Int, Int, Integer)
suffixPrimesCountUniqueLargest !b = explore (0∷Int) 0
  where explore !len !start =
          let (count0, ucount0, largest0) =
                foldl' (\(!count1, !ucount1, !largest1) (!count2, !ucount2, !largest2)
                        → (count1+count2, ucount1+ucount2, largest1 `max` largest2))
                (0,0,0)
                [explore (len+1) p | d←[1..b-1], let p = d ⋅ b^len + start, isPrime p]
          in
           if start ≡ 0
           then (count0, ucount0, largest0)
           else if count0 ≡ 0
                then (1, 1, start)
                else (count0+1, ucount0, largest0)

usage ∷ IO ()
usage = do
  progName ← getProgName
  putStrLn $ "Usage: " ⧺ progName ⧺ " [--all|--prime|base ...]"

main ∷ IO ()
main = do
  args ← getArgs
  case args of
    [] → main' [10]
    ["--all"] → main' [3..]
    ["--prime"] → main' (drop 1 primes)
    bss | all ((¬) ∘ null) bss && all (all isDigit) bss →
      main' (map read bss)
    _ →
      do usage
         exitFailure

main' ∷ [Integer] → IO ()
main' bs = do
  putStrLn "Suffix Primes"
  putStrLn "Base\tNumber\tUnique\tLargest"
  mapM_ (\b → let (!count, !ucount, !largest) = suffixPrimesCountUniqueLargest b
               in putStrLn $ concat $ intersperse "\t"
                  [show b, show count, show ucount, show largest])
    bs

{-
Efficiency Ideas

Keep track of the value mod k, for one or more k, then prune before
even constructing d ⋅ b^len + start when the result would be zero mod a
factor of k.  Do this for either some set of prime k, or some set of
products of primes.

d ⋅ b^len + start (mod k)
 = d ⋅ b ⋅ b^(len-1) + start (mod k)
 = d (mod k) ⋅ b (mod k) ⋅ b^(len-1) (mod k) + start (mod k)

-}

-- Notes and pointers.

-- The On-Line Encyclopedia of Integer Sequences
-- A103443, Largest left-truncatable prime in base n
-- https://oeis.org/A103443

-- http://www.johndcook.com/blog/2013/03/12/a-suffix-prime/

-- http://rosettacode.org/wiki/Find_largest_left_truncatable_prime_in_a_given_base

-- Comment from xayahrainie4793
--
-- These primes are called "left-truncatable primes", see
-- https://primes.utm.edu/glossary/page.php?sort=LeftTruncatablePrime
-- and https://oeis.org/A024785
--
-- I have a page for more data (in many bases) for left-truncatable
-- primes, right-truncatable primes, and minimal primes:
-- https://github.com/xayahrainie4793/minimal-primes-and-left-right-truncatable-primes

-- This code
-- https://gist.github.com/barak/5157964

{-
$ ./suffix-primes --prime
Base	Number	Largest
3	3	23
5	15	7817
7	22	817337
11	75	2276005673
13	100	812751503
17	362	13563641583101
19	685	546207129080421139
23	1372	116516557991412919458949
29	3926	4300289072819254369986567661
31	11314	645157007060845985903112107793191
37	31898	1227906370578703713220305155209183002341
41	62930	30936779512712351021129939907890870072806919
43	132356	103002723043391997634166979685453834058870099
47	196709	18846885009617750406310165885751950752131092800539
53	607427	1927313691687008616062631595723184712591552207338172661
59	1870270	4077508027639490147990095553206762226024137866036348886452326969
61	3680378	238089544539042409591603000630200657940967296327557776119980057913
-}

{-
Suffix Primes
Base	Number		Unique		Largest
3	3		1		23
4	16		5		4091
5	15		4		7817
6	454		148		4836525320399
7	22		7		817337
8	446		144		14005650767869
9	108		37		1676456897
10	4260		1442		357686312646216567629137
11	75		26		2276005673
12	170053		60026		13092430647736190817303130065827539
13	100		33		812751503
14	34393		12166		615419590422100474355767356763
15	9357		3266		34068645705927662447286191
16	27982		9925		1088303707153521644968345559987
17	362		138		13563641583101
18	14979714	5367877		571933398724668544269594979167602382822769202133808087
19	685		245		546207129080421139
20	3062899		1096991		1073289911449776273800623217566610940096241078373
21	59131		21014		391461911766647707547123429659688417
22	1599447		573339		33389741556593821170176571348673618833349516314271
23	1372		489		116516557991412919458949
--
25	10484		3746		8211352191239976819943978913
26	7028048		2531840		12399758424125504545829668298375903748028704243943878467
27	98336		35345		10681632250257028944950166363832301357693
28	69058060	24921547	720639908748666454129676051084863753107043032053999738835994276213
29	3926		1407		4300289072819254369986567661
--
31	11314		4075		645157007060845985903112107793191
32	35007483	12664826	1131569863270120248974817136287838489359936416046975582122661310411
33	2527304		910455		924039815258046855588818237912726885772934968646554431
34	240423316	86987571	982498935397824993800810311994840611581693708091339679644860318739434026149
35	607905		219133		448739985000415097566502155600731235704288431019152509
--
37	31898		11561		1227906370578703713220305155209183002341
--
39	12027247	4345818		63286636142126528910450194211546021902700646280145670068655463971
--
41	62930		22783		30936779512712351021129939907890870072806919
--
43	132356		47914		103002723043391997634166979685453834058870099
--
47	196709		71325		18846885009617750406310165885751950752131092800539
--
49	1902815		690111		123160384294554520419522161162422047328327908437480046644671
--
67	10769745	3919053		133287693324947940231446085140961171990912987629821695909722950802440523
--
-}
	{-# LANGUAGE UnicodeSyntax #-}
	{-# LANGUAGE ParallelListComp #-}
	{-# LANGUAGE BangPatterns #-}

	-- A suffix prime is a prime number all of whose suffixes are also
	-- prime. No zero digits allowed. Here we find all, in a given base.

	import Prelude.Unicode

	import Data.List (intersperse, foldl')
	import Data.Char (isDigit)
	import Data.Tuple.HT (mapPair) -- cabal install utility-ht
	import Math.NumberTheory.Primes (isPrime, primes) -- cabal install arithmoi
	import System.Environment (getArgs, getProgName)
	import System.Exit (exitFailure)

	-- Try to raise the abstraction bar a bit here
	-- import Control.Monad.SearchTree -- cabal install tree-monad
	-- import Control.Parallel.TreeSearch -- cabal install parallel-tree-search

	-- Depth-First Search should use less space when seeking the largest.
	-- Note that the time should be about the same, since it is dominated
	-- by isPrime which must be called on the same set of numbers, albeit
	-- in a different order.

	-- Breadth-First Search yields numerically sorted list
	suffixPrimes ∷ Integer → [Integer]
	suffixPrimes b = concat $ takeWhile ((¬) ∘ null) spl
	where spl = [filter isPrime [a+p \| a ← map (⋅b^n) [1..b-1], p ← ps]
	\| ps ← [0]:spl \| n ← [(0∷Int)..]]

	-- Depth-First Search yields alphabetically (by reversed base-b
	-- digits) sorted list.
	suffixPrimes' ∷ Integer → [Integer]
	suffixPrimes' b = loop (0∷Int) 0
	where loop len !start
	= concat [p:loop (len+1) p
	\| d←[1..b-1], let p = d ⋅ b^len + start, isPrime p]

	-- Depth-First Search for unique (non-extendable) suffix primes.
	uniqSuffixPrimes' ∷ Integer → [Integer]
	uniqSuffixPrimes' b = loop (0∷Int) 0
	where loop len !start
	= concat [if null ps then [p] else ps
	\| d←[1..b-1],
	let p = d ⋅ b^len + start,
	isPrime p,
	let ps=loop (len+1) p]

	-- Extract the length and last element of a long list in a fashion
	-- allowing the storage for the beginning of the list to be reclaimed
	-- as the computation proceeds. Unlike
	-- (length xs, last xs)
	-- which has poor space properties.
	lenLast ∷ [a] → (Int, Maybe a)
	lenLast = foldl' (\(n,_) x → (n+1, Just x)) (0, Nothing)

	lenMax ∷ Ord a ⇒ [a] → (Int, Maybe a)
	lenMax = foldl' (flip p1) (0, Nothing)
	where
	p1 x = mapPair (succ, return ∘ maybe x (max x))

	-- Breadth-First Search:
	suffixPrimesCountLargest ∷ Integer → (Int, Integer)
	suffixPrimesCountLargest b = (count, largest)
	where (count, Just largest) = lenLast (suffixPrimes b)

	-- Depth-First Search:
	suffixPrimesCountLargest' ∷ Integer → (Int, Integer)
	suffixPrimesCountLargest' b = (count, largest)
	where (count, Just largest) = lenMax (suffixPrimes' b)

	-- Depth-first Search for count, unique count, and largest.
	suffixPrimesCountUniqueLargest ∷ Integer → (Int, Int, Integer)
	suffixPrimesCountUniqueLargest !b = explore (0∷Int) 0
	where explore !len !start =
	let (count0, ucount0, largest0) =
	foldl' (\(!count1, !ucount1, !largest1) (!count2, !ucount2, !largest2)
	→ (count1+count2, ucount1+ucount2, largest1 `max` largest2))
	(0,0,0)
	[explore (len+1) p \| d←[1..b-1], let p = d ⋅ b^len + start, isPrime p]
	in
	if start ≡ 0
	then (count0, ucount0, largest0)
	else if count0 ≡ 0
	then (1, 1, start)
	else (count0+1, ucount0, largest0)

	usage ∷ IO ()
	usage = do
	progName ← getProgName
	putStrLn $ "Usage: " ⧺ progName ⧺ " [--all\|--prime\|base ...]"

	main ∷ IO ()
	main = do
	args ← getArgs
	case args of
	[] → main' [10]
	["--all"] → main' [3..]
	["--prime"] → main' (drop 1 primes)
	bss \| all ((¬) ∘ null) bss && all (all isDigit) bss →
	main' (map read bss)
	_ →
	do usage
	exitFailure

	main' ∷ [Integer] → IO ()
	main' bs = do
	putStrLn "Suffix Primes"
	putStrLn "Base\tNumber\tUnique\tLargest"
	mapM_ (\b → let (!count, !ucount, !largest) = suffixPrimesCountUniqueLargest b
	in putStrLn $ concat $ intersperse "\t"
	[show b, show count, show ucount, show largest])
	bs

	{-
	Efficiency Ideas

	Keep track of the value mod k, for one or more k, then prune before
	even constructing d ⋅ b^len + start when the result would be zero mod a
	factor of k. Do this for either some set of prime k, or some set of
	products of primes.

	d ⋅ b^len + start (mod k)
	= d ⋅ b ⋅ b^(len-1) + start (mod k)
	= d (mod k) ⋅ b (mod k) ⋅ b^(len-1) (mod k) + start (mod k)

	-}

	-- Notes and pointers.

	-- The On-Line Encyclopedia of Integer Sequences
	-- A103443, Largest left-truncatable prime in base n
	-- https://oeis.org/A103443

	-- http://www.johndcook.com/blog/2013/03/12/a-suffix-prime/

	-- http://rosettacode.org/wiki/Find_largest_left_truncatable_prime_in_a_given_base

	-- Comment from xayahrainie4793
	--
	-- These primes are called "left-truncatable primes", see
	-- https://primes.utm.edu/glossary/page.php?sort=LeftTruncatablePrime
	-- and https://oeis.org/A024785
	--
	-- I have a page for more data (in many bases) for left-truncatable
	-- primes, right-truncatable primes, and minimal primes:
	-- https://github.com/xayahrainie4793/minimal-primes-and-left-right-truncatable-primes

	-- This code
	-- https://gist.github.com/barak/5157964

	{-
	$ ./suffix-primes --prime
	Base Number Largest
	3 3 23
	5 15 7817
	7 22 817337
	11 75 2276005673
	13 100 812751503
	17 362 13563641583101
	19 685 546207129080421139
	23 1372 116516557991412919458949
	29 3926 4300289072819254369986567661
	31 11314 645157007060845985903112107793191
	37 31898 1227906370578703713220305155209183002341
	41 62930 30936779512712351021129939907890870072806919
	43 132356 103002723043391997634166979685453834058870099
	47 196709 18846885009617750406310165885751950752131092800539
	53 607427 1927313691687008616062631595723184712591552207338172661
	59 1870270 4077508027639490147990095553206762226024137866036348886452326969
	61 3680378 238089544539042409591603000630200657940967296327557776119980057913
	-}

	{-
	Suffix Primes
	Base Number Unique Largest
	3 3 1 23
	4 16 5 4091
	5 15 4 7817
	6 454 148 4836525320399
	7 22 7 817337
	8 446 144 14005650767869
	9 108 37 1676456897
	10 4260 1442 357686312646216567629137
	11 75 26 2276005673
	12 170053 60026 13092430647736190817303130065827539
	13 100 33 812751503
	14 34393 12166 615419590422100474355767356763
	15 9357 3266 34068645705927662447286191
	16 27982 9925 1088303707153521644968345559987
	17 362 138 13563641583101
	18 14979714 5367877 571933398724668544269594979167602382822769202133808087
	19 685 245 546207129080421139
	20 3062899 1096991 1073289911449776273800623217566610940096241078373
	21 59131 21014 391461911766647707547123429659688417
	22 1599447 573339 33389741556593821170176571348673618833349516314271
	23 1372 489 116516557991412919458949
	--
	25 10484 3746 8211352191239976819943978913
	26 7028048 2531840 12399758424125504545829668298375903748028704243943878467
	27 98336 35345 10681632250257028944950166363832301357693
	28 69058060 24921547 720639908748666454129676051084863753107043032053999738835994276213
	29 3926 1407 4300289072819254369986567661
	--
	31 11314 4075 645157007060845985903112107793191
	32 35007483 12664826 1131569863270120248974817136287838489359936416046975582122661310411
	33 2527304 910455 924039815258046855588818237912726885772934968646554431
	34 240423316 86987571 982498935397824993800810311994840611581693708091339679644860318739434026149
	35 607905 219133 448739985000415097566502155600731235704288431019152509
	--
	37 31898 11561 1227906370578703713220305155209183002341
	--
	39 12027247 4345818 63286636142126528910450194211546021902700646280145670068655463971
	--
	41 62930 22783 30936779512712351021129939907890870072806919
	--
	43 132356 47914 103002723043391997634166979685453834058870099
	--
	47 196709 71325 18846885009617750406310165885751950752131092800539
	--
	49 1902815 690111 123160384294554520419522161162422047328327908437480046644671
	--
	67 10769745 3919053 133287693324947940231446085140961171990912987629821695909722950802440523
	--
	-}