Homeに戻る  一覧に戻る 

Factorization of 4^473+1


[2021.07.30]4^473+1の素因数分解


■問題

       4^473+1を素因数分解せよ。


■小さな素因数を見つける。

       m = 4473+1 = 22*11*43+1
           = 594806763391113225119224999259960224052504080663757783622308743726376262864161749418067325798462540235919489516077189220181834098217962283116332232440957850313188336178983949577074563933719094748095678312940574882427099482751152035262839576139463233204818042181657565129506139525873665
は、285桁の合成数である。
実際に10000以下の素数で割ってみることにより、5, 173, 397, 2113, 9461などの小さな素因数を持つことは直ぐに分かる。

[pari/gpによる計算]
gp> m=4^473+1
%1 = 594806763391113225119224999259960224052504080663757783622308743726376262864161749418067325798462540235919489516077189220181834098217962283116332232440957850313188336178983949577074563933719094748095678312940574882427099482751152035262839576139463233204818042181657565129506139525873665
gp> log(m)/log(10)
%2 = 284.77437589812621067219699440937040332
gp> forprime(p=2,10000,if(m%p==0,print(p)))
5
173
397
2113
9461
time = 1 ms.

■多項式を使って、いくつかの因数に分解する。

このまま、直ちに楕円曲線法(ECM)などで、素因数を探しても良いが、m-1が2の(2*11*43)べき乗数であることをうまく活用すると、多項式の因数分解を使って、2つまたは4つの因数に分解することができる。

例えば、mがm1=22+1, m2=22*11+1, m3=22*43+1で割り切れることは直ぐに分かる。
m1=5, m2=4194305, m3=77371252455336267181195265は、以下のように、簡単に素因数分解できる。

[pari/gpによる計算]
gp> m1=2^2+1
%4 = 5
gp> m2=2^(2*11)+1
%5 = 4194305
gp> m3=2^(2*43)+1
%6 = 77371252455336267181195265
gp> factor(m1)
%7 =
[5 1]

gp> factor(m2)
%8 =
[   5 1]

[ 397 1]

[2113 1]

gp> factor(m3)
time = 2 ms.
%9 =
[            5 1]

[          173 1]

[       101653 1]

[       500177 1]

[1759217765581 1]


さらに、2*11*43=236*4+2より、x=2236とすると、
       m = 2236*4+2+1 = 22*(2236)4+1
           = 4*x4+1= (2*x2)2-(2*x)2 = (2*x2+2*x+1)*(2*x2-2*x+1)
           = 24388660549343689307668728357759111763660922989570087116087163747073216930385302004487932628095718617550962377821476350823057075195214890532865 * 24388660549343689307668728357759111763660922989570087116087163747073216488673535809891850232271343431821333420950502131918317544793664567377921
となり、2つの因数(143桁、143桁の合成数)に分解できる。

[pari/gpによる計算]
gp> factor(4*x^4+1)
time = 1 ms.
%10 =
[2*x^2 - 2*x + 1 1]

[2*x^2 + 2*x + 1 1]

gp> f1(x)=2*x^2 - 2*x + 1;
gp> f2(x)=2*x^2 + 2*x + 1;
gp> m4=f1(2^236)
time = -8 ms.
%13 = 24388660549343689307668728357759111763660922989570087116087163747073216488673535809891850232271343431821333420950502131918317544793664567377921
gp> m5=f2(2^236)
%14 = 24388660549343689307668728357759111763660922989570087116087163747073216930385302004487932628095718617550962377821476350823057075195214890532865
gp> m4*m5-m
%15 = 0
gp> log(m4)/log(10)
%16 = 142.38718794906310533609849720468520166
gp> log(m5)/log(10)
%17 = 142.38718794906310533609849720468520166
同様に、2*11*43=21*44+22より、y=221とすると、
       m = 221*4+22+1 = 222*(221)44+1 = 222*y44+1 = (2*y4)11+1
           = (2*y2 - 2*y + 1)*(2*y2 + 2*y + 1)*(1024*y20 - 1024*y19 + 512*y18 - 256*y16 + 256*y15 - 128*y14 + 64*y12 - 64*y11 + 32*y10 - 32*y9 + 16*y8 - 8*y6 + 8*y5 - 4*y4 + 2*y2 - 2*y + 1)*(1024*y20 + 1024*y19 + 512*y18 - 256*y16 - 256*y15 - 128*y14 + 64*y12 + 64*y11 + 32*y10 + 32*y9 + 16*y8 - 8*y6 - 8*y5 - 4*y4 + 2*y2 + 2*y + 1)
           = 8796088827905 * 8796097216513 * 2772668372009192495130787400310267083262504493885805446596344843418484971869580234903437835283456674146985706429133952914257084417 * 2772671016233067656125183207971465038554089072943218983374364980946027072519597181173133823689844885090955523250102076564865810433
のように、4つの因数(13桁、13桁、130桁、130桁の合成数)に分解できる。
13桁の合成数8796088827905, 8796097216513は簡単に素因数分解できて、
       8796088827905 = 5 * 1759217765581
       8796097216513 = 173 * 101653 * 500177
