Skip to content

Instantly share code, notes, and snippets.

@foxqstm
Created September 4, 2020 23:44
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save foxqstm/8f5734eccc1f5dafa8ed8ec7d04855da to your computer and use it in GitHub Desktop.
Save foxqstm/8f5734eccc1f5dafa8ed8ec7d04855da to your computer and use it in GitHub Desktop.
factorization ver1.ipynb
Display the source blob
Display the rendered blob
Raw
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
@myaosato
Copy link

myaosato commented Sep 5, 2020

複数プロセス試したコード貼り付けておきます。

import math
import time
import sys
from multiprocessing import Process, Value

plist = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103,
107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223,
227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347,
349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607,
613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743,
751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883,
887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031,
1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279,
1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423,
1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523,
1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627,
1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777,
1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907,
1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039,
2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,
2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309
]

cplist0 = [1, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107,109,
113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 169, 173, 179, 181, 191, 193, 197, 199, 211, 221, 223,
227, 229, 233, 239, 241, 247, 251, 257, 263, 269, 271, 277, 281, 283, 289, 293, 299, 307, 311, 313, 317,
323, 331, 337, 347, 349, 353, 359, 361, 367, 373, 377, 379, 383, 389, 391, 397, 401, 403, 409, 419, 421,
431, 433, 437, 439, 443, 449, 457, 461, 463, 467, 479, 481, 487, 491, 493, 499, 503, 509, 521, 523, 527,
529, 533, 541, 547, 551, 557, 559, 563, 569, 571, 577, 587, 589, 593, 599, 601, 607, 611, 613, 617, 619,
629, 631, 641, 643, 647, 653, 659, 661, 667, 673, 677, 683, 689, 691, 697, 701, 703, 709, 713, 719, 727,
731, 733, 739, 743, 751, 757, 761, 767, 769, 773, 779, 787, 793, 797, 799, 809, 811, 817, 821, 823, 827,
829, 839, 841, 851, 853, 857, 859, 863, 871, 877, 881, 883, 887, 893, 899, 901, 907, 911, 919, 923, 929,
937, 941, 943, 947, 949, 953, 961, 967, 971, 977, 983, 989, 991, 997, 1003, 1007, 1009, 1013, 1019, 1021,
1027, 1031, 1033, 1037, 1039, 1049, 1051, 1061, 1063, 1069, 1073, 1079, 1081, 1087, 1091, 1093, 1097,
1103, 1109, 1117, 1121, 1123, 1129, 1139, 1147, 1151, 1153, 1157, 1159, 1163, 1171, 1181, 1187, 1189,
1193, 1201, 1207, 1213, 1217, 1219, 1223, 1229, 1231, 1237, 1241, 1247, 1249, 1259, 1261, 1271, 1273,
1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1313, 1319, 1321, 1327, 1333, 1339, 1343, 1349,
1357, 1361, 1363, 1367, 1369, 1373, 1381, 1387, 1391, 1399, 1403, 1409, 1411, 1417, 1423, 1427, 1429,
1433, 1439, 1447, 1451, 1453, 1457, 1459, 1469, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1501, 1511,
1513, 1517, 1523, 1531, 1537, 1541, 1543, 1549, 1553, 1559, 1567, 1571, 1577, 1579, 1583, 1591, 1597,
1601, 1607, 1609, 1613, 1619, 1621, 1627, 1633, 1637, 1643, 1649, 1651, 1657, 1663, 1667, 1669, 1679,
1681, 1691, 1693, 1697, 1699, 1703, 1709, 1711, 1717, 1721, 1723, 1733, 1739, 1741, 1747, 1751, 1753,
1759, 1763, 1769, 1777, 1781, 1783, 1787, 1789, 1801, 1807, 1811, 1817, 1819, 1823, 1829, 1831, 1843,
1847, 1849, 1853, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1891, 1901, 1907, 1909, 1913, 1919, 1921,
1927, 1931, 1933, 1937, 1943, 1949, 1951, 1957, 1961, 1963, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
2011, 2017, 2021, 2027, 2029, 2033, 2039, 2041, 2047, 2053, 2059, 2063, 2069, 2071, 2077, 2081, 2083,
2087, 2089, 2099, 2111, 2113, 2117, 2119, 2129, 2131, 2137, 2141, 2143, 2147, 2153, 2159, 2161, 2171,
2173, 2179, 2183, 2197, 2201, 2203, 2207, 2209, 2213, 2221, 2227, 2231, 2237, 2239, 2243, 2249, 2251,
2257, 2263, 2267, 2269, 2273, 2279, 2281, 2287, 2291, 2293, 2297, 2309
]

