gcd_ doesn't give the correct result

[tueda]

In the following code, FORM is expected to give gcd(a*g,b*g) = g,

#-
Format 255;
S x1,...,x5;
#$a = 34*x2^2*x5 + x1^2*x2*x4*x5 + x1^5;
#$b = x4^5 + x3^5 + x2*x3*x5^3;
#$g = x3*x4^4 + x2^3*x4 + x1*x3;
#$p = $a * $g;
#$q = $b * $g;
#$gcd = gcd_($p,$q);
#message a = `$a'
#message b = `$b'
#message g = `$g'
#message a*g = `$p'
#message b*g = `$q'
#message gcd(a*g,b*g) = `$gcd'
.end

but gives 1:

FORM 4.2.0 (Dec 15 2017, v4.2.0-29-g26793e4) 64-bits  Run: Mon Jan 22 22:33:45 2018
    #-
~~~a = 34*x2^2*x5+x1^2*x2*x4*x5+x1^5
~~~b = x4^5+x3^5+x2*x3*x5^3
~~~g = x3*x4^4+x2^3*x4+x1*x3
~~~a*g = 34*x2^2*x3*x4^4*x5+34*x2^5*x4*x5+34*x1*x2^2*x3*x5+x1^2*x2*x3*x4^5*x5+x1^2*x2^4*x4^2*x5+x1^3*x2*x3*x4*x5+x1^5*x3*x4^4+x1^5*x2^3*x4+x1^6*x3
~~~b*g = x3*x4^9+x3^6*x4^4+x2*x3^2*x4^4*x5^3+x2^3*x4^6+x2^3*x3^5*x4+x2^4*x3*x4*x5^3+x1*x3*x4^5+x1*x3^6+x1*x2*x3^2*x5^3
~~~gcd(a*g,b*g) = 1

Maybe related to #258 (if the problem is, for example, in the C wrapper). If the problem is some bug in the C++ code, then rational polynomial arithmetic may also be unreliable.

[tueda]

Another (shorter) example:

#-
Format 255;
S x1,...,x4;
#$a = 79*x2 + x2^4 + x1*x3*x4;
#$b = x4^5 + x1*x3^4 + x1^5;
#$g = x2^4*x3 + 84*x1^5;
#$p = $a * $g;
#$q = $b * $g;
#$gcd = gcd_($p,$q);
#message a = `$a'
#message b = `$b'
#message g = `$g'
#message a*g = `$p'
#message b*g = `$q'
#message gcd(a*g,b*g) = `$gcd'
.end

FORM 4.2.0 (Dec 15 2017, v4.2.0-29-g26793e4) 64-bits  Run: Tue Jan 23 13:55:49 2018
    #-
~~~a = 79*x2+x2^4+x1*x3*x4
~~~b = x4^5+x1*x3^4+x1^5
~~~g = x2^4*x3+84*x1^5
~~~a*g = 79*x2^5*x3+x2^8*x3+x1*x2^4*x3^2*x4+6636*x1^5*x2+84*x1^5*x2^4+84*x1^6*x3*x4
~~~b*g = x2^4*x3*x4^5+x1*x2^4*x3^5+84*x1^5*x4^5+x1^5*x2^4*x3+84*x1^6*x3^4+84*x1^10
~~~gcd(a*g,b*g) = 1

[tueda]

This must be a bug in the C++ code, because pa and pb in poly_gcd() get the correct polynomials.

[benruijl]

Yes, it is. I'll go after it. Some preliminary info:

