[seqfan] A069105

[Don Reble]

Seqfans:

    Now that A069105 is corrected (Thanks, Ray.), note that one can
    compute such sequences quickly using factorizations of 10^n-1.
    See below.

-- 
Don Reble  djr at nk.ca


#!/usr/bin/python
# Print numbers X for which 1/X (in decimal) has a particular period.
# Terminators have period 1 because of the trailing zeroes.

import sys

# prime factors of (10^n - 1)
PrFacs = (
    (), # 0
    ( 3, ),
    ( 3, 11 ),
    ( 3, 37 ),
    ( 3, 11, 101 ),
    ( 3, 41, 271 ),
    ( 3, 7, 11, 13, 37 ),
    ( 3, 239, 4649 ),
    ( 3, 11, 73, 101, 137 ),
    ( 3, 37, 333667 ),
    ( 3, 11, 41, 271, 9091 ), # 10
    ( 3, 21649, 513239 ),
    ( 3, 7, 11, 13, 37, 101, 9901 ),
    ( 3, 53, 79, 265371653 ),
    ( 3, 11, 239, 4649, 909091 ),
    ( 3, 31, 37, 41, 271, 2906161 ),
    ( 3, 11, 17, 73, 101, 137, 5882353 ),
    ( 3, 2071723, 5363222357 ),
    ( 3, 7, 11, 13, 19, 37, 52579, 333667 ),
    ( 3, 1111111111111111111 ),
    ( 3, 11, 41, 101, 271, 3541, 9091, 27961 ), # 20
    ( 3, 37, 43, 239, 1933, 4649, 10838689 ),
    ( 3, 11, 23, 4093, 8779, 21649, 513239 ),
    ( 3, 11111111111111111111111 ),
    ( 3, 7, 11, 13, 37, 73, 101, 137, 9901, 99990001 ),
    ( 3, 41, 271, 21401, 25601, 182521213001 ),
    ( 3, 11, 53, 79, 859, 265371653, 1058313049 ),
    ( 3, 37, 757, 333667, 440334654777631 ),
    ( 3, 11, 29, 101, 239, 281, 4649, 909091, 121499449 ),
    ( 3, 3191, 16763, 43037, 62003, 77843839397 ),
    ( 3, 7, 11, 13, 31, 37, 41, 211, 241, 271, 2161, 9091, 2906161 ), # 30
    ( 3, 2791, 6943319, 57336415063790604359 ),
    ( 3, 11, 17, 73, 101, 137, 353, 449, 641, 1409, 69857, 5882353 ),
)

def generate(prod, facseq):
    global SmallPows, LargePow
    if len(facseq) == 0:
        for sp in SmallPows:
            if (sp % prod) == 0: return
        yield prod
    else:
        while (LargePow % prod) == 0:
            for val in generate(prod, facseq[1:]): yield val
            prod *= facseq[0]


# check period parameter
if len(sys.argv) != 2:
    print "usage:", sys.argv[0], "period"
    sys.exit(1)
Period = int(sys.argv[1])
if (Period < 1) or (Period >= len(PrFacs)):
    print "0 < period <", len(PrFacs), "please"
    sys.exit(1)

# compute numbers which divide 10^period-1 but no smaller 10^n-1:

# compute those 10^n-1 values
SmallPows = []
tenpow = 10
for exp in range(1,Period):
    if (Period % exp) == 0: SmallPows += [tenpow-1]
    tenpow *= 10
LargePow = tenpow - 1

# generate suitable products of PrFacs[Period]
Vals = set(generate(1, PrFacs[Period]))

# do powers of 2 and 5
while True:
    val = min(Vals)
    print val
    Vals.remove(val)
    Vals.add(val*2)
    Vals.add(val*5)