Refactor 8 sequences (#390)

* Remove while-break from A006577

* Refactor A000045 to remove `yield 0`

* Refactor A000217 (Triangle numbers) to avoid factorials

The previous code appears to use binomial(n+1, 2) rather than the faster n(n+1)/2 form. Reduced from O(n) multiplications to O(1)

* Refactor A000203 to keep running sum

Reduces space use (single int is much smaller than a list)

* Refactor A133058 to remove extra if branch

* Change range in A000005 to remove zero-check

* Square duplicated factorial from A000108

2*i-i == i, so just square i!

* Refactor A001622 to be correct for all n

The prior version would be inaccurate for values of n over 99999 (but it took a long time to evaluate that far)
Test validates 2000 decimal places (from goldennumber.net), to check that precision scales at the right time

---------

Co-authored-by: Julien Palard <[email protected]>

oeis.py

@@ -288,17 +288,13 @@ def A181391() -> Iterable[int]:
 @oeis.from_function(offset=1)
 def A006577(n: int) -> int:
     """Give the number of halving and tripling steps to reach 1 in '3x+1' problem."""
-    if n == 1:
-        return 0
     x = 0
-    while True:
+    while n > 1:
         if n % 2 == 0:
             n //= 2
         else:
             n = 3 * n + 1
         x += 1
-        if n < 2:
-            break
     return x
 
 
@@ -329,10 +325,9 @@ def A001221(n: int) -> int:
 def A000045() -> Iterable[int]:
     """Fibonacci numbers: F(n) = F(n-1) + F(n-2) with F(0) = 0 and F(1) = 1."""
     a, b = (0, 1)
-    yield 0
     while True:
-        a, b = b, a + b
         yield a
+        a, b = b, a + b
 
 
 @oeis.from_generator()
@@ -416,11 +411,9 @@ def A000142(n: int) -> int:
 
 
 @oeis.from_function()
-def A000217(i: int):
+def A000217(n: int):
     """Triangular numbers: a(n) = binomial(n+1,2) = n(n+1)/2 = 0 + 1 + 2 + ... + n."""
-    if i < 1:
-        return 0
-    return math.factorial(i + 1) // math.factorial(2) // math.factorial((i + 1) - 2)
+    return n * (n + 1) // 2
 
 
 @oeis.from_function()
@@ -487,17 +480,13 @@ def A000203(i: int) -> int:
 
     a(n) = sigma(n). Also called sigma_1(n).
     """
-    divisors = []
-    for j in range(int(math.sqrt(i)) + 1):
-        if j == 0:
-            continue
+    divisor_sum = 0
+    for j in range(1, int(math.sqrt(i)) + 1):
         if i % j == 0:
-            if i / j == j:
-                divisors.append(j)
-            else:
-                divisors.append(j)
-                divisors.append(i // j)
-    return int(sum(divisors))
+            divisor_sum += j
+            if i // j != j:
+                divisor_sum += i // j
+    return divisor_sum
 
 
 @oeis.from_function()
@@ -526,24 +515,21 @@ def A133058() -> Iterable[int]:
     otherwise a(n) = a(n-1)/gcd(a(n-1),n).
     """
     last = 1
-    for i in count():
-        if i in (0, 1):
-            yield 1
-        elif (math.gcd(i, last)) == 1:
+    yield 1
+    yield 1
+    for i in count(2):
+        if (math.gcd(i, last)) == 1:
             last = last + i + 1
-            yield last
         else:
             last = int(last / math.gcd(last, i))
-            yield last
+        yield last
 
 
 @oeis.from_function(offset=1)
 def A000005(i: int) -> int:
     """d(n) (also called tau(n) or sigma_0(n)), the number of divisors of n."""
     divisors = 0
-    for j in range(int(math.sqrt(i)) + 1):
-        if j == 0:
-            continue
+    for j in range(1, int(math.sqrt(i)) + 1):
         if i % j == 0:
             if i / j == j:
                 divisors += 1
@@ -558,12 +544,7 @@ def A000108(i: int) -> int:
 
     Also called Segner numbers.
     """
-    return (
-        math.factorial(2 * i)
-        // math.factorial(i)
-        // math.factorial(2 * i - i)
-        // (i + 1)
-    )
+    return math.factorial(2 * i) // math.factorial(i) ** 2 // (i + 1)
 
 
 @oeis.from_function()
@@ -596,10 +577,14 @@ def A000120(n: int) -> int:
 def A001622() -> Iterable[int]:
     """Decimal expansion of golden ratio phi (or tau) = (1 + sqrt(5))/2."""
     with localcontext() as ctx:
-        ctx.prec = 99999
-        tau = (1 + Decimal(5).sqrt()) / 2
-        for n in count():
-            yield math.floor(tau * 10**n) % 10
+        ctx.prec = 10
+        start = 0
+        while True:
+            ctx.prec *= 10
+            phi = (1 + Decimal(5).sqrt()) / 2
+            for n in range(start, ctx.prec - 1):
+                yield math.floor(phi * 10**n) % 10
+            start = n + 1
 
 
 @oeis.from_function(offset=1)

tests/test_A001622.py

@@ -10,3 +10,48 @@ def test_sequence():
         1 1 7 7 2 0 3 0 9 1 7 9 8 0 5 7 6 2 8 6 2 1 3 5 4 4 8 6 2 2 7 0 5 2 6 0 4 6 2 8
         1 8 9 0 2 4 4 9 7 0 7 2 0 7 2 0 4 1 8 9 3 9 1 1 3 7 4 8 4 7 5""".split()
     ]