*** [0]  CALL: gcd(84*a^10+84*a^6*c^4+a^5*b^4*c+84*a^5*d^5+a*b^4*c^5+b^4*c*d^5,84*a^6*c*d+84*a^5*b^4+6636*a^5*b+a*b^4*c^2*d+b^8*c+79*b^5*c)
*** [3.8e-05]  CALL: content_univar(84*a^10+84*a^6*c^4+a^5*b^4*c+84*a^5*d^5+a*b^4*c^5+b^4*c*d^5,0)
*** [5.5e-05]  CALL: gcd(0,84)
*** [6.2e-05]  RES : content_univar(84*a^10+84*a^6*c^4+a^5*b^4*c+84*a^5*d^5+a*b^4*c^5+b^4*c*d^5,0) = 1
*** [6.8e-05]  CALL: content_univar(84*a^6*c*d+84*a^5*b^4+6636*a^5*b+a*b^4*c^2*d+b^8*c+79*b^5*c,0)
*** [7.6e-05]  CALL: gcd(0,84*c*d)
*** [8e-05]  CALL: gcd(84*c*d,84*b^4+6636*b)
*** [9.1e-05]  RES : content_univar(84*a^6*c*d+84*a^5*b^4+6636*a^5*b+a*b^4*c^2*d+b^8*c+79*b^5*c,0) = 1
*** [9.7e-05]  CALL: gcd(1,1)
*** [0.000112]  CALL: gcd_heuristic(84*a^10+84*a^6*c^4+a^5*b^4*c+84*a^5*d^5+a*b^4*c^5+b^4*c*d^5,84*a^6*c*d+84*a^5*b^4+6636*a^5*b+a*b^4*c^2*d+b^8*c+79*b^5*c,{0,1,2,3})
*** [0.000155]  CALL: gcd_heuristic(13282*b^4*c^5+b^4*c*d^5+413349464675634214432*b^4*c+461169037545028985431209216*c^4+34721355032753274012288*d^5+14352053515601203688793202131000243594940416,b^8*c+79*b^5*c+13282*b^4*c^2*d+34721355032753274012288*b^4+2742987047587508646970752*b+461169037545028985431209216*c*d,{1,2,3})
*** [0.000209]  CALL: gcd_heuristic(9016534295484106205282010622254171751216731761128231857835499535295364381879750009711975296555102669890809050807102756928437138119680090269470956286679868427338538102686231997472*c^5+461169037545028985431209216*c^4+678853658747485785670984085397844582985750019660309581225380178835669656819737239098929023983970988547719398494737445936488265179918693741113609116599899746072770524219713296*c*d^5+280603796436368919296866715908630939686610017276182275057097179519080829795683673408444356009079019359684677597518994445735900025802575231090330878860970538851202556969050099372264166511225487872*c+34721355032753274012288*d^5+14352053515601203688793202131000243594940416,9016534295484106205282010622254171751216731761128231857835499535295364381879750009711975296555102669890809050807102756928437138119680090269470956286679868427338538102686231997472*c^2*d+461169037545028985431209216*c*d+460842289994847884370238061444390717993171359139330139905779197189241189002896573652510887031413323633387927920738461101534720430726380382877816302204067709050889371804353543210615322749451929049489623594672446431704253408091261917989530208090476756815267625785763671109898614694564916315758229592973268099888164759348815675065725563308886242126944*c+23570718900654989220936804136324998933675241451199311104796163079602789702837428566309325904762637626213512918191595533441815602246151313210741533180663645696573943937516080017136499312254080114432,{2,3})
*** [0.000281]  CALL: gcd_heuristic(625688949337144697331022987222542603160195495764857371397918695966983725666819212680084711748962632519288994747732234134257811599167706514248059501757733086702725864357688849776356426846389467623770050681601208976761896622654770858204668144200560103715061552886329532177102793370110469280693009743713307662242867605784123247742134738876432821416594366869376290167887138185659653293823085024482456477843760773175848924652070554898717564141938405663622631034562451649246790149309524215121747982080821829687095146740425586912*d^5+5997237492869676437001288044200227812981058466078151785180762125016971684028185917832555261640637030977181644721653683389853486337277550831298415983434815204395112928749195725765275982656553823023790260802729205805142553063868509538257166039739587174685844562724376974697528429639331754110064526033266752133284339239548675722976849166103505835136741665333865058732702994820016469475137507934077533492694840623396204523807908663745777761893508160322624097073668189894168262732534602099799498732192958189741026788824528005777824620263092211817544095116108588345268438127693354744624458822173245695917841673644124899288083605195206598950471829718654812798547597908580929379233801482705951282466187831412835597974890487847665729192920398710611346054429752084738880369571694533409878904446203548102538553517777717413137123760961506420571352540994665610287836791904850234919631070015029566414810613247496498319353672132362061826647227298375920592009253673674008080679145276443785225319995329737822323151300420539860802571717696631142672608144878458165742549739842288424010692664634532576142622871964207586873157027381084596468841934199357147169522697459128827143949334200818855644632174894561232063241852941536163575497218647709444909396185444984777173439344132210385554926922764566755308237819364898189937360221247264752102124720335584359528720649860316472622939672221411343017471433839252418338815269791071833382132971236431021379533141592661364092832925839490225623135930733534545794101560589121370678477717879393911098761578430713494382389145294290144131417868110827490443103921353873321416506768413524012382660665466708569133027164640886258772398847756005730857991177335887323859324541106986726513774620681756699696956570867516330377270858787966085230459321761676115824576092055385636294640083021457553521107127537682311468067653744521305519317336970887157247820428447060569292343023454258491970090807646977983604226048,7659568109687671255793716475420527259491775667437809734166869049062373508021853156948742757233608135550671599836538240852110124860462267863787823115045295152190408576767049429243844379671397656667694570045639627507543870129836565041663339393005964201426891158504111930903752724351966863112026234315729621137766970714366127067872190965070921240282105391372352998929495245554712829322759124629851723562842532298146529041267611677661275307075655044485162244792048066350147162805552727479923670269309280197884331510998161266228110139105362789211050230647014894032292961204376134581507827919195723890002528515958470745703880353770295777485111166811343039938529195541674059706619398060834120360640842989151452243634463401964804362552241048888019307670718704940242269774564592713912949424966738568454339197186410264695712154015297074177524921862982932274711530001998627459404026496*d+424751232495390948939395952915699779052315104203229772179869243166736266522689359274835666463980267470006289323649751947386034418850966691934968081518707001718872478230409114159433328868147606421501399863521878528174155675637569588327190009340156041166308111648083545161556757757892501868425648232289565585326935092483800396954301524135132033828922929843041062193043015071394164060935582292489820995189724153127083982612396132846979539869598457989071366971334410562650651335118468930908239073907306358638935798746400115463153698575324743966391499541897391958856197076944143731454395780293456846395554978083294261819341678585946442177361622883681478628457737755918485456742680118377314021554434368,{3})
*** [0.000558]  RES : gcd_heuristic(625688949337144697331022987222542603160195495764857371397918695966983725666819212680084711748962632519288994747732234134257811599167706514248059501757733086702725864357688849776356426846389467623770050681601208976761896622654770858204668144200560103715061552886329532177102793370110469280693009743713307662242867605784123247742134738876432821416594366869376290167887138185659653293823085024482456477843760773175848924652070554898717564141938405663622631034562451649246790149309524215121747982080821829687095146740425586912*d^5+5997237492869676437001288044200227812981058466078151785180762125016971684028185917832555261640637030977181644721653683389853486337277550831298415983434815204395112928749195725765275982656553823023790260802729205805142553063868509538257166039739587174685844562724376974697528429639331754110064526033266752133284339239548675722976849166103505835136741665333865058732702994820016469475137507934077533492694840623396204523807908663745777761893508160322624097073668189894168262732534602099799498732192958189741026788824528005777824620263092211817544095116108588345268438127693354744624458822173245695917841673644124899288083605195206598950471829718654812798547597908580929379233801482705951282466187831412835597974890487847665729192920398710611346054429752084738880369571694533409878904446203548102538553517777717413137123760961506420571352540994665610287836791904850234919631070015029566414810613247496498319353672132362061826647227298375920592009253673674008080679145276443785225319995329737822323151300420539860802571717696631142672608144878458165742549739842288424010692664634532576142622871964207586873157027381084596468841934199357147169522697459128827143949334200818855644632174894561232063241852941536163575497218647709444909396185444984777173439344132210385554926922764566755308237819364898189937360221247264752102124720335584359528720649860316472622939672221411343017471433839252418338815269791071833382132971236431021379533141592661364092832925839490225623135930733534545794101560589121370678477717879393911098761578430713494382389145294290144131417868110827490443103921353873321416506768413524012382660665466708569133027164640886258772398847756005730857991177335887323859324541106986726513774620681756699696956570867516330377270858787966085230459321761676115824576092055385636294640083021457553521107127537682311468067653744521305519317336970887157247820428447060569292343023454258491970090807646977983604226048,7659568109687671255793716475420527259491775667437809734166869049062373508021853156948742757233608135550671599836538240852110124860462267863787823115045295152190408576767049429243844379671397656667694570045639627507543870129836565041663339393005964201426891158504111930903752724351966863112026234315729621137766970714366127067872190965070921240282105391372352998929495245554712829322759124629851723562842532298146529041267611677661275307075655044485162244792048066350147162805552727479923670269309280197884331510998161266228110139105362789211050230647014894032292961204376134581507827919195723890002528515958470745703880353770295777485111166811343039938529195541674059706619398060834120360640842989151452243634463401964804362552241048888019307670718704940242269774564592713912949424966738568454339197186410264695712154015297074177524921862982932274711530001998627459404026496*d+424751232495390948939395952915699779052315104203229772179869243166736266522689359274835666463980267470006289323649751947386034418850966691934968081518707001718872478230409114159433328868147606421501399863521878528174155675637569588327190009340156041166308111648083545161556757757892501868425648232289565585326935092483800396954301524135132033828922929843041062193043015071394164060935582292489820995189724153127083982612396132846979539869598457989071366971334410562650651335118468930908239073907306358638935798746400115463153698575324743966391499541897391958856197076944143731454395780293456846395554978083294261819341678585946442177361622883681478628457737755918485456742680118377314021554434368,{3}) = overflow
*** [0.000875]  CALL: gcd_modular(84*a^5*b^4+6636*a^5*b,b^8+79*b^5,{0,1})
*** [0.001874]  CALL: gcd_modular_dense_interpolation(a^5*b^4+79*a^5*b (mod 2147483629),b^8+79*b^5 (mod 2147483629),{0,1},,0 (mod 2147483629))
*** [0.001906]  CALL: gcd_Euclidean(0 (mod 2147483629),b^4+79*b (mod 2147483629))
*** [0.001912]  CALL: gcd_Euclidean(0 (mod 2147483629),b^8+79*b^5 (mod 2147483629))
*** [0.001929]  CALL: gcd_Euclidean(b^4+79*b (mod 2147483629),b^8+79*b^5 (mod 2147483629))
*** [0.001941]  RES : gcd_Euclidean(b^4+79*b (mod 2147483629),b^8+79*b^5 (mod 2147483629)) = b^4+79*b (mod 2147483629)
*** [0.001948]  CALL: gcd_Euclidean(1 (mod 2147483629),1 (mod 2147483629))
*** [0.00208]  CALL: gcd_modular_dense_interpolation(a^5 (mod 2147483629),1 (mod 2147483629),{0},,0)
*** [0.002089]  CALL: gcd_Euclidean(a^5 (mod 2147483629),1 (mod 2147483629))
*** [0.002099]  CALL: gcd_modular_dense_interpolation(a^5 (mod 2147483629),1 (mod 2147483629),{0},,1 (mod 2147483629))
*** [0.002105]  CALL: gcd_Euclidean(a^5 (mod 2147483629),1 (mod 2147483629))
*** [0.002113]  CALL: gcd_Euclidean(0 (mod 2147483629),1 (mod 2147483629))
*** [0.002118]  RES : gcd_modular_dense_interpolation(a^5*b^4+79*a^5*b (mod 2147483629),b^8+79*b^5 (mod 2147483629),{0,1},,0 (mod 2147483629)) = b^4+79*b (mod 2147483629)
*** [0.00213]  CALL: content_univar(b^4+79*b,0)
*** [0.002134]  CALL: gcd(0,b^4+79*b)
*** [0.002138]  RES : content_univar(b^4+79*b,0) = b^4+79*b
*** [0.002142]  RES : gcd_modular(84*a^5*b^4+6636*a^5*b,b^8+79*b^5,{0,1}) = 1
*** [0.002158]  RES : gcd(84*a^10+84*a^6*c^4+a^5*b^4*c+84*a^5*d^5+a*b^4*c^5+b^4*c*d^5,84*a^6*c*d+84*a^5*b^4+6636*a^5*b+a*b^4*c^2*d+b^8*c+79*b^5*c) = 1

