Primal Pi

[theStack]

Hi folks,

I'm very curious about your way of solving this challenge, since all I could think of was just brute-forcing witht the fast Miller-Rabin prime test. I was too lazy to implement the algorithm myself in any language, so I just took Java and used BigInteger, who has a isProbablePrime() method:

Code: Select all

import java.Math.BigInteger;

public class PrimalPi {
    private static final String sequence = "14159 [.................]"
    
    public static void main(String[] args)
    {
        for (int i=0; i<sequence.length(); i++) {
            System.out.println("Testing with startindex " + i);
            for (int j=sequence.length(); j>i; j--) {
                String totest_str = sequence.substring(i, j);
                char lastchar = totest_str.charAt(totest_str.length()-1);
                if (lastchar == '0' ||
                    lastchar == '2' ||
                    lastchar == '4' ||
                    lastchar == '5' ||
                    lastchar == '6' ||
                    lastchar == '8')
                    continue;

                BigInteger totest = new BigInteger(totest_str);
                if (totest.isProbablePrime(4)) {
                    System.out.println("Probable prime: " + totest_str);
                }
            }
        }
    }
}

So, my only optimization here is to check for the last digit, and skip the test if it's 0/2/4/5/6/8 because numbers ending with one of these digits are never primes (at least numbers > 9).

The other idea I had was searching in a prime database for all primes between let's say 500 and 2048 digits, but unfortunately the only site I found on the net had a strange format and you couldn't see the full prime numbers itself but only the way they would be calculated.
See http://primes.utm.edu/primes/search.php

So, is there any other way than brute-forcing? Or any other optimization I've not thought of?

[MerickOWA]

Thats basically the method I chose, only I had to write a fermat primality test, which I enjoyed
I switched the loops around, trying all possibilities of a given length, but basically the same idea.

[gfoot]

Hmm, I might have got the names wrong, but I thought the fermats-little-theorem method just boiled down to this in Python:

Code: Select all

pow(2,n-1,n) == 1

for n a suspect prime, and that's what I used. Python's pow() with the optional modulo seems pretty well optimised for long integers.

You can run it again with 3, 5, or some other prime instead of 2, to reinforce the result; in any case though, it's just a statement of probability and you ought to do some more checking.

However, submitting it as the answer seems a workable shortcut test.

So I too just ran a program that looped down from 2048 testing subsequences ending in odd digits for being 2-PRP probable primes. I did pre-screen for divisibility by small primes (up to 2048, calculated with a simple sieve) because that can be done almost instantaneously, and it's worth avoiding the pow calculation.

I got most of this from a site on primes which I didn't keep a link to. It suggested the pre-screening, then the 2-PRP testing via this method, followed (if pow(...)==1), and then either using classical methods or modern (elliptic curve) methods, depending on whether you can find any useful factors of n-1 and n+1 (which are required by the quicker classical methods).

I was on the point of downloading ecpp, but it's a binary-only distribution and I was put off for long enough (30 minutes?) for my existing script to find the answer.

I can't see any other way than brute-forcing here, because you're essentially given a random pile of digits to take sequences from, and being asked for a specific 'longest' sequence, you need to actually find the longest one, not just heuristically search for something long enough. Dunno if that argument makes any sense though!

[TheHiveMind]

Basically the same as everyone else, using GMP with C.

For brevity I cut out the prescreening, that would lower the 1.5 minutes it took to run a little bit ( mpz_probab_prime_p() performs some trial divisions first as well).

Code: Select all

#include <stdio.h>
#include <gmp.h>

int main()
{

	mpz_t candidate;
	char pi[2049];
	
	mpz_init(candidate);
	scanf("%2048s",pi);
	
	int length=2049;
	int offset=0;
	char tempchar;
	
	while(--length)
		for(offset=0;offset<2049-length;offset++)
		{
			tempchar = pi[offset+length];
			pi[offset+length]=0;
			mpz_set_str(candidate,pi+offset,10);
			if(mpz_probab_prime_p(candidate,10)) return puts(pi+offset);
			pi[offset+length]=tempchar;
		}
	return 0;
}