def is_prime(N):
    if N<2310:
        if N in plist:
            return True
        else:
            return False
    else:
        for p in plist:
            if N%p==0:
                return False
        
        bound=math.floor(N**0.5)
        for cp in cplist:
            maxk=math.floor((bound-cp)/2310)
            for k in range(1,maxk+1):
                q=2310*k+cp
                if q>bound:
                    break
                if N%q==0:
                    return False

        return True

def makecp2310(n,cplist0):
    cplist=[]
    maxl=n//2310
    for l in range(0,maxl+1):
        l2310=2310*l
        for cp0 in cplist0:
            cp=l2310+cp0
            if cp>=n:
                break
            if math.gcd(cp,n)==1:
                cplist.append(cp)
    return cplist
    
    
def makecp(n):
    cpnlist=[]
    for i in range(1,n):
        if math.gcd(i,n)==1:
            cpnlist.append(i)
    return cpnlist
  
def intsqrt(n):
    sqnf=math.sqrt(n)
    dig=math.ceil(math.log10(sqnf))
    sq=0
    for k in reversed(range(0,dig+1)):
        for l in range(0,10):
            sq+=10**k
            if sq*sq==n:
                break
            if sq*sq>n:
                sq-=10**k
                break
    return sq
    
    
textn=input("n=")
n=int(textn)
rsa100=1522605027922533360535618378132637429718068114961380688657908494580122963258952897654000350692006139
ntest=12709189
ntest2=391961320150251294647
ntest3=20028641525215261
#n=rsa100
n=ntest2
p100=37975227936943673922808872755445627854565536638199

sf=time.time()

Nmax=30
countp=0
countmax=10**10
#countmax=5*10**6
pplist=[]

for N in range(2,min(2310,Nmax)+1):
    if is_prime(N):
        pplist.append(N)
        countp+=1
        if countp==countmax:
            break         

print('N=',N)
print('π(N)=',countp)
plist=pplist
pplist=[]

sqn=intsqrt(n)
sqn4=intsqrt(sqn)
dig=math.log10(sqn4)
primorial=1
countp=-1
for p in plist:
    countp+=1
    primorial*=p
    if primorial>max(2,10**(dig)):
        break
print('primo,count',primorial,countp)
m=primorial
print('m',m)
#for k in range(1,p):
#    m=primorial*k
#    if m>max(2,10**dig):
#        break


sqn4=intsqrt(sqn)
if m>=2310:
    cplist=makecp2310(m,cplist0)
else:
    cplist=makecp(m)
#print(cplist)

#percent=0
#delta=0.01
minc=0
c=minc
print('minc',minc)
maxc=sqn//m
print('maxc',maxc)
BOOL=False
for a in cplist:
    if a!=1 and n%a==0:
        BOOL=True
        break

#for c in range(1,maxc+1):
#    cm=c*m
#    if c/maxc>percent:
#        print(round(percent*100,2),'%',c)
#        percent+=delta
para = 8
def proc(b, flg):
    for d in range(0, math.ceil((maxc + 1) / para)):
        c = b + d * para
        cm=c*m
        for a in cplist:
            if flg.value == 1:
                return
            if n%(cm+a)==0:
                p=cm+a
                ff=time.time()
                print('full time',ff-sf,'[s]')            
                print('p',p)
                if n==1:
                    print('n=1 is first natural number.')
                else:
                    flg.value = 1
                    q=n//p
                    if q!=p:
                        print('q-p,log10(q-p)',q-p,math.log10(q-p))
                    textp=str(p)
                    textq=str(q)
                    print('product:',textn+'='+textp+'*'+textq)
                    if p*q==n:
                        print('Ans:',n,'=',p*q,'=',p,'*',q)                    
                    return
flg = Value('B', 0, lock = False)
ps = []
for b in range(1, para + 1):
    p = Process(target=proc, args=(b, flg))
    p.start()
    ps.append(p)
for p in ps:
    p.join()
if flg.value == 1:
   sys.exit()
    #for a in cplist:
    #    if n%(cm+a)==0:
    #        p=cm+a
    #        ff=time.time()
    #        print('full time',ff-sf,'[s]')
    #        
    #        print('p',p)
    #
    #        if n==1:
    #            print('n=1 is first natural number.')
    #        else:
    #            path='./'
    #            filename='RSA.txt'
    #            q=n//p
    #            if q!=p:
    #                print('q-p,log10(q-p)',q-p,math.log10(q-p))
    #            with open(path+filename,mode='w') as f:
    #                textp=str(p)
    #                textq=str(q)
    #                print('product:',textn+'='+textp+'*'+textq)
    #                if p*q==n:
    #                    print('Ans:',n,'=',p*q,'=',p,'*',q)
    #                f.write(textn+'='+textp+'*'+textq)
    #                sys.exit()

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment