RosettaCodeData/Task/Blum-integer/ALGOL-68/blum-integer.alg

91 lines
3.7 KiB
Plaintext

BEGIN # find Blum integers, semi-primes where both factors are 3 mod 4 #
# and the factors are distinct #
INT max number = 10 000 000; # maximum possible Blum we will consider #
[ 1 : max number ]INT upfc; # table of unique prime factor counts #
[ 1 : max number ]INT lpf; # table of last prime factors #
FOR i TO UPB upfc DO upfc[ i ] := 0; lpf[ i ] := 1 OD;
FOR i FROM 2 TO UPB upfc OVER 2 DO
IF upfc[ i ] = 0 THEN
FOR j FROM i + i BY i TO UPB upfc DO
upfc[ j ] +:= 1;
lpf[ j ] := i
OD
FI
OD;
# returns TRUE if n is a Blum integer, FALSE otherwise #
PROC is blum = ( INT n )BOOL:
IF n < 21 THEN FALSE # the lowest primes modulo 4 that are 3 are #
# 3 & 7, so 21 is the lowest possible number #
ELIF NOT ODD n THEN FALSE
ELSE
INT v := n;
INT f count := 0;
INT f := 3;
WHILE f <= v AND f count < 3 DO
IF v MOD f = 0 THEN
IF f MOD 4 /= 3 THEN
f count := 9 # the prime factor mod 4 isn't 3 #
ELSE
v OVERAB f;
f count +:= 1
FI;
IF v MOD f = 0 THEN
f count := 9 # f is not a distinct factor of n #
FI
FI;
f +:= 2
OD;
f count = 2
FI # is blum # ;
# show various Blum integers #
print( ( "The first 50 Blum integers:", newline ) );
INT b count := 0;
[ 0 : 9 ]INT last count := []INT( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 )[ AT 0 ];
FOR i FROM 21 BY 4 WHILE b count < 400 000 DO
# Blum integers must be 1 modulo 4 #
IF b count < 50 THEN
IF is blum( i ) THEN
b count +:= 1;
last count[ i MOD 10 ] +:= 1;
print( ( whole( i, -4 ) ) );
IF b count MOD 10 = 0 THEN print( ( newline ) ) FI
FI
ELIF upfc[ i ] = 2 THEN
# two prime factors - could be a Blum integer #
IF lpf[ i ] MOD 4 = 3 THEN
# the last prime factor mod 4 is three #
IF upfc[ i OVER lpf[ i ] ] = 0 THEN
# and the other prime factor is unique (e.g. not 3^3) #
b count +:= 1;
last count[ i MOD 10 ] +:= 1;
IF b count = 26 828
OR b count = 100 000
OR b count = 200 000
OR b count = 300 000
OR b count = 400 000
THEN
print( ( "The ", whole( b count, -6 )
, "th Blum integer is ", whole( i, -11 )
, newline
)
)
FI
FI
FI
FI
OD;
# show some statistics of the last digits #
print( ( "For Blum integers up to ", whole( b count, 0 ), ":", newline ) );
FOR i FROM LWB last count TO UPB last count DO
IF last count[ i ] /= 0 THEN
print( ( " ", fixed( ( last count[ i ] * 100 ) / b count, -8, 3 )
, "% end with ", whole( i, 0 )
, newline
)
)
FI
OD
END