Find largest left truncatable prime in a given base: Difference between revisions

Content added Content deleted

Inline

Revision as of 08:41, 8 October 2012

A truncatable prime is one where all non-empty substrings that finish at the end of the number (right-substrings) are also primes when understood as numbers in a particular base. The largest such prime in a given (integer) base is therefore computable, provided the base is larger than 2.

Let's consider what happens in base 10. Obviously the right most digit must be prime, so in base 10 candidates are 2,3,5,7. Putting a digit in the range 1 to base-1 in front of each candidate must result in a prime. So 2 and 5 like the whale and the carnations in Hitchickers come into existance only to be extinguished before they have time to realize it. 13,17, ... 83,97 are candidates. Again Putting a digit in the range 1 to base-1 in front of each candidate must be a prime. Repeating until there are no larger candidates finds the largest left truncatable prime.

Let's work base 3 by hand:

0 and 1 are not prime so the last digit must be 2. 12₃ = 5₁₀ which is prime, 22₃ = 8₁₀ which is not so 12₃ is the only candidate. 112₃ = 14₁₀ which is not prime, 212₃ = 23₁₀ which is, so 212₃ is the only candidate. 1212₃ = 50₁₀ which is not prime, 2212₃ = 77₁₀ which also is not prime. So there are no more candidates, therefore 23 is the largest left truncatable prime in base 3.

The task is to reconstruct as much, and possibly more, of the table in the OEIS as you are able.

Related Tasks:

Miller-Rabin primality test

C

<lang c>#include <stdio.h>

include <gmp.h>

typedef unsigned long ulong;

ulong small_primes[] = {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};

define MAX_STACK 128

mpz_t tens[MAX_STACK], value[MAX_STACK], answer;

ulong base, seen_depth;

void add_digit(ulong i) { ulong d; for (d = 1; d < base; d++) { mpz_set(value[i], value[i-1]); mpz_addmul_ui(value[i], tens[i], d); if (!mpz_probab_prime_p(value[i], 1)) continue;

if (i > seen_depth || (i == seen_depth && mpz_cmp(value[i], answer) == 1)) { if (!mpz_probab_prime_p(value[i], 50)) continue;

mpz_set(answer, value[i]); seen_depth = i; gmp_fprintf(stderr, "\tb=%lu d=%2lu | %Zd\n", base, i, answer); }

add_digit(i+1); } }

void do_base() { ulong i; mpz_set_ui(answer, 0); mpz_set_ui(tens[0], 1); for (i = 1; i < MAX_STACK; i++) mpz_mul_ui(tens[i], tens[i-1], base);

for (seen_depth = i = 0; small_primes[i] < base; i++) { fprintf(stderr, "\tb=%lu digit %lu\n", base, small_primes[i]); mpz_set_ui(value[0], small_primes[i]); add_digit(1); } gmp_printf("%d: %Zd\n", base, answer); }

int main(void) { ulong i; for (i = 0; i < MAX_STACK; i++) { mpz_init_set_ui(tens[i], 0); mpz_init_set_ui(value[i], 0); } mpz_init_set_ui(answer, 0);

for (base = 22; base < 30; base++) do_base();

return 0; }</lang>

Output:

$ ./a.out 2&>/dev/null
22: 33389741556593821170176571348673618833349516314271
23: 116516557991412919458949
24: 10594160686143126162708955915379656211582267119948391137176997290182218433
25: 8211352191239976819943978913
26: 12399758424125504545829668298375903748028704243943878467
27: 10681632250257028944950166363832301357693
28: 720639908748666454129676051084863753107043032053999738835994276213
29: 4300289072819254369986567661
...

Haskell

Miller-Rabin test code from HaskellWiki, with modifications. <lang haskell>primesTo100 = [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]

-- (eq. to) find2km (2^k * n) = (k,n) find2km :: Integral a => a -> (Int,a) find2km n = f 0 n where f k m

           | r == 1 = (k,m)
           | otherwise = f (k+1) q
           where (q,r) = quotRem m 2

-- n is the number to test; a is the (presumably randomly chosen) witness millerRabinPrimality :: Integer -> Integer -> Bool millerRabinPrimality n a

   | a >= n_ = True
   | b0 == 1 || b0 == n_ = True
   | otherwise = iter (tail b)
   where
       n_ = n-1
       (k,m) = find2km n_
       b0 = powMod n a m
       b = take k $ iterate (squareMod n) b0
       iter [] = False
       iter (x:xs)
           | x == 1 = False
           | x == n_ = True
           | otherwise = iter xs