となる。
残りの130桁の合成数は、それぞれ、2113*9641, 397を因数に持つことが分かる。

[pari/gpによる計算]
gp> 2*11*43-(21*44+22)
%18 = 0
gp> factor(2^22*x^44+1)
time = 4 ms.
%19 =
[2*x^2 - 2*x + 1 1]

[2*x^2 + 2*x + 1 1]

[1024*x^20 - 1024*x^19 + 512*x^18 - 256*x^16 + 256*x^15 - 128*x^14 + 64*x^12 - 64*x^11 + 32*x^10 - 32*x^9 + 16*x^8 - 8*x^6 + 8*x^5 - 4*x^4 + 2*x^2 - 2*x + 1 1]

[1024*x^20 + 1024*x^19 + 512*x^18 - 256*x^16 - 256*x^15 - 128*x^14 + 64*x^12 + 64*x^11 + 32*x^10 + 32*x^9 + 16*x^8 - 8*x^6 - 8*x^5 - 4*x^4 + 2*x^2 + 2*x + 1 1]

gp> f3(x)=1024*x^20 - 1024*x^19 + 512*x^18 - 256*x^16 + 256*x^15 - 128*x^14 + 64*x^12 - 64*x^11 + 32*x^10 - 32*x^9 + 16*x^8 - 8*x^6 + 8*x^5 - 4*x^4 + 2*x^2 - 2*x + 1;
gp> f4(x)=1024*x^20 + 1024*x^19 + 512*x^18 - 256*x^16 - 256*x^15 - 128*x^14 + 64*x^12 + 64*x^11 + 32*x^10 + 32*x^9 + 16*x^8 - 8*x^6 - 8*x^5 - 4*x^4 + 2*x^2 + 2*x + 1;
gp> m6=f1(2^21)
%22 = 8796088827905
gp> m7=f2(2^21)
%23 = 8796097216513
gp> m8=f3(2^21)
%24 = 2772668372009192495130787400310267083262504493885805446596344843418484971869580234903437835283456674146985706429133952914257084417
gp> m9=f4(2^21)
%25 = 2772671016233067656125183207971465038554089072943218983374364980946027072519597181173133823689844885090955523250102076564865810433
gp> m6*m7*m8*m9-m
%26 = 0
gp> factor(m6)
%27 =
[            5 1]

[1759217765581 1]

gp> factor(m7)
%28 =
[   173 1]

[101653 1]

[500177 1]

gp> forprime(p=2,10000,if(m8%p==0,print(p)))
2113
9461
time = 1 ms.
gp> forprime(p=2,10000,if(m9%p==0,print(p)))
397
time = 1 ms.
gp> c123=m8/(2113*9461)
%31 = 138695186501768187218717225731993097289002882127845908505170019639170553199346340637975013936629511660367229867278089943069
gp> c127=m9/397
%32 = 6984057975398155305101217148542733094594682803383423131925352596841378016422159146531823233475679811312230537153909512757848389
gp> [log(c123)/log(10),log(c127)/log(10)]
%33 = [122.14206138888874337745580319343503281, 126.84410783583654545006812786311199802]

これまでの結果をまとめると、
       m = 5 * 173 * 397 * 2113 * 9461 * 101653 * 500177 * 1759217765581 * c123 * c127
となる。ここで、c123, c127はそれぞれ、123桁、127桁の合成数
       c123 = 138695186501768187218717225731993097289002882127845908505170019639170553199346340637975013936629511660367229867278089943069
       c127 = 6984057975398155305101217148542733094594682803383423131925352596841378016422159146531823233475679811312230537153909512757848389
である。

■楕円曲線を使った素因数分解(ECM)で合成数c123を分解する。

kを小さな素因数を多く持つ大きな整数、ここでは、k=lcm(1,2,3,...,1000000)とする。
Z/(c123)Zは体ではないが、体と見なして、適当な楕円曲線の有理点Pのk倍点[k]Pを計算する。k倍点が運良く無限遠点Oになると、その計算が途中で失敗して、結果的に合成数c123の素因数が得られる。

楕円曲線y^2=x^2+4*x-15の有理点(2,1)のk倍点を計算することにより、素因数386299835479975297が求まった。
残りは105桁の合成数
       c105 = 359035064898332729835946471279657493150229553669855988585404007668454362020863096397816812546005247211677
である。

