Wolstenholme numbers: Difference between revisions
→{{header|jq}}
| (14 intermediate revisions by 12 users not shown) | |||
Line 24:
;* [[oeis:A123751|OEIS:A123751 - Prime Wolstenholme numbers]]
=={{header|Ada}}==
<syntaxhighlight lang="ada">
-- Wolstenholme numbers
-- J. Carter 2023 May
with Ada.Numerics.Generic_Elementary_Functions;
with Ada.Text_IO;
with System;
procedure Wolstenholme is
type Number is mod System.Max_Binary_Modulus;
type Real is digits System.Max_Digits;
package Math is new Ada.Numerics.Generic_Elementary_Functions (Float_Type => Real);
function GCD (Left : in Number; Right : in Number) return Number;
-- Greatest Common Divisor
function Prime (Value : in Number) return Boolean;
-- Returns True if Value is prime; False otherwise
function GCD (Left : in Number; Right : in Number) return Number is
Min : Number := Number'Min (Left, Right);
Max : Number := Number'Max (Left, Right);
Remainder : Number;
begin -- GCD
Reduce : loop
if Min = 0 then
return Max;
end if;
Remainder := Max rem Min;
Max := Min;
Min := Remainder;
end loop Reduce;
end GCD;
function Prime (Value : in Number) return Boolean is
Last : constant Number := Number (Real'Truncation (Math.Sqrt (Real (Value) ) ) );
Div : Number := 3;
begin -- Prime
All_Divisors : loop
exit All_Divisors when Div > Last;
if Value rem Div = 0 then
return False;
end if;
Div := Div + 2;
end loop All_Divisors;
return True;
end Prime;
Num : Number := 1;
Den : Number := 1;
Wrk : Number;
begin -- Wolstenholme
Ada.Text_IO.Put_Line (Item => (1 .. 17 => ' ') & '1');
All_Numbers : for K in Number range 2 .. 20 loop
Wrk := K ** 2;
Num := Num * Wrk + Den;
Den := Den * Wrk;
Wrk := GCD (Num, Den);
Num := Num / Wrk;
Den := Den / Wrk;
Put : declare
Image : constant String := Num'Image;
begin -- Put
Ada.Text_IO.Put (Item => (1 .. 18 - Image'Length => ' ') & Image);
end Put;
if Prime (Num) then
Ada.Text_IO.Put (Item => '*');
end if;
Ada.Text_IO.New_Line;
end loop All_Numbers;
Ada.Text_IO.Put_Line (Item => "* Number is prime");
end Wolstenholme;
</syntaxhighlight>
{{out}}
<pre>
1
5*
49
205
5269
5369
266681*
1077749
9778141
1968329
239437889
240505109
40799043101*
40931552621
205234915681
822968714749
238357395880861
238820721143261
86364397717734821*
17299975731542641
* Number is prime
</pre>
=={{header|ALGOL 68}}==
Line 29 ⟶ 139:
Uses Miller-Rabin for primality testing.
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
{{libheader|ALGOL 68-primes}}
Note, the source of primes.incl.a68 is on a page on Rosetta Code - see the above link.
<syntaxhighlight lang="algol68">
BEGIN # find some Wolstenholme numbers and indicate which are primes #
Line 43 ⟶ 155:
a
END # gcd # ;
# returns the numerator of the nth second order harmonic number
PROC
BEGIN
LONG LONG INT v := 0, f := 1;
Line 59 ⟶ 171:
print( ( "First ", whole( UPB wols, 0 ), " Wolstenholme numbers:", newline ) );
FOR i TO UPB wols DO
wols[ i ] :=
print( ( whole( wols[ i ], -18 ) ) );
IF i MOD 4 = 0 THEN print( ( newline ) ) FI
Line 84 ⟶ 196:
... of the above, the following are prime: 5 266681 40799043101 86364397717734821
</pre>
=={{header|Arturo}}==
<syntaxhighlight lang="arturo">wolstenholme: function [n][
numerator fold 1..n [s k] -> s + to :rational @[1 k*k]
]
print "First 20 Wolstenholme numbers:"
loop 1..20 => [print wolstenholme &]
print "\nFirst 4 Wolstenholme primes:"
select.first:4 1..∞ =>
[prime? wolstenholme &] | loop => [print wolstenholme &]</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
First 4 Wolstenholme primes:
5
266681
40799043101
86364397717734821</pre>
=={{header|C}}==
Line 316 ⟶ 471:
935 984 1202 1518 1539
1770 1811 2556 2762 2769</syntaxhighlight>
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
public final class WolstenholmeNumbers {
public static void main(String[] args) {
WolstenholmeIterator iterator = new WolstenholmeIterator();
List<BigInteger> wolstenholmePrimes = new ArrayList<BigInteger>();
System.out.println("Wolstenholme numbers:");
for ( int index = 1; index <= 10_000; index++ ) {
BigInteger wolstenholmeNumber = iterator.next();
if ( wolstenholmePrimes.size() < 15 && wolstenholmeNumber.isProbablePrime(15) ) {
wolstenholmePrimes.add(wolstenholmeNumber);
}
if ( index <= 20 || PRINT_POINTS.contains(index) ) {
System.out.println(String.format("%5d%s%s", index, ": ", abbreviate(wolstenholmeNumber)));
}
}
System.out.println();
System.out.println("Prime Wolstenholme numbers:");
for ( int index = 0; index < wolstenholmePrimes.size(); index++ ) {
System.out.println(String.format("%5d%s%s", index + 1, ": ", abbreviate(wolstenholmePrimes.get(index))));
}
}
private static String abbreviate(BigInteger number) {
String text = number.toString();
int digits = text.length();
return digits <= 20 ?
text : text.substring(0, 20) + " ... " + text.substring(digits - 20) + " (digits: " + digits + ")";
}
private static final List<Integer> PRINT_POINTS = List.of( 500, 1_000, 2_500, 5_000, 10_000 );
}
final class WolstenholmeIterator {
public BigInteger next() {
BigInteger result = secondHarmonic.numerator();
k += 1;
secondHarmonic = secondHarmonic.add( new Rational(BigInteger.ONE, BigInteger.valueOf(k * k)) );
return result;
}
private int k = 1;
private Rational secondHarmonic = new Rational(BigInteger.ONE, BigInteger.ONE);
}
final class Rational {
public Rational(BigInteger aNumerator, BigInteger aDenominator) {
numerator = aNumerator;
denominator = aDenominator;
BigInteger gcd = numerator.gcd(denominator);
numerator = numerator.divide(gcd);
denominator = denominator.divide(gcd);
}
public Rational add(Rational aRational) {
BigInteger numer = numerator.multiply(aRational.denominator).add(aRational.numerator.multiply(denominator));
BigInteger denom = aRational.denominator.multiply(denominator);
return new Rational(numer, denom);
}
public BigInteger numerator() {
return numerator;
}
private BigInteger numerator;
private BigInteger denominator;
}
</syntaxhighlight>
{{ out }}
<pre>
Wolstenholme numbers:
1: 1
2: 5
3: 49
4: 205
5: 5269
6: 5369
7: 266681
8: 1077749
9: 9778141
10: 1968329
11: 239437889
12: 240505109
13: 40799043101
14: 40931552621
15: 205234915681
16: 822968714749
17: 238357395880861
18: 238820721143261
19: 86364397717734821
20: 17299975731542641
500: 40989667509417020364 ... 48084984597965892703 (digits: 434)
1000: 83545938483149689478 ... 99094240207812766449 (digits: 866)
2500: 64537911900230612090 ... 12785535902976933153 (digits: 2164)
5000: 34472086597488537716 ... 22525144829082590451 (digits: 4340)
10000: 54714423173933343999 ... 49175649667700005717 (digits: 8693)
Prime Wolstenholme numbers:
1: 5
2: 266681
3: 40799043101
4: 86364397717734821
5: 36190908596780862323 ... 79995976006474252029 (digits: 104)
6: 33427988094524601237 ... 48446489305085140033 (digits: 156)
7: 22812704758392002353 ... 84405125167217413149 (digits: 218)
8: 28347687473208792918 ... 45794572911130248059 (digits: 318)
9: 78440559440644426017 ... 30422337523878698419 (digits: 520)
10: 22706893975121925531 ... 02173859396183964989 (digits: 649)
11: 27310394808585898968 ... 86311385662644388271 (digits: 935)
12: 13001072736642048751 ... 08635832246554146071 (digits: 984)
13: 15086863305391456002 ... 05367804007944918693 (digits: 1202)
14: 23541935187269979100 ... 02324742766220468879 (digits: 1518)
15: 40306783143871607599 ... 58901192511859288941 (digits: 1539)
</pre>
=={{header|jq}}==
'''Works with gojq, the Go implementation of jq'''
gojq supports infinite-precision integer arithmetic and has no problem
computing the first 10,000 Wolstenholme numbers, but the algorithm
used here for testing primality is only suitable for identifying
the first four prime numbers in the sequence.
See [[Arithmetic/Rational#jq]] for a library of functions that supports
rational arithmetic and that is suitable for inclusion using jq's `include`
directive.
<syntaxhighlight lang="jq">
include "rational" {search: "."}; # see comment above
# Take advantage of gojq's support for infinite-precision integer arithmetic:
# If $j is 0, then an error condition is raised;
# otherwise, assuming infinite-precision integer arithmetic,
# if the input and $j are integers, then the result will be an integer.
def idivide($j):
(. % $j) as $mod
| (. - $mod) / $j ;
# To take advantage of gojq's arbitrary-precision integer arithmetic:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);
def properfactors:
. as $in
| [2, $in, false]
| recurse(
. as [$p, $q, $valid, $s]
| if $q == 1 then empty
elif $q % $p == 0 then [$p, ($q|idivide($p)), true]
elif $p == 2 then [3, $q, false, $s]
else ($s // ($q | sqrt)) as $s
| if ($p + 2) <= $s then [$p + 2, $q, false, $s]
else [$q, 1, true]
end
end )
| if .[2] and .[0] != $in then .[0] else empty end ;
def is_prime:
[limit(1; properfactors)] | length == 0;
# Use $list to specify which Wolstenholme numbers are to be displayed;
# use $primes to specify the number of prime Wolstenholme numbers to identify.
def wolstenholme($max; $list; $primes):
{primes: [], w: [], h: 0}
| foreach range (1; 1+$max) as $k (.;
.emit = null
| .h = radd(.h; r(1; $k * $k))
| .w += [.h]
| (.h | .n) as $n
| if (.primes|length) < $primes and ($n|is_prime) then .primes += [$n] end
| if $k <= 20
then .emit = [$k, $n]
elif ($k | IN($list[]))
then .emit ="\($k): \($n|tostring[0:20]) (digits: \($n|tostring|length))"
else .
end;
select(.emit).emit,
(if $k == $max then {primes} else empty end) );
"Wolstenholme numbers:", wolstenholme(10000; [500, 1000, 2500, 5000, 10000]; 4)
</syntaxhighlight>
{{output}}
<pre>
Wolstenholme numbers:
[1,1]
[2,5]
[3,49]
[4,205]
[5,5269]
[6,5369]
[7,266681]
[8,1077749]
[9,9778141]
[10,1968329]
[11,239437889]
[12,240505109]
[13,40799043101]
[14,40931552621]
[15,205234915681]
[16,822968714749]
[17,238357395880861]
[18,238820721143261]
[19,86364397717734821]
[20,17299975731542641]
500: 40989667509417020364 (digits: 434)
1000: 83545938483149689478 (digits: 866)
2500: 64537911900230612090 (digits: 2164)
5000: 34472086597488537716 (digits: 4340)
10000: 54714423173933343999 (digits: 8693)
{"primes":[1,5,266681,40799043101]}
</pre>
=={{header|Julia}}==
Line 390 ⟶ 772:
15: 40306783143871607599250..98658901192511859288941 (1539 digits)
</pre>
=={{header|Maxima}}==
<syntaxhighlight lang="maxima">
wolstenholme(n):=num(apply("+",1/makelist(i^2,i,n)))$
/* Test cases */
makelist(wolstenholme(i),i,20);
/* The previous list is enough to find the first 4 Wolstenholme primes */
sublist(%,primep);
</syntaxhighlight>
{{out}}
<pre>
[1,5,49,205,5269,5369,266681,1077749,9778141,1968329,239437889,240505109,40799043101,40931552621,205234915681,822968714749,238357395880861,238820721143261,86364397717734821,17299975731542641]
[5,266681,40799043101,86364397717734821]
</pre>
=={{header|Nim}}==
===Task===
Using only the standard library modules.
<syntaxhighlight lang="Nim">import std/rationals
func isPrime(n: Natural): bool =
if n < 2: return false
if (n and 1) == 0: return n == 2
if n mod 3 == 0: return n == 3
var k = 5
while k * k <= n:
if n mod k == 0 or n mod (k + 2) == 0: return false
inc k, 6
result = true
iterator wolstenholme(): (int, int) =
var count = 0
var k = 1
var s = 1.toRational
while true:
inc count
yield (count, s.num)
inc k
s += 1 // (k * k)
echo "First 20 Wolstenholme numbers:"
var wprimes: seq[int]
for (count, n) in wolstenholme():
echo n
if n.isPrime:
wprimes.add n
if count == 20: break
echo "\nFirst 4 Wolstenholme primes:"
for i in 0..3:
echo wprimes[i]
</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
First 4 Wolstenholme primes:
5
266681
40799043101
86364397717734821</pre>
===Stretch task===
{{libheader|bignum}}
<syntaxhighlight lang="Nim">import std/strformat
import bignum
iterator wolstenholme(): (int, Int) =
var count = 0
var k = 1
var s = newRat(1)
while true:
inc count
yield (count, s.num)
inc k
s += newRat(1, k * k)
var wprimes: seq[Int]
for (count, n) in wolstenholme():
if wprimes.len < 15 and probablyPrime(n, 25) != 0:
wprimes.add n
if count in [500, 1000, 2500, 5000, 10000]:
echo &"The {count}th Wolstenholme number has {len($n)} digits."
if count == 10000: break
echo "\nDigit count of the first 15 Wolstenholme primes:"
for i in 0..14:
echo &"{i + 1: >2}: {len($wprimes[i])}"
</syntaxhighlight>
{{out}}
<pre>The 500th Wolstenholme number has 434 digits.
The 1000th Wolstenholme number has 866 digits.
The 2500th Wolstenholme number has 2164 digits.
The 5000th Wolstenholme number has 4340 digits.
The 10000th Wolstenholme number has 8693 digits.
Digit count of the first 15 Wolstenholme primes:
1: 1
2: 6
3: 11
4: 17
5: 104
6: 156
7: 218
8: 318
9: 520
10: 649
11: 935
12: 984
13: 1202
14: 1518
15: 1539</pre>
=={{header|Perl}}==
Line 399 ⟶ 917:
use Math::BigRat try => 'GMP';
sub abbr
my @W = Math::BigRat->new('1/1');
Line 407 ⟶ 925:
printf "%5s: %s\n", $_, abbr $W[$_-1]->numerator() for 1..20, 5e2, 1e3, 2.5e3, 5e3, 1e4;
my($n,$c
do { printf "%5s: %s\n", ++$c, abbr $W[$n]->numerator() if is_prime $W[
</syntaxhighlight>
{{out}}
Line 441 ⟶ 959:
Prime Wolstenholme numbers:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
=={{header|Phix}}==
Line 629 ⟶ 1,146:
14th: 23541935187269979100..02324742766220468879 (digits: 1518)
15th: 40306783143871607599..58901192511859288941 (digits: 1539)</pre>
=={{header|RPL}}==
{{works with|HP|49g}}
Wolstenholme numbers displayed as negative are prime.
≪
SQ SWAP ARRY→ DROP
SWAP 3 PICK * OVER +
UNROT * DUP2 GCD
UNROT { 2 } →ARRY
SWAP /
≫ '<span style="color:blue">NEXTHARM</span>' STO <span style="color:grey">''( [ a b ] n → [ a' b' ] )'' where a'/b' = a/b + 1/n²</span>
≪ { } [ 0 1 ]
1 20 '''FOR''' n
n <span style="color:blue">NEXTHARM</span>
SWAP OVER 1 GET
'''IF''' DUP ISPRIME? '''THEN''' NEG '''END'''
+ SWAP
'''NEXT''' DROP
≫ '<span style="color:blue">TASK</span>' STO
{{out}}
<pre>
1: {1 -5 49 205 5269 5369 -266681 1077749 9778141 1968329 239437889 240505109 -40799043101 40931552621 205234915681 822968714749 238357395880861 238820721143261 -86364397717734821 17299975731542641}
</pre>
=={{header|Ruby}}==
<syntaxhighlight lang="ruby">require 'prime'
res = (1..20).each.inject([0]){|memo, n| memo << memo.last + (1r/(n*n))}
wolstenholmes = res.map(&:numerator)[1..]
wolstenholmes_primes = wolstenholmes.select(&:prime?)
[wolstenholmes, wolstenholmes_primes].each do |whs|
whs.each.with_index(1){|n, i| puts "%-3d: %d" % [i, n] }
puts
end
</syntaxhighlight>
{{out}}
<pre>1 : 1
2 : 5
3 : 49
4 : 205
5 : 5269
6 : 5369
7 : 266681
8 : 1077749
9 : 9778141
10 : 1968329
11 : 239437889
12 : 240505109
13 : 40799043101
14 : 40931552621
15 : 205234915681
16 : 822968714749
17 : 238357395880861
18 : 238820721143261
19 : 86364397717734821
20 : 17299975731542641
1 : 5
2 : 266681
3 : 40799043101
4 : 86364397717734821
</pre>
=={{header|Sidef}}==
<syntaxhighlight lang="ruby">say "Wolstenholme numbers:"
for n in (1..20) {
say ("#{'%5s'%n}: ", sum(1..n, {|k| 1/k**2 }).nu)
}
for n in (500, 1000, 2500, 5000, 10000) {
var w = Str(sum(1..n, {|k| 1/k**2 }).nu)
say ("#{'%5s'%n}: ", w.first(20), '...', w.last(20), " (#{w.len} digits)")
}
say "\nPrime Wolstenholme numbers:"
Enumerator({|callback|
var w = 0
for k in (1..Inf) {
w += 1/k**2
if (w.nu.is_prime) {
callback([k, w.nu])
}
}
}).first(15).each_2d {|n,p|
var w = Str(p)
if (w.len > 50) {
w = join('', w.first(20), '...', w.last(20), " (#{w.len} digits)")
}
say ("#{'%5s'%n}: ", w)
}</syntaxhighlight>
{{out}}
<pre>
Wolstenholme numbers:
1: 1
2: 5
3: 49
4: 205
5: 5269
6: 5369
7: 266681
8: 1077749
9: 9778141
10: 1968329
11: 239437889
12: 240505109
13: 40799043101
14: 40931552621
15: 205234915681
16: 822968714749
17: 238357395880861
18: 238820721143261
19: 86364397717734821
20: 17299975731542641
500: 40989667509417020364...48084984597965892703 (434 digits)
1000: 83545938483149689478...99094240207812766449 (866 digits)
2500: 64537911900230612090...12785535902976933153 (2164 digits)
5000: 34472086597488537716...22525144829082590451 (4340 digits)
10000: 54714423173933343999...49175649667700005717 (8693 digits)
Prime Wolstenholme numbers:
2: 5
7: 266681
13: 40799043101
19: 86364397717734821
121: 36190908596780862323...79995976006474252029 (104 digits)
188: 33427988094524601237...48446489305085140033 (156 digits)
252: 22812704758392002353...84405125167217413149 (218 digits)
368: 28347687473208792918...45794572911130248059 (318 digits)
605: 78440559440644426017...30422337523878698419 (520 digits)
745: 22706893975121925531...02173859396183964989 (649 digits)
1085: 27310394808585898968...86311385662644388271 (935 digits)
1127: 13001072736642048751...08635832246554146071 (984 digits)
1406: 15086863305391456002...05367804007944918693 (1202 digits)
1743: 23541935187269979100...02324742766220468879 (1518 digits)
1774: 40306783143871607599...58901192511859288941 (1539 digits)
</pre>
=={{header|Wren}}==
{{libheader|Wren-gmp}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="
import "./fmt" for Fmt
| |||