There's more than one way to do things PerlMonks

### Finding Primes

by Tommy (Chaplain)
 on Aug 14, 2003 at 00:28 UTC Need Help??
Tommy has asked for the wisdom of the Perl Monks concerning the following question:

Task: to find two prime numbers whose product is a 400 digit number. My math professor says such a number exists, and I'm out to find it. I turn to Perl.

So here is the code I am using to do it thus far. It can't possibly be the most efficient way to do this, but it's all I've got. I'm not very good at math. How can I do it better?

```#!/usr/bin/perl -w
use strict;

BEGIN { ++\$| and print "Starting...\n\n" }

my(@primes) = &gen_primes(eval(1 . '0' x 10));

PRIMESCAN: foreach (@primes) {
my(\$a) = \$_;
foreach (@primes) {
my(\$b) = \$_;
if (length(\$a * \$b) == 400) {
last PRIMESCAN and print <<__FOUND__
NUMBER DISCOVERED!
\$a * \$b = \${\\$a * \$b}
__FOUND__
} else { print <<__NOTYET__ }
\$a * \$b = \${\\$a * \$b}
__NOTYET__
}
}

sub gen_primes {
my(\$ceil) = shift || 0;
print "Getting primes...\n";
for (my(\$i) = 0; \$i < \$ceil; ++\$i) {
if (isa_prime(\$i)) {
print "   \$i is prime\n";
push(@primes,\$i);
}
}
return @primes
}

sub isa_prime {
my(\$candidate) = shift;
my(\$i) = 1;
while (++\$i < \$candidate) {
return undef if (\$candidate % \$i) == 0
}
return \$candidate
}
Tommy Butler, a.k.a. TOMMY

edited by ybiC: s/pre/code/ and s/tab/spacespace/

Replies are listed 'Best First'.
Re: Finding Primes
by Zaxo (Archbishop) on Aug 14, 2003 at 01:07 UTC

Any two 200 digit primes will do it, but your isa_prime() and gen_primes() subs are way too inefficient to find primes of that size. You are hoping to test primality by exhaustively checking every number less than the candidate as a divisor. That will take more than 10**200 divisions. If you can do one every 10 nanoseconds, that gives you 10**192 seconds, or around 10**175 ages of the universe to check a single number.

Has your math professor given you some reading on the subject of primality testing? If not, check The Prime Pages.

Math::Pari is very good for number theory kind of things.

After Compline,
Zaxo

Any two 200 digit primes will do it

Without knowing all the 200 digit primes, I can't say whether or not that is a true statement. I wouldn't assume it is sufficient, however. The square root of 10^399 is a little greater than this 200 digit number:

```3162277660168379331998893544432718533719555139325216826857504852792594
+438639238221344248108379300295187347284152840055148548856030453880014
+6905195967001539033449216571792599406591501534741133394841240
So, if there are two primes between 10^199 and that number, then there will be two 200 digit primes whose product is not a 400 digit number.

```-sauoq
"My two cents aren't worth a dime.";
```
First of all, there are an infinite number of primes and there is a theorem by Bertrand which states that there is a prime between any natural number n and 2*n. So there are two prime numbers which product is a 400 digit number.
Well, as 10^1 is a 2 digit number, 10^399 is a 400 digit number....:)

thor

If you can do one every 10 nanoseconds, That gives you 10**192 seconds, or around 10**175 ages of the universe. That to check a single number.

Yikes! And you know my professor said that the smallest possible starting point would be two numbers of 10**200. At least this is what I gathered from our conversation. Is that right?

Math::Pari is very good for number theory kind of things.

I'll have a look at that when I get home from school here. Thanks for the tip.

I'm still pretty lost here. I don't know where to start.

Tommy Butler, a.k.a. TOMMY
Re: Finding Primes
by simonm (Vicar) on Aug 14, 2003 at 01:30 UTC
In addition to the more general problems pointed out by Zaxo, there are some run of the mill coding improvements that could be made.
&gen_primes(eval(1 . '0' x 10));

The number 1-with-ten-zeros-after-it can be written as 10 ** 10.

foreach (@primes) { ... foreach (@primes) { ... } }

You can make a brute force search run faster by cutting off areas of the search space that you know can be skipped.

In this case, there are a number of ways you can shrink the search space; for example, next if ( \$a > \$b ); allows you to avoid checking twice for \$a * \$b and \$b * \$a, and next PRIMESCAN if ( length(\$a * \$b) > 400; would let you stop testing a sequence when the numbers got too large.

Even with these tricks, brute-force isn't going to be an efficient approach, for the reasons discussed elsewhere.

last PRIMESCAN and print <<__FOUND__

This invokes "last" before "print" -- meaning that the print never gets invoked. You probably want print <<__FOUND__ and last PRIMESCAN.

Thanks for the tips :-)