[ruby 3.0.1p64(Windows 10)による計算]
C:\Users\hisay>ruby ecm-c123.rb
c123=138695186501768187218717225731993097289002882127845908505170019639170553199346340637975013936629511660367229867278089943069
k=lcm{1,..,1000000}
gcd(k,c123)=1
time=234.6681946[s]
y^2=x^3+0x+(-7)
k(2 1)=(115561061577155966313608978083872842379499163178071523878858423297161610838232084984049490366978353702302041970708092535323 91399716991250720787089249790287438430639139905307042890156611934787782909758364002889186192979830532740144782852630873806)
time=1276.4495205[s]
y^2=x^3+1x+(-9)
k(2 1)=(66629339454325307350310229240612495025879201088595495559998093733872406275441224140436778371496783595007077746314305163531 29069565795019524019510270332436460127081440208283079415074421030886581263495864319653645054156637037632743656315545243260)
time=1582.9924357[s]
y^2=x^3+2x+(-11)
k(2 1)=(133278515380580928067770751970448701705256681429872248721870376186232074369378952689529388874796232918899872740168223802656 47310382084582998909515601460353602356938929658764941066204019059014099177925398667796089833961495821874126799588373215074)
time=1697.042612[s]
y^2=x^3+3x+(-13)
k(2 1)=(28270579550935010731047698093622866683842448175187745784161975764002107503988156481053239661406896960665436144232518480965 19503674126905102568484536242221347798036058482779250292904731511704721253098451787204991384737840462896900971965525388438)
time=1647.0948661[s]
y^2=x^3+4x+(-15)
k(2 1)=(386299835479975297 )
time=1715.6540649[s]
total time=8153.9244139[s]


■楕円曲線を使った素因数分解(ECM)で合成数c127を分解する。

ここでは、k=lcm(1,2,3,...,1000000)とする。
Z/(c127)Zは体ではないが、体と見なして、適当な楕円曲線の有理点Pのk倍点[k]Pを計算する。

楕円曲線y^2=x^2+x-9の有理点(2,1)のk倍点を計算することにより、素因数1096710113137が求まった。
残りは115桁の合成数
       c115 = 6368189635291266184454627784828723548091461501180570317457822569882380873375125854136653695158328730643828303109397
である。

[ruby 3.0.1p64(Windows 10)による計算]
C:\Users\hisay>ruby ecm-c123.rb
c123=138695186501768187218717225731993097289002882127845908505170019639170553199346340637975013936629511660367229867278089943069
k=lcm{1,..,1000000}
gcd(k,c123)=1
time=234.6681946[s]
y^2=x^3+0x+(-7)
k(2 1)=(115561061577155966313608978083872842379499163178071523878858423297161610838232084984049490366978353702302041970708092535323 91399716991250720787089249790287438430639139905307042890156611934787782909758364002889186192979830532740144782852630873806)
time=1276.4495205[s]
y^2=x^3+1x+(-9)
k(2 1)=(66629339454325307350310229240612495025879201088595495559998093733872406275441224140436778371496783595007077746314305163531 29069565795019524019510270332436460127081440208283079415074421030886581263495864319653645054156637037632743656315545243260)
time=1582.9924357[s]
y^2=x^3+2x+(-11)
k(2 1)=(133278515380580928067770751970448701705256681429872248721870376186232074369378952689529388874796232918899872740168223802656 47310382084582998909515601460353602356938929658764941066204019059014099177925398667796089833961495821874126799588373215074)
time=1697.042612[s]
y^2=x^3+3x+(-13)
k(2 1)=(28270579550935010731047698093622866683842448175187745784161975764002107503988156481053239661406896960665436144232518480965 19503674126905102568484536242221347798036058482779250292904731511704721253098451787204991384737840462896900971965525388438)
time=1647.0948661[s]
y^2=x^3+4x+(-15)
k(2 1)=(386299835479975297 )
time=1715.6540649[s]
total time=8153.9244139[s]

■楕円曲線を使った素因数分解(ECM)で合成数c105を分解する。

k=lcm(1,2,3,...,1000000)とする。
楕円曲線y^2=x^3+193*x+(-602)の有理点(3,2)のk倍点を計算することにより、c105の素因数588073600518602649955549が見つかった。
残りは81桁の合成数
       c81=610527431569299465158022167405255836855896092107776363873073943375109371498913473
である。

[ruby-2.7.4p191(NetBSD-8.0/amsd64)による計算]
-bash-3.1$ ruby num7-c105-p32.rb
c105=359035064898332729835946471279657493150229553669855988585404007668454362020863096397816812546005247211677
k=lcm{1,..,1000000}
gcd(k,c105)=1
time=590.237324275[s]
y^2=x^3+0x+(-23)
k(3 2)=(7099702687473996741128769378198168570605630566566122667736177434591977123797036995584769309476185435923 309363359651701004243288094167971790043682940937611520366246916733650418863461650924388758714599337398060)
time=1201.408261279[s]
y^2=x^3+1x+(-26)
k(3 2)=(294777817120162275572104663242920686952154403819744601028205095526355084070723919657079601232775062384864 225805846646277097119345328247166532623758818395419893388419037612307619121346091595407528539435949944656)
time=1117.497301089[s]
y^2=x^3+2x+(-29)
k(3 2)=(88681998061135065874769785655777176265277159122793698319927112514869676852348843127926339693214105669395 274979835345297939142733686057459851236518737390952229449077056694194803613817691512156874821827910317089)
time=738.26131451[s]
...(省略)...
y^2=x^3+191x+(-596)
k(3 2)=(135855250384393011144837583094017958563893485908289164328088854001377545641696248801610875808263793623781 110034261329822847030640187884204769938917896957232831895952476732899571752471795905833731298078777855931)
time=1079.37399331[s]
y^2=x^3+192x+(-599)
k(3 2)=(186517708984166503588564256685115080735053003810732029930121482924491155649862905911357210965934364452117 31341315090567182525033448582391837781063511956753972723069362506950479923613601509288478744029301770086)
time=776.953968356[s]
y^2=x^3+193x+(-602)
k(3 2)=(588073600518602649955549 )
time=860.214374358[s]

