- (The Forth program for listing small primes is here: https://joels-programming-fun.blogspot.com/2015/11/small-prime-sieve-in-forth.html.
- The original C version is here: https://joels-programming-fun.blogspot.com/2015/11/small-prime-sieve-in-c.html.
- The C program for counting small primes is here: https://joels-programming-fun.blogspot.com/2019/02/taking-sieve-too-far.html.)
You say you didn't know 20,000,003 was a small number, much less a small prime?
It's small enough to calculate by this simple method on modern personal computers, so, yes, it's small. In terms of cryptology and cryptography, also, it is small.
Well, it turns out the Forth version was directly usable (on modern hardware, if you set the terminal emulator to remember enough lines) to determine if dates written as decimal numbers like this:
20191231for yesterday, are prime. But you have to go searching through the list for the number yourself.
It would be easier to be able to use the computer to look for it for you. This is easy to do in Forth. (Taking a clue from the BASIC version.)
Here's the code I'm going to start with, from the post on small primes in Forth, modified to only count, not print, the primes:
( Archetypical implementation of the counting sieve of eratosthenes in FORTH -- fig, bif-c -- )
( using more of the FORTH idiom. )
( Copyright Joel Rees, 2015, 2020 )
( By Joel Rees, Amagasaki, Japan, 2015, 2020 )
( All rights reserved. )
( Permission granted by the author to use this code )
( for any purpose, )
( on condition that substantial use )
( shall retain this copyright and permission notice. )
VOCABULARY sieve-local ( Make a local symbol table. )
sieve-local DEFINITIONS ( Switch to the local vocabulary. )
30000000 CONSTANT MAXSIEVE
MAXSIEVE 1 - 2 / CONSTANT FINALPASS
5 CONSTANT DISPWIDTH ( enough digits to display MAXSIEVE )
0 VARIABLE sieve ( Some old FORTHs don't provide a default behavior for CREATE )
( gforth will leave the zero on the stack. Some old FORTHs need an initial value. )
HERE sieve - DUP ( Some old FORTHs don't provide a CELL width. )
MAXSIEVE SWAP - ALLOT ( Allocate the rest of the byte array. )
CONSTANT CELLWIDTH ( To show how this can be done. )
( Arrays are handled differently by different Forths. )
( If CELLWIDTH were just to save a wasted entry in the sieve array, )
( it wouldn't be worth it here. )
( I'm assuming there are likely to be other uses. )
: sieveInit ( -- adr )
0 sieve C! ( 0 is not prime. )
0 sieve 1+ C! ( 1 is not prime. )
sieve MAXSIEVE 2 DO ( set flags to true for 2 to MAXSIEVE. )
-1 OVER I + C! LOOP ( sieve pointer -- still on stack. )
;
: primePass ( adr prime -- adr )
MAXSIEVE OVER DUP + DO ( start at first multiple -- double. )
OVER I + 0 SWAP C! ( clear at this multiple. )
DUP +LOOP ( next multiple )
DROP ; ( sieve address still on stack. )
: findPrimes ( adr -- adr )
FINALPASS 2 DO ( clear flags at all multiples. )
DUP I + C@ IF ( don't bother if not prime. )
I primePass
ENDIF
LOOP ; ( sieve still on stack. )
: countPrimes ( adr -- )
0 SWAP
MAXSIEVE 0 DO
DUP I + C@ IF
SWAP 1+ SWAP
ENDIF
LOOP DROP . ." Primes less than " MAXSIEVE . ;
FORTH DEFINITIONS
: sieveMain ( -- )
[ sieve-local ] sieveInit
findPrimes
countPrimes ;
sieveMain
People who do a lot of Forth will have a few complaints about this code. It runs fast, and it leaves the prime sieve in memory for later querying, but the queries are just slightly arcane:
sieve 20000003 + c@ .will print "255", which tells you that 20,000,003 is prime. So, we define a query word:
: isprime? sieve + c@ ifAnd now you can type
." Yes."
else
." No."
endif
;
20000003 isprime?and it responds, Yes. And you can type
20000001 isprime?and it responds, No. (You know 20000002 is not prime, so no need to test that, right? Right?)
Very convenient. Let's compare that to doing something similar in C.
(Incidentally, a full functionality prime query and factoring command, the primes (source) utility command (manual), can be found in bsd games, available in most Linux distros as well as in BSD, and in Cygwin, as well, for MSWindows sufferers.)
We'll need to pick one of the existing prime sieves and add some code to read numbers from the input. The prime counting routine deals with the largest primes, so I'll borrow it:
/* Archetypical queriable implementation of the counting sieve of eratosthenes in C
** By Joel Rees, Amagasaki, Japan, 2015, 2018, 2020
** All rights reserved.
** Permission granted by the author to use this code
** for any purpose,
** on condition that substantial use
** shall retain this copyright and permission notice.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXSIEVE 30000000L /* Maximum testable of 29999999L */
#define FINALPASS ( ( MAXSIEVE + 1 ) / 2 )
char * sieve;
int main( int argc, char * argv[] )
{
unsigned long index, prime;
unsigned long count = 0;
unsigned long maxsieve = MAXSIEVE;
unsigned long finalpass = FINALPASS;
unsigned parameterX = 1;
char * unscanned;
if ( argc > 1 )
{
if ( ( strcmp( "--help", argv[ 1 ] ) == 0 )
|| ( ( argv[ 1 ][ 0 ] == '-' )
&& ( ( argv[ 1 ][ 1 ] == 'h' ) || ( argv[ 1 ][ 1 ] == '?' ) ) ) )
{
printf( "usage:\n%s [ -max MAX ] ( prime < MAX ) ...\n", argv[ 0 ] );
return EXIT_SUCCESS;
}
else if ( ( argc > 2 ) && ( strcmp( "-max", argv[ 1 ] ) == 0 ) )
{
maxsieve = strtoul( argv[ 2 ], &unscanned, 0 );
if ( ! ( unscanned > argv[ 2 ] ) )
{
printf( "max: %s invalid, using %lu\n", argv[ 2 ], MAXSIEVE );
maxsieve = MAXSIEVE;
}
parameterX = 3;
}
}
finalpass = ( maxsieve + 1 ) / 2;
sieve = malloc( maxsieve );
if ( sieve == NULL )
{
printf( "Can't allocate %lu byte array, quitting.\n", maxsieve );
return EXIT_FAILURE;
}
sieve[ 0 ] = 0;
sieve[ 1 ] = 0;
for ( index = 2; index < maxsieve; index = index + 1 )
{ sieve[ index ] = 1;
}
for ( prime = 2; prime < finalpass; prime = prime + 1)
{
if ( sieve[ prime ] == 1 )
{ for ( index = prime + prime; index < maxsieve; index = index + prime )
{ sieve[ index ] = 0;
}
}
}
for ( index = 2; index < maxsieve; index = index + 1 )
{
if ( sieve[ index ] != 0 )
{
++count;
}
}
printf( "%lu primes less than %lu\n", count, maxsieve );
for ( ; parameterX < argc; ++parameterX )
{
unsigned long parameter = strtoul( argv[ parameterX ], &unscanned, 0 );
if ( ! ( unscanned > argv[ parameterX ] ) || ( parameter >= maxsieve ) )
{
printf( "parameter %s invalid\n", argv[ parameterX ] );
}
else
{
printf( "%lu is%sprime\n", parameter,
( sieve[ parameter ] == 0 ) ? " not " : " " );
}
}
return EXIT_SUCCESS;
}
Hmm. But the C program also counts.
It's not really very interactive. But ...
I should note that the Forth program does not make use of isprime? except in the querying phase. If the definition of the sieve and of isprime? somehow get modified so they don't agree, who's to know?
This kind of hidden dependency is called implicit coupling, and should be avoided in good design. Therefore, while I'm adding the counting code, I'll also redefine the sieve construction to use the same code as the query. This will slow it down just slightly, unless your Forth interpreter has some built-in leaf-flattening optimizations. You probably won't notice.
I'll also remove printPrimes, since we don't want millions of primes scrolling off the screen:
( Forth-like archetypical implementation of the sieve of eratosthenes in FORTH -- fig, bif-c -- )
( using even more of the FORTH idiom. )
( Copyright Joel Rees, 2015, 2020 )
( By Joel Rees, Amagasaki, Japan, 2015, 2020 )
( All rights reserved. )
( Permission granted by the author to use this code )
( for any purpose, )
( on condition that substantial use )
( shall retain this copyright and permission notice. )
VOCABULARY sieve-local ( Make a local symbol table. )
sieve-local DEFINITIONS ( Switch to the local vocabulary. )
30000000 CONSTANT MAXSIEVE
MAXSIEVE 1 - 2 / CONSTANT FINALPASS
0 VARIABLE sieve ( Some old FORTHs don't provide a default behavior for CREATE )
( gforth will leave the zero there. Some old FORTHs need an initial value. )
HERE sieve - DUP ( Some old FORTHs don't provide a CELL width. )
MAXSIEVE SWAP - ALLOT ( Allocate the rest of the byte array. )
CONSTANT CELLWIDTH ( To show how this can be done. )
( Arrays are handled differently by different Forths. )
( If CELLWIDTH were just to save a wasted entry in the sieve array, )
( it wouldn't be worth it here. )
( I'm assuming there are likely to be other uses.
: sieveCell ( index -- adr )
sieve + ; ( This level of indirection can slow the loop down a lot. )
: initAsPrime ( index -- )
sieveCell -1 swap C! ; ( The swap can also slow things down. )
: makeNotPrime ( index -- )
sieveCell 0 swap C! ; ( The swap can also slow things down. )
: ?prime ( index -- flag )
sieveCell C@ ;
: sieveInit ( -- )
0 makeNotPrime ( 0 is not prime. )
1 makeNotPrime ( 1 is not prime. )
MAXSIEVE 2 DO ( set flags to true for 2 to MAXSIEVE. )
I initAsPrime LOOP
;
: primePass ( prime -- )
MAXSIEVE OVER DUP + DO ( start at first multiple -- double. )
I makeNotPrime
DUP +LOOP ( next multiple )
DROP ;
: findPrimes ( -- )
FINALPASS 2 DO ( clear flags at all multiples. )
I ?prime IF ( don't bother if not prime. )
I primePass
ENDIF
LOOP ;
: countPrimes ( -- count )
0
MAXSIEVE 0 DO
I ?prime IF 1+
ENDIF
LOOP
;
FORTH DEFINITIONS
: sieveMain ( -- )
[ sieve-local ] sieveInit
findPrimes
countPrimes . ." primes less than " MAXSIEVE . ;
FORTH DEFINITIONS
: isPrime? ( index -- ) ( Say whether it's prime or not. )
[ sieve-local ]
DUP MAXSIEVE < 0= IF
." Too big!" CR
ELSE
?prime IF
." Yes."
ELSE
." No."
ENDIF
ENDIF
;
FORTH DEFINITIONS
sieveMain
Running time for the C program on a more recent and much higher spec machine than the old one I did this on before:
$ time ./queriable 20191231 20200101Comparing that to the running time of the gforth code with all the dependencies explicit (telling gforth to use a big enough dictionary):
1857859 primes less than 30000000
20191231 is prime
20200101 is not prime
real 0m1.069s
user 0m1.044s
sys 0m0.024s
$ time gforth -m 400M queriable.fs -e "20191231 isPrime? 20200101 isPrime? bye"And comparing that to the time for the counting forth version at the top of this post, with all the implicit linkage:
1857859 primes less than 30000000 Yes.No.
real 0m2.938s
user 0m2.937s
sys 0m0.000s
$ time gforth -m 400M sievecount.fs -e "20191231 isprime? 20200101 isprime? bye"Someday, I want to write a Forth capable of the leaf flattening optimizations I mentioned above.
1857859 Primes less than 30000000 Yes.No.
real 0m1.848s
user 0m1.844s
sys 0m0.004s
(Note that this current machine seems to be about 4 times as fast at primes as the machine I developed the earlier code on. The processor and memory are that much faster. But, for another measure of processor speed, see this post from way back in 2012, https://reiisi.blogspot.com/2012/08/too-late-to-complain-performance-of.html, and compare those times with this for the pi test:
time echo "scale=4000;a(1)" | bc -l
...
real 0m11.219sSo, yes, Intel did "fulfill" their roadmap. For higher spec, more expensive machines.
user 0m11.218s
sys 0m0.000s
Just for the perspective, we'll note that IBM has kept ahead of Intel's pace with their PPC and Freescale has been bought out. Wasteland.
Incidentally, this machine is from about the last batch of CPUs that doesn't have the currently trending vulnerabilities in speculative execution. Speculative execution? Heh. And they didn't expect the vulnerabilities? Well, ...)
No comments:
Post a Comment