Wolstenholme numbers
Wolstenholme numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
A Wolstenholme number is a number that is the fully reduced numerator of a second order harmonic number:
Prime Wolstenholme numbers are Wolstenholme numbers that are prime.
(Oddly enough, there is also a class of numbers called Wolstenholme primes that are not Wolstenholme numbers that are prime. This is not them.)
- Task
- Find and display the first 20 Wolstenholme numbers. (Or as many as reasonably supported by your language if it is fewer.)
- Find and display the first 4 prime Wolstenholme numbers. (Or as many as reasonably supported by your language if it is fewer.)
- Stretch
- Find and display the digit count of the 500th, 1000th, 2500th, 5000th, 10000th Wolstenholme numbers.
- Find and display the digit count of the first 15 prime Wolstenholme numbers.
- See also
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;
- Output:
1
5*
49
205
5269
5369
266681*
1077749
9778141
1968329
239437889
240505109
40799043101*
40931552621
205234915681
822968714749
238357395880861
238820721143261
86364397717734821*
17299975731542641
* Number is prime
ALGOL 68
Based on
...but using Algol 68G's LONG LONG INT (with default precision) and an iterative gcd routine.
Uses Miller-Rabin for primality testing.
Note, the source of primes.incl.a68 is on a page on Rosetta Code - see the above link.
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- Output:
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
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 &]
- Output:
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
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;
}
- Output:
Same as Wren example.
FreeBASIC
#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- Output:
First 12 Wolstenholme numbers: 1 5 49 205 5269 5369 266681 1077749 9778141 1968329 239437889 240505109 Two Wolstenholme primes: 5 266681
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))
}
}
}
- Output:
Same as Wren example.
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
Stretch:
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
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;
}
- Output:
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)
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.
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)- Output:
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]}
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()
- Output:
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)
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);
- Output:
[1,5,49,205,5269,5369,266681,1077749,9778141,1968329,239437889,240505109,40799043101,40931552621,205234915681,822968714749,238357395880861,238820721143261,86364397717734821,17299975731542641] [5,266681,40799043101,86364397717734821]
Nim
Task
Using only the standard library modules.
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]
- Output:
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
Stretch task
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])}"
- Output:
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
Perl
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;
- Output:
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)
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)
- Output:
Same as C, but capped at 1000/8 primes when run in a browser, to keep it under 10s.
Capping it to the first 10 primes makes it 25* faster, and also keeps it under 10s on desktop/Phix.
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()
- Output:
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
Raku
use Lingua::EN::Numbers;
sub abbr ($_) { .chars < 41 ?? $_ !! .substr(0,20) ~ '..' ~ .substr(*-20) ~ " (digits: {.chars})" }
my @wolstenholme = lazy ([\+] (1..∞).hyper.map: {FatRat.new: 1, .²}).map: *.numerator;
say 'Wolstenholme numbers:';
printf "%8s: %s\n", .&ordinal-digit(:c), abbr @wolstenholme[$_-1] for flat 1..20, 5e2, 1e3, 2.5e3, 5e3, 1e4;
say "\nPrime Wolstenholme numbers:";
printf "%8s: %s\n", .&ordinal-digit(:c), @wolstenholme.grep(&is-prime)[$_-1]».&abbr for 1..15;
- Output:
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)
RPL
Wolstenholme numbers displayed as negative are prime.
≪
SQ SWAP ARRY→ DROP
SWAP 3 PICK * OVER +
UNROT * DUP2 GCD
UNROT { 2 } →ARRY
SWAP /
≫ 'NEXTHARM' STO ( [ a b ] n → [ a' b' ] ) where a'/b' = a/b + 1/n²
≪ { } [ 0 1 ]
1 20 FOR n
n NEXTHARM
SWAP OVER 1 GET
IF DUP ISPRIME? THEN NEG END
+ SWAP
NEXT DROP
≫ 'TASK' STO
- Output:
1: {1 -5 49 205 5269 5369 -266681 1077749 9778141 1968329 239437889 240505109 -40799043101 40931552621 205234915681 822968714749 238357395880861 238820721143261 -86364397717734821 17299975731542641}
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
- Output:
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
Sidef
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)
}
- Output:
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)
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)
}
}
- Output:
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)