Wolstenholme numbers

Wolstenholme numbers
You are encouraged to solve this task according to the task description, using any language you may know.

A Wolstenholme number is a number that is the fully reduced numerator of a second order harmonic number:

$$H_{n,2}=\sum _{k=1}^{n}{\frac {1}{k^{2}}}.$$

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

• Wikipedia: Wolstenholme number
• OEIS:A007406 - Wolstenholme numbers
• OEIS:A123751 - Prime Wolstenholme numbers

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

Translation of: FreeBASIC

...but using Algol 68G's LONG LONG INT (with default precision) and an iterative gcd routine.

Uses Miller-Rabin for primality testing.

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32

Library: ALGOL 68-primes

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

Translation of: Wren

Library: GMP

#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.

EasyLang

Translation of: ALGOL 68

fastfunc isprim num .
   if num < 2
      return 0
   .
   i = 2
   while i <= sqrt num
      if num mod i = 0
         return 0
      .
      i += 1
   .
   return 1
.
func gcd a b .
   while b <> 0
      h = b
      b = a mod b
      a = h
   .
   return a
.
func wolstenholme n .
   f = 1
   for i = 2 to n
      f *= i * i
   .
   for i = 1 to n
      v += f / (i * i)
   .
   return v / gcd v f
.
for i to 11
   w = wolstenholme i
   w[] &= w
   print w
.
print ""
for w in w[]
   if isprim w = 1
      print w
   .
.

Output:

1
5
49
205
5269
5369
266681
1077749
9778141
1968329
239437889

5
266681

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

Translation of: Wren

Library: GMP(Go wrapper)

Library: Go-rcu

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

Translation of: Python

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]

Mathematica / Wolfram Language

ClearAll[wolstenhomeNumber, primeWolstenhomeNumber, indexOfNextPrimeWolstenhomeNumber, 
         ns, first15PrimeWolstenhomeNumbers, digitCount];

wolstenhomeNumber[n_Integer?Positive] := Numerator@HarmonicNumber[n, 2];
primeWolstenhomeNumber[n_Integer?Positive] := 
    wolstenhomeNumber@Nest[indexOfNextPrimeWolstenhomeNumber, 1, n];

indexOfNextPrimeWolstenhomeNumber[n_Integer?Positive] :=
    NestWhile[# + 1 &, n + 1, ! (PrimeQ@wolstenhomeNumber@#) &];

digitCount[n_Integer] := Length@IntegerDigits@n;

SetAttributes[{wolstenhomeNumber, indexOfNextPrimeWolstenhomeNumber, digitCount}, 
    Listable];

ns = {500, 1000, 2500, 5000, 10000};
first15PrimeWolstenhomeNumbers = 
    wolstenhomeNumber@NestList[indexOfNextPrimeWolstenhomeNumber, 2, 14];

Print["The first 20 Wolstenhome numbers are: ", wolstenhomeNumber[Range[20]], "."];
Print["The first four prime Wolstenhome numbers are: ", first15PrimeWolstenhomeNumbers[[1 ;; 4]], "."];
Print["The ", First@#, "th Wolstenhome number has ", Last@#, " digits."] & /@ 
    Transpose[{ns, digitCount@wolstenhomeNumber@ns}];
Print["Digit counts of the first n prime Wolstenhome numbers:"];
Print[
    TableForm[Transpose[{Range@Length@first15PrimeWolstenhomeNumbers, digitCount@first15PrimeWolstenhomeNumbers}],
    TableHeadings -> {None, {"n", "# of digits"}}]
];

Output:

The first 20 Wolstenhome numbers are: {1,5,49,205,5269,5369,266681,1077749,9778141,1968329,239437889,240505109,40799043101,40931552621,205234915681,822968714749,238357395880861,238820721143261,86364397717734821,17299975731542641}.

The first four prime Wolstenhome numbers are: {5,266681,40799043101,86364397717734821}.

The 500th Wolstenhome number has 434 digits.
The 1000th Wolstenhome number has 866 digits.
The 2500th Wolstenhome number has 2164 digits.
The 5000th Wolstenhome number has 4340 digits.
The 10000th Wolstenhome number has 8693 digits.

Digit counts of the first n prime Wolstenhome numbers:

n	# of digits
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

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

Library: bignum

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

Library: ntheory

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

Translation of: C

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

Works with: HP version 49g

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

Library: Wren-gmp

Library: Wren-fmt

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)