Playing with CADO's individual polyselect binaries

[Dubslow]

Originally Posted by jasonp

1) The value of the integrals would be the same, but remember that we are optimizing very large polynomial expressions, and rectangular vs radial coordinates yield different polynomials that have different gradients and hessians. Our experience is that radial coordinates are better behaved numerically, presumably because polynomial size in the sieving region naturally has radial patterns in it

2) I believe you are correct, the size optimization is considered part of the initial hashtable-based searching.

3) My memory is hazy on what the difference was; the formula for the E-value does include the sieving bound and sieve area size as parameters, but those are pretty much always chosen to be the same numbers across the different implementations I know about. Actually I remember now they gave different answers because Dickman's function is computed in segments, and the RSA-896 poly search required a few more segments than that version of CADO had tabulated. So we should always have produced the same answers, but in that specific case there was a small bug that they fortunately found quickly

1) Again, I think "the" polynomials and their gradients/hessians, the mathematical objects, are "the same" amongst coordinate changes, but I understand you to mean that the numerical computation of the mathematical objects use different numbers (vacuously), and the different numbers tend to behave better from a floating point implementation standpoint, presumably due to some symmetry of the mathematical objects in question?

2) I've just discovered that the binaries produced by building CADO include polyselect, which does "stage 1" and size opt at the same time, as well as a separate sopt binary which takes "raw" polynomials (presumably the results of the hashtable search...?) for size opting as well, and the equivalent output of either binary can then be fed to polyselect_ropt whose function should be obvious in this context.

3) I've just started a run of the latter for my Aliquot 4788 efforts, and it seems that it prints the input values to the MurphyE score as well as the score itself, so I guess we can compare against whatever Msieve uses (though I don't really feel like digging through its source, maybe you can comment? Here's an example:)

### root-optimized polynomial 13 ###
# n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
# Y0: -36095620248719813997359851316167411851
# Y1: 287364767617774550963053
# c0: -56088615892505493935209426379951715142392160
# c1: -144890118607664376088446114937854122208
# c2: -8348445372213634471790816940388
# c3: 28963477043189504738748320
# c4: -11102776326974469
# c5: 1725480
# skew: 3267584.000
# # lognorm 63.38, alpha -5.35 (proj -1.53), E 58.03, 3 real roots
# # MurphyE(Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16)=3.27e-15
### Best MurphyE so far is 4.50e-15

Incidentally, I'm really a fan of taking the log of all these numbers with large exponents. Wouldn't it be nicer for everyone if we also used logMe (and lognorm obviously) in place of the raw MurphyE and norm?

[Dubslow]

Originally Posted by VBCurtis

Attached is my top 500 hits from 15 GPU-days. Worst norm is 1.46e26, best 1.33e25. Please post best poly or two from CADO root-opt on this file.

Err... I may be misinterpreting the .ms file (been years since I worked with one), but... were you searching for degree 6 polys...?

[VBCurtis]

Nope, definitely degree 5.

[Dubslow]

Originally Posted by VBCurtis

Nope, definitely degree 5.

Could you help me spot my error then?

Here's the code (msieve/gnfs/poly/poly_skew.c):

static void sizeopt_callback_log(uint32 deg, mpz_t *alg_coeffs, mpz_t *rat_coeffs, 
				double sizeopt_norm, double projective_alpha, 
				void *extra)
{
	uint32 i;
	FILE *mfile = (FILE *)extra;

	for (i = deg; (int32)i >= 0; i--)
		gmp_fprintf(mfile, "%Zd ", alg_coeffs[i]);

	gmp_fprintf(mfile, "%Zd %Zd %.2lf %le\n", rat_coeffs[1], 
				rat_coeffs[0], projective_alpha, 
				exp(projective_alpha) * sizeopt_norm);
	fflush(mfile);
}

And here's the first three lines of your file I downloaded:

971700 273849234706546 225962405583700796407442 -84054861767065123330509949866 -29517104017082095088929447630073626321 3262277732054109330095369390256411899572767 1660633240515241123 -40487932918725137061265492892721380626 -1.59 1.338783e+25
628488 25807572840266955 303999146710932834523652 -3119549682722644921895399900207 -14123903064615248232023501030205940006 48260972168255749244473149176534441863659851 29048885255261071969 -44174601963086515462086860672795847500 -1.77 1.614272e+25
546960 -5605325602623737 -1955381755454195745958 6621960285938958484871001763602 -1854857638980877962690869298795624892 -1032866398309654043779803193974135126928581602 4835900641980806243 -45419351810482876976420189466026262897 -1.37 1.842692e+25