It's been a while since I looked at this code, but the numbers in the heuristics function don't seem quite right. However, the real wrong result is caused by the gcd_modular.

[benruijl]

The correct result is obtained if the univariate content is not removed from the result at the end of gcd_modular. In this specific case the univariate content is computed on a constant term: content_univar(b^4+79*b,a). Perhaps this is an exception?

Edit: we should likely only divide by the integer content, as described in Algorithms for Computer Algebra by Geddes. In the paper Computing the Greatest Common Divisor of
Multivariate Polynomials over Finite Fields by Suling Yang there is a division by the univariate content, but that is compensated by taking the univariate content for other constants too. I don't know which solution is better.

[tueda]

@benruijl Thanks. Your commit will fail to pass the tests on Travis CI, but this is due to some upgrade (against Meltdown) on the Travis side. I will fix it.

[tueda]

Should we suppress "New GCD attempt (unused vars)..." message?

[benruijl]

Yes, that should be removed. Sorry! Op di 23 jan. 2018 18:42 schreef Takahiro Ueda <[email protected]>:

…

Should we suppress "New GCD attempt (unused vars)..." message? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#260 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AARGGfVTBX9UeAZjfreJOSkrXEFi3DRdks5tNhnwgaJpZM4RorU3> .

[tueda]

Do you want to do that? Or I can do with adding test cases.

[benruijl]

You can do it :) Op di 23 jan. 2018 om 18:57 schreef Takahiro Ueda <[email protected]

…

: Do you want to do that? Or I can do with adding test cases. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#260 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AARGGcrFniACflZMIMHZzWq09EOvq0Iwks5tNh2SgaJpZM4RorU3> .

[tueda]

This fix surely affects on PolyRatFun with many symbols:

S x1,...,x4;
CF rat;
PolyRatFun rat;
#$a = 79*x2 + x2^4 + x1*x3*x4;
#$b = x4^5 + x1*x3^4 + x1^5;
#$g = x2^4*x3 + 84*x1^5;
#$p = $a * $g;
#$q = $b * $g;
L F = rat(1,$p) + rat(1,$q);
P;
.end

In the above case, the reduction of the fraction was not perfect with the previous version.

[tueda]

Thanks to this fix, now I can make some benchmark plot for GCD of random polynomials, inspired by plots of Rings. In case you are interested, the slowest gcd (the right-most point; FORM: 140.5s, Mathematica: 3.3s) is here.

[spj101]

Hey @tueda, this benchmarking script is really interesting. I added support for Fermat which is frequently used for its fast GCD in integral reduction codes (e.g. Reduze and Fire). You can find the extended script here.

Surprisingly the FORM GCD performs much better than Fermat on relatively simple polynomials (the default settings of your script). Perhaps this indicates that I implemented the link to Fermat poorly?

