Advertisement
Kaelygon

Terrible python code divisibility test

Nov 22nd, 2024 (edited)
396
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 4.53 KB | None | 0 0
  1. import math
  2. import sys
  3. from tabulate import tabulate
  4.  
  5. # This terrible python code is provided by Kaelygon
  6. #
  7. # Find divisibility rule by prime P: P * N = a * 10^b + c
  8. # such that an integer can be split in two parts A and B
  9. # and reduced to a smaller number by simulating 10^b division and modulus
  10. # using string manipulation and addition: B*a-A*c
  11. # Matt Parker explanation https://youtu.be/6pLz8wEQYkA
  12.  
  13. #Find divisibility rules P * N = a * 10^b + c
  14. def findDivRules(P, aRange, bRange, cRange):
  15.     rules = []
  16.    
  17.     bMin=int(math.log10(abs(P))+1)
  18.     for b in range(bMin, bRange):
  19.         bPowTen=10**b # 10^b
  20.         for a in range(1, aRange):
  21.             ab = a*bPowTen # a*b^10
  22.             if ab<P: continue
  23.             for c in range(-cRange, cRange):
  24.                 PN = ab + c # P * N = a * 10^b + c
  25.                 if PN<=P or PN%P!=0 or c==0: continue # skip trivial solutions
  26.                 N = PN//P
  27.                 rules.append((P, N, a, b, c))
  28.    
  29.     # Sort by smallest b and c
  30.     rules.sort(key=lambda x: (x[3], abs(x[4])))
  31.    
  32.     return rules
  33.  
  34. #Output formatting
  35. def printRuleTables(rules):
  36.     headers = ["P", "N", "a", "b", "c", "P*N"]
  37.     table = [[P, N, a, b, c, P * N] for P, N, a, b, c in rules]
  38.     print(tabulate(table, headers=headers, tablefmt="plain"))  # Use "plain", "grid", or others
  39.  
  40. def printDivIters(iterTable):
  41.     headers = ["A.B", "B*a-A*c", "Intermed"]
  42.     table = [[f"{A}.{B}", f"{B}*{a}-{A}*{c}", iterNum] for A, B, iterNum, a, c in iterTable]
  43.     print(tabulate(table, headers=headers, tablefmt="plain"))
  44.    
  45.  
  46. # Test function that reduces a number by using the given rules
  47. def divByRule(rules,num):
  48.     P,N,a,b,c=rules
  49.     print(f"Using rules P={P} a={a} b={b} c={c}")
  50.     print(f"Testing {num} divisibility by {P}\n")
  51.    
  52.     iterNum=num
  53.     iterHistory=[num]
  54.     iterTable=[]
  55.     maxIters=100
  56.     iterCount=0
  57.     iterSmallest=num
  58.    
  59.     while(iterCount<maxIters): #Iterate till 10s table has been reached
  60.         iterCount+=1
  61.         if iterNum<iterSmallest: iterSmallest=iterNum
  62.        
  63.         powTen=10**b
  64.         A=iterNum//powTen
  65.         B=iterNum%powTen
  66.         if A==0 or B==0: break #Break if invalid A or B
  67.         iterNum=B*a-A*c
  68.         iterTable.append([A,B,iterNum,a,c])
  69.        
  70.         #Keep history to check whether the iteration is in a loop
  71.         if iterNum in iterHistory: break
  72.         iterHistory.append(iterNum)
  73.    
  74.     printDivIters(iterTable)
  75.     print("\n")
  76.    
  77.     factor=iterSmallest/P
  78.     print(f"Smallest iteration {iterSmallest} = {P}*{factor:g}")
  79.     wholeNum=math.floor(factor)
  80.     if factor != wholeNum:
  81.         print(f"Not divisible by {P}")
  82.         print(f"Remainder of {num}/{P} is {num%P}")
  83.     else:
  84.         print(f"{num} is divisible by {P}")
  85.        
  86.     return iterNum
  87.    
  88. #Simple rotate right LCG for deterministic results
  89. def rorLCG(seed):
  90.     modulus=0xFFFFFFFFFFFFFFFF #2^64-1
  91.     bits=modulus.bit_length()
  92.     shift=int((bits+1)/3)
  93.  
  94.     randNum=(seed>>shift)|((seed<<(bits-shift))&modulus) #RORR
  95.     randNum=randNum*3343+11770513 #LCG
  96.     randNum&=modulus
  97.     return randNum
  98.  
  99.  
  100. def randDivRuleTest(rules,testProductsP):
  101.     print("\nExample:")
  102.    
  103.     P=rules[0]
  104.     randNum=rorLCG(P) #use the prime as seed
  105.    
  106.     minSize=P+17
  107.     maxSize=P+P//2+minSize
  108.    
  109.     randPN=1
  110.     if testProductsP:
  111.         mult=randNum%maxSize+minSize
  112.         randPN=P*mult # random product of P
  113.     else:
  114.         randPN=randNum%(maxSize**2)+minSize**2
  115.         if randPN%P==0: randPN+=1 # ensure the test number is not product of P
  116.    
  117.     num=divByRule(rules,randPN)
  118.  
  119. def main():
  120.     defaultP = 313
  121.     P = int(sys.argv[1]) if len(sys.argv) > 1 else defaultP
  122.     if P==0: P=defaultP
  123.    
  124.     digitCount=int(math.log10(abs(P))+1)
  125.  
  126.     aRange = 10 # Second split 'B' multiplier
  127.     bRange = digitCount*2+1 # 'b' sets the tested number split right to left
  128.     cRange = abs(P//20) # First split 'A' multiplier
  129.  
  130.     #minimums
  131.     if bRange<5 : bRange=5
  132.     if cRange<16 : cRange=16
  133.    
  134.     maxDisplay=10
  135.     testProductsP = True # Whether we should test divisibility rules with known products of P
  136.  
  137.     #Computation starts
  138.     rules = findDivRules(P, aRange, bRange, cRange)
  139.  
  140.     if rules:
  141.         print(f"Found {len(rules)} rules for {P}, showing first {maxDisplay}:\n")
  142.         printRuleTables(rules[:maxDisplay])  # Display first 10 rules
  143.     else:
  144.         print(f"No rules found for {P}")
  145.         return 0
  146.  
  147.     randDivRuleTest(rules[0],testProductsP)
  148.    
  149.  
  150. if __name__=="__main__":
  151.     main()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement