Better divisibility rule finder

import math
import sys
from tabulate import tabulate

# 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
# Matt Parker explanation https://youtu.be/6pLz8wEQYkA


def countDigits(n): #floor(log10(n))+1
    return int( n.bit_length() * math.log10(2) ) + 1

def getNthReciprocal(n,pos):
    digits=countDigits(n)
    powTen=(10**(digits+pos)) #shift decimal
    digit=(powTen//n)%10 #get last digit
    return digit

def getReciprocalDigits(n,count): #get non-zero digits of 1/n
    array=[]
    for i in range(count):
        digit=getNthReciprocal(n,i)
        array.append(digit)

    return array

def getNthRound(array,n): #round element if next elem is >=5
    if n>=len(array)-1: return array[n] #last element
    isRoundUp = int(array[n+1]>=5)
    return array[n]+isRoundUp

def getRoundReciprocal(array,n0,n1):
    n1-=1 #compensate for the last number
    num=0
    for i in range(n0,n1):
       num=num*10+array[i]
    num=num*10+getNthRound(array,n1)
    return num

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

#Find divisibility rules P * N = a * 10^b + c
def findDivRulesByInverse(P, aRange, bRange, cRange):
    rules = []

    reciList=getReciprocalDigits(P,bRange[1]+1) #non-zero reciprocal digit list
    magP=countDigits(P)-1
    for b in range(bRange[0], bRange[1]):
        bPowTen=10**b # 10^b
        magReci=b-magP+1 #magnitude such that P*reci magnitude will be same as 10^(b+1) # P+r=b+1
        reci=getRoundReciprocal(reciList,0,magReci)

        scalerMax=reci #next b will cover numbers when scaler>scalerMax
        scaler=1
        for i in range(aRange[0],aRange[1]):
            scaler=intDiv(i*bPowTen,P) #Find nearest multiples that results in a*P=a*10^n+c where c is a small and n is near b
            if(scaler>scalerMax):break
            PN = P*scaler

            #Calculate rules
            a=intDiv(PN,bPowTen)
            c=PN-a*bPowTen
            N=scaler
            if(c>cRange[0] and c<cRange[1]):
                rules.append((P, N, a, b, c))

    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"))  # Use "plain", "grid", or others

def printDivIters(iterTable):
    headers = ["A.B", "B*a-A*c", "Intermed"]
    table = [[f"{A}.{B}", f"{B}*{a}-{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(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
        if iterNum<iterSmallest: iterSmallest=iterNum
        #Simulate splitting number at b:th digit
        powTen=10**b
        A=iterNum//powTen
        B=iterNum%powTen
        if A==0 or B==0: break #Break if invalid A or B
        iterNum=B*a-A*c
        iterTable.append([A,B,iterNum,a,c])

        #Keep history to check whether the iteration is in a loop
        if iterNum in iterHistory: break
        iterHistory.append(iterNum)

    printDivIters(iterTable)
    print("\n")

    factor=iterSmallest/P
    print(f"Smallest iteration {iterSmallest} = {P}*{factor:g}")
    wholeNum=math.floor(factor)
    if factor != wholeNum:
        print(f"Not divisible by {P}")
        print(f"Remainder of {num}/{P} is {num%P}")
    else:
        print(f"{num} is divisible by {P}")

    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,testProductsP):
    print("\nExample:")

    P=rules[0]
    randNum=rorLCG(P) #use the prime as seed

    minSize=P+17
    maxSize=P+P//2+minSize

    randPN=1
    if testProductsP:
        mult=randNum%maxSize+minSize
        randPN=P*mult # random product of P
    else:
        randPN=randNum%(maxSize**2)+minSize**2
        if randPN%P==0: randPN+=1 # ensure the test number is not product of P

    num=divByRule(rules,randPN)

def main():
    #Configuration
    defaultP = 313 # number which divisibility rule is searched
    P = int(sys.argv[1]) if len(sys.argv) > 1 else defaultP
    if P==0: P=defaultP

    testBestRule  = True # Verify rule by iterating B*a-A*c
    testProductsP = True # Whether we should test divisibility rules with known products of P
    maxDisplay=10 # Limit shown rule count

    # '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,
    # which sets the viable magnitude of the tested numbers
    # e.g. b=3 works well with numbers 10^4 to 10^7
    #
    # Chosen a,b,c reasoning:
    # checking a>9 is redundant these values will be checked by next b (a*10^b)
    # b 2-4x the magnitude of P are doable by hand. 'b' sets the split 'B' magnitude
    # increase c when numbers grows as it becomes harder to find divisibility rules
    #
    magP=countDigits(P)
    aRange = [ 1, 10             ]
    bRange = [ magP, magP*2+1 ]
    cRange = [ -abs(P//20), abs(P//20) ]
    #upper bound minimums
    if bRange[1]<5  : bRange[1]=5
    if cRange[1]<16 : cRange[1]=16

    #Computation starts
    rules = findDivRulesByInverse(P, aRange, bRange, cRange)
    rules.sort(key=lambda x: (x[3], abs(x[4])))

    #Print
    if rules:
        print(f"Found {len(rules)} rules for {P}, showing first {maxDisplay}:\n")
        printRuleTables(rules[:maxDisplay])  # Display first 10 rules
    else:
        print(f"No rules found for {P}")
        return 0

    if testBestRule:
        randDivRuleTest(rules[0],testProductsP)


if __name__=="__main__":
    main()