Prime reciprocal sum: Difference between revisions
Content deleted Content added
m →{{header|ALGOL 68}}: layout |
Thundergnat (talk | contribs) m Promote to full task. Over a year as a draft, multiple implementations, little controversy. |
||
| (18 intermediate revisions by 9 users not shown) | |||
Line 1:
{{
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'''.
Line 15:
;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.
If any of the tests for primality used in
your program are probabilistic please so indicate.
Line 32 ⟶ 34:
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
Uses Algol 68G's LONG LONG INT which has programmer defined precision.<br>
Re-uses code from the [[Arithmetic/Rational]] task, modified to use LONG LONG INT.<br/>
Also uses Miller-Rabin (probabilistic) primality testing.
{{libheader|ALGOL 68-primes}}
<syntaxhighlight lang="algol68">
Line 39 ⟶ 42:
# previous primes and the sum remain below 1 #
PR precision
PR read "primes.incl.a68" PR # include prime utilities #
Line 87 ⟶ 90:
# 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
LONG LONG INT next := IF num OF sum + 1 = den OF sum
THEN den OF sum + 1
Line 129 ⟶ 132:
13: 12748246592672078196...20708715953110886963 592 digits
14: 46749025165138838243...65355869250350888941 1180 digits
15: 11390125639471674628...31060548964273180103 2358 digits
</pre>
Line 244 ⟶ 248:
Tested in J9.4
=={{header|jq}}==
'''Works with jq, the C implementation of jq''' (*)
'''Works with gojq, the Go implementation of jq''' (**)
(*) Using the two programs presented here, the C implementation
cannot get beyond the 7th prime in the series because it lacks
support for both "big integers" and "big floats".
(**) The `nextPrime` algorithm effectively prevents the Go
implementation from generating more than eight primes in the sequence.
gojq does support infinite-precision integer arithmetic, which in
principle should allow the first program, which uses rational number
arithmetic, to proceed indefinitely but for the slowness of
`nextPrime`.
==={{header|Using rational.jq}}===
The "rational" module referenced here is available at
[[:Category:jq/rational.jq]].
<syntaxhighlight lang="jq">
include "rational" {search:"."}; # see above
# itrunc/0 attempts to compute the "trunc" of the input number in such
# a way that gojq will recognize the result as an integer, while
# leveraging the support that both the C and Go implementations have
# for integer literals.
# It is mainly intended for numbers for which the tostring value does NOT contain an exponent.
# The details are complicated because the goal is to provide portability between jq implementations.
#
# Input: a number or a string matching ^[0-9]+([.][0-9]*)$
# Output: as specified by the implementation.
def itrunc:
if type == "number" and . < 0 then - ((-.) | itrunc)
else . as $in
| tostring as $s
| ($s | test("[Ee]")) as $exp
| ($s | test("[.]")) as $decimal
| if ($exp|not) and ($decimal|not) then $s | tonumber # do not simply return $in
elif ($exp|not) and $decimal then $s | sub("[.].*";"") | tonumber
else trunc
| tostring
| if test("[Ee]")
then if $exp then "itrunc cannot be used on \($in)" | error end
else tonumber
end
end
end;
# rtrunc is like trunc but for Rational numbers, though the input may be
# either a number or a Rational.
def rtrunc:
if type == "number" then r(itrunc;1)
elif 0 == .n or (0 < .n and .n < .d) then r(0;1)
elif 0 < .n or (.n % .d == 0) then .d as $d | r(.n | idivide($d); 1)
else rminus( r( - .n; .d) | rtrunc | rminus; 1)
end;
def is_prime:
. as $n
| if ($n < 2) then false
elif ($n % 2 == 0) then $n == 2
elif ($n % 3 == 0) then $n == 3
elif ($n % 5 == 0) then $n == 5
elif ($n % 7 == 0) then $n == 7
elif ($n % 11 == 0) then $n == 11
elif ($n % 13 == 0) then $n == 13
elif ($n % 17 == 0) then $n == 17
elif ($n % 19 == 0) then $n == 19
else sqrt as $s
| 23
| until( . > $s or ($n % . == 0); . + 2)
| . > $s
end;
# Input: any number
# Output: the next prime
def nextPrime:
if . < 2 then 2
else itrunc # for gojq (to ensure result is a Go integer)
| (if .%2 == 1 then .+2 else .+1 end) as $n
| first( range($n; infinite; 2) | select(is_prime))
end;
# Input: a Rational
# Output: an integer
def rnextPrime:
if rlessthan(r(2;1)) then 2
else rtrunc
| .n
| nextPrime
end ;
def synopsis:
tostring as $i
| ($i|length) as $n
| if $n > 40
then "\($i[:20]) ... \($i[20:]) (\($n))"
else $i
end;
def prs_sequence:
r(1;1) as $one
| {p: 0, s: r(0;1)}
| recurse(
.emit = null
| .count += 1
| .p = (rminus($one; .s) | rinv | rnextPrime)
| .s = radd(.s; r(1; .p))
| .emit = .p )
| select(.emit).emit ;
16 as $n
| "First \($n) elements of the sequence:",
limit($n; prs_sequence)
</syntaxhighlight>
{{output}}
The following output was generated using gojq.
<pre>
First 16 elements of the sequence:
2
3
7
43
1811
654149
27082315109
153694141992520880899
(execution terminated)
</pre>
==={{header|Using floats}}===
The program is essentially as above, except that the directive for
including the "rational" module can be omitted, and the `prs_sequence`
function is replaced by the following definition:
<syntaxhighlight lang="jq">
def prs_sequence:
{p: 0, s: 0}
| recurse(
.emit = null
| .count += 1
| .p = ((1 / (1 - .s)) | nextPrime)
| .s += (1 / .p)
| if .s >= 1 then "Insufficient precision to proceed"|error end
| .emit = .p )
| select(.emit).emit;
</syntaxhighlight>
<pre>
"First 16 elements of the sequence:"
2
3
7
43
1811
654149
27082315109
(program halted)
</pre>
=={{header|Julia}}==
Line 294 ⟶ 462:
16: 36961763505630520555..02467094377885929191 (4711 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>
Line 510 ⟶ 723:
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 ''probably'' prime - the standard prime tests are only deterministic to 3,317,044,064,679,887,385,961,981 but at least the 9th and later results agree with others on this page.</small>
<!--(phixonline)-->
<syntaxhighlight lang="phix">
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
</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>
Line 581 ⟶ 887:
15: 11390125639471674628..31060548964273180103 (2358 digits)
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}}==
Line 586 ⟶ 937:
{{libheader|Wren-fmt}}
Even with GMP takes about 4½ minutes to find the first 16.
GMP’s mpz_nextprime function (which Wren is calling under the hood) uses a probabilistic algorithm to identify primes but, in practice, the chance of a composite number passing is extremely small and should be nil for numbers less than 2^64.
<syntaxhighlight lang="wren">import "./gmp" for Mpz, Mpq
import "./fmt" for Fmt
| |||