■楕円曲線を使った素因数分解(ECM)で合成数c115を分解する。

k=lcm(1,2,3,...,1000000)とする。楕円曲線y^2=x^3+88*x+(-287)の有理点(3,2)のk倍点を計算することにより、c115の素因数32272871788802004477006473が見つかった。

[ruby-2.7.4p191(NetBSD-8.0/amsd64)による計算]
-bash-3.1$ ruby num7-c115-p32.rb
c115=6368189635291266184454627784828723548091461501180570317457822569882380873375125854136653695158328730643828303109397
k=lcm{1,..,1000000}
gcd(k,c115)=1
time=556.147905592[s]
y^2=x^3+0x+(-23)
k(3 2)=(4791738625221646910095601410060544056301386504068296451959637506008074245123128492660984840811160236879609933170785 2863152332448155680687960067548247040819052137562587182051872348240685324083998941740706967927470871749315291174663)
time=838.254252362[s]
y^2=x^3+1x+(-26)
k(3 2)=(302502495318981654652259117134771266046790214142646250453910506981559774684370995546751372269685921390363807299828 5567171427803968515093186778093994849700254426588049198821041189131005725420050244498381562165175845097020353319048)
time=837.229152439[s]
y^2=x^3+2x+(-29)
k(3 2)=(1586944552203508745150997512403738576375569761509318075608515837387272637547659253219319520687755631059365259520925 2074794469473656163567202933409683723283695136369521017475324435661777229832073950707383281023114610727262283125875)
time=838.99791345[s]
y^2=x^3+3x+(-32)
k(3 2)=(2376139828197027077701728431347533466480127682505301947802554129036600766686618666675882074435807808347338932309781 5552106478762324533991655927437911131547382235126329638899312941352405885268836924189720436454897153548468066920152)
time=840.492821541[s]
...(省略)...
y^2=x^3+86x+(-281)
k(3 2)=(4283787121977966956530299717541155443169583791574976204988124854475352528303854196147617853166671454368833794058503 311322946977003251477269753083357897955527495330733413365374345859350743201345280476622294834406699234069242873316)
time=1161.296805405[s]
y^2=x^3+87x+(-284)
k(3 2)=(4562202951553084340852737961168480943271141374827591328211829615545430094990138303077365572375276867387741045482908 5123121946725180699548982985815570138730215366443257107229221753207579466066323333667455434173447400082968933763040)
time=1374.155029665[s]
y^2=x^3+88x+(-287)
k(3 2)=(32272871788802004477006473 )
time=1368.50509546[s]

並行して、楕円曲線y^2=x^3+91*x+(-83)の有理点(1,3)のk倍点を計算することにより、c115の別の素因数35303956875887883008321が見つかった。

[ruby-2.7.4p191(NetBSD-8.0/amsd64)による計算]
-bash-3.1$ ruby num7-c115-p13.rb
c115=6368189635291266184454627784828723548091461501180570317457822569882380873375125854136653695158328730643828303109397
k=lcm{1,..,1000000}
gcd(k,c115)=1
time=806.731861545[s]
y^2=x^3+0x+(8)
k(1 3)=(2739877078968013861708658491762285539807153462804043452911287508272671952429508534496806113819210300228720745104511 4502485554470859839244183153593109666810649191700751379488846801740170350203874650307436917592005268681058868606502)
time=1417.168185849[s]
y^2=x^3+1x+(7)
k(1 3)=(4149362662021261457659633273586250471319491087653714160903613551517555108741767437363909337754220259227629234203076 4887131674894981012509681680791843807081776675188027429337306131585101324575371991837374858272678953796988518802728)
time=1415.223096713[s]
y^2=x^3+2x+(6)
k(1 3)=(4953053593718158736356038669720397256756162598546056728349391638282901496964727479262619090286357851588557220043637 1261178794253973762681402498156549364755129049694792919978746803509220811695640233067654667549798385217444904182449)
time=1415.897923371[s]
...(省略)...
y^2=x^3+89x+(-81)
k(1 3)=(4051234985853779522109118660123017051939426065480952621134019423162689265423567285371465186097146539258342741741907 3972462659932294579158862411953981233358868818459459069591033863086673510825829533384239670540035563851010792402274)
time=944.348369466[s]
y^2=x^3+90x+(-82)
k(1 3)=(4606502766375207002196052245708105081150543910199760868214053680541038576594082548488637244035418609826090199896047 2467736246969210000435352515087571218893108427865898761442348687321013279557835793335459522233604727280895781386153)
time=951.034601548[s]
y^2=x^3+91x+(-83)
k(1 3)=(2674782178548971060435453123002848946252636083378460040030850456412490204581951493435737296460692590916682782706612 5433074863344394606018541984346720747994394077819540287713917576239151521713974145229476384723271492939544571530502)
time=1004.05304071[s]
y^2=x^3+92x+(-84)
k(1 3)=(35303956875887883008321 )
time=1425.705084757[s]
残りは67桁の素数5589268731817733787231040982654760345321756061873730705031153091309である(ECPPにより確認)ので、合成数c115は、
       c115 = 35303956875887883008321 * 32272871788802004477006473 * 5589268731817733787231040982654760345321756061873730705031153091309