However, as I increase the complexity (especially the number of terms in the polynomial) Fermat becomes much more competitive and soon greatly outperforms FORM for the types of examples we typically encounter during our reduction. I added to the README some example parameters that more closely mirror what we encounter in our reduction.

[tueda]

Actually the updated script now has options to run Fermat and Rings, via libFermat and Ammonite REPL, respectively :-) I got results that FORM wins against them by the default benchmark parameters, but I didn't upload the plots because these may be totally dependent on the specific examples (and I'm not sure about the license for Fermat at Nikhef).

libFermat is written in C++ and just uses stdin and stdout (redirected by pipes), so it is almost what we can do best, I think, and still I got a bit disappointed result. In this sense, the slowness of Fermat for the default benchmark parameters via pipes in your case is not due to your Python implementation (to some extent; certainly running Fermat for each problem has some overhead. I run Fermat once for the whole problem set).

On the other hand, it is claimed that Rings is faster than Mathematica (for their examples), but this is not the case in my default examples. I'm not sure the REPL would have heavy overhead, or just because of difference in examples.

Again, these results may totally depend on the benchmark parameters. Even for the default parameters, FORM has some unlucky cases where FORM is somehow very slow (this would be a clue to improve the FORM's GCD performance for more complicated problems). I haven't tried for more complicated inputs with more terms so much, where Fermat (and the other programs) would outperform FORM.

[tueda]

I just tried ./run.py --nvars 3 --degree 50 --nterms 5000 --coeffpow 100 --nwarmups 0 --nproblems 10 --form --mathematica --rings --fermat on a machine with AMD Opteron 8439 SE and enough memory. The printed timings were:

FORM: 18845.189s
Mathematica: 872.407s
Rings: 126594.310s
Fermat: 3142.147s

[spj101]

Hi @tueda, another code which is used for GCD by Reduze is GiNaC, though its performance is usually very much worse than Fermat (see Fig 3 of the Reduze 2 paper). I wish it was legal/possible for you to make your full testing script public so that there was a single public tool for GCD testing of all interesting packages.

One thing that I find odd with your above numbers is that Mathematica is very much faster than Fermat. I wonder: why does FIRE bother to implement a link to Fermat from a mostly Mathematica code if it has a fast GCD available within Mathematica?

[tueda]

The script also has --ginac option which makes a C++ program using GiNaC. I think I tried once for the default setting but it was too slow and I gave it up.

I'm not sure but when FIRE was firstly programmed maybe the situation was different (the paper for the public FIRE 3.0 appeared in 2008). Or, the benchmark just tests polynomial GCD; performance for rational function arithmetic would be different, possibly.

[PoslavskySV]

Dear @tueda,

I just looked your benchmark script -- the part where you call Rings library. You got a dramatically low performance for Rings because of one reason: you run a separate Java instance for each GCD example (via RingsSolver method) -- this is not how Rings (and Java in general) should be used:) Java's just-in-time compiler will do not have time to make optimizations and actually you obtain the result produced by interpreted (not compiled) code, which is 10-100 (or even 1000) times slower than the compiled one.

