Rosetta Code

Wolstenholme numbers

Revision as of 21:41, 3 May 2023 by Tigerofdarkness (talk | contribs) (→‎{{header|ALGOL 68}}: Tweak)

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

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.

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


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
BEGIN # find some Wolstenholme numbers and indicate which are primes         #
    PR read "primes.incl.a68" PR                   # include prime utilities #
    # iterative Greatest Common Divisor routine, returns the gcd of m and n  #
    PROC gcd = ( LONG LONG INT m, n )LONG LONG INT:
         BEGIN
            LONG LONG INT a := ABS m, b := ABS n;
            WHILE b /= 0 DO
                LONG LONG INT new a = b;
                b        := a MOD b;
                a        := new a
            OD;
            a
        END # gcd # ;
    # returns the nth second order harmonic number                           #
    PROC harmonic2 = ( 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 ] := harmonic2( i );
        print( ( whole( wols[ i ], -18 ) ) );
        IF i MOD 4 = 0 THEN print( ( newline ) ) FI
    OD;
    print( ( newline, "... of the above, the following are prime:" ) );
    FOR i TO UPB wols DO
        IF is probably prime( wols[ i ] ) THEN
            print( ( " ", whole( wols[ i ], 0 ) ) )
        FI
    OD;
    print( ( newline ) )

END
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

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.

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

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 sympy import isprime


def wolstenholme(k):
    """ Wolstenholme numbers """
    return sum(Fraction(1, i*i) for i in range(1, k+1)).numerator

def process_wolstenholmes(nmax):
    """ Run the tasks at rosettacode.org/wiki/Wolstenholme_numbers """
    wols = [wolstenholme(n) for n in range(1, nmax+1)]
    print('First 20 Wolstenholme numbers:')
    for i in wols[:20]:
        print(i)
    print('\nFour Wolstenholme primes:')
    for j in filter(isprime, wols):
        print(j)


process_wolstenholmes(20)
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

Four Wolstenholme primes:
5
266681
40799043101
86364397717734821

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)

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)
Retrieved from "https://rosettacode.org/wiki/Wolstenholme_numbers?oldid=341498"