Afterwards, it took 13 hours to generate an Atkin-Goldwasser-Kilian-Morain certificate using PRIMO, so rest assured that the number here is really a prime. (Message me if you're interested in the certificate )

[klavierspieler21]

I used a simple nested loop in Maple, which has the function isprime(n). And the first prime I found, in decreasing string length, was outputted.

[nighthalk]

i did it in c# using a biginteger library that had a "isprobablyprime" function built in... started from the max length and worked backwords. first checked last digit like above, then used "first 2000 primes" to quick check, added all of those that passed to a list, then did the full blown "is probably prime" from there, even with a confidence of one the first match was the right answer. took 10 minutes to correctly code, and like 20 to solve it. did anyone actually check if its truly prime?

[Stack]

My python code ran in about 3 hours. It used the Miller-Rabin test.

OPTIMIZATIONS:
In response to user theStack (Not to be confused with me, Stack) your optimization really wasn't that great. If you want, Wikipedia "Prime Sieve" or "Sieve of Erastosthenes." Basically, if it's modulus by a prime number turns out to be zero, then you've got yourself a composite number. I just sieved using the first few primes. I was too lazy to type a lot of code.

Also, I made my code so that if it's found a prime of length 1024, it will only check substrings longer than 1024. (Clever, eh? ... Not really.)

CODE:

Code: Select all

#This code is for the Primal Pi challenge on hacker(D0T)org.
#Written by stack.  (With a bit of leeched code)
#Runs in about 3-4 hours on my relatively fast computer.
#If you want to run this:
#Make a file with the first 2048 digits of pi and call it pi2048.txt in the same directory

import sys
import random

##### Start leeched Miller Rabin code:

def toBinary(n):
  r = []
  while (n > 0):
    r.append(n % 2)
    n = n / 2
  return r

def test(a, n):
  b = toBinary(n - 1)
  d = 1
  for i in xrange(len(b) - 1, -1, -1):
    x = d
    d = (d * d) % n
    if d == 1 and x != 1 and x != n - 1:
      return True # Complex
    if b[i] == 1:
      d = (d * a) % n
  if d != 1:
    return True # Complex
  return False # Prime

def MillerRabin(n, s = 50):
  for j in xrange(1, s + 1):
    a = random.randint(1, n - 1)
    if (test(a, n)):
      return False # n is complex
  return True # n is prime

##### The ABOVE modified Miller-Rabin code leeched from:
##### http://snippets.dzone.com/posts/show/4200
##### Why reinvent the wheel when you're too lazy to? ;)
##### I won't be commenting the Miller Rabin code for
##### you to understand.  That's what Wikipedia is for.

file = open("pi2048.txt")#Inside of pi2048.txt is the first 2048 digits of pi
pi=file.read()#Reads pi2048.txt into pi
file.close()#Frees up a bit of memory
pilen=len(pi)#Length of pi
glen=2 #Greatest length found (Changes while the code runs)

l=0
status=l-1
sys.stdout.write("Starting with ")
while l<pilen-200:
  m=glen
  while m<=2048:
    
    
#####Too lazy to make a longer sieve:
    if ((int(pi[l:m+l])%2!=0) and (int(pi[l:m+l])%3!=0) and (int(pi[l:m+l])%5!=0) and (int(pi[l:m+l])%7!=0) and (int(pi[l:m+l])%11!=0) and (int(pi[l:m+l])%13!=0) and (int(pi[l:m+l])%17!=0) and (int(pi[l:m+l])%19!=0) and (int(pi[l:m+l])%23!=0) and (int(pi[l:m+l])%29!=0) and (int(pi[l:m+l])%31!=0) and (int(pi[l:m+l])%37!=0) and (int(pi[l:m+l])%41!=0) and (int(pi[l:m+l])%43!=0) and (int(pi[l:m+l])%47!=0) and (int(pi[l:m+l])%53!=0) and (int(pi[l:m+l])%59!=0) and (int(pi[l:m+l])%61!=0) and (int(pi[l:m+l])%67!=0) and (int(pi[l:m+l])%71!=0) and (int(pi[l:m+l])%73!=0) and (int(pi[l:m+l])%79!=0) and (int(pi[l:m+l])%83!=0) and (int(pi[l:m+l])%89!=0) and (int(pi[l:m+l])%97!=0) and (int(pi[l:m+l])%101!=0)):
      if MillerRabin(int(pi[l:m+l]), 2):#First param is the number to test
                                         #Second param is the number of times
                                         #to run the test (the more times, the more accurate)
                                         #Because it only gives the probability of the number being a prime
                                         #Again, read more in Wikipedia
        print "-----"
        print "Loop: ", l
        print "Len: ", len(pi[l:m+l])
        print pi[l:m+l]
        glen=len(pi[l:m+l])-1
    if l>status:#Prints the status of the code
      print "Loop", l
      status+=100#Change this for the status to update more frequently
    m+=1
  l+=1

[/size]

inb4 tl;dr

[Hippo]

As in the first post ... Biginteger with MillerRabin test were used, may be my order of tests was better.
At first attempt I have started shortening the substrings, but it was not giving subresults ... so I could not guess when it finishes so I used another approach. It's advantage is eliminating lot of tests by checking the last digit. (I have probably used a subresult for initialisation in next run ...)

Code: Select all

	/**
	 * finds longest prime number in first 2048 decimal digits of pi
	 */
	public static void primes() {
		int bestlen=1977;int bestend=2003;
		for(int e=2050;e>bestlen+2;e--) {
			System.out.printf("\nending "+e+" on "+piDigits.substring(e-1,e));
			if ("024568".indexOf(piDigits.substring(e-1,e))<0) {
				for(int l=e-2;l>bestlen;l--) {
					BigInteger maybePrime = new BigInteger(piDigits.substring(e-l,e));  
					if (maybePrime.isProbablePrime(20)) {
						System.out.printf("match at "+e+" of length "+l);
						bestlen=l;bestend=e;
					} else {
						System.out.printf(".");
					}
				}
			}
		}
		System.out.printf("best match of length "+bestlen+" ending on "+bestend);
		writeToClipboard(piDigits.substring(bestend-bestlen,bestend));
	}

[eulerscheZahl]

I used pari/gp for this, as it can handle huge numbers and has built-in mathematical functions (it could even solve Bigger Fib directly).

Code: Select all

{
	pidigits = 1415926535897932[...];
	found = 0;
	forstep(len = 2048,1,-1,
		for (start = 0,2048-len,
			part = (pidigits \ 10^start) % 10^len;
			if (ispseudoprime(part) && found == 0,print(part); found++);
		);
		if (found > 0, break);
	);
}

took about 20 seconds.