Wolstenholme numbers: Difference between revisions
Content deleted Content added
Thundergnat (talk | contribs) m →{{header|Perl}}: Also, won't run under 5.30 |
SqrtNegInf (talk | contribs) →{{header|Perl}}: now runs correctly, and on older Perl |
||
Line 392:
=={{header|Perl}}==
{{libheader|ntheory}}
<syntaxhighlight lang="perl" line>
Line 400 ⟶ 399:
use Math::BigRat try => 'GMP';
sub abbr
my @W = Math::BigRat->new('1/1');
Line 409 ⟶ 408:
push @res, "\nPrime Wolstenholme numbers:";
my($n,$c
do { printf "%5s: %s\n", ++$c, abbr $W[$n]->numerator() if is_prime $W[
</syntaxhighlight>
{{out}}
Line 442 ⟶ 441:
Prime Wolstenholme numbers:
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
=={{header|Phix}}==
|
Revision as of 12:05, 4 May 2023
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
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.
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
#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
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)
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;
push @res, "\nPrime Wolstenholme numbers:";
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)
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)