-- (eq. to) pow_ (*) (^2) n k = n^k pow_ :: (Num a, Integral b) => (a->a->a) -> (a->a) -> a -> b -> a pow_ _ _ _ 0 = 1 pow_ mul sq x_ n_ = f x_ n_ 1

   where
       f x n y
           | n == 1 = x `mul` y
           | r == 0 = f x2 q y
           | otherwise = f x2 q (x `mul` y)
           where
               (q,r) = quotRem n 2
               x2 = sq x

mulMod :: Integral a => a -> a -> a -> a mulMod a b c = (b * c) `mod` a squareMod :: Integral a => a -> a -> a squareMod a b = (b * b) `rem` a

-- (eq. to) powMod m n k = n^k `mod` m powMod :: Integral a => a -> a -> a -> a powMod m = pow_ (mulMod m) (squareMod m)

-- Caller supplies a witness list w, which may be used for MR test. -- Use faster trial division against a small primes list first, to -- weed out more obvious composites. is_prime w n | n < 100 = n `elem` primesTo100 | any ((==0).(n`mod`)) primesTo100 = False | otherwise = all (millerRabinPrimality n) w

-- final result gets a more thorough Miller-Rabin left_trunc base = head $ filter (is_prime primesTo100) (reverse hopeful) where hopeful = extend base $ takeWhile (<base) primesTo100 where extend b x = if null d then x else extend (b*base) d where d = concatMap addDigit [1..base-1] -- we do *one* prime test, which seems good enough in practice addDigit a = filter (is_prime [3]) $ map (a*b+) x

main = mapM_ print $ map (\x->(x, left_trunc x)) [3..21]</lang>

(3,23)
(4,4091)
(5,7817)
(6,4836525320399)
(7,817337)
(8,14005650767869)
(9,1676456897)
(10,357686312646216567629137)
(11,2276005673)
(12,13092430647736190817303130065827539)
(13,812751503)
(14,615419590422100474355767356763)
(15,34068645705927662447286191)
(16,1088303707153521644968345559987)
(17,13563641583101)
(18,571933398724668544269594979167602382822769202133808087)
(19,546207129080421139)
(20,1073289911449776273800623217566610940096241078373)
(21,391461911766647707547123429659688417)

Perl 6

Using the Miller-Rabin definition of is-prime from [1], or the built-in, if available. For speed we do a single try, which allows a smattering of composites in, but this does not appear to damage the algorithm much, and the correct answer appears to be within a few candidates of the end all the time. We keep the last few hundred candidates, sort them into descending order, then pick the largest that passes Miller-Rabin with 100 tries.

I believe the composites don't hurt the algorithm because they do not confer any selective advantage on their "children", so their average length is no longer than the correct solution. In addition, the candidates in the finalist set all appear to be largely independent of each other, so any given composite tends to contribute only a single candidate, if any. I have no proof of this; I prefer my math to have more vigor than rigor. :-) <lang perl6>for 3..* -> $base {

   say "Starting base $base...";
   my @stems = grep { is-prime($_, 100)}, ^$base;
   my @runoff;
   for 1 .. * -> $digits {
     print ' ', @stems.elems;
     my @new;
     my $place = $base ** $digits;
     my $tries = 1;
     for @stems -> $stem {

for 1 ..^ $base -> $digit { my $new = $digit * $place + $stem; @new.push($new) if is-prime($new, $tries);

       }
     }
     last unless @new;
     push @runoff, @stems if @new < @stems and @new < 100;
     @stems = @new;
   }
   push @runoff, @stems;
   say "\n  Finalists: ", +@runoff;

   for @runoff.sort(-*) -> $finalist {

my $b = $finalist.base($base); say " Checking :$base\<", $b, '>'; my $ok = True; for $base ** 2, $base ** 3, $base ** 4 ... $base ** $b.chars -> $place { my $f = $finalist % $place; if not is-prime($f, 100) { say " Oops, :$base\<", $f.base($base), '> is not a prime!!!'; $ok = False; last; } } next unless $ok;

say " Largest ltp in base $base = $finalist"; last;

}</lang>

