Wolstenholme numbers: Difference between revisions
Content deleted Content added
m →{{header|Python}}: simplify output functions |
m →{{header|Phix}}: use pygments |
||
(31 intermediate revisions by 17 users not shown) | |||
Line 1:
{{
A [[wp:Wolstenholme number|Wolstenholme number]] is a number that is the fully reduced numerator of a second order harmonic number:
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}}==
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.
{{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 #
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 numerator of the nth second order harmonic number #
PROC harmonic2numerator = ( INT n )LONG LONG INT:
BEGIN
LONG LONG INT v := 0, f := 1;
FOR i TO n DO
f *:= i * i
OD;
FOR i TO 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 ] := harmonic2numerator( 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|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
Line 72 ⟶ 234:
17299975731542641
5
266681
40799043101
86364397717734821</pre>
=={{header|C}}==
{{trans|Wren}}
{{libheader|GMP}}
<syntaxhighlight lang="c">#include <stdio.h>
#include <string.h>
#include <gmp.h>
#include <locale.h>
const char *ord(unsigned long c) {
int m = c % 100;
if (m >= 4 && m <= 20) return "th";
m %= 10;
return (m == 1) ? "st" :
(m == 2) ? "nd" :
(m == 3) ? "rd" : "th";
}
void abbreviate(char a[], const char *s) {
size_t len = strlen(s);
if (len < 40) {
strcpy(a, s);
return;
}
strncpy(a, s, 20);
strcpy(a + 20, "...");
strncpy(a + 23, s + len - 20, 21);
}
int main() {
int i, pc = 0, si = 0;
unsigned long k, l[5] = {500, 1000, 2500, 5000, 10000};
char *s, a[44];
mpq_t w, h;
mpz_t n, primes[15];
mpq_init(w);
mpq_init(h);
mpz_init(n);
for (i = 0; i < 15; ++i) mpz_init(primes[i]);
printf("Wolstenholme numbers:\n");
setlocale(LC_NUMERIC, "");
for (k = 1; k <= 10000; ++k) {
mpq_set_ui(h, 1, k * k);
mpq_add(w, w, h);
mpq_get_num(n, w);
if (pc < 15 && mpz_probab_prime_p(n, 15) > 0) mpz_set(primes[pc++], n);
if (k <= 20) {
s = mpz_get_str(NULL, 10, n);
printf("%6ld%s: %s\n", k, ord(k), s);
} else if (k == l[si]) {
s = mpz_get_str(NULL, 10, n);
abbreviate(a, s);
printf("%'6ld%s: %s (digits: %ld)\n", k, ord(k), a, strlen(s));
++si;
}
}
printf("\nPrime Wolstenholme numbers:\n");
for (i = 0; i < 15; ++i) {
s = mpz_get_str(NULL, 10, primes[i]);
if (i < 4) {
printf("%6d%s: %s\n", i+1, ord(i+1), s);
} else {
abbreviate(a, s);
printf("%'6d%s: %s (digits: %ld)\n", i+1, ord(i+1), a, strlen(s));
}
}
mpq_clear(w);
mpq_clear(h);
mpz_clear(n);
for (i = 0; i < 15; ++i) mpz_clear(primes[i]);
return 0;
}</syntaxhighlight>
{{out}}
<pre>
Same as Wren example.
</pre>
=={{header|FreeBASIC}}==
<syntaxhighlight lang="vb">#include "isprime.bas"
Function GCD(n As Uinteger, d As Uinteger) As Uinteger
Return Iif(d = 0, n, GCD(d, n Mod d))
End Function
Function numArmonico(n As Uinteger, m As Uinteger) As Uinteger
Dim As Integer i, v = 0, f = 1
For i = 1 To n
f *= i ^ m
Next i
For i = 1 To n
v += f / i ^ m
Next i
v /= GCD(v, f)
Return v
End Function
Dim As Integer i, wols(12)
Print "First 12 Wolstenholme numbers:"
For i = 1 To 12
wols(i) = numArmonico (i,2)
Print wols(i)
Next i
Print !"\nTwo Wolstenholme primes:"
For i = 1 To 12
If isPrime(wols(i)) Then Print wols(i)
Next i
Sleep</syntaxhighlight>
{{out}}
<pre>First 12 Wolstenholme numbers:
1
5
49
205
5269
5369
266681
1077749
9778141
1968329
239437889
240505109
Two Wolstenholme primes:
5
266681</pre>
=={{header|Go}}==
{{trans|Wren}}
{{libheader|GMP(Go wrapper)}}
{{libheader|Go-rcu}}
<syntaxhighlight lang="go">package main
import (
"fmt"
big "github.com/ncw/gmp"
"golang.org/x/exp/constraints"
"rcu"
)
type Int = constraints.Integer
func ord[T Int](c T) string {
m := c % 100
if m >= 4 && m <= 20 {
return "th"
}
m %= 10
if m == 1 {
return "st"
} else if m == 2 {
return "nd"
} else if m == 3 {
return "rd"
} else {
return "th"
}
}
func abbreviate(s string) string {
le := len(s)
if le < 40 {
return s
}
return s[:20] + "..." + s[le-20:]
}
func main() {
pc, si := 0, 0
l := [5]int64{500, 1000, 2500, 5000, 10000}
w, h := new(big.Rat), new(big.Rat)
primes := make([]*big.Int, 15)
fmt.Println("Wolstenholme numbers:")
for k := int64(1); k <= 10000; k++ {
h.SetFrac64(1, k*k)
w.Add(w, h)
n := w.Num()
if pc < 15 && n.ProbablyPrime(15) {
primes[pc] = n
pc++
}
if k <= 20 {
fmt.Printf("%6d%s: %s\n", k, ord(k), n.String())
} else if k == l[si] {
s := n.String()
a := abbreviate(s)
ks := rcu.Commatize(k)
fmt.Printf("%6s%s: %s (digits: %d)\n", ks, ord(k), a, len(s))
si++
}
}
fmt.Println("\nPrime Wolstenholme numbers:")
for i, p := range primes {
s := p.String()
if i < 4 {
fmt.Printf("%6d%s: %s\n", i+1, ord(i+1), s)
} else {
a := abbreviate(s)
fmt.Printf("%6d%s: %s (digits: %d)\n", i+1, ord(i+1), a, len(s))
}
}
}</syntaxhighlight>
{{out}}
<pre>
Same as Wren example.
</pre>
=={{header|J}}==
<syntaxhighlight lang=J>NB. the first y members of the Wholstenhome sequence
WolstenholmeSeq=: {{ 2 {."1@x:+/\%*:1x+i.y }}
WolstenholmeSeq 20
1 5 49 205 5269 5369 266681 1077749 9778141 1968329 239437889 240505109 40799043101 40931552621 205234915681 822968714749 238357395880861 238820721143261 86364397717734821 17299975731542641
(#~ 1 p: ])WolstenholmeSeq 20
5 266681 40799043101 86364397717734821
</syntaxhighlight>
Stretch:
<syntaxhighlight lang="J"> W=: WolstenholmeSeq 1e4
#@":@> (<:500 1000 2500 5000 10000) { W
434 866 2164 4340 8693
#WP=: W#~ 1 p: W
23
#@":@> 4 5$WP
1 6 11 17 104
156 218 318 520 649
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}}==
{{trans|Python}}
<syntaxhighlight lang="julia">using Primes
""" Get N numbers in the series of Wolstenholme numbers """
wolstenholme(N) = map(numerator, accumulate(+, (1 // (i*i) for i in big"1":N)))
""" Abbreviate a large string by showing beginning / end and number of chars """
function abbreviate(s, term = "digits", thresh = 50, idx = thresh ÷ 2 - 3)
w = length(s)
return w < thresh ? s : s[begin:begin+idx] * ".." * s[end-idx:end] * " ($w $term)"
end
""" Run the tasks at rosettacode.org/wiki/Wolstenholme_numbers """
function process_wolstenholmes(nmax = 10000, filterto = 2000)
wols = wolstenholme(nmax)
println("Wolstenholme numbers 1:20, 500, 1000, 2500, 5000, 10000:")
for i in [1:20; [500, 1000, 2500, 5000, 10000]]
println(rpad(i, 5), ": ", abbreviate(string(wols[i])))
end
println("\nFifteen Wolstenholme primes:")
for (i, n) in enumerate(filter(isprime, @view wols[begin:filterto]))
println(rpad(i, 2), ": ", abbreviate(string(n)))
end
end
process_wolstenholmes()
</syntaxhighlight>{{out}}
<pre>
Wolstenholme numbers 1:20, 500, 1000, 2500, 5000, 10000:
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 : 40989667509417020364501..12248084984597965892703 (434 digits)
1000 : 83545938483149689478187..58699094240207812766449 (866 digits)
2500 : 64537911900230612090849..91212785535902976933153 (2164 digits)
5000 : 34472086597488537716198..20022525144829082590451 (4340 digits)
10000: 54714423173933343999582..37149175649667700005717 (8693 digits)
Fifteen Wolstenholme primes:
1 : 5
2 : 266681
3 : 40799043101
4 : 86364397717734821
5 : 36190908596780862323291..68379995976006474252029 (104 digits)
6 : 33427988094524601237303..12048446489305085140033 (156 digits)
7 : 22812704758392002353774..54384405125167217413149 (218 digits)
8 : 28347687473208792918550..41045794572911130248059 (318 digits)
9 : 78440559440644426017290..55630422337523878698419 (520 digits)
10: 22706893975121925531372..03702173859396183964989 (649 digits)
11: 27310394808585898968805..03886311385662644388271 (935 digits)
12: 13001072736642048751114..78008635832246554146071 (984 digits)
13: 15086863305391456002934..94405367804007944918693 (1202 digits)
14: 23541935187269979100228..81502324742766220468879 (1518 digits)
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}}==
{{libheader|ntheory}}
<syntaxhighlight lang="perl" line>
use strict;
use warnings;
use ntheory 'is_prime';
use Math::BigRat try => 'GMP';
sub abbr { my $d = shift; my $l = length $d; $l < 41 ? $d : substr($d,0,20) . '..' . substr($d,-20) . " ($l digits)" }
my @W = Math::BigRat->new('1/1');
push @W, $W[-1] + Math::BigRat->new(join '/', 1, $_**2) for 2..10000;
print "Wolstenholme numbers:\n";
printf "%5s: %s\n", $_, abbr $W[$_-1]->numerator() for 1..20, 5e2, 1e3, 2.5e3, 5e3, 1e4;
print "\nPrime Wolstenholme numbers:\n";
my($n,$c);
do { printf "%5s: %s\n", ++$c, abbr $W[$n]->numerator() if is_prime $W[++$n]->numerator() } until $c == 15;
</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:
1: 5
2: 266681
3: 40799043101
4: 86364397717734821
5: 36190908596780862323..79995976006474252029 (104 digits)
6: 33427988094524601237..48446489305085140033 (156 digits)
7: 22812704758392002353..84405125167217413149 (218 digits)
8: 28347687473208792918..45794572911130248059 (318 digits)
9: 78440559440644426017..30422337523878698419 (520 digits)
10: 22706893975121925531..02173859396183964989 (649 digits)
11: 27310394808585898968..86311385662644388271 (935 digits)
12: 13001072736642048751..08635832246554146071 (984 digits)
13: 15086863305391456002..05367804007944918693 (1202 digits)
14: 23541935187269979100..02324742766220468879 (1518 digits)
15: 40306783143871607599..58901192511859288941 (1539 digits)</pre>
=={{header|Phix}}==
{{trans|C}}
<!--(phixonline)-->
<syntaxhighlight lang="phix">
include mpfr.e
atom t0 = time(), t1 = t0+1
sequence primes = {}, l5 = {500, 1000, 2500, 5000, 10000}
mpq {h,w} = mpq_inits(2)
mpz n = mpz_init()
printf(1,"Wolstenholme numbers:\n");
for k=1 to iff(platform()=JS?1000:10000) do
mpq_set_si(h, 1, k*k)
mpq_add(w, w, h)
mpq_get_num(n, w)
-- if length(primes)<15 and mpz_prime(n) then
if length(primes)<iff(platform()=JS?8:10) and mpz_prime(n) then
primes = append(primes,mpz_get_short_str(n))
end if
if k<=20 or k=l5[1] then
printf(1,"%6d%s: %s\n", {k, ord(k), mpz_get_short_str(n)})
if k>20 then l5 = l5[2..$] end if
elsif time()>t1 then
progress("checking k=%d, length(primes)=%d...\r",{k,length(primes)})
t1 = time()+1
end if
end for
progress("")
printf(1,"\nPrime Wolstenholme numbers:\n");
for i,pw in primes do
printf(1,"%6d%s: %s\n", {i, ord(i), pw})
end for
?elapsed(time()-t0)
</syntaxhighlight>
{{out}}
Same as C, but capped at 1000/8 primes when run in a browser, to keep it under 10s.<br>
Capping it to the first 10 primes makes it 25* faster, and also keeps it under 10s on desktop/Phix.
=={{header|Python}}==
<syntaxhighlight lang="python">""" rosettacode.orgwiki/Wolstenholme_numbers """
from fractions import Fraction
from itertools import accumulate
from gmpy2 import is_prime
def wolstenholme(k):
""" Get the first k Wolstenholme numbers """
return [r.numerator
for r in accumulate((Fraction(1, i*i) for i in range(1, k+1)), lambda x, y: x+y)]
def abbreviated(wstr, thresh=60, term='digits'):
""" return an abbreviated string with beginning / end and actual number of chars """
i, wlen = max(thresh // 2, 5), len(wstr)
return wstr if wlen < thresh else ' '.join([wstr[:i], '...', wstr[-i-1:], str(wlen), term])
def process_wolstenholmes():
""" Run the tasks at rosettacode.org/wiki/Wolstenholme_numbers """
wols = wolstenholme(10000)
print('Wolstenholme numbers 1 through 20, 500, 1000, 2500, 5000, 10000:')
for i in list(range(1, 21)) + [500, 1000, 2500, 5000, 10000]:
print(f'{i:5d}: {abbreviated(str(wols[i-1]))}')
print('\nFifteen Wolstenholme primes:')
for i, num in enumerate(filter(is_prime, wols[:2000])):
print(f'{i+1:2d}: {abbreviated(str(num))}')
process_wolstenholmes()
</syntaxhighlight>{{out}}
<pre>
Wolstenholme numbers 1 through 20, 500, 1000, 2500, 5000, 10000:
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: 409896675094170203645010936925 ... 0081053712248084984597965892703 434 digits
1000: 835459384831496894781878542648 ... 2396236858699094240207812766449 866 digits
2500: 645379119002306120908498982314 ... 6779664291212785535902976933153 2164 digits
5000: 344720865974885377161982285737 ... 0347952420022525144829082590451 4340 digits
10000: 547144231739333439995821771601 ... 0312824537149175649667700005717 8693 digits
Fifteen Wolstenholme primes:
1: 5
2: 266681
3: 40799043101
4: 86364397717734821
5: 361909085967808623232911476131 ... 2280300268379995976006474252029 104 digits
6: 334279880945246012373031736295 ... 5361841612048446489305085140033 156 digits
7: 228127047583920023537742255081 ... 7225367354384405125167217413149 218 digits
8: 283476874732087929185503884347 ... 1887987141045794572911130248059 318 digits
9: 784405594406444260172907289698 ... 7380014755630422337523878698419 520 digits
10: 227068939751219255313726661483 ... 9454324603702173859396183964989 649 digits
11: 273103948085858989688057737445 ... 3536939403886311385662644388271 935 digits
12: 130010727366420487511145995120 ... 9791520578008635832246554146071 984 digits
13: 150868633053914560029346562867 ... 5398810094405367804007944918693 1202 digits
14: 235419351872699791002285017174 ... 2993590081502324742766220468879 1518 digits
15: 403067831438716075992501031099 ... 8022698298658901192511859288941 1539 digits
</pre>
Line 135 ⟶ 1,147:
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="wren">import "./gmp" for Mpq
import "./fmt" for Fmt
var w = Mpq.new()
var h = Mpq.new()
var primes = []
var l = [500, 1000, 2500, 5000, 10000]
System.print("Wolstenholme numbers:")
for (k in 1..10000) {
h.set(1, k * k)
w.add(h)
var n = w.num
if (primes.count < 15 && n.probPrime(15) > 0) primes.add(n)
if (k <= 20) {
Fmt.print("$8r: $i", k, n)
} else if (l.contains(k)) {
var ns = n.toString
Fmt.print("$,8r: $20a (digits: $d)", k, ns, ns.count)
}
}
System.print("\nPrime Wolstenholme numbers:")
for (i in 0...primes.count) {
if (i < 4) {
Fmt.print("$8r: $i", i+1, primes[i])
} else {
var ps = primes[i].toString
Fmt.print("$8r: $20a (digits: $d)", i+1, ps, ps.count)
}
}</syntaxhighlight>
{{out}}
<pre>
Wolstenholme numbers:
1st: 1
2nd: 5
3rd: 49
4th: 205
5th: 5269
6th: 5369
7th: 266681
8th: 1077749
9th: 9778141
10th: 1968329
11th: 239437889
12th: 240505109
13th: 40799043101
14th: 40931552621
15th: 205234915681
16th: 822968714749
17th: 238357395880861
18th: 238820721143261
19th: 86364397717734821
20th: 17299975731542641
500th: 40989667509417020364...48084984597965892703 (digits: 434)
1,000th: 83545938483149689478...99094240207812766449 (digits: 866)
2,500th: 64537911900230612090...12785535902976933153 (digits: 2164)
5,000th: 34472086597488537716...22525144829082590451 (digits: 4340)
10,000th: 54714423173933343999...49175649667700005717 (digits: 8693)
Prime Wolstenholme numbers:
1st: 5
2nd: 266681
3rd: 40799043101
4th: 86364397717734821
5th: 36190908596780862323...79995976006474252029 (digits: 104)
6th: 33427988094524601237...48446489305085140033 (digits: 156)
7th: 22812704758392002353...84405125167217413149 (digits: 218)
8th: 28347687473208792918...45794572911130248059 (digits: 318)
9th: 78440559440644426017...30422337523878698419 (digits: 520)
10th: 22706893975121925531...02173859396183964989 (digits: 649)
11th: 27310394808585898968...86311385662644388271 (digits: 935)
12th: 13001072736642048751...08635832246554146071 (digits: 984)
13th: 15086863305391456002...05367804007944918693 (digits: 1202)
14th: 23541935187269979100...02324742766220468879 (digits: 1518)
15th: 40306783143871607599...58901192511859288941 (digits: 1539)
</pre>
|