Prime reciprocal sum: Difference between revisions
Added Algol 68 |
Added Sidef |
||
(9 intermediate revisions by 6 users not shown) | |||
Line 39: | Line 39: | ||
# previous primes and the sum remain below 1 # |
# previous primes and the sum remain below 1 # |
||
PR precision |
PR precision 5000 PR # set the precision of LONG LONG INT # |
||
PR read "primes.incl.a68" PR # include prime utilities # |
PR read "primes.incl.a68" PR # include prime utilities # |
||
Line 87: | Line 87: | ||
# each succeeding prime is the next prime > the denominator of the sum # |
# each succeeding prime is the next prime > the denominator of the sum # |
||
# divided by the difference of the denominator and the numerator # |
# divided by the difference of the denominator and the numerator # |
||
FOR i FROM 2 TO |
FOR i FROM 2 TO 15 DO |
||
LONG LONG INT next := IF num OF sum + 1 = den OF sum |
LONG LONG INT next := IF num OF sum + 1 = den OF sum |
||
THEN den OF sum + 1 |
THEN den OF sum + 1 |
||
Line 129: | Line 129: | ||
13: 12748246592672078196...20708715953110886963 592 digits |
13: 12748246592672078196...20708715953110886963 592 digits |
||
14: 46749025165138838243...65355869250350888941 1180 digits |
14: 46749025165138838243...65355869250350888941 1180 digits |
||
15: 11390125639471674628...31060548964273180103 2358 digits |
|||
</pre> |
</pre> |
||
Line 294: | Line 295: | ||
16: 36961763505630520555..02467094377885929191 (4711 digits) |
16: 36961763505630520555..02467094377885929191 (4711 digits) |
||
17: 21043364724798439508..14594683820359204509 (9418 digits) |
17: 21043364724798439508..14594683820359204509 (9418 digits) |
||
</pre> |
|||
=={{header|Nim}}== |
|||
{{libheader|bignum}} |
|||
<syntaxhighlight lang="Nim">import std/strformat |
|||
import bignum |
|||
iterator a075442(): Int = |
|||
let One = newRat(1) |
|||
var sum = newRat(0) |
|||
var p = newInt(0) |
|||
while true: |
|||
let q = reciprocal(One - sum) |
|||
p = nextPrime(if q.isInt: q.num else: q.toInt + 1) |
|||
yield p |
|||
sum += newRat(1, p) |
|||
func compressed(str: string; size: int): string = |
|||
## Return a compressed value for long strings of digits. |
|||
if str.len <= 2 * size: str |
|||
else: &"{str[0..<size]}...{str[^size..^1]} ({str.len} digits)" |
|||
var count = 0 |
|||
for p in a075442(): |
|||
inc count |
|||
echo &"{count:2}: {compressed($p, 20)}" |
|||
if count == 15: break |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> 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) |
|||
</pre> |
</pre> |
||
Line 510: | Line 556: | ||
Real time: 21.024 s User time: 20.768 s Sys. time: 0.067 s CPU share: 99.10 % |
Real time: 21.024 s User time: 20.768 s Sys. time: 0.067 s CPU share: 99.10 % |
||
@home (4.4Ghz 5600G ) modified with primes up to 1E6 ( Uint16-> Uint32 ) |
|||
78498 999979 reducing factor 0.04059798 |
|||
... |
|||
15 10015 ms ,delta 3510 // first guess than searching for next prime |
|||
1139012563947167462...31060548964273180103 ( 2358 digits) |
|||
1731172 999979 |
|||
16 110030 ms ,delta 6493 |
|||
3696176350563052055...02467094377885929191 ( 4711 digits) |
|||
17 5098268 ms ,delta 55552 // first guess than searching for next prime |
|||
2104336472479843950...14594683820359204509 ( 9418 digits) |
|||
</pre> |
|||
=={{header|Perl}}== |
|||
<syntaxhighlight lang="perl" line> |
|||
use strict; use warnings; use feature 'state'; |
|||
use Math::AnyNum <next_prime ceil>; |
|||
sub abbr { my $d=shift; my $l = length $d; $l < 41 ? $d : substr($d,0,20) . '..' . substr($d,-20) . " ($l digits)" } |
|||
sub succ_prime { |
|||
state $sum = 0; |
|||
my $next = next_prime ceil( 1 / (1-$sum) ); |
|||
$sum += 1/$next; |
|||
$next |
|||
} |
|||
printf "%2d: %s\n", $_, abbr succ_prime for 1..14; |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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) |
|||
</pre> |
|||
=={{header|Phix}}== |
|||
The Julia entry took over 4 hours to get the 17th on this box, so I just went with "keep it all under 10s"...<br> |
|||
<small>(As others have alluded, it is all about getting the next prime. Even should you land on one straightaway, it still takes quite some time to prove an n-thousand digit number ''is'' prime.)</small> |
|||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">limit</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()=</span><span style="color: #004600;">JS</span><span style="color: #0000FF;">?</span><span style="color: #000000;">12</span><span style="color: #0000FF;">:</span><span style="color: #000000;">14</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> |
|||
<span style="color: #004080;">mpq</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">u</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpq_inits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">})</span> |
|||
<span style="color: #004080;">mpz</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_inits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"First %d elements of the sequence:\n"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">limit</span><span style="color: #0000FF;">);</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">count</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">limit</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #7060A8;">mpq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">u</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpq_inv</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">q</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpq_get_den</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">q</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpz_set_q</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">q</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #7060A8;">mpz_add_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
<span style="color: #7060A8;">mpz_nextprime</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%2d: %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">count</span><span style="color: #0000FF;">,</span> <span style="color: #7060A8;">mpz_get_short_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)})</span> |
|||
<span style="color: #7060A8;">mpq_set_z</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpq_inv</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #7060A8;">mpq_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
|||
<pre> |
|||
First 14 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 (76 digits) |
|||
11: 20300753813848234767...91313959045797597991 (150 digits) |
|||
12: 20323705381471272842...21649394434192763213 (297 digits) |
|||
13: 12748246592672078196...20708715953110886963 (592 digits) |
|||
14: 46749025165138838243...65355869250350888941 (1,180 digits) |
|||
</pre> |
|||
If you've really got nothing better to do, you can also get, in 10 mins on 64-bit, thrice that on 32-bit: |
|||
<pre> |
|||
15: 11390125639471674628...31060548964273180103 (digits: 2358) |
|||
16: 36961763505630520555...02467094377885929191 (digits: 4711) |
|||
</pre> |
</pre> |
||
Line 581: | Line 719: | ||
15: 11390125639471674628..31060548964273180103 (2358 digits) |
15: 11390125639471674628..31060548964273180103 (2358 digits) |
||
16: 36961763505630520555..02467094377885929191 (4711 digits)</pre> |
16: 36961763505630520555..02467094377885929191 (4711 digits)</pre> |
||
=={{header|RPL}}== |
|||
Limited floating-point precision prevents from finding the correct 7th term. Program starts with a non-empty sequence to avoid the <code>∑LIST</code> calculation bug that occurs when the input list has less than two items. |
|||
{{works with|HP|49g}} |
|||
≪ {2 3} |
|||
1 4 '''START''' 1 OVER INV ∑LIST - INV →NUM IP NEXTPRIME + '''NEXT''' |
|||
≫ '<span style="color:blue">∑INVPR</span>' STO |
|||
{{out}} |
|||
<pre> |
|||
1: {2 3 7 43 1811 654149} |
|||
</pre> |
|||
=={{header|Sidef}}== |
|||
<syntaxhighlight lang="ruby">var A075442 = Enumerator({|callback| |
|||
var sum = 0 |
|||
loop { |
|||
var p = next_prime(ceil(1/(1-sum))) |
|||
sum += 1/p |
|||
callback(p) |
|||
} |
|||
}) |
|||
A075442.first(15).each_kv {|k,n| |
|||
var s = Str(n) |
|||
s = "#{s.first(20)}..#{s.last(20)} (#{s.len} digits)" if (s.len > 50) |
|||
say "#{'%2d' % k+1}: #{s}" |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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) |
|||
</pre> |
|||
=={{header|Wren}}== |
=={{header|Wren}}== |
||
Line 586: | Line 769: | ||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
Even with GMP takes about 4½ minutes to find the first 16. |
Even with GMP takes about 4½ minutes to find the first 16. |
||
<syntaxhighlight lang=" |
<syntaxhighlight lang="wren">import "./gmp" for Mpz, Mpq |
||
import "./fmt" for Fmt |
import "./fmt" for Fmt |
||
Revision as of 22:20, 15 December 2023
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
ALGOL 68
Uses Algol 68G's LONG LONG INT which has programmer defined precision.
Re-uses code from the Arithmetic/Rational task, modified to use LONG LONG INT.
BEGIN # find a sequence of primes whose members are the smallest prime whose #
# reciprocal can be added to the sum or the reciprocals of the #
# previous primes and the sum remain below 1 #
PR precision 5000 PR # set the precision of LONG LONG INT #
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 # ;
# code from the Arithmetic/Rational task modified to use LONG LONG INT #
MODE FRAC = STRUCT( LONG LONG INT num #erator#, den #ominator# );
PROC lcm = ( LONG LONG INT a, b )LONG LONG INT: # least common multiple #
a OVER gcd(a, b) * b;
PRIO // = 9; # higher then the ** operator #
OP // = ( LONG LONG INT num, den )FRAC: ( # initialise and normalise #
LONG LONG INT common = gcd( num, den );
IF den < 0 THEN
( -num OVER common, -den OVER common )
ELSE
( num OVER common, den OVER common )
FI
);
OP + = (FRAC a, b)FRAC: (
LONG LONG INT common = lcm( den OF a, den OF b );
FRAC result := ( common OVER den OF a * num OF a + common OVER den OF b * num OF b, common );
num OF result//den OF result
);
OP +:= = (REF FRAC a, FRAC b)REF FRAC: ( a := a + b );
# end code from the Arithmetic/Rational task modified to use LONG LONG INT #
# the sequence starts with the reciprocal of the first prime > 0, i.e. 2 #
LONG LONG INT one = 1;
FRAC sum := one // LONG LONG 2;
print( ( " 1: 2", newline ) );
# each succeeding prime is the next prime > the denominator of the sum #
# divided by the difference of the denominator and the numerator #
FOR i FROM 2 TO 15 DO
LONG LONG INT next := IF num OF sum + 1 = den OF sum
THEN den OF sum + 1
ELSE ( den OF sum OVER ( den OF sum - num OF sum ) ) + 1
FI;
IF NOT ODD next THEN next +:= 1 FI;
WHILE NOT is probably prime( next ) DO next +:= 2 OD;
print( ( whole( i, -2 ), ": " ) );
STRING prime text = whole( next, 0 );
IF INT digits = ( UPB prime text - LWB prime text ) + 1;
digits <= 40
THEN
print( ( prime text ) )
ELSE
print( ( prime text[ : LWB prime text + 19 ], "..."
, prime text[ UPB prime text - 19 : ], " "
, whole( digits, -6 ), " digits"
)
)
FI;
print( ( newline ) );
sum +:= one // next
OD
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
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)
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)
Nim
import std/strformat
import bignum
iterator a075442(): Int =
let One = newRat(1)
var sum = newRat(0)
var p = newInt(0)
while true:
let q = reciprocal(One - sum)
p = nextPrime(if q.isInt: q.num else: q.toInt + 1)
yield p
sum += newRat(1, p)
func compressed(str: string; size: int): string =
## Return a compressed value for long strings of digits.
if str.len <= 2 * size: str
else: &"{str[0..<size]}...{str[^size..^1]} ({str.len} digits)"
var count = 0
for p in a075442():
inc count
echo &"{count:2}: {compressed($p, 20)}"
if count == 15: break
- 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)
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 % @home (4.4Ghz 5600G ) modified with primes up to 1E6 ( Uint16-> Uint32 ) 78498 999979 reducing factor 0.04059798 ... 15 10015 ms ,delta 3510 // first guess than searching for next prime 1139012563947167462...31060548964273180103 ( 2358 digits) 1731172 999979 16 110030 ms ,delta 6493 3696176350563052055...02467094377885929191 ( 4711 digits) 17 5098268 ms ,delta 55552 // first guess than searching for next prime 2104336472479843950...14594683820359204509 ( 9418 digits)
Perl
use strict; use warnings; use feature 'state';
use Math::AnyNum <next_prime ceil>;
sub abbr { my $d=shift; my $l = length $d; $l < 41 ? $d : substr($d,0,20) . '..' . substr($d,-20) . " ($l digits)" }
sub succ_prime {
state $sum = 0;
my $next = next_prime ceil( 1 / (1-$sum) );
$sum += 1/$next;
$next
}
printf "%2d: %s\n", $_, abbr succ_prime for 1..14;
- 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)
Phix
The Julia entry took over 4 hours to get the 17th on this box, so I just went with "keep it all under 10s"...
(As others have alluded, it is all about getting the next prime. Even should you land on one straightaway, it still takes quite some time to prove an n-thousand digit number is prime.)
with javascript_semantics constant limit = iff(platform()=JS?12:14) include mpfr.e mpq {q, r, s, u} = mpq_inits(4,{0,0,0,1}) mpz {p, t} = mpz_inits(2) printf(1,"First %d elements of the sequence:\n", limit); for count=1 to limit do mpq_sub(q, u, s) mpq_inv(q, q) mpq_get_den(t, q) mpz_set_q(p, q) if mpz_cmp_si(t, 1) then mpz_add_si(p, p, 1) end if mpz_nextprime(p, p) printf(1,"%2d: %s\n", {count, mpz_get_short_str(p)}) mpq_set_z(r, p) mpq_inv(r, r) mpq_add(s, s, r) end for
- Output:
First 14 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 (76 digits) 11: 20300753813848234767...91313959045797597991 (150 digits) 12: 20323705381471272842...21649394434192763213 (297 digits) 13: 12748246592672078196...20708715953110886963 (592 digits) 14: 46749025165138838243...65355869250350888941 (1,180 digits)
If you've really got nothing better to do, you can also get, in 10 mins on 64-bit, thrice that on 32-bit:
15: 11390125639471674628...31060548964273180103 (digits: 2358) 16: 36961763505630520555...02467094377885929191 (digits: 4711)
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)
RPL
Limited floating-point precision prevents from finding the correct 7th term. Program starts with a non-empty sequence to avoid the ∑LIST
calculation bug that occurs when the input list has less than two items.
≪ {2 3}
1 4 START 1 OVER INV ∑LIST - INV →NUM IP NEXTPRIME + NEXT
≫ '∑INVPR' STO
- Output:
1: {2 3 7 43 1811 654149}
Sidef
var A075442 = Enumerator({|callback|
var sum = 0
loop {
var p = next_prime(ceil(1/(1-sum)))
sum += 1/p
callback(p)
}
})
A075442.first(15).each_kv {|k,n|
var s = Str(n)
s = "#{s.first(20)}..#{s.last(20)} (#{s.len} digits)" if (s.len > 50)
say "#{'%2d' % k+1}: #{s}"
}
- 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)
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)