こちらのqiita記事でも言及されております。フロイドの循環検出法に基づく素因数発見のための乱択式アルゴリズムとかいう意味不明な単語が並びますが、いくら読んでもまったく理解できないのでそういうもんなんだと割り切って先に進みましょう。素因数分解と割り切るをかけた秀逸な表現です。
""" ポラード・ロー法で素因数分解 """ target = 36610051291281 seed = 0 # 最大公約数 def gcd(a, b): if a<b: return gcd(b, a) while b>0: a, b = b, a%b return a # 素数判定 def isprime(n): if n<2: return False for m in range(2, int(n**0.5+1)): if n%m==0: return False return True # 疑似乱数発生器(らしい) def f(n): ret = n**2+seed return ret%target # 素数になるまで割る? def factoring(n, _seed = 3): if isprime(n): return n x, y, d = 2, 2, 1 tmp = n seed = _seed while d==1: x = f(x) y = f(f(y)) dif = x-y if x>y else y-x d = gcd(dif, tmp) seed+=1 seed = 0 if(d%2==0): return 2 if isprime(d) is False: return factoring(d, seed+1) return factoring(d, 1) # メイン def factors(): dic = {} tmp = target if(isprime(tmp)): dic[tmp] = 1 return dic ret = factoring(tmp) while tmp/ret is not 0: if ret in dic.keys(): dic[ret]+=1 else: dic[ret]=1 tmp = int(tmp/ret) if tmp==1: break ret = factoring(tmp) return dic def main(): print(factors()) if __name__=='__main__': main()
>>> {653: 1, 593783: 1, 3: 3, 13: 1, 269: 1}
ウルフラム先生に聞いてみるとたしかにあっているみたいです。
あと25桁くらいの数値に対してやると、間違った結果が出力されました。具体的には3661005129128136162001001227
とかを渡してやると、先生と異なる答えが出ました。このへんは割り算してるところの精度っぽい気がするんですが現時点ではわからないです。。。