Sisyphus sequence
Sisyphus sequence 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.
The Sisyphus sequence is an infinite sequence of positive integers that was devised in 2022 by Eric Angelini and Carole Dubois.
The first term is 1. Subsequent terms are found by applying the following rule:
- If the previous term was even, then halve it.
- If the previous term was odd, then add the smallest prime number that has not yet been added.
1 is odd and so the second term is: 1 + 2 = 3, because 2 is the smallest prime not yet added.
3 is odd and so the third term is: 3 + 3 = 6, because 3 is the smallest prime not yet added.
6 is even and so the fourth term is : 6 ÷ 2 = 3, and so on.
- Task
Find and show on this page (in 10 lines of 10 terms), the first 100 terms of the sequence.
What are the 1,000th, 10,000th, 100,000th and 1,000,000th terms of the sequence and, in each case, what is the highest prime needed to reach them?
If it is difficult or impossible for your language or output device to meet all of these requirements, then just do what you reasonably can.
- Stretch
What are the 10 millionth and 100 millionth terms and the highest prime needed to reach each one?
By the time the 100 millionth term is reached, which number(s) under 250:
- Have not yet occurred in the sequence.
- Have occurred the most times and their number of occurrences.
- Extreme stretch
What is the number of the first term to equal 36?
This was originally set as a challenge by Neil Sloane who was worried by its non-appearance and found by Russ Cox.
- References
- OEIS sequence A350877: The Sisyphus sequence
ALGOL 68
Using a "generous" estimate of the number of primes required to reach the sequence element 100 000 000 (see the Wren sample).
Note that a sieve size of 1 000 000 000 is not possible with Algol 68G (so you'll need to use another compiler), though a sieve size of 100 000 000 is allowed - specify e.g.: -heap 768M on the command line. It took in the region of 2-3 minutes for Algol 68G to find the elements up to 10 000 000 on a Windows 11 laptop.
This sample also shows first the positions of 1..100 in the sequence - apart from the values which don't appear up to element 100 000 000, some of them are quite a long way into the sequence before they appear.
BEGIN # generate elements of the Sysiphus Sequence: see OEIS A350877 #
# returns the largest element in a #
OP MAX = ( []INT a )INT:
IF LWB a >UPB a THEN 0
ELSE
INT result := a[ LWB a ];
FOR i FROM LWB a + 1 TO UPB a DO
IF result < a[ i ] THEN result := a[ i ] FI
OD;
result
FI # MAX # ;
# sieve the primes #
INT sieve max = 1 000 000 000; ### CHANGE TO E.G.: 100 000 000 FOR ALGOL 68G ###
BOOL is odd := TRUE;
[ sieve max ]BOOL sieve; FOR i TO UPB sieve DO sieve[ i ] := is odd; is odd := NOT is odd OD;
sieve[ 1 ] := FALSE;
sieve[ 2 ] := TRUE;
FOR s FROM 3 BY 2 TO ENTIER sqrt( sieve max ) DO
IF sieve[ s ] THEN
FOR p FROM s * s BY s TO sieve max DO sieve[ p ] := FALSE OD
FI
OD;
[ 1 : 250 ]INT pos; # positions of 1..250 in the sequence #
[ 1 : 250 ]INT occurs; # occurances of 1..250 in the sequence #
FOR i TO UPB pos DO pos[ i ] := occurs[ i ] := 0 OD;
INT max count = sieve max OVER 10; # highest element required #
INT s := 1; # the first element is defined as 1 #
INT count := 1; # count of elements found so far #
print( ( "Sysiphus sequence - first 100 elements:", newline ) );
print( ( whole( s, -4 ) ) );
pos[ s ] := count;
INT next to show := 1000; # next power-of-10 element to show #
INT last used prime := 0; # latest prime from the list #
INT p pos := 0; # current position in the sieve #
WHILE count < max count DO
# calculate the next element #
IF NOT ODD s THEN
# the previous element was even - halve it #
s OVERAB 2
ELSE
# the previous element was odd: add the next prime from the list #
WHILE p pos +:= 1;
NOT sieve[ p pos ]
DO SKIP OD;
s +:= ( last used prime := p pos )
FI;
count +:= 1;
IF count <= 100 THEN # have one of the first 100 elements #
print( ( whole( s, -4 ) ) );
IF count MOD 10 = 0 THEN print( ( newline ) ) FI;
IF count = 100 THEN print( ( newline ) ) FI
ELIF count = next to show THEN
# reached a power of ten count #
print( ( "sequence element ", whole( count, -10 )
, " is ", whole( s, -10 )
, ", highest used prime is ", whole( last used prime, -10 )
, newline
)
);
next to show *:= 10
FI;
IF s < UPB pos THEN
IF pos[ s ] = 0 THEN
# have the first appearence of s in the sequence #
pos[ s ] := count
FI;
occurs[ s ] +:= 1
FI
OD;
print( ( newline ) );
print( ( "Integers in 1..", whole( UPB pos, 0 )
, " not found in the sequence up to element ", whole( max count, 0 )
, ":", newline
)
);
FOR i TO UPB pos DO
IF pos[ i ] = 0 THEN print( ( " ", whole( i, 0 ) ) ) FI
OD;
print( ( newline ) );
INT max occurs = MAX occurs;
print( ( "Integers in 1..", whole( UPB pos, 0 )
, " that occur most often ( ", whole( max occurs, 0 )
, " times ) up to element ", whole( max count, 0 )
, ":", newline
)
);
FOR i TO UPB occurs DO
IF occurs[ i ] = max occurs THEN print( ( " ", whole( i, 0 ) ) ) FI
OD;
print( ( newline, newline ) );
print( ( "Position in the sequence of 1..100 up to element ", whole( max count, 0 )
, ":", newline
)
);
FOR i TO 100 DO
print( ( IF pos[ i ] = 0 THEN " unknown" ELSE whole( pos[ i ], -8 ) FI ) );
IF i MOD 8 = 0 THEN print( ( newline ) ) FI
OD
END- Output:
Sysiphus sequence - first 100 elements:
1 3 6 3 8 4 2 1 8 4
2 1 12 6 3 16 8 4 2 1
18 9 28 14 7 30 15 44 22 11
42 21 58 29 70 35 78 39 86 43
96 48 24 12 6 3 62 31 92 46
23 90 45 116 58 29 102 51 130 65
148 74 37 126 63 160 80 40 20 10
5 106 53 156 78 39 146 73 182 91
204 102 51 178 89 220 110 55 192 96
48 24 12 6 3 142 71 220 110 55
sequence element 1000 is 990, highest used prime is 2273
sequence element 10000 is 24975, highest used prime is 30713
sequence element 100000 is 265781, highest used prime is 392111
sequence element 1000000 is 8820834, highest used prime is 4761697
sequence element 10000000 is 41369713, highest used prime is 55900829
sequence element 100000000 is 1179614168, highest used prime is 640692323
Integers in 1..250 not found in the sequence up to element 100000000:
36 72 97 107 115 127 144 167 194 211 214 230 232 250
Integers in 1..250 that occur most often ( 7 times ) up to element 100000000:
7 14 28
Position in the sequence of 1..100 up to element 100000000:
1 7 2 6 71 3 25 5
22 70 30 13 345 24 27 16
161 21 148 69 32 29 51 43
1154 344 161336 23 34 26 48 737
156 160 36 unknown 63 147 38 68
234 31 40 28 53 50 126 42
639 1153 58 343 73 161335 88 111
108 33 135 614667 192 47 65 736
60 155 454 159 186 35 97 unknown
78 62 2340 146 143 37 24841 67
476 233 433 10579 140 39 359 169
85 52 80 49 195 125 166 41
unknown 638 204 1152
Using Algol 68G with max sieve set to 100 000 000 (sequence elements up to 10 000 000) is as above, except for the missing elements and counts of occurrences (and obviously, the line showing the 100 000 000th element is not present):
Integers in 1..250 not found in the sequence up to element 10000000: 36 72 97 107 115 127 144 167 194 211 214 223 230 232 250 Integers in 1..250 that occur most often ( 6 times ) up to element 10000000: 3 57 65 85 114 125 130 170 228
Phix
atom t0 = time() constant limit = 1e8 sequence sisyphus = {1}, under250 = reinstate(repeat(0,250),1,1) integer np = 0, specific = 1000, count = 1, m atom next = 1 while true do if even(next) then next /= 2 else np += 1 next += get_prime(np) end if if next <= 250 then under250[next] += 1 end if count += 1 if count <= 100 then sisyphus &= next if count == 100 then printf(1,"The first 100 members of the Sisyphus sequence are:\n%s\n", {join_by(sisyphus,1,10," ",fmt:="%3d")}) end if elsif count == specific then printf(1,"%,13d%s member is: %,13d and highest prime needed: %,11d\n", {count, ord(count), next, get_prime(np)}) if count == limit then m = largest(under250,true); exit end if specific *= 10 end if end while printf(1,"\nThese numbers under 250 do not occur in the first %,d terms:\n",count) printf(1," %v\n",{find_all(0,under250)}) printf(1,"\nThese numbers under 250 occur the most in the first %,d terms:\n",count) printf(1," %v all occur %d times.\n",{find_all(m,under250),m}) ?elapsed(time()-t0)
- Output:
The first 100 members of the Sisyphus sequence are:
1 3 6 3 8 4 2 1 8 4
2 1 12 6 3 16 8 4 2 1
18 9 28 14 7 30 15 44 22 11
42 21 58 29 70 35 78 39 86 43
96 48 24 12 6 3 62 31 92 46
23 90 45 116 58 29 102 51 130 65
148 74 37 126 63 160 80 40 20 10
5 106 53 156 78 39 146 73 182 91
204 102 51 178 89 220 110 55 192 96
48 24 12 6 3 142 71 220 110 55
1,000th member is: 990 and highest prime needed: 2,273
10,000th member is: 24,975 and highest prime needed: 30,713
100,000th member is: 265,781 and highest prime needed: 392,111
1,000,000th member is: 8,820,834 and highest prime needed: 4,761,697
10,000,000th member is: 41,369,713 and highest prime needed: 55,900,829
100,000,000th member is: 1,179,614,168 and highest prime needed: 640,692,323
These numbers under 250 do not occur in the first 100,000,000 terms:
{36,72,97,107,115,127,144,167,194,211,214,230,232}
These numbers under 250 occur the most in the first 100,000,000 terms:
{7,14,28} all occur 7 times.
"12.8s"
Extreme stretch (just a note about)
Knowing that s(77,534,485,877)=36, a crude extrapolation of the above figures suggests it would need primes up to 570 billion, or about 30 billion of them, or at least 480GB of RAM, which is a tad more than I have available on this machine... If, on the other hand, I could figure out how to not need quite so much ram (I believe the theoretical minimum to be 36GB), it might be achievable in as little as between 4 and 28 hours... Update: as per the talk page, having a prebuilt 36GB odds-only bitsieve knocking about on your hard drive would certainly help, but as it happens I don't.
Python
from collections import Counter
from itertools import islice
from typing import Iterable
from typing import Iterator
from typing import Tuple
from typing import TypeVar
import primesieve
def primes() -> Iterator[int]:
it = primesieve.Iterator()
while True:
yield it.next_prime()
def sisyphus() -> Iterator[Tuple[int, int]]:
prime = primes()
n = 1
p = 0
yield n, p
while True:
if n % 2:
p = next(prime)
n = n + p
else:
n = n // 2
yield n, p
def consume(it: Iterator[Tuple[int, int]], n) -> Tuple[int, int]:
next(islice(it, n - 1, n - 1), None)
return next(it)
T = TypeVar("T")
def batched(it: Iterable[T], n: int) -> Iterable[Tuple[T, ...]]:
_it = iter(it)
batch = tuple(islice(_it, n))
while batch:
yield batch
batch = tuple(islice(_it, n))
if __name__ == "__main__":
it = sisyphus()
first_100 = list(islice(it, 100))
print("The first 100 members of the Sisyphus sequence are:")
for row in batched(first_100, 10):
print(" ".join(str(n).rjust(3) for n, _ in row))
print("")
for interval in [10**x for x in range(3, 9)]:
n, prime = consume(it, interval - (interval // 10))
print(f"{interval:11,}th number: {n:13,} highest prime needed: {prime:11,}")
print("")
sisyphus_lt_250 = Counter(n for n, _ in islice(sisyphus(), 10**8) if n < 250)
print("These numbers under 250 do not occur in the first 100,000,000 terms:")
print(" ", [n for n in range(1, 250) if n not in sisyphus_lt_250])
print("")
most_common = sisyphus_lt_250.most_common(1)[0][1]
print("These numbers under 250 occur the most in the first 100,000,000 terms:")
print(
f" {[n for n, c in sisyphus_lt_250.items() if c == most_common]} "
f"all occur {most_common} times."
)
- Output:
The first 100 members of the Sisyphus sequence are:
1 3 6 3 8 4 2 1 8 4
2 1 12 6 3 16 8 4 2 1
18 9 28 14 7 30 15 44 22 11
42 21 58 29 70 35 78 39 86 43
96 48 24 12 6 3 62 31 92 46
23 90 45 116 58 29 102 51 130 65
148 74 37 126 63 160 80 40 20 10
5 106 53 156 78 39 146 73 182 91
204 102 51 178 89 220 110 55 192 96
48 24 12 6 3 142 71 220 110 55
1,000th number: 990 highest prime needed: 2,273
10,000th number: 24,975 highest prime needed: 30,713
100,000th number: 265,781 highest prime needed: 392,111
1,000,000th number: 8,820,834 highest prime needed: 4,761,697
10,000,000th number: 41,369,713 highest prime needed: 55,900,829
100,000,000th number: 1,179,614,168 highest prime needed: 640,692,323
These numbers under 250 do not occur in the first 100,000,000 terms:
[36, 72, 97, 107, 115, 127, 144, 167, 194, 211, 214, 230, 232]
These numbers under 250 occur the most in the first 100,000,000 terms:
[28, 14, 7] all occur 7 times.
Extreme stretch
Using Python bindings for the primesieve C++ library and PyPy 3.9 instead of CPython (PyPy is nearly twice as fast up to s(100,000,000)), this completed in about 1 hour and 53 minutes on a AMD Ryzen 5 1500X. Memory usage was minimal thanks to primesieve.Iterator().
import primesieve
def sisyphus36():
primes = primesieve.Iterator()
n = 1
p = 0
i = 1
while True:
i += 1
if n % 2:
p = primes.next_prime()
n = n + p
else:
n = n // 2
if n == 36:
print(f"{i:,}, {n:,}, {p:,}")
break
if __name__ == "__main__":
sisyphus36()
- Output:
$ time pypy3.9 sisyphus_36.py 77,534,485,877, 36, 677,121,348,413 real 113m2.636s user 112m56.973s sys 0m0.060s
Raku
use Math::Primesieve;
use Lingua::EN::Numbers;
my ($exp1, $exp2, $limit1, $limit2) = 3, 8, 100, 250;
my ($n, $s0, $s1, $p, @S1, %S) = 1, 1, Any, Any, 1;
my $iterator = Math::Primesieve::iterator.new;
my @Nth = ($exp1..$exp2)».exp(10);
my $S2 = BagHash.new;
repeat {
$n++;
$s1 = $s0 %% 2 ?? $s0 div 2 !! $s0 + ($p = $iterator.next);
@S1.push: $s1 if $n ≤ $limit1;
$S2.add: $s1 if $s1 ≤ $limit2;
%S{$n}{'value', 'prime'} = $s1, $p if $n ∈ @Nth;
$s0 = $s1;
} until $n == @Nth[*-1];
say "The first $limit1 members of the Sisyphus sequence are:";
say @S1.batch(10)».fmt('%4d').join("\n") ~ "\n";
printf "%12sth member is: %13s with prime: %11s\n", ($_, %S{$_}{'value'}, %S{$_}{'prime'})».&comma for @Nth;
say "\nNumbers under $limit2 that do not occur in the first {comma @Nth[*-1]} terms:";
say (1..$limit2).grep: * ∉ $S2.keys;
say "\nNumbers under $limit2 that occur the most ({$S2.values.max} times) in the first {comma @Nth[*-1]} terms:";
say $S2.keys.grep({ $S2{$_} == $S2.values.max}).sort;
- Output:
The first 100 members of the Sisyphus sequence are:
1 3 6 3 8 4 2 1 8 4
2 1 12 6 3 16 8 4 2 1
18 9 28 14 7 30 15 44 22 11
42 21 58 29 70 35 78 39 86 43
96 48 24 12 6 3 62 31 92 46
23 90 45 116 58 29 102 51 130 65
148 74 37 126 63 160 80 40 20 10
5 106 53 156 78 39 146 73 182 91
204 102 51 178 89 220 110 55 192 96
48 24 12 6 3 142 71 220 110 55
1,000th member is: 990 with prime: 2,273
10,000th member is: 24,975 with prime: 30,713
100,000th member is: 265,781 with prime: 392,111
1,000,000th member is: 8,820,834 with prime: 4,761,697
10,000,000th member is: 41,369,713 with prime: 55,900,829
100,000,000th member is: 1,179,614,168 with prime: 640,692,323
Numbers under 250 that do not occur in the first 100,000,000 terms:
36 72 97 107 115 127 144 167 194 211 214 230 232
Numbers under 250 that occur the most (7 times) in the first 100,000,000 terms:
7 14 28
RPL
This is a very slow version (15 minutes to get the 10,000th term on a HP-50g) , but I'm not sure to have enough memory to generate a large enough sieve to find the millionth term in a reasonable amount of time. Anyway, even the first requirement (displaying 10 lines of 10 items) is out of range with a 22-character screen.
≪ 1 CF IF DUP 0 < THEN 1 SF NEG END @ A negative argument requires the sequence to be stored then displayed 0 { 1 } → n cnt seq ≪ 2 1 WHILE 'cnt' INCR n < REPEAT IF DUP 2 MOD THEN OVER + SWAP NEXTPRIME SWAP ELSE 2 / END IF 1 FS? THEN 'seq' OVER STO+ END END IF 1 FS? THEN DROP2 seq ELSE SWAP PREVPRIME R→C END ≫ ≫ 'SISYPH' STO
-100 SISYPH 1000 SISYPH 10000 SISYPH
- Output:
3: {1 3 6 3 8 4 2 1 8 4 2 1 12 6 3 16 8 4 2 1 18 9 28 14 7 30 15 44 22 11 42 21 58 29 70 35 78 39 86 43 96 48 24 12 6 3 62 31 92 46 23 90 45 116 58 29 102 51 130 65 148 74 37 126 63 160 80 40 20 10 5 106 53 156 78 39 146 73 182 91 204 102 51 178 89 220 110 55 192 96 48 24 12 6 3 142 71 220 110 55 }
2: (990.,2273.)
1: (24975.,30713.)
Wren
No option here but to use a sieve as relying on the 'nextPrime' method would be far too slow to achieve the stretch goal in a reasonable time using Wren. Sieve limit found by experimentation.
Extreme stretch not attempted and probably out of the question for Wren.
import "./math" for Int, Nums
import "./fmt" for Fmt
var limit = 1e8
var primes = Int.primeSieve(7 * limit)
var under250 = List.filled(250, 0)
var sisyphus = [1]
under250[1] = 1
var prev = 1
var nextPrimeIndex = 0
var specific = 1000
var count = 1
while (true) {
var next
if (prev % 2 == 0) {
next = prev / 2
} else {
next = prev + primes[nextPrimeIndex]
nextPrimeIndex = nextPrimeIndex + 1
}
count = count + 1
if (count <= 100) sisyphus.add(next)
if (next < 250) under250[next] = under250[next] + 1
if (count == 100) {
System.print("The first 100 members of the Sisyphus sequence are:")
Fmt.tprint("$3d ", sisyphus, 10)
System.print()
} else if (count == specific) {
var prime = primes[nextPrimeIndex-1]
Fmt.print("$,13r member is: $,13d and highest prime needed: $,11d", count, next, prime)
if (count == limit) {
var notFound = (1..249).where { |i| under250[i] == 0 }.toList
var max = Nums.max(under250)
var maxFound = (1..249).where { |i| under250[i] == max }.toList
Fmt.print("\nThese numbers under 250 do not occur in the first $,d terms:", count)
Fmt.print(" $n", notFound)
Fmt.print("\nThese numbers under 250 occur the most in the first $,d terms:", count)
Fmt.print(" $n all occur $d times.", maxFound, max)
return
}
specific = 10 * specific
}
prev = next
}
- Output:
The first 100 members of the Sisyphus sequence are:
1 3 6 3 8 4 2 1 8 4
2 1 12 6 3 16 8 4 2 1
18 9 28 14 7 30 15 44 22 11
42 21 58 29 70 35 78 39 86 43
96 48 24 12 6 3 62 31 92 46
23 90 45 116 58 29 102 51 130 65
148 74 37 126 63 160 80 40 20 10
5 106 53 156 78 39 146 73 182 91
204 102 51 178 89 220 110 55 192 96
48 24 12 6 3 142 71 220 110 55
1,000th member is: 990 and highest prime needed: 2,273
10,000th member is: 24,975 and highest prime needed: 30,713
100,000th member is: 265,781 and highest prime needed: 392,111
1,000,000th member is: 8,820,834 and highest prime needed: 4,761,697
10,000,000th member is: 41,369,713 and highest prime needed: 55,900,829
100,000,000th member is: 1,179,614,168 and highest prime needed: 640,692,323
These numbers under 250 do not occur in the first 100,000,000 terms:
[36, 72, 97, 107, 115, 127, 144, 167, 194, 211, 214, 230, 232]
These numbers under 250 occur the most in the first 100,000,000 terms:
[7, 14, 28] all occur 7 times.