Output:

Starting base 3...
 1 1 1
  Finalists: 1
  Checking :3<212>
  Largest ltp in base 3 = 23
Starting base 4...
 2 3 5 5 6 4 2
  Finalists: 12
  Checking :4<3323233>
    Oops, :4<33> is not a prime!!!
  Checking :4<1323233>
    Oops, :4<33> is not a prime!!!
  Checking :4<333323>
  Largest ltp in base 4 = 4091
Starting base 5...
 2 4 4 3 1 1
  Finalists: 8
  Checking :5<222232>
  Largest ltp in base 5 = 7817
Starting base 6...
 3 4 12 26 45 56 65 67 66 57 40 25 14 9 4 2 1
  Finalists: 285
  Checking :6<14141511414451435>
  Largest ltp in base 6 = 4836525320399
Starting base 7...
 3 6 6 4 1 1 1
  Finalists: 11
  Checking :7<6642623>
  Largest ltp in base 7 = 817337
Starting base 8...
 4 12 29 50 67 79 65 52 39 28 17 8 3 2 1
  Finalists: 294
  Checking :8<313636165537775>
  Largest ltp in base 8 = 14005650767869
Starting base 9...
 4 9 15 17 24 16 9 6 5 3
  Finalists: 63
  Checking :9<4284484465>
  Largest ltp in base 9 = 1676456897
Starting base 10...
 4 11 40 101 197 335 439 536 556 528 456 358 278 212 117 72 42 24 13 6 5 4 3 1
  Finalists: 287
  Checking :10<357686312646216567629137>
  Largest ltp in base 10 = 357686312646216567629137
Starting base 11...
 4 8 15 18 15 8 4 2 1
  Finalists: 48
  Checking :11<A68822827>
  Largest ltp in base 11 = 2276005673
Starting base 12...
 5 24 124 431 1179 2616 4948 8054 11732 15460 18100 19546 19777 18280 15771 12574 9636 6866 4625 2990 1867 1152 627 357 189 107 50 20 9 3 1 1
  Finalists: 190
  Checking :12<471A34A164259BA16B324AB8A32B7817>
  Largest ltp in base 12 = 13092430647736190817303130065827539
Starting base 13...
 5 13 21 23 17 11 7 4
  Finalists: 62
  Checking :13<CC4C8C65>
  Largest ltp in base 13 = 812751503
Starting base 14...
 6 28 103 308 694 1329 2148 3089 3844 4290 4311 3923 3290 2609 1948 1298 809 529 332 167 97 55 27 13 5 2
  Finalists: 366
  Checking :14<D967CCD63388522619883A7D23>
  Largest ltp in base 14 = 615419590422100474355767356763
Starting base 15...
 6 22 80 213 401 658 955 1208 1307 1209 1075 841 626 434 268 144 75 46 27 14 7 1
  Finalists: 314
  Checking :15<6C6C2CE2CEEEA4826E642B>
  Largest ltp in base 15 = 34068645705927662447286191
Starting base 16...
 6 32 131 355 792 1369 2093 2862 3405 3693 3519 3140 2660 1930 1342 910 571 335 180 95 50 25 9 5 1
  Finalists: 365
  Checking :16<DBC7FBA24FE6AEC462ABF63B3>
  Largest ltp in base 16 = 1088303707153521644968345559987
Starting base 17...
 6 23 47 64 80 62 43 32 23 8 1
  Finalists: 249
  Checking :17<6C66CC4CC83>
  Largest ltp in base 17 = 13563641583101
Starting base 18...
 7 49 311 1396 5117 15243 38818 85684 167133 293518 468456 680171 911723 1133959 1313343 1423597 1449405 1395514 1270222 1097353 902477 707896 529887 381239 263275 174684 111046 67969 40704 23201 12793 6722 3444 1714 859 422 205 98 39 14 7 1 1
  Finalists: 364
  Checking :18<AF93E41A586HE75A7HHAAB7HE12FG79992GA7741B3D>
  Largest ltp in base 18 = 571933398724668544269594979167602382822769202133808087
Starting base 19...
 7 29 61 106 122 117 93 66 36 18 10 10 6 4
  Finalists: 350
  Checking :19<CIEG86GCEA2C6H>
  Largest ltp in base 19 = 546207129080421139
