Skip to content

Instantly share code, notes, and snippets.

@barak barak/suffix-primes.hs
Last active Dec 14, 2015

Embed
What would you like to do?
suffix primes in haskell
{-# 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 [(0Int)..]]
-- Depth-First Search yields alphabetically (by reversed base-b
-- digits) sorted list.
suffixPrimes' Integer [Integer]
suffixPrimes' b = loop (0Int) 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 (0Int) 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 (0Int) 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
-- 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
--
-}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.