I count 8 coefficients per line (plus the two scores), so it definitely looks like degree 6...

[VBCurtis]

No idea, sir- but the polys in msieve.log are degree 5, and I am certain I didn't invoke the degree=6 flag in the command line. Msieve doesn't pick degree 6 on its own...

Try running -npr on the file from msieve to convince yourself?

[jcrombie]

Originally Posted by Dubslow

Could you help me spot my error then?

Sure, deg 5 polys have 6 coefficients. (I'm sure you'll be kicking yourself for that )

[Dubslow]

Originally Posted by jcrombie

Sure, deg 5 polys have 6 coefficients. (I'm sure you'll be kicking yourself for that )

Aww crud. I usually have a feeling that I'm doing something stupid, hence all the ellipses, and... well I wasn't wrong, about the feeling stupid part at least... ugh.

Thanks Obviously I'm rusty.

[Dubslow]

General note -- before realizing that CADO can and will happily read Msieve format size opted polys, I wrote a short Python script to convert formats. Unfortunately, it took me a few days (only 20 minutes of actual effort) to discover that having more than one blank line separating polys in the file is beyond CADO's ability to parse. Just a heads up, and for future reference.

[Dubslow]

Okay, here's the results from my cluster of messing around with CADO.

A search of LCs from 0-10M (only those divisible by 60, per CADO default) (roughly 3 CPU days), then passing the top 500 size norm hits to CADO rootopt, produced the following top 5 polys (I have all 500 available):

n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
skew: 62078976.0
c0: 2078861482606939262732348788772234729500066505
c1: -81452413435682187177974315743754413126
c2: -1256739981486858947548697657268
c3: -22396930453004670461906
c4: -263922878237933
c5: 2547840
Y0: -33388556319178604525067212126801994606
Y1: 96874809441465515717599
# MurphyE (Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16) = 1.16e-14
# lognorm 59.06

n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
skew: 57720832.0
c0: -17160072336623445558486491447680793046654059460
c1: -1654348007928800032775401504767487803936
c2: 28947429876755966986844833643017
c3: 946553263211289715531399
c4: -3924318619058088
c5: 5516640
Y0: -47025518976311715704856348352373464079
Y1: 176154525925723359784759
# MurphyE (Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16) = 8.85e-15
# lognorm 61.5

n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
skew: 25583616.0
c0: -685992159184539512613548044513984252041883140
c1: 134610898781364897801610290286626955999
c2: 7109903060621761009580559757488
c3: -238819588636731179920283
c4: -8746008810620432
c5: 53587200
Y0: -27519388373819676833051410839329305921
Y1: 1839671560743535632496133
# MurphyE (Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16) = 8.51e-15
# lognorm 60.22

n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
skew: 14503936.0
c0: -65698494574585654869665676718988331163562112
c1: 97332536632788185915561208310178169236
c2: -371931585676372953901140963249
c3: -870401356387883046717037
c4: -910339736204828
c5: 3553440
Y0: -41220600122861670237999244493991322477
Y1: 118715258484298437008683
# MurphyE (Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16) = 7.93e-15
# lognorm 60.58

n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
skew: 58310656.0
c0: -18280501391805263953245129313762844281481248662
c1: 646040699763406262391803981659489087237
c2: 36712368286900883073060051628284
c3: -529564887634112218615187
c4: -6061192475321092
c5: 10035120
Y0: -33492176102045778534763735786623971315
Y1: 2449174464852706655304347
# MurphyE (Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16) = 7.81e-15
# lognorm 61.1

Note that I had originally searched 0-50M (roughly two weeks CPU time), but lost the results in one of the greater ups of recent times, and no doubt the top 5 would have been substantially better using those results, if not a better top 1. (Edit: 0-50M produced 1 hit below lognorm 59, and 33 total below 60; 0-10M produced the same hit below 59, and only 6 others below 60.)

And secondly, I also ran an identical rootsieve of VBCurtis' top 500 Msieve GPU stage 1 and sizeopt hits (~2 weeks GPU+CPU time?), which produced the following top 5 (again, all 500 available upon request):

n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
skew: 146997248.0
c0: -131414304044478185366394231747951586868400235157
c1: 4029601527701148189310418069976213395221
c2: 30149105678605166244788285142337
c3: -268424081661710226278041
c4: -1338222885809036
c5: 1505028
Y0: -37095662305095086253404478445430777250
Y1: 9762700048268674463
# MurphyE (Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16) = 9.54e-15
# lognorm 61.02

