Find largest left truncatable prime in a given base
You are encouraged to solve this task according to the task description, using any language you may know.
- ref 1: https://oeis.org/A103443
- related task 1: http://rosettacode.org/wiki/Truncatable_primes
Obviously the right most digit must be prime, so in base 10 candidate 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 ref 1 as you are able.
C
<lang c>#include <stdlib.h>
- 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(i);
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
<lang 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!