So, to use Rings correctly in your benchmark, you need to have a single instance of rings.repl process, and communicate to it without creating and destroying the system process for each GCD example. That would give you a completely different timings.

I'm not sure, but I suppose you can do something like this in python:

# this should be a single instance for all GCD computations
ringsProcess = subprocess.Popen(['rings.repl'], stdin=PIPE, stdout=PIPE, bufsize=1)

# do some GCD example
def solveWithRings(some_data):
    # write input
    print >>ringsProcess.stdin, "rings_code_for_gcd"
    #  flush process input
    ringsProcess.stdin.flush()
    result = ""
    for i in range(how_much_lines_you_expect_in_output):
        # read result
        result = result + ringsProcess.stdout.readline()
    return result

# first invocation is quite slow (code is just interpreted)
result_1 = solveWithRings(some_example_1)
# this will be much faster (JIT will compile "hot" methods)
result_2 = solveWithRings(some_example_2)
# this will be even more faster (JIT will compile more "hot" methods)
result_3 = solveWithRings(some_example_3)
# etc...

I'll try to incorporate this in your python script later this week, but if you do this faster, I'll really appreciate it since it is very interesting)

[tueda]

@PoslavskySV, thanks for your comment. Actually I worried that maybe I was doing something wrong. But I use only one rings.repl process. In RingsSolver._solve() the Python code generates a Scala script that processes all the gcd problems, iterated by

for (i <- 1 to N) {
    ...
}

Then this script is put into the standard input of rings.repl as

os.system('rings.repl <rings.tmp.prog >/dev/null 2>&1')

though I'm not sure if JIT is guaranteed to work in REPL.

[PoslavskySV]

@tueda That does make sense! Sorry, I didn't figure that promptly:) I am worried about the numbers