n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
skew: 115310592.0
c0: 72654783078207317663016087360899076774225974549
c1: -13183316839111760341950789781721971964575
c2: -29276900156268889268041437680369
c3: 1780141248353293122105769
c4: 1020057436955814
c5: 912600
Y0: -40999253791772676218846852869952650362
Y1: 4848218606269847041
# MurphyE (Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16) = 8.51e-15
# lognorm 62.36

n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
skew: 81494016.0
c0: -98942718328354657571921135515000654102809275280
c1: -1094123364329131040120216467254391237420
c2: 65859755925338247464662996618032
c3: -1618738878212462977815579
c4: -4830711805039576
c5: 687420
Y0: -43389796499263425737859091822712579063
Y1: 47084851870118530061
# MurphyE (Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16) = 8.25e-15
# lognorm 62.48

n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
skew: 81166336.0
c0: 45557711067352835790870412530520906044046352000
c1: 538182327666734598053665140544142650228
c2: -35247498407001430812508122651020
c3: -85598128594418788150756
c4: 4204857918193797
c5: 1081860
Y0: -39627605646451711638391781283296477421
Y1: 9440158732784987191
# MurphyE (Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16) = 8.21e-15
# lognorm 61.04

n: 105720747827131650775565137946594727648048676926428801477497713261333062158444658783837181718187127016255169032147325982158719006483898971407998273736975091062494070213867530300194317862973608499
skew: 54214656.0
c0: 21724609728591871248210729947586274756512302379
c1: -317862276216193685941919711583925203153
c2: -10170954010451741993810971872497
c3: 13022944940440380507497
c4: -8295877607702306
c5: 1175280
Y0: -38976584207204966069784298467971861060
Y1: 559932214212309221
# MurphyE (Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16) = 8.17e-15
# lognorm 61.54

I haven't done any test sieving on any poly, CADO or otherwise, for this search.

Note that CADO's default "rootsieve effort" is 5, on a scale of 1-10; one line of inquiry for the next poly search is to see if effort 10 over say, 50 or 100 hits is better than effort 5 over 500 hits.

And btw, I definitely vote for taking the log of Murphy E in much the same way CADO takes the log of the sizeopt norm.

[Max0526]

Originally Posted by Dubslow

2) I've just discovered that the binaries produced by building CADO include polyselect, which does "stage 1" and size opt at the same time, as well as a separate sopt binary which takes "raw" polynomials (presumably the results of the hashtable search...?) for size opting as well, and the equivalent output of either binary can then be fed to polyselect_ropt whose function should be obvious in this context.

On your machine, what is the location of the separate sopt binary? I use cado-nfs-2.2.0 and can't find sopt binary anywhere.