に素因数分解できた。

■2次ふるい法で、合成数c81を分解する。
素因数が見つかるかどうかは、楕円曲線法(ECM)は完全に運任せであるが、2次ふるい法は、時間が掛かるものの、ほぼ確実に見つけることができる。

ここでは、参考文献[5]のmpq-0.1.tar.gz(2005-04-29)を利用する。

最初に、sieveで2次ふるいにかける。ここに一番時間がかかる(18663.7[s])。

[NetBSD-8.0/amd64上での実行結果]
-bash-3.1$ cat c81
610527431569299465158022167405255836855896092107776363873073943375109371498913473
500000
-bash-3.1$ ./sieve -p c81 -n 50000
Welcome to Multiple Polynomial Quadratic Sieve
using THRESH1=50 ERROR2=4
trying to factor:
  610527431569299465158022167405255836855896092107776363873073943375109371498913473
using M=600000 (10 blocks of length 60000)
factor base bound 500000
large prime bound 64000000
sieve threshold is 176776695296636
stop sieving after 50000 smooths
factor base primes will be put in file 'c81.fbase'
smooths will be put in file 'c81.fulls'
partials will be put in file 'c81.partials'
pps will be put in file 'c81.pps'
using multiplier k = 1
k*N = 610527431569299465158022167405255836855896092107776363873073943375109371498913473
sieve thresholds: 50 69
modified: 14 64 83
creating factor base primes file!
sieving begins at 223
Fri Jul 30 05:19:22 GMT 2021
PANTHER
time=0.2s: 0 ffs, 0 fps, 0 pps, est. 0 cycles (20892 wanted)
time=1.1s: 1 ffs, 36 fps, 142 pps, est. 0 cycles (20892 wanted)
time=2.0s: 6 ffs, 71 fps, 312 pps, est. 0 cycles (20892 wanted)
time=2.9s: 7 ffs, 104 fps, 454 pps, est. 0 cycles (20892 wanted)
time=3.8s: 10 ffs, 142 fps, 615 pps, est. 0 cycles (20892 wanted)
time=4.7s: 10 ffs, 177 fps, 778 pps, est. 0 cycles (20892 wanted)
time=5.7s: 13 ffs, 202 fps, 942 pps, est. 0 cycles (20892 wanted)
time=6.6s: 16 ffs, 226 fps, 1126 pps, est. 0 cycles (20892 wanted)
time=7.3s: 19 ffs, 252 fps, 1298 pps, est. 0 cycles (20892 wanted)
time=8.1s: 22 ffs, 289 fps, 1460 pps, est. 0 cycles (20892 wanted)
time=9.0s: 23 ffs, 319 fps, 1610 pps, est. 0 cycles (20892 wanted)
time=9.9s: 24 ffs, 353 fps, 1797 pps, est. 0 cycles (20892 wanted)
time=10.8s: 27 ffs, 385 fps, 1964 pps, est. 0 cycles (20892 wanted)
...(省略)...
time=18658.4s: 49982 ffs, 793986 fps, 3830467 pps, est. 482534 cycles (20892 wanted)
time=18659.1s: 49984 ffs, 794028 fps, 3830617 pps, est. 482585 cycles (20892 wanted)
time=18659.8s: 49985 ffs, 794066 fps, 3830758 pps, est. 482631 cycles (20892 wanted)
time=18660.5s: 49987 ffs, 794092 fps, 3830906 pps, est. 482663 cycles (20892 wanted)
time=18661.2s: 49991 ffs, 794120 fps, 3831064 pps, est. 482697 cycles (20892 wanted)
time=18661.9s: 49994 ffs, 794152 fps, 3831221 pps, est. 482736 cycles (20892 wanted)
time=18662.6s: 49996 ffs, 794187 fps, 3831380 pps, est. 482778 cycles (20892 wanted)
time=18663.3s: 49998 ffs, 794218 fps, 3831522 pps, est. 482816 cycles (20892 wanted)
sieving took 18663.7s (factor 5361.0s)
50000 smooth equations
794231 semi-smooth equations
3831601 partial partials

次に、combineにより、partial relationsを結合する(96.8[s])。

[NetBSD-8.0/amd64上での実行結果]
-bash-3.1$ ./combine -p c81
read 794231 fps, 3831601 pps
found 1338841 large primes appearing at least twice
found 255174 cycles from fps
kept 753259 fps and 3581769 pps
627 connected component(s)
found 2744152 cycles from pps, total 2999326
combining partials took 96.8s