+    assert A001622[2:2002] == list(
+        map(
+            int,
+            "61803398874989484820458683436563811772030917980576"
+            "28621354486227052604628189024497072072041893911374"
+            "84754088075386891752126633862223536931793180060766"
+            "72635443338908659593958290563832266131992829026788"
+            "06752087668925017116962070322210432162695486262963"
+            "13614438149758701220340805887954454749246185695364"
+            "86444924104432077134494704956584678850987433944221"
+            "25448770664780915884607499887124007652170575179788"
+            "34166256249407589069704000281210427621771117778053"
+            "15317141011704666599146697987317613560067087480710"
+            "13179523689427521948435305678300228785699782977834"
+            "78458782289110976250030269615617002504643382437764"
+            "86102838312683303724292675263116533924731671112115"
+            "88186385133162038400522216579128667529465490681131"
+            "71599343235973494985090409476213222981017261070596"
+            "11645629909816290555208524790352406020172799747175"
+            "34277759277862561943208275051312181562855122248093"
+            "94712341451702237358057727861600868838295230459264"
+            "78780178899219902707769038953219681986151437803149"
+            "97411069260886742962267575605231727775203536139362"
+            "10767389376455606060592165894667595519004005559089"
+            "50229530942312482355212212415444006470340565734797"
+            "66397239494994658457887303962309037503399385621024"
+            "23690251386804145779956981224457471780341731264532"
+            "20416397232134044449487302315417676893752103068737"
+            "88034417009395440962795589867872320951242689355730"
+            "97045095956844017555198819218020640529055189349475"
+            "92600734852282101088194644544222318891319294689622"
+            "00230144377026992300780308526118075451928877050210"
+            "96842493627135925187607778846658361502389134933331"
+            "22310533923213624319263728910670503399282265263556"
+            "20902979864247275977256550861548754357482647181414"
+            "51270006023890162077732244994353088999095016803281"
+            "12194320481964387675863314798571911397815397807476"
+            "15077221175082694586393204565209896985556781410696"
+            "83728840587461033781054443909436835835813811311689"
+            "93855576975484149144534150912954070050194775486163"
+            "07542264172939468036731980586183391832859913039607"
+            "20144559504497792120761247856459161608370594987860"
+            "06970189409886400764436170933417270919143365013715",
+        )
+    )