I need help size-optimizing the following msieve poly (http://stdkmd.com/nrr/c.cgi?q=68881_216) in cado:

# Murphy_E = 4.296e-12
# expecting poly E from 4.51e-012 to > 5.19e-012; sieved all c5 < 564
n: 6964469412395724414400887219282737626566096536914345017855771078382532118955354745772117936678678317022004114419569353913029298176668228168002524415033
Y0: -719974750345759133972843201212
Y1: 2509148016863593
c0: -645947407201258180676614400938960976175
c1: 123446188733445544772289511592948
c2: 18319373780518695347898595
c3: -2567126999321415502
c4: -18908689182
c5: 36
skew: 32970491.17

I don't think I can recreate this polynomial in cado because Y1=2509148016863593=7.31.10223.11447.98809 and 98809>2*11447. Maybe somebody knows the way to recreate it?

When I pass the msieve polynomial to polyselect_ropt, the best output is:

# Best polynomial found:
n: 6964469412395724414400887219282737626566096536914345017855771078382532118955354745772117936678678317022004114419569353913029298176668228168002524415033
Y0: -719974797818859686216157669956
Y1: 2509148016863593
c0: -74098902767640540677544413831582460335175
c1: 3710094808934559412128037232637690
c2: 120979270663498122139675645
c3: -1007248783970378638
c4: -22314290622
c5: 36
skew: 63127552.000
# lognorm 48.59, alpha -7.30 (proj -1.24), E 41.28, 5 real roots
# MurphyE(Bf=1.00e+07,Bg=5.00e+06,area=1.00e+16)=4.12e-12

After that, I can change the skew to 83758158.21476 (that increases the lognorm I assume) which gives Murphy_E=4.14415907e-12. (http://myfactors.mooo.com/, Optimal Skew)

I didn't do actual sieving, just trying to compare Murphy_E scores.

[Dubslow]

bill@Gravemind⌚2258 ~/cado/build/Gravemind/polyselect ∰∂ ls | grep -v "[Cc]1"
cadops.bash
CMakeFiles
cmake_install.cmake
CTestTestfile.cmake
dlpolyselect
libpolyselect_common.a
Makefile
ms_to_cado.py
num
polyselect
polyselect_ropt
rotate
rotate_all
sopt
sort_polys.py
sys

Importing other polys can be a bit difficult, you'll have to play around with the individual binaries to figure out how they work (or read code or view the -help output).

For instance I wrote this small bash script to glue things together for a start-to-stop polyselect run:

bill@Gravemind⌚2300 ~/cado/build/Gravemind/polyselect ∰∂ cat cadops.bash
#! /usr/bin/env bash

n=51703345090678007603283094883953312763757798669035713906213433384993369348322376058745609326491127794881496108288948261956762322783390589392048698391

name=C149_3408_1631
threads=8

# Typically you should copy these from cado/parameters/factor/params.cxxx
deg=5
P=500000
admin=0
admax=20000000
incr=60
nq=1000
keep=900

##############################################################################

# http://stackoverflow.com/a/9194117/1497645
# Round admin up to nearest multiple of incr
admin=$(( ((admin + incr - 1) / incr) * incr ))

f1="$name.sizeopt"
f2="$name.rootopt"

nice -n 19 ./polyselect -degree "$deg" -P "$P" -admax "$admax" -nq "$nq" -keep "$keep" -t "$threads" -N "$n" | tee "$f1"

cmd="./sort_polys.py $f1 $keep"

$cmd && nice -n 19 ./polyselect_ropt -inputpolys "$f1.sorted" -t "$threads" | tee "$f2"

./sort_polys.py "$f2"

It also makes use of the following small Python script, which is largely lifted from CADO's own Python glue code (which I otherwise found too inflexible for my taste):

bill@Gravemind⌚2300 ~/cado/build/Gravemind/polyselect ∰∂ cat sort_polys.py 
#! /usr/bin/env python3

import sys

sys.path.append('/home/bill/cado/scripts/cadofactor')

from cadotask import Polynomials, PolynomialParseException

def read_blocks(input):
        """ Return blocks of consecutive non-empty lines from input
        
        Whitespace is stripped; a line containing only whitespace is
        considered empty. An empty block is never returned.
        
        >>> list(Polysel1Task.read_blocks(['', 'a', 'b', '', 'c', '', '', 'd', 'e', '']))
        [['a', 'b'], ['c'], ['d', 'e']]
        """
        block = []
        for line in input:
            line = line.strip()
            if line:
                block.append(line)
            else:
                if block:
                    yield block
                block = []
        if block:
            yield block

def parse_block(block):
     try:
          poly = Polynomials(block)
     except:
          print(block)
          raise
     if not poly:
          raise ValueError("useless poly:\n{}".format(block))
     return poly

def sort(filename, keep=None):
     lst = []
     with open(filename, 'r') as f:
          for block in read_blocks(f):
               if '### root-optimized polynomial' in block[0]:
                    # For whatever reason, all output from `polyselect_ropt` is prefaced by "# ",
                    # so remove it
                    block = [line[2:] for line in block if '# Stat:' not in line]
                    # The very last output poly has no blank line between it and stat output, so
                    # knock that out here
                    lst.append(parse_block(block))
               elif any('# Size-optimized polynomial:' in line for line in block):
                    lst.append(parse_block(block))
          print("Read {} polys from {}".format(len(lst), filename))
          have_murphyE = bool(lst[0].MurphyE)
          if have_murphyE:
               if not all(poly.MurphyE for poly in lst):
                    raise ValueError("Some polys have MurphyE, some don't")
               key, reverse = lambda poly: poly.MurphyE, True
          else:
               if any(poly.MurphyE for poly in lst):
                    raise ValueError("Some polys have MurphyE, some don't")
               key, reverse = lambda poly: poly.lognorm, False
          lst.sort(key=key, reverse=reverse)
          if keep:
               lst = lst[:keep]
          with open(filename+'.sorted', 'w') as g:
               for poly in lst:
                    poly = '\n'.join(l for l in str(poly).splitlines() if '(x) =' not in l)
                    g.write(poly+'\n\n')

if __name__ == '__main__':
     filename = sys.argv[1]
     keep = None
     if len(sys.argv) > 2:
          keep = int(sys.argv[2])
     sort(filename, keep)

Note that ./polyselect does both first stage and size opt at the same time. You'll again have to play with the binary/review code to figure out how to separate the two stages (I believe it is possible, but I might be remembering incorrectly).