The number 1-with-ten-zeros-after-it can be written as 10 ** 10.

Thanks! That's what I was looking for there. Yup, the eval(1 . '0' x 10) was a hack.

This invokes "last" before "print" -- meaning that the print never gets invoked. You probably want print <<__FOUND__ and last PRIMESCAN.

Yes, quite. Thanks.

-- Tommy Butler, a.k.a. TOMMY
Re: Finding Primes
by Abigail-II (Bishop) on Aug 14, 2003 at 09:56 UTC
You might want to start reading at http://mathworld.wolfram.com/PrimalityTest.html. It's only recent that someone discovered a "efficient" algorithm to determine the primality of a number. But it's still an algorithm that's bounded by the 12th power of the number of bits the number is written in. There are faster algorithms, but they either may give false positives, or false negatives. You may want to look into the Lucas-Lehmer test: http://mathworld.wolfram.com/Lucas-LehmerTest.html

You should realize that this is not a trivial problem.

Abigail

The Rabin-Miller test is a fast method to determine either a) if a number is composite or b) a (arbitrarily high) probability that a number is prime. "Introduction to Algorithms" by Cormen, et al has a good section on this (they call it Miller-Rabin). Here's my (very lightly documented) implementation--
```use warnings;
use strict;

use Getopt::Long;
use Math::BigInt;

my %pars = ( verify => 0 );
GetOptions ( \%pars, 'verify', 'help', 'usage' )
or die "Bad GetOptions";

die "primes number [ more numbers ] \nreturns the number's primality s
+tatus"
if \$pars{help} || \$pars{usage};

my @inputs = @ARGV
or die "no numbers given for primality test\n";

foreach ( @inputs ) {

my \$number = Math::BigInt->new ( \$_ );
my ( \$primality_prob, \$witness ) = miller_rabin_prime ( \$number );

print "\$number is ",
( \$primality_prob )
? "prime with prob \$primality_prob\n"
: "not prime, witness \$witness\n";
}

sub miller_rabin_prime {

my \$n       = shift () or die "no n in MillerRabin";
my \$s       = shift () || 25;

my \$a = Math::BigInt->new ( 0 );

for ( my \$j = 1; \$j <= \$s; ++\$j ) {

# \$a is a random int in [1,n-1]
#
my \$ranlim = \$n-1;
do {
\$a = Math::BigInt->new ( sprintf "%d", rand ( \$ranlim ) +
+1 );
\$ranlim /= 2;
} while ( 0 > \$a );

if ( 1 == Witness ( \$a, \$n ) ) {
return 0, \$a;
}
}
my \$prob = 1 - (3/4)**\$s;
return \$prob, undef;
}

sub Witness {

my \$a = shift () or die "no a in Witness";
my \$n = shift () or die "no n in Witness";

# n-1 = ( 2**t ) * u
#
my ( \$t, \$u ) = get_t_u ( \$n - 1 );

my @x;
\$x[0] = mod_exp ( \$a, \$u, \$n );

for ( my \$i = 1; \$i <= \$t; ++\$i ) {
\$x[\$i] = mod_exp ( \$x[\$i-1], 2, \$n );
return 1
if (1 == \$x[\$i]) && (1 != \$x[\$i-1]) && (( \$n - 1 ) != \$x[\$
+i-1]);
}
return 1
if 1 != \$x[\$t];

return 0;
}

sub get_t_u {

my \$m = shift () or die "no m in get_t_u";

my ( \$t, \$u ) = ( 0, 0 );

while ( 0 == \$m % 2 ) {
++\$t;
\$m /= 2;
}
\$u = \$m;
return ( \$t, \$u );
}

sub mod_exp {

# i**j mod n
#
# works with Math::BigInt in addition to regular scalars
#
my \$i = shift () or die "no i in mod_exp";
my \$j = shift () or die "no j in mod_exp";
my \$n = shift () or die "no n in mod_exp";

# return 1 for zero exponents
#
my \$result = \$i - \$i + 1;
return \$result unless \$j;

my \$pow2 = \$i;

while (1) {
if ( \$j % 2 ) {
\$result = ( \$pow2 * \$result ) % \$n;
return \$result unless --\$j;
}
\$j /= 2;
\$pow2 = ( \$pow2 * \$pow2 ) % \$n;
}
}
Re: Finding Primes
by Foggy Bottoms (Monk) on Aug 14, 2003 at 09:58 UTC
There are a lot of algorithms out there to help you deal with prime numbers. Here are a couple ordered by efficiency, the first being the most efficient. If you run the code in C, you'll notice there's a sharp drop in time consumption between the first and 3rd examples but that the change is very slight between the 3rd and 4th...
\$n is the number we're trying to figure out whether it's prime or not...
• divide \$n by all the numbers preceding it until you find a multiple of \$n or until you run out of numbers
• Same method as above, but instead of going all the way to \$n, you can solely consider all the numbers till 1+trunc(sqrt(\$n))...
• Same method as above but exclude all the pair numbers (easily done by having a loop where the index is increased twice)
• Same method as above plus store all the previous prime numbers you've already found and exclude their product...
None of them are at all feasible when trying to determine the primeness of a 200 digit number. The square root of a 200 digit number is a 100 digit number. Even if the first 200 digit number you are going to try is prime, and take your last suggestion, you'd need to keep track of all the prime numbers less than 100 digits long.

Considering there are about n / ln n prime numbers less than n, you'd need to keep track of about 4 * 10**98 prime numbers.

Many people in this thread just don't seem to realize how big a 200 digit number is.

Abigail

Hi Abigail,

You're quite right but I never actually meant to give a solution to the 200-digit long number, I was only recapping the typical algorithms one can come up with, pretenselessly ;-)...
There must be an option with n! numbers because they grow so quickly - I thought that n!+1 might a prime number but an easy counter-example is n=4 (which'd give n!=24 hence n!+1=25 = 5*5, barely what you can call a prime number...)
Re: Finding Primes
by mildside (Friar) on Aug 14, 2003 at 03:11 UTC
I think you will need to use something like Math::BigInt so that you can handle integers with lots of digits of precision. Good luck!

Cheers!

Math::BigInt is way too slow for such a task. You want to do the number crunching in C or FORTRAN.

Abigail

Re: Finding Primes
by BrowserUk (Pope) on Aug 14, 2003 at 13:13 UTC

Here's a 400-digit number that is probably the product of 2 primes. Disproving it is left as an exercise for the teacher:)

If he succeeds, then multiply any pairing of numbers from the following lists and try again. That gives you 100 possible candidates. If they all turn out bogus, generating another few thousand possibles doesn't take too long:)

```10000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000188
70000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000054621