Starting base 20...
 8 56 321 1311 4156 10963 24589 47737 83011 129098 181707 234685 278792 306852 315113 302684 273080 233070 188331 145016 105557 73276 48819 31244 19237 11209 6209 3383 1689 917 430 196 80 44 20 7 2
  Finalists: 349
  Checking :20<FC777G3CG1FIDI9I31IE5FFB379C7A3F6EFID>
  Largest ltp in base 20 = 1073289911449776273800623217566610940096241078373
Starting base 21...
 8 41 165 457 1079 2072 3316 4727 6003 6801 7051 6732 5862 4721 3505 2474 1662 1039 628 369 219 112 52 17 13 4 2
  Finalists: 200
  Checking :21<G8AGG2GCA8CAK4K68GEA4G2K22H>
  Largest ltp in base 21 = 391461911766647707547123429659688417
Starting base 22...
 8 48 261 1035 3259 8247 17727 33303 55355 82380 111919 137859 158048 167447 165789 153003 132516 108086 83974 61950 43453 29212 18493 11352 6693 3738 2053 1125 594 293 145 70 31 13 6 2 1
  Finalists: 268
  Checking :22<FFHALC8JFB9JKA2AH9FAB4I9L5I9L3GF8D5L5>
  Largest ltp in base 22 = 33389741556593821170176571348673618833349516314271
Starting base 23...
 8 30 78 137 181 200 186 171 121 100 67 41 24 16 9 2 1
  Finalists: 260
  Checking :23<IMMGM6C6IMCI66A4H>
  Largest ltp in base 23 = 116516557991412919458949

Ruby

Ruby Ruby

Compute the largest left truncatable prime
Nigel_Galloway
September 15th., 2012.

require 'prime' BASE = 3 MAX = 500 stems = Prime.each(BASE-1).to_a (1..MAX-1).each {|i|

 print "#{stems.length} "
 t = []
 b = BASE ** i
 stems.each {|z|
   (1..BASE-1).each {|n|
     c = n*b+z
     t.push(c) if c.prime?
 }}
 break if t.empty?
 stems = t

} puts "The largest left truncatable prime #{"less than #{BASE ** MAX} " if MAX < 500}in base #{BASE} is #{stems.max}" </lang> By changing BASE from 3 to 14 this produces the solutions in 'Number of left truncatable primes in a given base' on the Discussion Page for bases except 10, 12 and 14.

The maximum left truncatable prime in bases 10 , 12, and 14 are very large. By changing MAX to 6 and BASE to 10 solves related task 1:

The largest left truncatable prime less than 1000000 in base 10 is 998443

JRuby

I require a fast probably prime test. Java has one, is it any good? Let's find out. Ruby can borrow from Java using JRuby. Modifying the Ruby solution: <lang Ruby>

Compute the largest left truncatable prime
Nigel_Galloway
September 15th., 2012.

require 'prime' require 'java' BASE = 10 MAX = 500 stems = Prime.each(BASE-1).to_a (1..MAX-1).each {|i|

 print "#{stems.length} "
 t = []
 b = BASE ** i
 stems.each {|z|
   (1..BASE-1).each {|n|
     c = n*b+z
     t.push(c) if java.math.BigInteger.new(c.to_s).isProbablePrime(100)
 }}
 break if t.empty?
 stems = t

} puts "\nThe largest left truncatable prime #{"less than #{BASE ** MAX} " if MAX < 500}in base #{BASE} is #{stems.max}" </lang> Produces all the reults in 'Number of left truncatable primes in a given base' on the discussion page. For bases 18, 20, and 22 I changed the confidence level from 100 to 5 and checked the final answer. Even so base 18 takes a while. For base 24:

9 87 677 3808 17096 63509 199432 545332 1319708 2863180 Error: Your application
used more memory than the safety cap of 500m.
Specify -J-Xmx####m to increase it (#### = cap size in MB).
Specify -w for full OutOfMemoryError stack trace

That is going to be big!

Tcl

<lang tcl>package require Tcl 8.5

proc tcl::mathfunc::modexp {a b n} {

   for {set c 1} {$b} {set a [expr {$a*$a%$n}]} {
       if {$b & 1} {
           set c [expr {$c*$a%$n}]
       }
       set b [expr {$b >> 1}]
   }
   return $c

}