さらに、buildmatrixにより、行列を構築する(183.0[s])。

[NetBSD-8.0/amd64上での実行結果]
-bash-3.1$ ./buildmatrix -p c81
factor base has 20892 primes
2999326 cycles
50000 full relations
3049326 total relations
remains 3049326 rels on 20893 primes (excess 3028433)
removed 3028423 heavy relations
remains 20903 rels on 20893 primes (excess 10)
discarded 1709 singleton(s)
remains 19194 rels on 18483 primes (excess 711)
discarded 404 singleton(s)
remains 18790 rels on 18061 primes (excess 729)
discarded 101 singleton(s)
remains 18689 rels on 17959 primes (excess 730)
discarded 19 singleton(s)
remains 18670 rels on 17940 primes (excess 730)
discarded 8 singleton(s)
remains 18662 rels on 17932 primes (excess 730)
discarded 1 singleton(s)
remains 18661 rels on 17931 primes (excess 730)
removed 720 heavy relations
remains 17941 rels on 17931 primes (excess 10)
discarded 247 singleton(s)
remains 17694 rels on 17678 primes (excess 16)
discarded 79 singleton(s)
remains 17615 rels on 17599 primes (excess 16)
discarded 27 singleton(s)
remains 17588 rels on 17572 primes (excess 16)
discarded 7 singleton(s)
remains 17581 rels on 17565 primes (excess 16)
removed 6 heavy relations
remains 17575 rels on 17565 primes (excess 10)
discarded 3 singleton(s)
remains 17572 rels on 17562 primes (excess 10)
discarded 1 singleton(s)
remains 17571 rels on 17561 primes (excess 10)
building matrix took 183.0s

最後に、gaussにより、行列のGauss消去法を実行し、dependenciesを解決する(42.5[s])。
[NetBSD-8.0/amd64上での実行結果]
-bash-3.1$ ./gauss -p c81
17571 relations
17561 primes with odd exponent
...............................................................................
Found 11 dependencies
Process dependency 1:
gcd(X+Y,N)=610527431569299465158022167405255836855896092107776363873073943375109371498913473
gcd(X-Y,N)=1
Process dependency 2:
gcd(X+Y,N)=610527431569299465158022167405255836855896092107776363873073943375109371498913473
gcd(X-Y,N)=1
Process dependency 3:
gcd(X+Y,N)=1
gcd(X-Y,N)=610527431569299465158022167405255836855896092107776363873073943375109371498913473
Process dependency 4:
gcd(X+Y,N)=610527431569299465158022167405255836855896092107776363873073943375109371498913473
gcd(X-Y,N)=1
Process dependency 5:
gcd(X+Y,N)=1
gcd(X-Y,N)=610527431569299465158022167405255836855896092107776363873073943375109371498913473
Process dependency 6:
gcd(X+Y,N)=1
gcd(X-Y,N)=610527431569299465158022167405255836855896092107776363873073943375109371498913473
Process dependency 7:
gcd(X+Y,N)=368596695116186434714529041561
gcd(X-Y,N)=1656356227982058657394041821129733144353123888780393
Gaussian elimination took 42.5s

2次ふるい法により、30桁の因数368596695116186434714529041561と52桁の因数1656356227982058657394041821129733144353123888780393が見つかり、いずれも素数である(ECPPにより確認)ので、81桁の合成数は、
       c81 = 368596695116186434714529041561 * 1656356227982058657394041821129733144353123888780393
のように素因数分解できた。

■以上をまとめると、285桁の合成数m = 4473+1は、以下のように、素因数分解できた。

       m = 5 * 173 * 397 * 2113 * 9461 * 101653 * 500177 * 1096710113137 * 1759217765581 * 386299835479975297 * 35303956875887883008321 * 32272871788802004477006473 * 368596695116186434714529041561 * 1656356227982058657394041821129733144353123888780393 * 5589268731817733787231040982654760345321756061873730705031153091309
となる。
各因子が素数であることは、ECPPなどの方法により、簡単に確認できる。

[pari/gpによる計算]
gp> 5*173*397*2113*9461*101653*500177*1096710113137*1759217765581*386299835479975297*35303956875887883008321*588073600518602649955549*32272871788802004477006473*368596695116186434714529041561*1656356227982058657394041821129733144353123888780393*5589268731817733787231040982654760345321756061873730705031153091309
%36 = 594806763391113225119224999259960224052504080663757783622308743726376262864161749418067325798462540235919489516077189220181834098217962283116332232440957850313188336178983949577074563933719094748095678312940574882427099482751152035262839576139463233204818042181657565129506139525873665
gp> 4^473+1
%37 = 594806763391113225119224999259960224052504080663757783622308743726376262864161749418067325798462540235919489516077189220181834098217962283116332232440957850313188336178983949577074563933719094748095678312940574882427099482751152035262839576139463233204818042181657565129506139525873665
gp> %36-%37
%38 = 0
gp> isprime(101653,3)
%39 = 1
gp> isprime(500177,3)
%40 = 1
gp> isprime(1096710113137,3)
%41 = 1
gp> isprime(1759217765581,3)
%42 = 1
gp> isprime(386299835479975297,3)
%43 = 1
gp> isprime(35303956875887883008321,3)
time = 55 ms.
%44 = 1
gp> isprime(588073600518602649955549,3)
time = 43 ms.
%45 = 1
gp> isprime(32272871788802004477006473,3)
time = 57 ms.
%46 = 1
gp> isprime(368596695116186434714529041561,3)
time = 36 ms.
%47 = 1
gp> isprime(1656356227982058657394041821129733144353123888780393,3)
time = 66 ms.
%48 = 1
gp> isprime(5589268731817733787231040982654760345321756061873730705031153091309,3)
time = 181 ms.
%49 = 1