Mathematica: 872.407s
Rings: 126594.310s

because that is completely (even extremely!) different from what I had earlier. I'll look closely what is the reason and let you know. Thanks for sharing your script!

[spj101]

@tueda On my machine (Intel Core i7-7700 CPU @ 4.0 GHz, openSUSE Tumbleweed 20180219, 64 GB RAM) the time that your script reports for Mathematica 11 is not very accurate.
Running:

time python run.py --nvars 3 --degree 50 --nterms 5000 --coeffpow 100 --nwarmups 0 --nproblems 10

I get:

// No Solver (estimate setup time)
real	0m35.462s
user	0m32.601s
sys	0m3.284s

// Mathematica 11.2.0
Mathematica: 340.467s
real	37m26.707s
user	22m46.853s
sys	22m36.008s

Measuring this with a wall clock suggests that the "real" time reported is accurate. Can you reproduce this?

[spj101]

On the same machine for Mathematica 8.0 - 10.4.1 the timing is also inaccurate.
Running:

time python run.py --nvars 3 --degree 50 --nterms 5000 --coeffpow 100 --nwarmups 0 --nproblems 10

I get:

// Mathematica 10.4.1
Mathematica: 335.254s
real	37m8.622s
user	21m52.255s
sys	23m16.229s

// Mathematica 9.0
Mathematica: 300.038s
real	36m56.102s
user	21m13.997s
sys	23m41.568s

// Mathematica 8.0
Mathematica: 230.160s
real	35m37.581s
user	19m46.900s
sys	23m46.013s

However, older versions of Mathematica are very much faster:

// Mathematica 7.0
Mathematica: 254.641s
real	5m32.534s
user	5m28.988s
sys	0m4.024s

// Mathematica 6.0
Mathematica: 242.787s
real	5m18.879s
user	5m16.200s
sys	0m4.075s

For comparison:

// Fermat 6.19 (January 5, 2018)
Fermat: 645.884s
real	11m25.884s
user	11m23.150s
sys	0m4.889s

The results of Mathematica 6.0 and Fermat 6.19 agree.

[tueda]

I got the following numbers on my Desktop PC:

$ time python run.py --nvars 3 --degree 50 --nterms 5000 --coeffpow 100 --nwarmups 0 --nproblems 10 --mathematica
Mathematica: 462.271s

real    9m18.317s
user    8m59.136s
sys     0m8.729s

The setup is less than 1 minute, so, yes there is some gap. This would be attributed to I/O: the number of Mathematica: 462.271s comes from timings only for gcd. Reading from the input file and writing to the output file etc. are not considered. The main Python program also needs to collect results from the output file.

Actually the FORM benchmark has a disadvantage for this point: it utilizes _sover_for() and the measured timing includes also input/output via pipe (well, but via pipes, not via files and should be rather fast).

Your observation that older Mathematica is faster than newer one is very interesting!

[tueda]

Wait, you got 37m wall clock time for 340.467s of gcds...?? hmm... Apparently I didn't look in your real time numbers carefully, and they have really big gaps.

[spj101]

Yes, around 30 of 37 minutes was lost with Mathematica v8-11. I find it surprising too. If you can't reproduce this then probably it is a quirk of my machine and it is best to ignore it for this discussion. I will look into it.

[PoslavskySV]

@tueda Takahiro, I also noticed that random polynomials generated with random_poly have special form: they are monic in some variables. This condition meets not very often in practice (but it is also not quite rare), and, what is important, is that Zippel algorithm for GCD used by Mathematica and Rings (and as far as I understand by FORM) is very sensitive to this special monic case.

So, to make random polys generated closer to those occur in practice, I recommend to change method as follows:

def random_poly(nvars, degree, nterms, coeffpow):

    ...

    while True:
        polynomial = ''
        for i in range(nterms):
            coeff = random_coeff(m)
            monomial = ''
            sum_degree = degree + random.randint(0, degree)
            
            # if i == 0:
            #     if flag1 and nvars <= sum_degree:
            #         for j in range(nvars):
            #             monomial += '*x{0}'.format(j + 1)
            #         sum_degree -= nvars

            xx = list(range(nvars))
            random.shuffle(xx)
            for j in xx:
                # these 2 lines caused polys to be monic in some variables
                # p = random.randint(0, sum_degree)
                # sum_degree -= p

                p = random.randint(0, sum_degree / nvars)
                if p != 0:
                    monomial += '*x{0}^{1}'.format(j + 1, p)
            polynomial += '+({0}){1}'.format(coeff, monomial)
        
        if random.randint(0, 3) == 0:
            polynomial += "+1"


        form.write('#$a={0};'.format(polynomial))
        polynomial = form.read('$a')
        if polynomial != '0':
            return polynomial

With this modification the timings distribution become different. Also, I noticed that FORM hangs in some cases (I digged into FORM's code, and probably the reason is that it misses LinZip algorithm for non-monic case GCD ? but I'm not sure).

I tried the following benchmark:

./run.py --nvars 5 --degree 40 --nterms 100 --coeffpow 1 --nwarmups 30 --nproblems 10 --mathematica --rings --form

and got the following result:

Mathematica: 9.425s
Rings: 25.992s
FORM: 500.759s

BTW: I've also found the reason of performance dip in Rings; I already fixed it locally, will release in version 2.3.1 on the next week.

[tueda]

Thanks. I was trying to make some uneven distribution of powers of variables with some probability, but I didn't know that the GCD algorithm is so sensitive for such cases. Maybe I would put an command line option to select random polynomial generator based on your suggestion.

[benruijl]

@PoslavskySV Could you file a bug report for cases where the gcd computation hangs? The non-monic linzip should be in there, but there is likely an error in the algorithm (I suspect an improper extraction of the content or something like that).

[PoslavskySV]

@benruijl I just re-run the benchmark, waited more time (I used a large --nwarmups so the the result didn't appear for nearly an hour), and found that FORM actually not hangs, but just takes too long -- I've updated my comment. I've now re-run the benchmark again -- I'll pick the slowest FORM example and make it a separate issue.

[PoslavskySV]

@benruijl I've added several debugging flags to the script, and found that FORM is constantly slow with random_poly function I pasted before and nvariables = 5 . This is the table of timings for generated problems (see attached problems.txt, two lines span one problem as generated by script run.py):

Problem	FORM	Rings (Mathematica is +- similar)
problem 0	10.5s	538ms
problem 1	34.2s	723ms
problem 2	44.4s	493ms
problem 3	45.5s	724ms
problem 4 (LinZip)	43.8s	1440ms
problem 5	59.7s	954ms
problem 6 (LinZip)	47.3s	1009ms
problem 7	39.8s	566ms
problem 8	37.0s	971ms
problem 9	36.3s	736ms
problem 10 (LinZip)	82.9s	1040ms
problem 11	56.1s	933ms
problem 12	36.9s	635ms
problem 13	108.0s	814ms
problem 14 (LinZip)	33.2s	526ms
problem 15	23.4s	707ms
problem 16	50.8s	1050ms
problem 17	120.0s	836ms
...	...	...

(Rings timings are not very accurate -- they are from my develop, which contains some fixes, but they are still relevant , Mathematica timings are close to Rings on these examples, so I didn't paste them here).

All these examples are non-monic case GCD, however in many of them, still the monic-case algorithm with some "scaling corrections" may be applied (as described in [dKMW05]); problems where the use of LinZip is unavoidable (i.e. where Rings switched to LinZip) are marked in bold.

I'm not sure whether the reason of slowness is just the non-monic case, since GCD over Z requires some more work to do than plain Zippel/heuristic algorithm. To detect this certainly, it would be useful to check GCD over some (sufficiently large) finite field, say Z_{2^31 - 1}. I'm not quite familiar with FORM -- is that possible to do and could you extend the script with this functionality?

[dKMW05] J de Kleine, M B Monagan, A D Wittkopf. Algorithms for the Non-monic Case of the Sparse Modular GCD Algorithm. Proceeding of ISSAC ‘05, ACM Press, pp. 124-131 , 2005.