Based on Miller-Rabin primality testing, but with small prime check first

proc is_prime {n {count 10}} {

   # fast check against small primes
   foreach p {

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

} {

if {$n == $p} {return true} if {$n % $p == 0} {return false}

   # write n-1 as 2^s·d with d odd by factoring powers of 2 from n-1
   set d [expr {$n - 1}]
   for {set s 0} {$d & 1 == 0} {incr s} {
       set d [expr {$d >> 1}]
   }

   for {} {$count > 0} {incr count -1} {
       set a [expr {2 + int(rand()*($n - 4))}]
       set x [expr {modexp($a, $d, $n)}]
       if {$x == 1 || $x == $n - 1} continue
       for {set r 1} {$r < $s} {incr r} {
           set x [expr {modexp($x, 2, $n)}]
           if {$x == 1} {return false}
           if {$x == $n - 1} break
       }

if {$x != $n-1} {return false}

   }
   return true

}

proc max_left_truncatable_prime {base} {

   set stems {}
   for {set i 2} {$i < $base} {incr i} {

if {[is_prime $i]} { lappend stems $i }

   }
   set primes $stems
   set size 0
   for {set b $base} {[llength $stems]} {set b [expr {$b * $base}]} {

# Progress monitoring is nice once we get to 10 and beyond... if {$base > 9} { puts "\t[llength $stems] candidates at length [incr size]" } set primes $stems set stems {} foreach s $primes { for {set i 1} {$i < $base} {incr i} { set n [expr {$b*$i + $s}] if {[is_prime $n 1]} { lappend stems $n } } }

   }
   # Could be several at same length; choose largest
   return [tcl::mathfunc::max {*}$primes]

}

for {set i 3} {$i <= 20} {incr i} {

   puts "$i: [max_left_truncatable_prime $i]"

}</lang>

Output up to base 12 (tab-indented parts are progress messages):

3: 23
4: 4091
5: 7817
6: 4836525320399
7: 817337
8: 14005650767869
9: 1676456897
	4 candidates at length 1
	11 candidates at length 2
	39 candidates at length 3
	99 candidates at length 4
	192 candidates at length 5
	326 candidates at length 6
	429 candidates at length 7
	521 candidates at length 8
	545 candidates at length 9
	517 candidates at length 10
	448 candidates at length 11
	354 candidates at length 12
	276 candidates at length 13
	212 candidates at length 14
	117 candidates at length 15
	72 candidates at length 16
	42 candidates at length 17
	24 candidates at length 18
	13 candidates at length 19
	6 candidates at length 20
	5 candidates at length 21
	4 candidates at length 22
	3 candidates at length 23
	1 candidates at length 24
10: 357686312646216567629137
	4 candidates at length 1
	8 candidates at length 2
	15 candidates at length 3
	18 candidates at length 4
	15 candidates at length 5
	8 candidates at length 6
	4 candidates at length 7
	2 candidates at length 8
	1 candidates at length 9
11: 2276005673
	5 candidates at length 1
	23 candidates at length 2
	119 candidates at length 3
	409 candidates at length 4
	1126 candidates at length 5
	2504 candidates at length 6
	4746 candidates at length 7
	7727 candidates at length 8
	11257 candidates at length 9
	14860 candidates at length 10
	17375 candidates at length 11
	18817 candidates at length 12
	19027 candidates at length 13
	17594 candidates at length 14
	15192 candidates at length 15
	12106 candidates at length 16
	9292 candidates at length 17
	6621 candidates at length 18
	4466 candidates at length 19
	2889 candidates at length 20
	1799 candidates at length 21
	1109 candidates at length 22
	601 candidates at length 23
	346 candidates at length 24
	181 candidates at length 25
	103 candidates at length 26
	49 candidates at length 27
	19 candidates at length 28
	8 candidates at length 29
	2 candidates at length 30
	1 candidates at length 31
	1 candidates at length 32
12: 13092430647736190817303130065827539
	5 candidates at length 1
	13 candidates at length 2
	20 candidates at length 3
	23 candidates at length 4
	17 candidates at length 5
	11 candidates at length 6
	7 candidates at length 7
	4 candidates at length 8
13: 812751503

I think I'll need to find a faster computer to calculate much more of the sequence, but memory consumption is currently negligible so there's no reason to expect there to be any major problems.