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.
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 petunias in The Hitchhiker's Guide to the Galaxy, come into existence only to be extinguished before they have time to realize it, because 2 and 5 preceded by any digit in the range 1 to base-1 is not prime. Some numbers formed by preceding 3 or 7 by a digit in the range 1 to base-1 are prime. So 13,17,23,37,43,47,53,67,73,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. 123 = 510 which is prime, 223 = 810 which is not so 123 is the only candidate. 1123 = 1410 which is not prime, 2123 = 2310 which is, so 2123 is the only candidate. 12123 = 5010 which is not prime, 22123 = 7710 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:
BBC BASIC
Uses the Huge Integer Math & Encryption library from http://devotechs.com/ <lang bbcbasic> HIMEM = PAGE + 3000000
INSTALL @lib$+"HIMELIB" PROC_himeinit("HIMEkey") DIM old$(20000), new$(20000) h1% = 1 : h2% = 2 : h3% = 3 : h4% = 4 FOR base% = 3 TO 17 PRINT "Base "; base% " : " FN_largest_left_truncated_prime(base%) NEXT END DEF FN_largest_left_truncated_prime(base%) LOCAL digit%, i%, new%, old%, prime%, fast%, slow% fast% = 1 : slow% = 50 old$() = "" PROC_hiputdec(1, STR$(base%)) PROC_hiputdec(2, "1") REPEAT new% = 0 : new$() = "" PROC_hiputdec(3, "0") FOR digit% = 1 TO base%-1 SYS `hi_Add`, ^h2%, ^h3%, ^h3% FOR i% = 0 TO old%-1 PROC_hiputdec(4, old$(i%)) SYS `hi_Add`, ^h3%, ^h4%, ^h4% IF old% OR digit% > 1 THEN IF old% > 100 THEN SYS `hi_IsPrime_RB`, ^fast%, ^h4% TO prime% ELSE SYS `hi_IsPrime_RB`, ^slow%, ^h4% TO prime% ENDIF IF prime% THEN new$(new%) = FN_higetdec(4) : new% += 1 ENDIF NEXT NEXT SYS `hi_Mul`, ^h1%, ^h2%, ^h2% SWAP old$(), new$() SWAP old%, new% UNTIL old% = 0 = new$(new%-1)</lang>
Output:
Base 3 : 23 Base 4 : 4091 Base 5 : 7817 Base 6 : 4836525320399 Base 7 : 817337 Base 8 : 14005650767869 Base 9 : 1676456897 Base 10 : 357686312646216567629137 Base 11 : 2276005673 Base 12 : 13092430647736190817303130065827539 Base 13 : 812751503 Base 14 : 615419590422100474355767356763 Base 15 : 34068645705927662447286191 Base 16 : 1088303707153521644968345559987 Base 17 : 13563641583101
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 ...
Eiffel
As there is currently no implementation for arbitrary precision integers this example only works for base 3 to base 9. <lang Eiffel> class LARGEST_LEFT_TRUNCABLE_PRIME
create make
feature
make -- Tests find_prime for different bases. local i: INTEGER decimal: INTEGER_64 do from i := 3 until i = 10 loop largest := 0 find_prime ("", i) decimal := convert_to_decimal (largest, i) io.put_string (i.out + ":%T" + decimal.out) io.new_line i := i + 1 end end
find_prime (right_part: STRING; base: INTEGER) -- Largest left truncable prime for a given 'base'. local i, larger, larger_dec: INTEGER_64 right: STRING prime: BOOLEAN do from i := 1 until i = base loop create right.make_empty right.deep_copy (right_part) right.prepend (i.out) larger := right.to_integer_64 if base /= 10 then larger_dec := convert_to_decimal (larger, base) if larger_dec < 0 then io.put_string ("overflow") prime := False else prime := is_prime (larger_dec) end else prime := is_prime (larger) end if prime = TRUE then find_prime (larger.out, base) else if right_part.count > 0 and right_part.to_integer_64 > largest then largest := right_part.to_integer_64 end end i := i + 1 end end
largest: INTEGER_64
convert_to_decimal (given, base: INTEGER_64): INTEGER_64 -- 'given' converted to base ten. require local n, i: INTEGER st_digits: STRING dec: REAL_64 do n := given.out.count dec := 0 st_digits := given.out from i := 1 until n < 0 or i > given.out.count loop n := n - 1 dec := dec + st_digits.at (i).out.to_integer * base ^ n i := i + 1 end Result := dec.truncated_to_integer_64 end
is_prime (n: INTEGER_64): BOOLEAN --Is 'n' a prime number? require positiv_input: n > 0 local i: INTEGER max: REAL_64 math: DOUBLE_MATH do create math if n = 2 then Result := True elseif n <= 1 or n \\ 2 = 0 then Result := False else Result := True max := math.sqrt (n) from i := 3 until i > max loop if n \\ i = 0 then Result := False end i := i + 2 end end end
end </lang>
- Output:
3: 23 4: 4091 5: 7817 6: 4836525320399 7: 817337 8: 14005650767869 9: 1676456897
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)
J
I'm fairly sure this is correct but it sure has long run time for base 12. <lang J>Filter=: (#~`)(`:6) Filter=: (#~`)(`:6)
f=: 3 : 0 BASE=. y primes=. (1&p:)Filter DIGITS=. }. i. x: BASE R=. 0 A=. ,0 while. #A do.
R=. >: R B=. A A=. primes , , A +/ BASE #. (_,R) {. ,. DIGITS
end. >./ B ) l NB. f N where N is the base.
(,. f"0) (3+i.8) 3 23 4 4091 5 7817 6 4836525320399 7 817337 8 14005650767869 9 1676456897
10 357686312646216567629137
f 11
2276005673
</lang> What can I say? J is an executable mathematical notation better than algebraic notation. That is the point, not that the execution time of a particular implementation is such and such. The current implementation of J uses custom arbitrary length integer code instead of gmp (j is not optimized to the extent gmp is optimized with its assembly code for various platforms. J is now open source. Join the project and install gmp!). J is terser than algebra. J handles text. J is simpler than algebra. J's precedence rules of noun-verb sequences evaluate strictly right to left. Versus the 17 precedence rules of C++. Versus algebraic confusion. "What is the intent of "force/area*time"? Surprise! It probably means physical power, force/(area*time). Different precedence, confusing, from the algebraic evaluation of, in c or mathematica, 8/3*7 == 7*8/3. Consider the simplifications of electrodynamics equations from Maxwell's original 20 to the 4 in modern notation. J benefits thought. Think in J. www.jsoftware.com Clearly, the various ways to express a thought in J perform differently upon execution. Languages like forth, python, c were developed AS PROGRAMMING LANGUAGES. J was developed as a TOOL OF THOUGHT. To untrained eye J looks like line noise with bit 8 clear. Get over it. Colon (:) and dot (.) are inflections, creating related concepts to the uninflected idea. Dyadic + means add. Inflected dyadic +, means "greatest common denominator" which for binary operands is "or". Am I wrong? Does the j solution not have 3 times fewer lines of code than the other languages have? Are my programs typically longer than those programs of others? Honeycombs/Python Honeycombs#C
Java
Code:
<lang java>import java.math.BigInteger; import java.util.*;
class LeftTruncatablePrime {
private static List<BigInteger> getNextLeftTruncatablePrimes(BigInteger n, int radix, int millerRabinCertainty) { List<BigInteger> probablePrimes = new ArrayList<BigInteger>(); String baseString = n.equals(BigInteger.ZERO) ? "" : n.toString(radix); for (int i = 1; i < radix; i++) { BigInteger p = new BigInteger(Integer.toString(i, radix) + baseString, radix); if (p.isProbablePrime(millerRabinCertainty)) probablePrimes.add(p); } return probablePrimes; } public static BigInteger getLargestLeftTruncatablePrime(int radix, int millerRabinCertainty) { List<BigInteger> lastList = null; List<BigInteger> list = getNextLeftTruncatablePrimes(BigInteger.ZERO, radix, millerRabinCertainty); while (!list.isEmpty()) { lastList = list; list = new ArrayList<BigInteger>(); for (BigInteger n : lastList) list.addAll(getNextLeftTruncatablePrimes(n, radix, millerRabinCertainty)); } if (lastList == null) return null; Collections.sort(lastList); return lastList.get(lastList.size() - 1); } public static void main(String[] args) { int maxRadix = Integer.parseInt(args[0]); int millerRabinCertainty = Integer.parseInt(args[1]); for (int radix = 3; radix <= maxRadix; radix++) { BigInteger largest = getLargestLeftTruncatablePrime(radix, millerRabinCertainty); System.out.print("n=" + radix + ": "); if (largest == null) System.out.println("No left-truncatable prime"); else System.out.println(largest + " (in base " + radix + "): " + largest.toString(radix)); } }
}</lang>
Example:
java LeftTruncatablePrime 17 1000000 n=3: 23 (in base 3): 212 n=4: 4091 (in base 4): 333323 n=5: 7817 (in base 5): 222232 n=6: 4836525320399 (in base 6): 14141511414451435 n=7: 817337 (in base 7): 6642623 n=8: 14005650767869 (in base 8): 313636165537775 n=9: 1676456897 (in base 9): 4284484465 n=10: 357686312646216567629137 (in base 10): 357686312646216567629137 n=11: 2276005673 (in base 11): a68822827 n=12: 13092430647736190817303130065827539 (in base 12): 471a34a164259ba16b324ab8a32b7817 n=13: 812751503 (in base 13): cc4c8c65 n=14: 615419590422100474355767356763 (in base 14): d967ccd63388522619883a7d23 n=15: 34068645705927662447286191 (in base 15): 6c6c2ce2ceeea4826e642b n=16: 1088303707153521644968345559987 (in base 16): dbc7fba24fe6aec462abf63b3 n=17: 13563641583101 (in base 17): 6c66cc4cc83
Mathematica
<lang>LargestLeftTruncatablePrimeInBase[n_] :=
Max[NestWhile[{Select[ Flatten@Outer[Function[{a, b}, #2 a + b], Range[1, n - 1], #1], PrimeQ], n #2} &, {{0}, 1}, #1 != {} &, 1, Infinity, -1]1]</lang>
Example:
<lang>Do[Print[n, "\t", LargestLeftTruncatablePrimeInBase@n], {n, 3, 17}]</lang>
Output:
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
PARI/GP
Takes half a second to find the terms up to 10, with proofs of primality. The time can be halved without proofs (use ispseudoprime
in place of isprime
).
<lang parigp>a(n)=my(v=primes(primepi(n-1)),u,t,b=1,best); while(#v, best=vecmax(v); b*=n; u=List(); for(i=1,#v,for(k=1,n-1,if(isprime(t=v[i]+k*b), listput(u,t)))); v=Vec(u)); best</lang>
Perl
Similar to the Pari solution. Uses ntheory for primality tests and Math::GMPz for bigints that aren't dog slow. We use is_prob_prime in the loop which does the BPSW test, then generate a proof for the selected result.
a(18) has a max candidate list of 1,449,405 entries and takes a bit over 20 minutes to solve.
<lang perl>use ntheory qw/:all/; use Math::GMPz;
sub lltp {
my($n, $b, $best) = (shift, Math::GMPz->new(1)); my @v = map { Math::GMPz->new($_) } @{primes($n-1)}; while (@v) { $best = vecmax(@v); $b *= $n; my @u; foreach my $vi (@v) { push @u, grep { is_prob_prime($_) } map { $vi + $_*$b } 1 .. $n-1; } @v = @u; } die unless is_provable_prime($best); $best;
}
printf "%2d %s\n", $_, lltp($_) for 3 .. 17;</lang>
- Output:
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
Perl 6
Using the Miller-Rabin definition of is-prime from Miller-Rabin primality test#Perl 6, 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
Python
<lang python>import random
def is_probable_prime(n,k):
#this uses the miller-rabin primality test found from rosetta code if n==0 or n==1: return False if n==2: return True if n % 2 == 0: return False s = 0 d = n-1
while True: quotient, remainder = divmod(d, 2) if remainder == 1: break s += 1 d = quotient
def try_composite(a): if pow(a, d, n) == 1: return False for i in range(s): if pow(a, 2**i * d, n) == n-1: return False return True # n is definitely composite for i in range(k): a = random.randrange(2, n) if try_composite(a): return False return True # no base tested showed n as composite
def largest_left_truncatable_prime(base):
radix = 0 candidates = [0] while True: new_candidates=[] multiplier = base**radix for i in range(1,base): new_candidates += [x+i*multiplier for x in candidates if is_probable_prime(x+i*multiplier,30)] if len(new_candidates)==0: return max(candidates) candidates = new_candidates radix += 1
for b in range(3,24):
print("%d:%d\n" % (b,largest_left_truncatable_prime(b)))
</lang>
Output:
3:23 4:4091 5:7817 6:4836525320399 7:817337 8:14005650767869 9:1676456897 10:357686312646216567629137 11:2276005673 12:13092430647736190817303130065827539 13:812751503
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!
Racket
<lang racket>
- lang racket
(require math/number-theory)
(define (prepend-digit b d i n)
(+ (* d (expt b i)) n))
(define (extend b i ts)
(define ts* (for/list ([t (in-set ts)]) (for/set ([d (in-range 1 b)] #:when (prime? (prepend-digit b d i t))) (prepend-digit b d i t)))) (apply set-union ts*))
(define (truncables b n)
; return set of truncables of length n in base b (if (= n 1) (for/set ([d (in-range 1 b)] #:when (prime? d)) d) (extend b (- n 1) (truncables b (- n 1)))))
(define (largest b)
(let loop ([ts (truncables b 1)] [n 1]) (define ts* (extend b n ts)) (if (set-empty? ts*) (apply max (set->list ts)) (loop ts* (+ n 1)))))
(for/list ([b (in-range 3 18)])
(define l (largest b)) ; (displayln (list b l)) (list b l))
- Output
'((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))
</lang>
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 certainty [expr {[llength $primes] > 100 ? 1 : 5}] set stems {} foreach s $primes { for {set i 1} {$i < $base} {incr i} { set n [expr {$b*$i + $s}] if {[is_prime $n $certainty]} { 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.