[2023.03.08追記]
数体篩法(Number Field Sieve method)のプログラムmsieve-1.46を使って、c81を素因数分解してみた。6313335個の関係データを集める処理に約3日掛かったが、問題なく30桁の素因数368596695116186434714529041561と52桁の素因数1656356227982058657394041821129733144353123888780393の積に分解できた。

[mesieve 1.46による計算結果]
-bash-3.1$ ./msieve -h

Msieve v. 1.46

usage: ./msieve [options] [one_number]

numbers starting with '0' are treated as octal,
numbers starting with '0x' are treated as hexadecimal

options:
   -s  save intermediate results to 
             instead of the default msieve.dat
   -l  append log information to 
             instead of the default msieve.log
   -i  read one or more integers to factor from
              (default worktodo.ini) instead of
             from the command line
   -m        manual mode: enter numbers via standard input
   -mb  hint for number of megabytes of memory for
             postprocessing (set automatically if unspec-
             ified or zero)
   -q        quiet: do not generate any log information,
             only print any factors found
   -d   deadline: if still sieving after 
             minutes, shut down gracefully (default off)
   -r   stop sieving after finding  relations
   -p        run at idle priority
   -v        verbose: write log information to screen
             as well as to logfile
   -t   use at most  threads

 elliptic curve options:
   -e        perform 'deep' ECM, seek factors > 15 digits

 quadratic sieve options:
   -c        client: only perform sieving

 number field sieve options:
   -n        use the number field sieve (80+ digits only;
             performs all NFS tasks in order)
   -nf  read from / write to NFS factor base file
              instead of the default msieve.fb
   -np [X,Y] perform only NFS polynomial selection; if
             specified, only cover leading coefficients
             in the range from X to Y inclusive
   -np1 [X,Y] perform stage 1 of NFS polynomial selection; if
             specified, only cover leading coefficients
             in the range from X to Y inclusive
   -np2      perform stage 2 of NFS polynomial selection
   -ns [X,Y] perform only NFS sieving; if specified,
             handle sieve lines X to Y inclusive
   -nc       perform only NFS combining (all phases)
   -nc1 [X,Y] perform only NFS filtering. Filtering will
             track ideals >= X (determined automatically
             if 0 or unspecified) and will only use the
             first Y relations (or all relations, if 0
             or unspecified)
   -nc2      perform only NFS linear algebra
   -ncr      perform only NFS linear algebra, restarting
             from a previous checkpoint
   -nc3 [X,Y] perform only NFS square root (compute
             dependency numbers X through Y, 1<=X<=Y<=64)
-bash-3.1$ ./msieve -s c81.log -n -v 6105274315692994651580221674052558368558960
92107776363873073943375109371498913473


Msieve v. 1.46
Sun Mar  5 17:03:22 2023
random seeds: 74adf541 c609a0a3
factoring 610527431569299465158022167405255836855896092107776363873073943375109371498913473 (81 digits)
no P-1/P+1/ECM available, skipping
commencing number field sieve (81-digit input)
R0: -56478943616282822222
R1:  6897069839
A0: -896756383370093653912725
A1: -9268795085226922257
A2: -3262201794610
A3:  8679282
A4:  60
skew 495302.19, size 3.987e-11, alpha -4.588, combined = 1.572e-07 rroots = 2
factor base loaded:
137317 rational ideals (max prime = 1832029)
137471 algebraic ideals (max prime = 1832029)
a range: [-4256410, 4256410]
b range: [1, 4294967295]
number of hash buckets: 28
sieve block size: 65536

maximum RFB prime: 1832029
RFB entries: 137317
medium RFB entries: 6542
resieved RFB entries: 6374
small RFB prime powers: 26
projective RFB roots: 4
RFB trial factoring cutoff: 54 or 81 bits
single large prime RFB range: 21 - 26 bits
double large prime RFB range: 42 - 50 bits
triple large prime RFB range: 65 - 76 bits

maximum AFB prime: 1832029
AFB entries: 137471
medium AFB entries: 6592
resieved AFB entries: 6384
small AFB prime powers: 67
projective AFB roots: 3
AFB trial factoring cutoff: 54 or 81 bits
single large prime AFB range: 21 - 26 bits
double large prime AFB range: 42 - 50 bits
triple large prime AFB range: 65 - 76 bits

multiplying 940555 primes from 1832029 to 16777216
multiply complete, product has 21556204 bits

