Wolstenholme numbers: Difference between revisions
Added Algol 68
(→{{header|J}}: stretch) |
(Added Algol 68) |
||
Line 24:
;* [[oeis:A123751|OEIS:A123751 - Prime Wolstenholme numbers]]
=={{header|ALGOL 68}}==
Based on {{Trans|FreeBASIC}}...but using Algol 68G's LONG LONG INT (with default precision) and an iterative gcd routine.<br>
Uses Miller-Rabin for primality testing, shows the first 20 numbers.
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
<syntaxhighlight lang="algol68">
BEGIN # find some Wolstenholme numbers and indicate which are primes #
PR read "primes.incl.a68" PR # include prime utilities #
# iterative Greatest Common Divisor routine, returns the gcd of m and n #
PROC gcd = ( LONG LONG INT m, n )LONG LONG INT:
BEGIN
LONG LONG INT a := ABS m, b := ABS n;
WHILE b /= 0 DO
LONG LONG INT new a = b;
b := a MOD b;
a := new a
OD;
a
END # gcd # ;
# returns the nth second order harmonic number #
PROC harmonic2 = ( LONG LONG INT n )LONG LONG INT:
BEGIN
LONG LONG INT v := 0, f := 1;
LONG LONG INT i;
i := 0;
WHILE ( i +:= 1 ) <= n DO
f *:= i * i
OD;
i := 0;
WHILE ( i +:= 1 ) <= n DO
v +:= f OVER ( i * i )
OD;
v OVER gcd( v, f )
END # harmonic2 # ;
[ 20 ]LONG LONG INT wols;
print( ( "First ", whole( UPB wols, 0 ), " Wolstenholme numbers:", newline ) );
FOR i TO UPB wols DO
wols[ i ] := harmonic2( i );
print( ( whole( wols[ i ], -18 ) ) );
IF i MOD 4 = 0 THEN print( ( newline ) ) FI
OD;
print( ( newline, "... of the above, the following are prime:" ) );
FOR i TO UPB wols DO
IF is probably prime( wols[ i ] ) THEN
print( ( " ", whole( wols[ i ], 0 ) ) )
FI
OD;
print( ( newline ) )
END
</syntaxhighlight>
{{out}}
<pre>
First 20 Wolstenholme numbers:
1 5 49 205
5269 5369 266681 1077749
9778141 1968329 239437889 240505109
40799043101 40931552621 205234915681 822968714749
238357395880861 238820721143261 86364397717734821 17299975731542641
... of the above, the following are prime: 5 266681 40799043101 86364397717734821
</pre>
=={{header|C}}==
| |||