This is the product of the first of each of the following two groups of probable primes:

10 x 200-digit probable primes:

10 x 201-digit probable primes

Credit goes to Olivier Langlois ( olanglois @ sympatico . ca ), John Moyer ( jrm @ rsok . com ) (and Shakuntala Devi).

Examine what is said, not who speaks.
"Efficiency is intelligent laziness." -David Dunham
"When I'm working on a problem, I never think about beauty. I think only how to solve the problem. But when I have finished, if the solution is not beautiful, I know it is wrong." -Richard Buckminster Fuller
If I understand your problem, I can solve it! Of course, the same can be said for you.

ROTFLH!!! Oh dear! BrowserUK ++ except I've used all mine up for today. I could have spent them all in this thread!

I'm no scientist (but I suspect a significant number of Monks in this thread really are). However, as I understand it, the following is true for some arbitrary approximation of truth

10**79 <= fpitu <= 10**81

`fpitu' = fundamental particles in the universe

So a 200 digit number is rather more than (111 orders of magnitude than?) is needed to allocate to each of these particles a unique number.

The estimation of the secondary storage space (!) requirement for this ultimate registry is also left to the numerically challenged.

hagen
Re: Finding Primes
by rinceWind (Monsignor) on Aug 14, 2003 at 11:54 UTC
Others have pointed out how long the computation could potentially take. I think that your math professor is reasonably confident that the problem will not be solved by brute force, in your lifetime.

However, discovery of a new method of factorizing products of large primes could make you a rich man indeed. This is the mathematics that is underlying the RSA algorithm - a general prime factorization algorithm that delivered results quickly, would render RSA crackable. Besides the money you could get from licensing the software, think how much the NSA would pay you to shut up ;).

think how much the NSA would pay you to shut up ;)

An automobile wreck with a drunk driver doesn't cost much. A gangland turf battle at an area McDonald's costs a little more. National Security has nothing to do with individual security.

--
[ e d @ h a l l e y . c c ]

National Security is globally scoped since your is only lexical. :)

Re: Finding Primes
by Anonymous Monk on Aug 14, 2003 at 10:55 UTC
The Following number is an example of a number is a product of two primes and has 400 digits

1000000000000000000000000000000000000000
0000000000000000000000000000000000000000
0000000000000000000000000000000000000000
0000000000000000000000000000000000000000
0000000000000000000000000000000000000666
0000000000000000000000000000000000000000
0000000000000000000000000000000000000000
0000000000000000000000000000000000000000
0000000000000000000000000000000000000000
000000000000000000000000000000000078489

Is it the product of two prime one of them being 10**199 + 153 the other 1**199 + 513.

note that I didn't found these numbers by using Perl

...nor did you use your fingers and toes to count the digits?
I make it 399, not 400.... and isn't (1**199)=1 ??

arithmetically_challenged
Re: Finding Primes
by thor (Priest) on Aug 14, 2003 at 12:47 UTC
If you're looking for 200-digit primes, here are a couple...

thor

Re: Finding Primes (A Solution)
by sauoq (Abbot) on Aug 15, 2003 at 00:08 UTC
Task: to find two prime numbers whose product is a 400 digit number.

This

```1022731116849984128382970140461725893199
2883532546140280576806704317960648119367
2377124699979471607606314465365589539394
1805302550414778332749072722843441732197
4205545703859751798082548809071278533442
1766932914859818046781739244057376689689
2255723329228491768773454858783052995068
4787708109345332331133139141202295689475
2588959029409151843111542915460444247984
1841531569215551416798910788156968134891
is the product of a very well known 386 digit prime, 2**1279-1, and a somewhat less famous 14 digit prime, 98264582985493, one of the known 14 digit prime divisors of Googolplex + 10.