sieving in progress (press Ctrl-C to pause)
b = 1079034, 6308031 complete / 266651 batched relations (need 6308029)
completed b = 1079034, found 6313335 relations

commencing relation filtering
estimated available RAM is 16268.5 MB
commencing duplicate removal, pass 1
found 279649 hash collisions in 6313335 relations
added 163500 free relations
commencing duplicate removal, pass 2
found 0 duplicates and 6476835 unique relations
memory use: 16.3 MB
reading ideals above 100000
commencing singleton removal, initial pass
memory use: 172.2 MB
reading all ideals from disk
memory use: 196.9 MB
keeping 5902059 ideals with weight <= 200, target excess is 35014
commencing in-memory singleton removal
begin with 6476835 relations and 5902059 unique ideals
reduce to 3886175 relations and 2961446 ideals in 11 passes
max relations containing the same ideal: 146
removing 1220006 relations and 820006 ideals in 400000 cliques
commencing in-memory singleton removal
begin with 2666169 relations and 2961446 unique ideals
reduce to 2457746 relations and 1906185 ideals in 7 passes
max relations containing the same ideal: 104
removing 972437 relations and 572437 ideals in 400000 cliques
commencing in-memory singleton removal
begin with 1485309 relations and 1906185 unique ideals
reduce to 1289874 relations and 1108715 ideals in 8 passes
max relations containing the same ideal: 67
removing 405562 relations and 265020 ideals in 140542 cliques
commencing in-memory singleton removal
begin with 884312 relations and 1108715 unique ideals
reduce to 780532 relations and 726808 ideals in 9 passes
max relations containing the same ideal: 46
removing 71528 relations and 58421 ideals in 13107 cliques
commencing in-memory singleton removal
begin with 709004 relations and 726808 unique ideals
reduce to 704126 relations and 663398 ideals in 7 passes
max relations containing the same ideal: 46
relations with 0 large ideals: 1296
relations with 1 large ideals: 9667
relations with 2 large ideals: 43756
relations with 3 large ideals: 111476
relations with 4 large ideals: 172604
relations with 5 large ideals: 181170
relations with 6 large ideals: 115480
relations with 7+ large ideals: 68677
commencing 2-way merge
reduce to 459639 relation sets and 418911 unique ideals
commencing full merge
memory use: 42.5 MB
found 211950 cycles, need 207111
weight of 207111 cycles is about 14637453 (70.67/cycle)
distribution of cycle lengths:
1 relations: 15092
2 relations: 16667
3 relations: 19196
4 relations: 19658
5 relations: 19521
6 relations: 18009
7 relations: 16520
8 relations: 14531
9 relations: 12708
10+ relations: 55209
heaviest cycle: 24 relations
commencing cycle optimization
start with 1459611 relations
pruned 44114 relations
memory use: 43.4 MB
distribution of cycle lengths:
1 relations: 15092
2 relations: 16926
3 relations: 19856
4 relations: 20318
5 relations: 20143
6 relations: 18627
7 relations: 17004
8 relations: 14855
9 relations: 13019
10+ relations: 51271
heaviest cycle: 23 relations
RelProcTime: 137

commencing linear algebra
read 207111 cycles
cycles contain 679333 unique relations
read 679333 relations
using 20 quadratic characters above 67105454
building initial matrix
memory use: 83.8 MB
read 207111 cycles
matrix is 206933 x 207111 (61.9 MB) with weight 19820352 (95.70/col)
sparse part has weight 13958004 (67.39/col)
filtering completed in 2 passes
matrix is 206720 x 206898 (61.9 MB) with weight 19810849 (95.75/col)
sparse part has weight 13954741 (67.45/col)
read 206898 cycles
matrix is 206720 x 206898 (61.9 MB) with weight 19810849 (95.75/col)
sparse part has weight 13954741 (67.45/col)
saving the first 48 matrix rows for later
matrix is 206672 x 206898 (59.3 MB) with weight 15600165 (75.40/col)
sparse part has weight 13486882 (65.19/col)
matrix includes 64 packed rows
using block size 65536 for processor cache size 16384 kB
commencing Lanczos iteration
memory use: 55.0 MB
linear algebra at 2.9%, ETA 0h 3m206898 dimensions (2.9%, ETA 0h 3m)
linear algebra completed 206430 of 206898 dimensions (99.8%, ETA 0h 0m)
lanczos halted after 3269 iterations (dim = 206672)
recovered 34 nontrivial dependencies
BLanczosTime: 203

commencing square root phase
reading relations for dependency 1
read 103491 cycles
cycles contain 339620 unique relations
read 339620 relations
multiplying 339620 relations
multiply complete, coefficients have about 12.90 million bits
initial square root is modulo 25975123
sqrtTime: 20
prp30 factor: 368596695116186434714529041561
prp52 factor: 1656356227982058657394041821129733144353123888780393
elapsed time 58:24:03


[参考文献]


Last Update: 2023.03.08
H.Nakao

Homeに戻る[Homeに戻る]  一覧に戻る[一覧に戻る]