Prime reciprocal sum

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.

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.

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

• OEIS:A075442 - Slowest-growing sequence of primes whose reciprocals sum to 1

C

Library: GMP

Translation of: Wren

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

Here, +/@:% is the sum of reciprocals, so 1%1-+/@:% is the reciprocal of the amount remaining, and 4 p:1%1-+/@:% is the smallest prime which is larger than that value.

Tested in J9.4

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.
Now pre-sieving with the primes til 65535, to reduce tests.
That is the same as checking all numbers in sieve as divisible by the small primes.Since the prime are in big distance in that region, that's an improvement.

program primeRezSum;
{$IFDEF FPC}  {$MODE DELPHI}{$Optimization ON,ALL} {$ENDIF}
{$IFDEF WINDOWS}{$APPTYPE CONSOLE}{$ENDIF}
uses
  sysutils,
  gmp;
const
  PrimeCount =   6542;
  PRIMEMAXVAL = 65535;
  SIEVESIZE = 65536;
type
  Tsieveprime = record
                  prime,
                  offset : Uint16;
                end;
  tSievePrimes = array[0..PrimeCount-1] of Tsieveprime;
  tSieve       = array[0..SIEVESIZE-1] of byte;
var
  s : AnsiString;
  MyPrimes : tSievePrimes;
  sieve : tSieve;
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;

function InitPrimes:Uint32;
var
  f : extended;
  idx,p,pr_cnt : Uint32;
Begin
  fillchar(sieve,Sizeof(sieve),#0);
  pr_cnt := 0;
  p := 2;
  f := 1.0;
  repeat
     while Sieve[p]<> 0 do
       inc(p);
     MyPrimes[pr_cnt].prime := p;
     f := f*(p-1)/p;
     inc(pr_cnt);
     idx := p*p;
     if idx > PRIMEMAXVAL then
       Break;
     repeat
       sieve[idx] := 1;
       inc(idx,p);
     until idx > high(sieve);
     inc(p);
  until sqr(p)>PRIMEMAXVAL;

  while (pr_cnt<= High(MyPrimes)) AND (p<PRIMEMAXVAL)  do
  begin
    while Sieve[p]<> 0 do
      inc(p);
    MyPrimes[pr_cnt].prime := p;
    f := f*(p-1)/p;
    inc(p);
    inc(pr_cnt);
  end;
  Writeln ('reducing factor ',f:10:8);
  result := pr_cnt-1;
end;

procedure DoSieveOffsetInit(var tmp:mpz_t);
var
  dummy :mpz_t;
  i,j,p : Uint32;
Begin
  mpz_init(dummy);
  for i := 0 to High(MyPrimes) do
    with MyPrimes[i] do
    Begin
      if prime = 0 then  Begin  writeln(i);halt;end;
      offset := prime-mpz_tdiv_q_ui(dummy,tmp,prime);
    end;
  mpz_set(dummy,tmp);
  repeat
    //one sieve
    fillchar(sieve,Sizeof(sieve),#0);
    //sieve
    For i := 0 to High(MyPrimes) do
    begin
      with MyPrimes[i] do
      begin
        p := prime;
        j := offset;
      end;
      repeat
        sieve[j] := 1;
        j += p;
      until j >= SIEVESIZE;
      MyPrimes[i].offset := j-SIEVESIZE;
    end;

    j := 0;
    For i := 0 to High(sieve) do
    begin
//  function mpz_probab_prime_p(var n: mpz_t; reps: longint): longint;1 = prime
      if (sieve[i]= 0) then
      begin
        mpz_add_ui(dummy,dummy,i-j);
        j := i;
        if (mpz_probab_prime_p(dummy,1) >0)  then
        begin
          mpz_set(tmp,dummy);
          mpz_clear(dummy);
          EXIT;
        end;
      end;
    end;
  until false;
end;

var
  nominator,denominator,tmp,tmpDemP,p : mpz_t;
  T1,T0:Int64;
  cnt : NativeUint;
begin
  InitPrimes;
  setlength(s,100000);
  cnt := 1;
  mpz_init(nominator);
  mpz_init(tmp);
  mpz_init(tmpDemP);
  mpz_init_set_ui(denominator,1);
  mpz_init_set_ui(p,1);

  repeat
    mpz_set(tmpDemP,p);
    T0 := GetTickCount64;
    if cnt > 9 then
      DoSieveOffsetInit(p)
    else
      mpz_nextprime(p,p);
    T1 := GetTickCount64;
    write(cnt:3,'  ',T1-T0,' ms ,delta ');
    mpz_sub(tmpDemP,p,tmpDemP);
    mpz_get_str(pChar(@s[1]),10, tmpDemP); OutStr(s);
    mpz_get_str(pChar(@s[1]),10, p); OutStr(s);
    if cnt >=15 then
      Break;
    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_add_ui(p,p,1);
    until false;
    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> 17;
end.

@TIO.RUN:

reducing factor 0.05041709
  1  0 ms ,delta 1
2
  2  0 ms ,delta 1
3
  3  0 ms ,delta 1
7
  4  0 ms ,delta 1
43
  5  0 ms ,delta 5
1811
  6  0 ms ,delta 16
654149
  7  0 ms ,delta 5
27082315109
  8  0 ms ,delta 71
153694141992520880899
  9  1 ms ,delta 14
337110658273917297268061074384231117039
 10  1 ms ,delta 350
8424197597064114319...13803869133407474043 (    76 digits)
 11  1 ms ,delta 203
2030075381384823476...91313959045797597991 (   150 digits)
 12  3 ms ,delta 33
2032370538147127284...21649394434192763213 (   297 digits)
 13  69 ms ,delta 348
1274824659267207819...20708715953110886963 (   592 digits)
 14  234 ms ,delta 192
4674902516513883824...65355869250350888941 (  1180 digits)
 15  20391 ms ,delta 3510
1139012563947167462...31060548964273180103 (  2358 digits)

Real time: 21.024 s User time: 20.768 s Sys. time: 0.067 s CPU share: 99.10 %

Python

''' rosettacode.org/wiki/Prime_reciprocal_sum '''
from fractions import Fraction
from sympy import nextprime


def abbreviated(bigstr, label="digits", maxlen=40, j=20):
    ''' Abbreviate string by showing beginning / end and number of chars '''
    wid = len(bigstr)
    return bigstr if wid <= maxlen else bigstr[:j] + '..' + bigstr[-j:] + f' ({wid} {label})'


def ceil(rat):
    ''' ceil function for Fractions '''
    return rat.numerator if rat.denominator == 1 else rat.numerator // rat.denominator + 1


psum = Fraction(0, 1)
for i in range(1, 15):  # get first 14 in sequence
    next_in_seq = nextprime(ceil(Fraction(1, 1 - psum)))
    psum += Fraction(1, next_in_seq)
    print(f'{i:2}:', abbreviated(str(next_in_seq)))

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)

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

Library: Wren-gmp

Library: Wren-fmt

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)