Prime reciprocal sum
Generate the sequence of primes where each term is the smallest prime whose reciprocal can be added to the cumulative sum and remain smaller than 1.
Prime reciprocal sum 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.
- E.G.
- The cumulative reciprocal sum with zero terms is 0. The smallest prime whose reciprocal can be added to the sum without reaching or exceeding 1 is 2, so the first term is 2 and the new cumulative reciprocal sum is 1/2.
- The smallest prime whose reciprocal can be added to the sum without reaching or exceeding 1 is 3. (2 would cause the cumulative reciprocal sum to reach 1.) So the next term is 3, and the new cumulative sum is 5/6.
- and so on...
- Task
- Find and display the first 10 terms of the sequence. (Or as many as reasonably supported by your language if it is less.) For any values with more than 40 digits, show the first and last 20 digits and the overall digit count.
- Stretch
- Find and display the next 5 terms of the sequence. (Or as many as you have the patience for.) Show only the first and last 20 digits and the overall digit count.
- See also
C
#include <stdio.h>
#include <stdbool.h>
#include <string.h>
#include <gmp.h>
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() {
mpq_t q, r, s, u;
mpz_t p, t;
int count = 0, limit = 16;
size_t len;
bool isInteger;
char *ps, a[44];
mpq_inits(q, r, s, u, NULL);
mpq_set_ui(u, 1, 1);
mpz_inits(p, t, NULL);
printf("First %d elements of the sequence:\n", limit);
while (count < limit) {
mpq_sub(q, u, s);
mpq_inv(q, q);
mpq_get_den(t, q);
isInteger = !mpz_cmp_ui(t, 1);
mpz_set_q(p, q);
if (!isInteger) mpz_add_ui(p, p, 1);
mpz_nextprime(p, p);
++count;
ps = mpz_get_str(NULL, 10, p);
len = strlen(ps);
if (len <= 40) {
printf("%2d: %s\n", count, ps);
} else {
abbreviate(a, ps);
printf("%2d: %s (digits: %ld)\n", count, a, len);
}
mpq_set_z(r, p);
mpq_inv(r, r);
mpq_add(s, s, r);
}
mpq_clears(q, r, s, u, NULL);
mpz_clears(p, t, NULL);
return 0;
}
- Output:
First 16 elements of the sequence: 1: 2 2: 3 3: 7 4: 43 5: 1811 6: 654149 7: 27082315109 8: 153694141992520880899 9: 337110658273917297268061074384231117039 10: 84241975970641143191...13803869133407474043 (digits: 76) 11: 20300753813848234767...91313959045797597991 (digits: 150) 12: 20323705381471272842...21649394434192763213 (digits: 297) 13: 12748246592672078196...20708715953110886963 (digits: 592) 14: 46749025165138838243...65355869250350888941 (digits: 1180) 15: 11390125639471674628...31060548964273180103 (digits: 2358) 16: 36961763505630520555...02467094377885929191 (digits: 4711)
J
Given:
taskfmt=: {{
if. 40<#t=. ":y do.
d=. ":#t
(20{.t),'..',(_20{.t),' (',d,' digits)'
else.
t
end.
}}@>
The task and stretch could be:
taskfmt (, 4 p:1%1-+/@:%)^:(15)i.0
2
3
7
43
1811
654149
27082315109
153694141992520880899
337110658273917297268061074384231117039
84241975970641143191..13803869133407474043 (76 digits)
20300753813848234767..91313959045797597991 (150 digits)
20323705381471272842..21649394434192763213 (297 digits)
12748246592672078196..20708715953110886963 (592 digits)
46749025165138838243..65355869250350888941 (1180 digits)
11390125639471674628..31060548964273180103 (2358 digits)
Julia
""" rosettacode.org/wiki/Prime_reciprocal_sum """
using Primes
using ResumableFunctions
""" Abbreviate a large string by showing beginning / end and number of chars """
function abbreviate(s; ds = "digits", t = 40, x = (t - 1) ÷ 2)
wid = length(s)
return wid < t ? s : s[begin:begin+x] * ".." * s[end-x:end] * " ($wid $ds)"
end
@resumable function generate_oeis75442()
psum = big"0" // big"1"
while true
n = BigInt(ceil(big"1" // (1 - psum)))
while true
n = nextprime(n + 1)
if psum + 1 // n < 1
psum += 1 // n
@yield n
break
end
end
end
end
for (i, n) in enumerate(Iterators.take(generate_oeis75442(), 17))
println(lpad(i, 2), ": ", abbreviate(string(n)))
end
- Output:
1: 2 2: 3 3: 7 4: 43 5: 1811 6: 654149 7: 27082315109 8: 153694141992520880899 9: 337110658273917297268061074384231117039 10: 84241975970641143191..13803869133407474043 (76 digits) 11: 20300753813848234767..91313959045797597991 (150 digits) 12: 20323705381471272842..21649394434192763213 (297 digits) 13: 12748246592672078196..20708715953110886963 (592 digits) 14: 46749025165138838243..65355869250350888941 (1180 digits) 15: 11390125639471674628..31060548964273180103 (2358 digits) 16: 36961763505630520555..02467094377885929191 (4711 digits) 17: 21043364724798439508..14594683820359204509 (9418 digits)
Pascal
Free Pascal
Most time consuming is finding the next prime.
program primeRezSum;
{$IFDEF FPC} {$MODE DELPHI}{$Optimization ON,ALL} {$ENDIF}
{$IFDEF WINDOWS}{$APPTYPE CONSOLE}{$ENDIF}
uses
sysutils,
gmp;
procedure OutStr(const s: AnsiString);
var
myString : AnsiString;
l : Integer;
begin
myString := copy(s,1,length(s));
l :=length(pChar(@s[1]));
setlength(myString,l);
IF l< 41 then
writeln(myString)
else
begin
myString[20]:= #0;
myString[21]:= #0;
writeln(pChar(@myString[1]),'...',pChar(@myString[l-19]),' (',l:6,' digits)');
end;
end;
var
nominator,denominator,tmp,tmpDemP,p : mpz_t;
s : AnsiString;
cnt : NativeUint;
begin
setlength(s,10000);
cnt := 1;
mpz_init(nominator);
mpz_init(tmp);
mpz_init(tmpDemP);
mpz_init_set_ui(denominator,1);
mpz_init_set_ui(p,2);
repeat
mpz_sub_ui(p,p,1);
mpz_nextprime(p,p);
write(cnt:3,' ');
repeat
mpz_mul(tmp,nominator,p);
mpz_add(tmp,tmp,denominator);
mpz_mul(tmpDemP,denominator,p);
if mpz_cmp(tmp,tmpDemP)< 0 then
BREAK;
mpz_nextprime (p,p);
until false;
mpz_get_str(pChar(@s[1]),10, p);
OutStr(s);
mpz_set(nominator,tmp);
mpz_mul(denominator,denominator,p);
//next smallest possible number denominator/delta
mpz_sub(tmp,denominator,nominator);
mpz_fdiv_q(p,denominator,tmp);
inc(cnt);
until cnt> 14;
end.
- @TIO.RUN:
1 2 2 3 3 7 4 43 5 1811 6 654149 7 27082315109 8 153694141992520880899 9 337110658273917297268061074384231117039 10 8424197597064114319...13803869133407474043 ( 76 digits) 11 2030075381384823476...91313959045797597991 ( 150 digits) 12 2032370538147127284...21649394434192763213 ( 297 digits) 13 1274824659267207819...20708715953110886963 ( 592 digits) 14 4674902516513883824...65355869250350888941 ( 1180 digits) Real time: 1.147 s User time: 1.093 s Sys. time: 0.046 s CPU share: 99.28 %
Raku
The sixteenth is very slow to emerge. Didn't bother trying for the seventeenth. OEIS only lists the first fourteen values.
sub abbr ($_) { .chars < 41 ?? $_ !! .substr(0,20) ~ '..' ~ .substr(*-20) ~ " ({.chars} digits)" }
sub next-prime {
state $sum = 0;
my $next = ((1 / (1 - $sum)).ceiling+1..*).hyper(:2batch).grep(&is-prime)[0];
$sum += FatRat.new(1,$next);
$next;
}
printf "%2d: %s\n", $_, abbr next-prime for 1..15;
- Output:
1: 2 2: 3 3: 7 4: 43 5: 1811 6: 654149 7: 27082315109 8: 153694141992520880899 9: 337110658273917297268061074384231117039 10: 84241975970641143191..13803869133407474043 (76 digits) 11: 20300753813848234767..91313959045797597991 (150 digits) 12: 20323705381471272842..21649394434192763213 (297 digits) 13: 12748246592672078196..20708715953110886963 (592 digits) 14: 46749025165138838243..65355869250350888941 (1180 digits) 15: 11390125639471674628..31060548964273180103 (2358 digits) 16: 36961763505630520555..02467094377885929191 (4711 digits)
Wren
Even with GMP takes about 4½ minutes to find the first 16.
import "./gmp" for Mpz, Mpq
import "./fmt" for Fmt
var q = Mpq.new()
var p = Mpz.new()
var r = Mpq.new()
var s = Mpq.new()
var count = 0
var limit = 16
Fmt.print("First $d elements of the sequence:", limit)
while (count < limit) {
q.set(Mpq.one.sub(s)).inv
var isInteger = (q.den == Mpz.one)
if (isInteger) p.set(q.toMpz) else p.set(q.toMpz.inc)
p.nextPrime(p)
count = count + 1
var ps = p.toString
Fmt.write("$2d: $20a", count, p)
if (ps.count > 40) {
Fmt.print(" (digits: $d)", ps.count)
} else {
System.print()
}
r.set(p).inv
s.add(r)
}- Output:
First 16 elements of the sequence: 1: 2 2: 3 3: 7 4: 43 5: 1811 6: 654149 7: 27082315109 8: 153694141992520880899 9: 337110658273917297268061074384231117039 10: 84241975970641143191...13803869133407474043 (digits: 76) 11: 20300753813848234767...91313959045797597991 (digits: 150) 12: 20323705381471272842...21649394434192763213 (digits: 297) 13: 12748246592672078196...20708715953110886963 (digits: 592) 14: 46749025165138838243...65355869250350888941 (digits: 1180) 15: 11390125639471674628...31060548964273180103 (digits: 2358) 16: 36961763505630520555...02467094377885929191 (digits: 4711)