Sometimes perl is the right tool for the job... and sometimes Google is.

```-sauoq
"My two cents aren't worth a dime.";
```

Amazing! Wow! That's awesome! ++sauoq :-)

This thread has been very educational to say the least. It has opened up a door to a whole realm of thought that I didn't even know existed. I certainly didn't have any idea how big a 400 digit number was, and I didn't know that finding prime numbers was such an important quest.

Thanks everyone for your very helpful input.

--
Tommy Butler, a.k.a. TOMMY

Re: Finding Primes
by johndageek (Hermit) on Aug 14, 2003 at 15:59 UTC
A quick question:

100 * 100 = 10000 (5 digits)
900 * 900 = 810000 (6 digits)

not all sets of numbers with the same mubers of digits will result in a product with exactly double the number of digits. So is exactly 400 digits the requirement?

Just for fun.
John

It's easy to determine how many digits you'll have in the resulting number.
In your example : 100*100 is actually 102*102 hence 104 which turns out to be 5 digits. The rule, hence is :
Take any number, divide repeatedly by ten until you reach a number greater than or equal to 1 and strictly lower than 10. The number of times you divided that number is actually the number of digits minus 1 that number has.
Example : 84212 divided by 10 is 8421.2 (1)
8421.2 divided by 10 is 842.12 (2)
842.12 divided by 10 is 84.212 (3)
84.212 divided by 10 is 8.4212 (4)
Hence the number has 4+1=5 digits !!!
It's easy to determine how many digits you'll have in the resulting number.

I'm a bit puzzled: sure, logarithms let you replace multiplication with addition -- part of the magic of slide rules -- but isn't taking the base-10 log of something typically as much or more labor intensive then doing the multiplication?

For example: 321 * 311 produces five digits, while 321 * 312 produces six digits.

Sure, you could work out that 10 ** 2.5065 * 10 ** 2.493 = 10 ** 4.9995 and thus five digits, while 10 ** 2.5065 * 10 ** 2.494 = 10 ** 5.0005 and thus six, but is that really easier?

Why not just use something like this:

```#!/usr/bin/perl -wl

sub logb10 {
return log(+shift)/log(10);
}

print int(logb10(100) + logb10(100) + 1);
print int(logb10(321) + logb10(311) + 1);
print int(logb10(321) + logb10(312) + 1);
__DATA__
output:
5
5
6

Or maybe we could just pull out the magic length method and stop trying to focus on a problem that isn't really a problem?

Re: Finding Primes
by jarich (Curate) on Aug 20, 2003 at 08:29 UTC

I don't know if you've seen it (it seems to get way more hits than I'd expect it to) but I did a final year project on Prime numbers and their generation (with source code for Pari gp) that you might be interested in.

If I remember correctly I was able to generate probable primes of 200 digits very rapidly with the code I provided. So you should be able to get lots of 400 composite numbers if you want to. You'll have to download Pari gp though, but that shouldn't be too hard. ;)

Have a look over at: http://perltraining.com.au/~jarich/Primes/index.html

I hope that helps.

jarich.

Create A New User
Node Status?
node history
Node Type: perlquestion [id://283725]
Approved by fglock
Front-paged by broquaint
help
Chatterbox?
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others examining the Monastery: (2)
As of 2018-01-21 06:05 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
How did you see in the new year?

Results (227 votes). Check out past polls.

Notices?