Better divisibility rule finder

import math
import sys
from tabulate import tabulate
from dataclasses import dataclass

# This terrible python code is provided by Kaelygon
#
# Find divisibility rule by prime P: P * N = a * 10^b + c
# such that an integer can be split in two parts A and B
# and reduced to a smaller number by simulating 10^b division and modulus
# using string manipulation and addition: B*a-A*c
#
# This doesn't always hold for composites if a and c share factors
#
# Matt Parker explanation https://youtu.be/6pLz8wEQYkA

# Define the custom countDigits function
def countDigits(n): return int(math.log10(abs(n))) + 1

def intDiv(a,b): return (a+b//2)//b #Rounded integer division

#Find divisibility rules using the fact P*N = a*10^b+c
def findDivRulesByInverse(P, aRange, bRange, cRange, recurseDepth=0, maxRecursion=2):
    rules = []
    recurseOffset=0
    if recurseDepth: # nudge N for reattempts, +1, -1, +2, -2, +3...
        recurseOffset = recurseDepth+1
        recurseOffset = -recurseOffset//2 if recurseOffset%2==0 else recurseOffset//2


    for b in range(bRange[0], bRange[1]):
        bPowTen=10**b # 10^b
        for a in range(1,aRange):
            # Find nearest multiples by reciprocal*10^b
            abPowTen = a*bPowTen # a*10^b
            N = intDiv(abPowTen,P)
            if recurseDepth: #recursion offset
                N = max(N+recurseOffset,1)

            PN = P*N
            c  = PN-abPowTen # c = P*N - a*10^b

            if ( # Verify conditions
                ( a==0 or c==0 ) or # a or c is zero
                ( (c<cRange[0] or c>cRange[1]) and (not recurseDepth) ) or  # c outside range
                ( math.gcd(a, c)>1 ) # a and c have common factor
            ): continue

            rules.append((P, N, a, b, c))

    # If no rules, reattempt. This rarely happens with default values
    # recurseDepth>0 activates recurseOffset and ignores cRange
    if not rules and recurseDepth<maxRecursion:
        recurseDepth+=1
        rules = findDivRulesByInverse(P, aRange, bRange, cRange, recurseDepth)
    return rules

#Output formatting
def printRuleTables(rules):
    headers = ["P", "N", "a", "b", "c", "P*N"]
    table = [[P, N, a, b, c, P * N] for P, N, a, b, c in rules]
    print(tabulate(table, headers=headers, tablefmt="plain"))

def printDivIters(iterTable):
    def getEquation(A, B, a, c):
        BSign = ""  if B*a  >=0 else "-"
        ASign = "+" if A*-c >=0 else "-"
        BPart=f"{abs(B)}*{abs(a)}"
        APart=f"{abs(A)}*{abs(c)}"
        return f"{BSign}{BPart}{ASign}{APart}"
    headers = ["A|B", "B*a-A*c", "Intermed"]
    table = [[f"{A}|{B}", getEquation(A, B, a, c), iterNum] for A, B, iterNum, a, c in iterTable]
    print(tabulate(table, headers=headers, tablefmt="plain"))

# Test function that reduces a number by using the given rules
def divByRule(rules,num):
    P,N,a,b,c=rules
    print("Example:")
    print(f"Using rules P={P} a={a} b={b} c={c}")
    print(f"Testing {num} divisibility by {P}\n")

    iterNum=num
    iterHistory=[num]
    iterTable=[]
    maxIters=100 #Very unlikely to reach +100 loop
    iterCount=0
    iterSmallest=num #best valid iteration

    while(iterCount<maxIters):
        iterCount+=1
        #Simulate splitting number at b:th digit
        powTen=10**b
        A=iterNum//powTen
        B=iterNum%powTen
        iterNum=B*a-A*c
        if abs(iterNum)<abs(iterSmallest): iterSmallest=iterNum
        #Keep history to check whether the iteration is in a loop
        if iterNum in iterHistory: break
        iterHistory.append(iterNum)
        iterTable.append([A,B,iterNum,a,c])
        if A==0 or B==0: break #Break if invalid A or B

    printDivIters(iterTable)
    print(f"\nSmallest iteration {iterSmallest} = {P}*{iterSmallest//P}")

    falsePositive=0
    floorIter=(iterSmallest//P)*P
    if floorIter != iterSmallest:
        print(f"Not divisible by {P}")
        print(f"Remainder of {num}/{P} is {num%P}")
        if num%P==0: falsePositive=1
    else:
        print(f"{num} is divisible by {P}")
        if num%P!=0: falsePositive=1

    if falsePositive:
        print("Something's not right")
        print("Bad rule perhaps, P*N=a*10^b+c ")
        print("doesn't always hold for highly ")
        print("composite numbers ")

    return iterNum

#Simple rotate right LCG for deterministic results
def rorLCG(seed):
    modulus=0xFFFFFFFFFFFFFFFF #2^64-1
    bits=modulus.bit_length()
    shift=int((bits+1)/3)
    randNum=(seed>>shift)|((seed<<(bits-shift))&modulus) #RORR
    randNum=(randNum*3343+11770513)&modulus #LCG
    return randNum

def randDivRuleTest(rules,isDivisible):
    minSize=10**(rules[3]+1)+17 # at least 10^b
    maxSize=10**(2*rules[3])+minSize
    P=rules[0]
    randPN=rorLCG(P)
    while(randPN<maxSize):
        randPN=randPN<<64|rorLCG(randPN) #cocat random values till maxSize
    randPN=randPN%maxSize+minSize
    if isDivisible:
        mult=(randPN%maxSize+minSize )//P
        randPN=P*mult # random product of P
    else:
        while abs(P)>1: # 1 would always be divisible by itself
            if randPN%P!=0: break # ensure the test number is not product of P
            randPN+=rorLCG(randPN)%(maxSize)+minSize

    num=divByRule(rules,randPN)

def parkerPrinting(rules):
    def ordinal(num):
        if 10 <= num % 100 <= 20:  # special cases 11th, 12th...
            return f"{num}th"
        return f"{num}{['th', 'st', 'nd', 'rd', 'th'][min(num % 10, 4)]}"
    P,N,a,b,c=rules
    bStr = ordinal(b)

    APart="A"
    BPart="B"
    print(f"{P} has following divisibility rule using B*a-A*c")
    print(f"Split the tested number into {APart} and {BPart} after {bStr} digit.")

    aStr = f"{abs(a)}" if abs(a) != 1 else ""
    cStr = f"{abs(c)}" if abs(c) != 1 else ""

    if cStr:
        cStr=f"{cStr}" # add only non 1 mult strings
        print(f"Multiply {APart} by {cStr} ",end="") # has c mult

    if aStr and cStr:
        print("and multiply ",end="") # has c and a mult
    elif aStr:
        print("Multiply ",end="") # has only a mult

    if aStr:
        print(f"{BPart} by {aStr}",end="")

    if aStr or cStr: #Don't print new line if both multipliers are 1
        print("",end="\n")

    #B*a-A*c sign
    BNum,ANum=a,-c

    if (ANum<0 and BNum<0): # if both are negative, make them positive
        ANum,BNum = abs(ANum),abs(BNum)

    esign=""
    if ANum>=0 and BNum>=0: #if added
        esign="+"
        print(f"Add {APart} and {BPart} together",end="")
    else: #if subtracted
        esign="-"
        if ANum>0:
            addend,subbed= (APart,BPart) # A minus B
        else:
            addend,subbed= (BPart,APart) # B minus A
        print(f"Subtract {subbed} from {addend}",end="")

    # add non 1 mult strings
    APart+= f"*{cStr}" if cStr else ""
    BPart+= f"*{aStr}" if aStr else ""

    if ANum>0: # first part is positive, it goes first
        equation=f"{APart}{esign}{BPart}"
    else:
        equation=f"{BPart}{esign}{APart}"

    print(f" = {equation}")


def main():
    @dataclass
    class Config:
        P: (int,int)
        PDefault: (int,int)
        useExample: bool
        isDivisible: bool
        parkerStyle: bool
        ruleCount: int
        ruleIndex: int

    #Configuration
    cfg = Config(
        P           = [None, None],
        PDefault    = [313 , None],   # [start,end) Any non-zero number which divisibility rule is
        useExample  = True,  # Verify rule by iterating B*a-A*c
        isDivisible = True,  # If tested example divided num is products of P
        parkerStyle = True,  # Parker style printing
        ruleCount   = 10,    # Limit shown rule count
        ruleIndex   = 0      # Index of used rule
    )
    if (not cfg.PDefault[1]): cfg.PDefault[1]=cfg.PDefault[0]+1

    # print help
    print("")
    if len(sys.argv) > 1:
        arg = str(sys.argv[1])
        if arg == "-h" or arg == "-help":
            spacerCount=len(str(f"python {sys.argv[0]} "))-1
            spacer=" "*spacerCount
            print(f"python {sys.argv[0]} [Integer or \"(Start,End)\"] [Show example? (0,1)] ")
            print(spacer,f"[Example is divisible? (0,1)] [Parker style? (0,1)] [Rule count] [Rule index]")
            print(f"Default command : python {sys.argv[0]} {cfg.PDefault[0]} {cfg.useExample} {cfg.isDivisible} {cfg.parkerStyle} {cfg.ruleCount} {cfg.ruleIndex}")
            print(f"Range example   : python {sys.argv[0]} \"[1,11]\" 0 0 1 2 0\n")
            return 0

    # Argument parsing
    # Detect if input is a single number or range tuple
    if len(sys.argv) > 1:
        arg1 = eval(sys.argv[1])
        try:
            if len(arg1)>1:
                cfg.P=arg1
        except:
            cfg.P=[arg1,arg1+1]

    # set unset values
    if (not cfg.P[0]): cfg.P[0]=cfg.PDefault[0]
    if (not cfg.P[1]): cfg.P[1]=cfg.P[0]+1

    if (cfg.P[1]-cfg.P[0]<=0): #invalid range
        cfg.P=cfg.PDefault

    cfg.useExample  = eval(sys.argv[2]) if len(sys.argv) > 2 else cfg.useExample
    cfg.isDivisible = eval(sys.argv[3]) if len(sys.argv) > 3 else cfg.isDivisible
    cfg.parkerStyle = eval(sys.argv[4]) if len(sys.argv) > 4 else cfg.parkerStyle
    cfg.ruleCount   = eval(sys.argv[5]) if len(sys.argv) > 5 else cfg.ruleCount
    cfg.ruleIndex   = eval(sys.argv[6]) if len(sys.argv) > 6 else cfg.ruleIndex

    for P in range(cfg.P[0],cfg.P[1]):
        if P==0: continue
        # 'a' and 'c' are multipliers in B*a-A*c, the smaller they are
        # the easier they will be to calculate by hand
        # 'b' determines where the tested number split right to left,
        #
        # Chosen a,b,c reasoning:
        # checking a>9 is redundant, these values will be checked by next b (a*10^b).
        # b 1-3x the magnitude of (small) P are doable by hand. 'b' sets the split 'B' magnitude.
        # increase c when numbers grows as it becomes harder to find divisibility rules.
        # if P has factors 2 or 5
        #
        # Compute time is limited by the number of digits in P due to division,
        # In terms of P, O( log(P)^1.466 * log(log(P)) * bRange^aRange )
        #
        digitsP=countDigits(P)
        bBase = digitsP # search b around P digit count
        bOffset = 10
        cBase = abs(P)//20+10
        aRange = 10
        bRange = [ bBase-bOffset, bBase+bOffset ]
        cRange = [ -cBase, cBase ]
        #limits
        bRange[0] = max(bRange[0],digitsP)
        bRange[1] = min(bRange[1],2*bBase)

        #Computation starts
        rules = findDivRulesByInverse(P, aRange, bRange, cRange)
        rules.sort(key=lambda x: (x[3], x[2]+abs(x[4]))) # sort by b, then a+c

        #Print
        if rules:
            rulesLen=len(rules)
            cfg.ruleIndex=max(0,min(cfg.ruleIndex,rulesLen-1))

            ruleWord = "rules" if rulesLen!=1 else "rule"
            print(f"Found {rulesLen} {ruleWord} for {P}, showing first {cfg.ruleCount}:")
            printRuleTables(rules[:cfg.ruleCount])  # Display first 10 rules
            print("")
            if cfg.parkerStyle:
                parkerPrinting(rules[cfg.ruleIndex])
                print("")
            if cfg.useExample: #random test
                randDivRuleTest(rules[cfg.ruleIndex],cfg.isDivisible)
        else:
            print(f"No rules found for {P}",end="")

        print("")

if __name__=="__main__":
    main()