From a1fefc07cc9f3616b29e8879751401226b5f4ea5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Foll=C3=A1th=20J=C3=A1nos?= Date: Fri, 13 Jun 2014 12:28:08 +0800 Subject: [PATCH] Boundary function optimization working. --- README | 77 +++++- TODO | 72 ++--- app/CMakeLists.txt | 4 + app/boundtool.cpp | 82 ++++-- app/eprune.cpp | 202 ++++++++++++++ app/genlattice.cpp | 32 ++- app/timing.cpp | 94 ++++++- debug.pdf | Bin 10830 -> 10669 bytes include/cleanbkz/cjloss.hpp | 12 +- include/cleanbkz/tools.hpp | 15 +- include/cleanbkz/version.hpp | 6 +- ratio.tmp | 2 +- src/boundary.cpp | 493 +++++++++++++++-------------------- src/cjloss.cpp | 54 ++-- src/enumeration.cpp | 135 +++++++++- src/tools.cpp | 20 -- test/CMakeLists.txt | 12 +- test/act_test.cpp | 168 ++++-------- test/unit_tests.cpp | 186 +++++-------- xaa | 460 ++++++++++++++++++++++++++++++++ 20 files changed, 1413 insertions(+), 713 deletions(-) create mode 100644 app/eprune.cpp create mode 100644 xaa diff --git a/README b/README index a0cbacf..f152c5d 100644 --- a/README +++ b/README @@ -2,10 +2,8 @@ cleanbkz is a library to find short vectors in a lattice. Once finished it will contain an improved enumeration algorithm using extreme pruning, the bkz2.0 and probably an improved random sampler algorithm. -To compile and test it, you can follow these steps: - ------------------------ -1. Compile cleanbkz sources +1. Compiling the cleanbkz sources ------------------------ cleanbkz requires the NTL library. Version 5.4.2 (or higher) is needed. @@ -15,11 +13,70 @@ $ cmake . $ make Now both the library and the tests are compiled. You can find the test -executables in the bin/ folder. +executables and the command line interface in the bin/ folder. + +The scripts plot.sh and debug.sh require gnuplot to function properly +and visualize the generated data. + + +------------------------ +2. The command line interface +------------------------ + +Every command line tool has a built-in help accessible with the -h switch. +Examples: + +0) Measuring the rate of reduction and enumeration: + +./bin/timing -d 60 -k 20 -n 960 + +Wait until it finishes. On a AMD Opteron(tm) 3380 it took 15 minutes and +gave the following result: + +Measuring dimension 60: +t_node= 6.26296e-08 +t_reduc= 0.855437 + + +1) Generating a lattice: + +./bin/genlattice -d 60 -c 0.94 -s 11 -k 20 > data/d60s11k20c0_94.lat + + +2) Generating a boundary function: + +./bin/boundtool -f data/d60s11k20c0_94.lat -n 6.26296e-08 -r 0.855437 -l 60 > data/d60s11k20c0_94.plt + +Wait until it finishes. On a AMD Opteron(tm) 3380 it took 20 minutes. You +also can generate a pdf chart of the bounding function: + +./plot.sh data/d60s11k20c0_94.plt + +This bounding function will not yet be optimal. To get a close to optimal +bounding function run the optimization a little longer: + +./bin/boundtool -f data/d60s11k20c0_94.lat -n 6.26296e-08 -r 0.855437 -l 60 -c 3000 > data/d60s11k20c0_94_tri.plt + + +3) Running the actual enumeration: + +The file produced in the previous step contains data for gnuplot and +other information like the bounding function in NTL's format, the +predicted success probability and the expected running time. + +The following program requires the boundary function in NTL's format, +so first we need to copy the line (without the starting '#' sign) in +a new text file. In this example the name of this new text file is +'d60s11k20c0_94.bnd'. Now issue the command + +./bin/eprune -l data/d60s11k20c0_94.lat -b data/d60s11k20c0_94.bnd + +and wait until it finishes. On a AMD Opteron(tm) 3380 it took around +89 seconds, fairly close to the promised 95 seconds. ------------------------ -2. Generate documentation +3. Generate API documentation ------------------------ The source files are annotated for doxygen documentation. To generate @@ -31,14 +88,16 @@ Now the documentation is in the doc/ folder. ------------------------ -3. File extensions and formats +4. File extensions and formats ------------------------ -.lat These are the data files generated by the lattice_gen tool and +.lat These are the data files generated by the genlattice tool and contain the lattice basis in NTL's output format. -.plt Boundary function generated by the boudtool application in a form - readable by gnuplot. +.plt Boundary function generated by the boundtool application in a + form readable by gnuplot. .bnd Boundary function in NTL's format. It is usually created by copying the corresponding line from a .plt file. + + diff --git a/TODO b/TODO index cbaa14a..ae002c8 100644 --- a/TODO +++ b/TODO @@ -1,65 +1,29 @@ -- Minimális pruning function keresését kikisérletezni: - - tesztelni a gyakorlatban a képletemet állitható kezdőértékkel (dupla és sima random lattice-el is) - - konstans pruninggal megkeresni a megfelelő arányokat az egyes dimenziókban (dupla és sima random lattice-el is) - - több különböző lattice-ra és összehasonlitani - -- a random lattice-s tesztprogramból automatikus teszteket irni - - node számlálás full (double) - - extreme pruning tesztet külön lefuttatni, hogy ha van még valami gond, akkor az itt kiderüljön - - node számlálás pruned (double) - - node számlálás full - - node számlálás pruned -- elküldeni a többieknek a linket a github projektre -- implementálni az adaptiv a keresést: kezd egy viszonylag nagy deltával, megy amig nem tud javitani, osztja a deltát 10-el majd szintén megy amig nem tud javitani és igy tovább (betenni a unit tesztbe) + +- elküldeni a többieknek a linket a github projektre és a generált teszt-adatokat + - Közvetlenül előtte persze commitolni + +************ Szünet ******************* +- optimalizálni a számitásokat: + - a négyzetreemeléseket kihozni amennyire csak lehet + - előre számitani és táblázatból olvasni + - az egységgömbök térfogatát + - a faktoriálisokat + - a normalizáló konstansokat + - kikisérletezni, hogy mi az a legkisebb RR pontosság amire a 120-150 dimenziós számitás stabil + - letekerni a pontosságot amennyire csak lehet +- implementálni az adaptiv a keresést: kezd egy viszonylag nagy deltával, megy amig nem tud javitani, osztja a deltát 10-el majd szintén megy amig nem tud javitani és igy tovább +- print scriptet egyesiteni és paraméterezhetőre megirni + +******** Itt lesz kész az extreme pruning ***** - Megnézni az early termination-ös cikkeket, hogy milyen blokkméretet lenne érdemes használni egy 128 dimenziós bázis előfeldolgozásához - Megirni a pruningos BKZ-t, enumerationt külön függvénybe tenni - Minden enumeration hivást és a vonatkozó unit teszteket átirni az új enumeration-re - kikisérletezni a paramétereket a bkz automatikus bounding function keresőjéhez - Implementálni a radius reductiont és az early terminationt - Implementálni az unexpected early terminationt (mindig kiirja lemezre, hogy hanyadik iterációnál jár és hogy mi az aktuális bázis) -Kérdések: -- miért használt a bkz 2.0-nál 10%-os sikeresélyűt, amikor nekem 95% körüliek sikerülnek és az is baromi gyorsra igéri magát? -- miért más a formája az én boundary functionjeimnek? -Opcionálisan: -- Megnézni, hogy milyen hatással van a mu korrigálása az enumerationra - - valszeg semmilyen, mert a főátlót eleve 1-nek veszi. Megkeresni a vonatkozó részt az algoritmusban -+ átirni, hogy mindig az R_i négyzetekkel dolgozzon és ne kelljen gyököt vonogatni. - Nem éri meg: ugyanolyan macerás, mert a gömbök sugaránál R_i kell. -- hermite faktoros optimalizálást implementálni -- refaktorálni a toolok paraméterfeldolgozását -- letisztázni és kibőviteni a unit teszteket -- a unit tesztekhez használt referenciaimplementációkat kiszedni a libből a teszt könyvtárba -- a boundtool.cpp -ek és a boundary.cpp -nek értelmesebb neveket találni -- Tesztet irni, hogy a double sqrt négyzete egyezik-e az eredetivel (különbség 0) + ******** Itt lesz kész a BKZ1.5 ********** - a rekurzive hivás paramétereit kikisérletezni : - Hermite faktorokat lemérni (early terminatiönnel és anélkül) - várható futási időket a bounding function számolóval összevetni -- random lattice-ekre is elvégezni a méréseket és összehasonlitani a cjlossosokkal (ehhez előtte megnézni Gama-Nguyent) -- saját namespace-t csinálni hasonlón dinamikusan mint az NTL-ben -- polytopvol unit tesztet kiegésziteni 15 dimenzióra - megirni a cikket és a megjelenés után átirni a dokumentációban a hivatkozásokat amiket lehet a saját cikkemre -- Általános BKZ-t implementálni: LLL -el redukálni bkz helyett és úgy lemérni az időket (10000000 node kell a viszonylagos stabilitáshoz és 45-nél már elég gyakran nem éri ezt el, szóval az az eredmény nem lesz olyan pontos) - + paraméterezni a programot a sample mérettel (hogy tudjam becsülni a magas mintás futási időt alacsony mintamérettel), futási időket megbecsülni - - 50 dimenzióval kikisérletezni a sample méretet ami kell a 3 stabil tizedesjegyhez - - lemérni 50, 55, 60, 65, 70, 75, 80 -ra (automatizálva) - - extrapolálni 40, 45 -re (automatizálni) -- Általános bounding function-öket találni (erősen opcionális: lehet róla irni, de kicsi a siker valószinűsége és nem egy nagy eredmény) - - GSA feltételes bounding function keresőt implementálni - - tesztelni, hogy hogyan változik a bounding function a gs hossz változtatásával - - tesztelni, hogy hogyan változik a futási idők változtatásával - -Problémák: -- az NTL szabványos függvényeit használva külön kell kiszámolni a GS együtthatókat, pedig azokat már a redukcióhoz kiszámolja. Lehetséges megoldások: - - igy hagyni és a számitásokat úgy végezni mintha nem lenne ott (ugyanis magas dimenziókban tényleg elhanyagolható a költsége) - - duplikálni és módositani a szükséges NTL függvényeket -- a cjloss bázisát publikussá kellett tennem, hogy működjenek a dolgok, lehet, hogy át kéne szerveznem, hogy a bázis ne csak adattag legyen, hanem a mat_ZZ kiterjesztése - vagy nem: az egész csak egy teszt, nem is kell a libbe, át fogom tenni a tests könyvtárba -- a méréshez egyéb utasitásokat is bele kell tennem az algoritmusba. Ez tényleges futás esetén lassithatja és olyan paramétereket is megkövetel, amik egy sima enumeration hiváskor nem kellenek. Lehetőségek - - bevállalom az ezzel járó lassulást (mennyivel is lesz ettől lassabb? (méréshez rögziteni a cjloss seedet)) - - duplikálom a kódot -- alacsony dimenzióban többnyire csak nagyon kevés node-ot kell megvizsgálnia, igy olyan kevés időt tölt vele, hogy nem tudja mérni. Lehetséges megoldások: - - ha megtalálta a legrövidebb vektort akkor resetelem és futtatom amig elég idő nem lesz -nagyon torzit - - magasabb dimenziós mérésekből approximálom -- bármilyen dimenzióban változó és bázistól függő a futási idő, kérdés, hogy hogyan rendeljek egy adott dimenzióhoz egy konkrét értéket. Tapasztalati várható érték (átlag tűnik a legjobbnak), de kérdés, - - hogy mennyi mintát vegyek alapul - - milyen (cjloss vagy teljesen random) lattice-t vegyek alapul (majdnem mindegy, mert az LLL redukció úgyis hasonló formára hozza) -- eddig b_0^* hosszát használtam kezdő sugárnak. Mi az amit a BKZ használ és mi az amit Phong használt a BKZ 2.0-ában? A méréseknek is úgy lenne értelme, ha onnan inditanám (és akkor a kicsik is lehet hogy végigmennek) diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt index e64d0f6..6a57689 100644 --- a/app/CMakeLists.txt +++ b/app/CMakeLists.txt @@ -1,6 +1,7 @@ SET(TARGET1 timing) SET(TARGET2 boundtool) SET(TARGET3 genlattice) +SET(TARGET4 eprune) INCLUDE_DIRECTORIES(${PROJECT_SOURCE_DIR}/include) LINK_DIRECTORIES(${LIBRARY_OUTPUT_PATH}) @@ -12,3 +13,6 @@ TARGET_LINK_LIBRARIES(${TARGET2} cleanbkz ntl) ADD_EXECUTABLE(${TARGET3} genlattice.cpp) TARGET_LINK_LIBRARIES(${TARGET3} cleanbkz ntl) + +ADD_EXECUTABLE(${TARGET4} eprune.cpp) +TARGET_LINK_LIBRARIES(${TARGET4} cleanbkz ntl) diff --git a/app/boundtool.cpp b/app/boundtool.cpp index 2ecee61..74b83ca 100644 --- a/app/boundtool.cpp +++ b/app/boundtool.cpp @@ -30,8 +30,7 @@ using namespace std; // Defined in boundary.cpp -// TODO: ehelyett majd egy get_p_succ nevűt irni, ami megcsinálja amit kell -double polytope_volume(double vols[], double bounds[], int dim); +extern double ball_vol(int k, double r); // Defined in tools.cpp char* get_cmd_option(char** begin, char** end, const string& option); @@ -39,7 +38,7 @@ bool cmd_option_exists(char** begin, char** end, const string& option); int main(int argc, char** argv) { - + double v= 0; double t_node= 0; double t_reduc= 0; double delta= 1e-1; @@ -56,8 +55,8 @@ int main(int argc, char** argv) { << "\t-n n\t\tThe (avarage) time the enumeration algorithm takes to process a single node. (required)" << endl << "\t-r n\t\tThe (avarege) running time of the preprocessing reduction algorithm. (required)" << endl << "\t-c n\t\tNumber of random changes to make during the numerical optimization. (default: 1000)" << endl - << "\t-t n\t\tThe number of unchanged rounds to require before exiting. It is ignored when -c is given." << endl - << "\t-d n\t\tThe size of the single random changes to meke during the optimization. (default: 0.1)" << endl; + << "\t-d n\t\tThe size of the single random changes to meke during the optimization. (default: 0.1)" << endl + << "\t-l n\t\tOptimize for shortest vector with length square root of n. (the Gaussian heuristic is default)" << endl; return 0; } @@ -120,6 +119,18 @@ int main(int argc, char** argv) { } } + act_arg= get_cmd_option(argv, argv + argc, "-l"); + if (act_arg) { + ss << act_arg; + ss >> v; + ss.clear(); + if(v <= 0) { + cerr << "ERROR: invalid shortest vector length. The shortest vector length should be greater than zero. Aborting." << endl; + return 1; + } + } + + act_arg= get_cmd_option(argv, argv + argc, "-f"); if (act_arg) { ifstream basis_file(act_arg); @@ -145,32 +156,49 @@ int main(int argc, char** argv) { int dim= mu1.NumRows(); double* boundary= new double[dim]; - //length of the shortest vector in the cjloss latticr - double R= sqrt(dim-1); + cout << "# Generated with cleanbkz " << CBKZ_VERSION << endl + << "# Copyright (C) 2014 Janos Follath" << endl + << "# This is free software with ABSOLUTELY NO WARRANTY." << endl << "#" << endl; + + //length of the shortest vector in the cjloss lattice + if(v>0) { + v= sqrt(v); + cout << "# Using supplied lambda : " << v << endl; + } + else { + double tmp,gh= 1; + double* gsghs= new double[dim]; + for(int i= 0; i < mu1.NumRows(); i++) { + conv(tmp, c1[i]); + gh*= tmp; + gsghs[i]= pow(gh/ball_vol(i+1, 1),1.0/(i+1)); + } + delete gsghs; + gh= pow(gh/ball_vol(dim, 1),1.0/dim); + cout << "# No lambda supplied, using Gaussian heuristic: " << gh << endl; + } + + // TODO: csinálni egy változatot, ahol nem iterationt hanem thressholdot adunk meg double p_succ; double t_enum; - generate_boundary(c, t_node, t_reduc, dim, boundary, R, delta, iterations, p_succ, t_enum, false); - - cout << "# Generated with cleanbkz " << CBKZ_VERSION << endl - << "# Copyright (C) 2014 Janos Follath" << endl - << "# This is free software with ABSOLUTELY NO WARRANTY." << endl << "#" << endl; - - cout << "# basis: '" << act_arg << "' " << endl - << "# estimated enumeration time: " << t_enum << endl - << "# success probability: " << p_succ << endl - << "# boudary function: " << endl; - - vec_RR out; - out.SetLength(dim); - for(int i= 0; i < dim; i++) - out[i]= boundary[i]; - cout << "# " << out << endl << endl; - - for(int i= 0; i < dim; i++) - cout << i << " " << boundary[i] << endl; - cout << endl; + generate_boundary(c, t_node, t_reduc, dim, boundary, v, delta, iterations, p_succ, t_enum, false); + + cout << "# basis: '" << act_arg << "' " << endl + << "# estimated enumeration time: " << t_enum << endl + << "# success probability: " << p_succ << endl + << "# boudary function: " << endl; + + vec_RR out; + out.SetLength(dim); + for(int i= 0; i < dim; i++) + out[i]= boundary[i]; + cout << "# " << out << endl << endl; + + for(int i= 0; i < dim; i++) + cout << i << " " << boundary[i] << endl; + cout << endl; return 0; } diff --git a/app/eprune.cpp b/app/eprune.cpp new file mode 100644 index 0000000..c7000fa --- /dev/null +++ b/app/eprune.cpp @@ -0,0 +1,202 @@ +/* + Copyright (C) 2014 Janos Follath. + + This file is part of cleanbkz. + + cleanbkz is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + cleanbkz is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with cleanbkz. If not, see . + +*/ + +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +extern void predict_nodes_RR(RR Rvec[], double b_star_norm[], int n); +extern void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result, unsigned long &termination, double &time); + +extern void profile_enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result); + +mat_ZZ unimod(int dim){ + mat_ZZ L, U; + L.SetDims(dim, dim); + U.SetDims(dim, dim); + + for(int i= 0; i < dim; i++) + L[i][i]= U[i][i]= 1; + + for(int i= 0; i < dim; i++) + for(int j= i+1; j < dim; j++) { + U[i][j]= RandomBnd(2); + L[j][i]= RandomBnd(2); + } + + return L*U; +} + +int main(int argc, char** argv) { + mat_ZZ basis; + vec_RR bnd; + char* act_arg; + int beta= 0; + + + if (argc==1 || cmd_option_exists(argv, argv+argc, "-h")) { + cout << "This program performs lattice enumeration with extreme pruning:" << endl + << "\t-h \t\tPrint this help." << endl + << "\t-v \t\tPerforms a single enumeration round and compares the nodes to the prediction." << endl + << "\t-l filename\tReads the lattice from the file filename. (This option is mandatory)" << endl + << "\t-b filename\tReads the boundary function from the file filename. (This option is mandatory)" << endl; + return 0; + } + + act_arg= get_cmd_option(argv, argv + argc, "-l"); + if (act_arg) { + ifstream basis_file(act_arg); + if (basis_file.is_open()) { + basis_file >> basis; + basis_file >> beta; + if(beta<2) { + cerr << "ERROR: the reduction block size given in the lattice file can't be less than 2. Aborting." << endl; + return 1; + } + } else { + cerr << "ERROR: can't open input file: '" << act_arg << "'. Aborting." << endl; + return 2; + } + } + else { + cerr << "ERROR: The option '-l' is mandatory. Aborting." << endl; + return 1; + } + + act_arg= get_cmd_option(argv, argv + argc, "-b"); + if (act_arg) { + ifstream bnd_file(act_arg); + if (bnd_file.is_open()) + bnd_file >> bnd; + else { + cerr << "ERROR: can't open input file: '" << act_arg << "'. Aborting." << endl; + return 2; + } + } + else { + cerr << "ERROR: The option '-b' is mandatory. Aborting." << endl; + return 1; + } + + mat_RR mu1; + vec_RR c1; + ComputeGS(basis,mu1,c1); + + double* c= new double[mu1.NumRows()]; + for(int i= 0; i < mu1.NumRows(); i++) + conv(c[i], c1[i]); + + double** mu= new double*[mu1.NumRows()]; + for(int i= 0; i < mu1.NumRows(); i++) + mu[i]= new double[mu1.NumCols()]; + for(int i= 0; i < mu1.NumRows(); i++) + for(int j= 0; j < mu1.NumCols(); j++) + conv(mu[i][j], mu1[i][j]); + + double* boundary= new double[basis.NumRows()]; + for(int i= 0; i < basis.NumRows(); i++){ + conv(boundary[i], bnd[i]); + boundary[i]*= boundary[i]; + } + + vec_RR solution; + + if (cmd_option_exists(argv, argv+argc, "-v")) { + profile_enumerate_epr(mu, c, boundary, basis.NumRows(), solution); + + for(int i= 0; i < mu1.NumRows(); i++) + conv(c[i], sqrt(c1[i])); + + predict_nodes_RR(bnd.elts(), c, basis.NumRows()); + + return 0; + } + + clock_t begin; + unsigned long termination= 0; + double t_enum, t_reduc, t_all; + t_all= t_enum= t_reduc= 0; + enumerate_epr(mu, c, boundary, basis.NumRows(), solution, termination, t_enum); + t_all+= t_enum; + + while(solution.length()==0) { + cout << "Enumeration failed. Randomizing lattice basis..." << endl; + + begin= clock(); + basis= basis*unimod(basis.NumRows()); + BKZ_QP1(basis, 0.99, beta); + + ComputeGS(basis,mu1,c1); + + for(int i= 0; i < mu1.NumRows(); i++) + conv(c[i], c1[i]); + + for(int i= 0; i < mu1.NumRows(); i++) + for(int j= 0; j < mu1.NumCols(); j++) + conv(mu[i][j], mu1[i][j]); + + t_reduc+= clock()-begin; + + cout << "Done." << endl; + cout << "Performing enumeration..." << endl; + + begin= clock(); + enumerate_epr(mu, c, boundary, basis.NumRows(), solution, termination, t_enum); + t_all+= clock()-begin; + cout << "Done." << endl; + } + + cout << "Enumeration successful. Shortest vector coordinates in the basis: " << endl << solution << endl; + + double* sol= new double[basis.NumRows()]; + for(int i= 0; i < basis.NumRows(); i++) + sol[i]= 0; + + double dbase, dsol; + for(int i= 0; i< basis.NumRows(); i++) + for(int j= 0; j< basis.NumCols(); j++) { + conv(dbase, basis[i][j]); + conv(dsol, solution[i]); + sol[j]+= dbase*dsol; + } + + cout << "Coordinates in the lattice: " << endl << "["; + for(int i= 0; i< basis.NumRows()-1; i++) + cout << sol[i] << " "; + cout << sol[basis.NumRows()-1] << "]" << endl; + + double sqrdLength= 0; + for(int i= 0; i< basis.NumRows(); i++) + sqrdLength+= sol[i]*sol[i]; + + cout << "Squared length of the shortest vector: " << sqrdLength << endl << endl; + + cout << "Time spent with enumeration: " << t_all/CLOCKS_PER_SEC << endl; + cout << "Time spent with reduction: " << t_reduc/CLOCKS_PER_SEC << endl; + cout << "Total running time: " << (t_all + t_reduc)/CLOCKS_PER_SEC << endl; + cout << "(Note that the basis was already reduced and thus the time of the first reduction is not included. Also the enumeration subrutine exits after finding the first acceptable vector and does not search through the whole pruned enumeration tree.)" << endl; +} + diff --git a/app/genlattice.cpp b/app/genlattice.cpp index 4347a14..a50ff55 100644 --- a/app/genlattice.cpp +++ b/app/genlattice.cpp @@ -23,16 +23,27 @@ #include #include #include +#include #include using namespace std; -// Defined in tools.cpp -char* get_cmd_option(char** begin, char** end, const string& option); -bool cmd_option_exists(char** begin, char** end, const string& option); - //TODO: refaktorálni a paraméterfeldolgozást: csinálni hozzá int és double változatot, hibaüzenetekkel és kivételkezeléssel +void gen_randlat(mat_ZZ& basis, ZZ determinant, int dim) { + basis.SetDims(dim, dim); + + ZZ s; + s= time(NULL); + SetSeed(s); + for(int i= 0; i < dim; i++) { + basis[i][0]= RandomBnd(determinant); + basis[i][i]= 1; + } + + basis[0][0]= determinant; +} + int main(int argc, char** argv) { int maxsize= 0; @@ -143,19 +154,18 @@ int main(int argc, char** argv) { mat_ZZ m; if (cmd_option_exists(argv, argv+argc, "-c")) { cjloss l(dim, density, seed); - m= l.basis; + m= l.get_basis(beta); } else if (cmd_option_exists(argv, argv+argc, "-m")) { - m.SetDims(dim, dim); - for(int i= 0; i < dim; i++) - for(int j= 0; j < dim; j++) - RandomBits(m[i][j], maxsize); + ZZ determinant; + GenPrime(determinant,maxsize); + gen_randlat(m,determinant,dim); } - if(beta > 1) + if(beta > 1 && !(cmd_option_exists(argv, argv+argc, "-c"))) BKZ_QP1(m, 0.99, beta); - cout << m; + cout << m << endl << beta; return 0; } diff --git a/app/timing.cpp b/app/timing.cpp index 3f827b9..755f1d7 100644 --- a/app/timing.cpp +++ b/app/timing.cpp @@ -21,20 +21,24 @@ #include #include #include -#include +#include +#include +#include #include using namespace std; +using namespace NTL; -// Defined in tools.cpp -char* get_cmd_option(char** begin, char** end, const string& option); -bool cmd_option_exists(char** begin, char** end, const string& option); +//TODO: random lattices + +void measure_epr(int dimension, int beta, unsigned long samples, unsigned long termination, + double& t_node, double& t_reduc, unsigned long& zeroes, unsigned long& bias); int main(int argc, char** argv) { int dimension= 0; int bkz_blocksize= 2; int nodes= 1000000; - int samples= 1000; + int samples= 100; int start= 0; int end= 0; stringstream ss; @@ -42,15 +46,15 @@ int main(int argc, char** argv) { clock_t begin, finish; if (argc==1 || cmd_option_exists(argv, argv+argc, "-h")) { - cout << "This program measures the parameters t_node and t_reduc required for the computation of boundary functions. These are the running time of the enumeration and reduction algorithms on the current platform. Program options:" << endl + cout << "This program uses cjloss lattices to measure the parameters t_node and t_reduc required for the computation of boundary functions. These are the running time of the enumeration and reduction algorithms on the current platform (in the case the cjloss lattices the dimension of the preprocessing is one higher than the enumeration). Program options:" << endl << "\t-h \t\tPrint this help." << endl << "\t-d n\t\tMeasure running times in dimension n." << endl << "\t-s n\t\tStart the measurements in dimension n. It is ignored when -d is given." << endl << "\t-e n\t\tContinue measuring in all dimensions until dimension n with a step five. It is ignored when -d is given." << endl << "\t-b n\t\tAbort enumeration after processing n nodes. (default: 10,000,000)" << endl - << "\t-n n\t\tNumber of experiments to make. (default: 1000)" << endl + << "\t-n n\t\tNumber of experiments to make. (default: 100)" << endl << "\t-k n\t\tThe blocksize of BKZ used for preprocessing. (default: 2)" << endl - << "\t-t n\t\tPredict the running time of the experiments." << endl; + << "\t-t\t\tPredict the running time of the experiments." << endl; return 0; } @@ -153,7 +157,8 @@ int main(int argc, char** argv) { cout << ((double)(finish - begin) / CLOCKS_PER_SEC)/10*samples << "s" << endl; total+= (unsigned long int) (((double)(finish - begin) / CLOCKS_PER_SEC)/10*samples); } - cout << "Expected total running time of the experiments: " << total/3600 << "h" << (total%3600)/60 << "m" << (total%36000)%60 << "s" << endl << endl; + cout << "Expected total running time of the experiments: " << total/3600 << "h" + << (total%3600)/60 << "m" << (total%36000)%60 << "s" << endl << endl; } for(int i= start; i <= end; i+=5) { @@ -165,7 +170,76 @@ int main(int argc, char** argv) { } } - cout << "(Sometimes the enumeration doesn't have to inspect the prescribed number of nodes or even nodes enough to make the runing time measureable, the number of these cases are in the \"biased/zeroes\" line. The first number measures only the biased cases with nonzero running time.)" << endl << endl; + cout << "(Sometimes the enumeration doesn't have to inspect the prescribed number of nodes or even nodes enough " + << "to make the runing time measureable, the number of these cases are in the \"biased/zeroes\" line." + << "The first number measures only the biased cases with nonzero running time.)" << endl << endl; return 0; } + +extern void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result, + unsigned long &termination, double &time); + +void measure_epr(int dimension, int beta, unsigned long samples, unsigned long termination, + double& t_node, double& t_reduc, unsigned long& zeroes, unsigned long& bias) { + double enu_time= 0; + cjloss* lattice; + mat_ZZ tmp,basis; + mat_RR mu1; + vec_RR c1; + clock_t begin, end; + vec_RR result; + unsigned long nodes; + + t_node= t_reduc= 0; + zeroes= bias= 0; + + double** mu= new double*[dimension]; + for(int i= 0; i < dimension; i++) + mu[i]= new double[dimension]; + + double* c= new double[dimension]; + + double* Rvec= new double[dimension]; + for(int i= 0; i < dimension; i++) + Rvec[i]= dimension; + + for(unsigned long i= 0; i < samples; i++) { + lattice= new cjloss(dimension,0.94,i); + + tmp= lattice->get_basis(0); + begin= clock(); + BKZ_QP1(tmp, 0.99, beta); + end= clock(); + t_reduc+= (double)(end - begin) / CLOCKS_PER_SEC; + + basis.SetDims(tmp.NumRows()-1, tmp.NumRows()-1); + for(int i= 0; iz(;9Gi&X=XMLDiGeef;7LSd9R=W=sTiq!ZUU5#|OW4A$0Y^^A7_*7XbRX+VO>HE<}UZr1#VV95N;zzo%mL4X%NX3xv zZ4tIS4h?@zeIi=ic51B*ed8`udA&{AF_J70!J@G-wzBjFf{MvIIK(-ow?*#plRNml z99{L0300hprVAZ>0}6hjgZ%2%?xTeM=!cBE8skQQ@d5FMR#)~>KoRN+ikw} zpRK6P;E2(1JC>PG_Zbj0c?@r}VzWB7<@;mm;-5BbSxLOdGs%+DPPHCCpgjflL9Qor zyAvZN8y{0?yCT`$PBcu<8IwMy((s2J`Bv;*xmV$u5rR4-ul>b zUR}}e=n-Zs5;2u_Y^t|k>^ml(Cws0O23lA~(8+cjVLe`9kJ2Tr>8kg%UOf^3Zgo)E zG}t8bU@tScPrEHn?uq-iX8}b$F)ZBYum1vrMWCu=`;+$je$0*7oj&I{70;t|p@iTR z>^~CpAzgcoDIh$Tf+j~xgHXr#|4}ynf)}~|S32!{jG~TS2bgAJgvL*5DJqoN8k@%} zgNq9C6T69Be99PkM?ACSfM{NC1x#8kL0Dh-mNwyf(`MWyF|)WBexcY}@-F|^Vfr-w zBi0r*Rqu7tuV3{wL$4;Zl&i%tM2mWV*7Vt(L0(b|->4NjnNMmA4{J5MU;{yt>rOH$|dQ)ejtf^H1hB|e4nTS+tKb_zY!%@TE?Er5jV{c>Jw zXB;Dnzo-}8Pj-VIk2x1ffi;~1kG=@;-2ntVk^hobyW z{H;5x!>lno3Q*y>Z+l=RohF}O+EJ8lFUJT>8>6oIzQ zzxGtG2B7cb>EgdD4#bj(;2uLvaGfPwO`w&#sFCyCwatR@p@5hb-)Tf=_*Hh9s;XdU zO!lQG#FNl9O!6Xh=~e;zj)Jy-V6%qyvKGl)d@v*Xl36uqQ&=nvD4BDx-x-PfsxCXd-pw zZu{{8%m#Dp^L=j1Y+p8AAGY|-*Ys0lMd|jj^o$Ja`H@YjUa0Op9-pR(PJz6Uf(FfV zaxb)x zq^?eTuCC80J=9K*wRnD{lqscW{klas$M&wAD}vIwpRksI9ZkngWhVVa;v~*Fq^~Ym z&THS6_GmdSCBAQB6Xbr;Inww%Fcz5m2#%&T|M*Q_ml^Dk3YjWnFvfJh#+4c z&}71c0lA~II8vIZNz)^t&!}t(X54I{>{agQMYlqE@qN*Oo=0gZr7)FI5j7Nd2M+gn z0X9yZ24U9?NK@Pg{iKvYH5rN=TN{~YIcUAVP=zr0!2!A|=!l%tK1SU;l zV{*M7AD;2NFtCA5IH@@>P*)fM=}b7JE2L$-5bNHd(wJGDoi7WNjM`osa${Tp5r({5 z)*f`f&P6phX9Bw(^XeUnRb!asV|cQlbdHomK`+w-+1$%Y_UbNST48Vlgu0xNaevOd zkmR)&AxVfa+nfc<0f&LlX*DspbL7Ye-PAC}^A!FHE!gH>YT?ya`jqXzj7Jtah{)>x zsUsfl><_K@o-57|EuP&Rw;yQy&YLmtKzpbuiaw6>wm#rjE$vjOh^|_8>7C^AUrV)K zMXxjOdK~NDqkQ*@1%{)2o{LS{j1We%vZk&VV zyORX|tOi99S#IOHR!q)KM#^x{5~J|}kxWS$W8?Ra(=if)`oUg57{kn4kv*Ob%{kR8 z^REw|;eo`OT`w1gH5KXHk0xivY0A(QPg{tFUbIiuSqZUJvN98>*?1u{Ut)`r47y2) z$eU9Xy`xX?@}|Zt)s9bJ=;wgLYVjV+YzK^X=k2p$8PKIZz86zCWFEMi8pXW<# z_g#u(5qQ#Icq73;QnF5jV#rNtpv^6Ls}71! zUXt*7-Qz2TnB!OJyus@7V)`~FJI?iaxz#V$M7daO&EbeiqgU0vQZxU1T<{9cvLA-m zbC~9W3dew%-Z#>H;rL$Ht9E~%IlDErz=q9dxl6o7J^rW&fe|{bWq9S8Yv%#Q*C|9U ztOY@Pze`BPu`SG%o@2NI%(w(=o_VvaP`Y{tT1d>Q#ez{>Dh|Z)w^SV^rAW=3R6K~o zP^!uo@zLQ;DY2=M|x3HHT1u2-lD0= zh1Y4CuFz?z-Qykd_tF@5J^Rw&=(G&q;yH9}CkalOrm4kcTL3+ue^&DwK|)Za59XR7 zLx7A=EO&Meq~a!7_*NnTyy~!?k^I_TY%2!-b*jqAmd3Jd{99+tgwrj>cHv=}qy9N2 zpWGNG$xy6CF4{X-ZqY-g>5wv2SX59GHKpI!E`FO`X|zR5yy~5>xGT5l70d4la{J}cJ}e9zjM>+u7M)#q}ZE)H%K8B7>pZ zPkC*w&c_7*9oe!(Q7N*l`@15{XdISc+u+=d;PH~=t*9e)CVY_pat9hUt60m64N^~f zVuJg|X7X7zkJ_q_vx~NV2zNk(BFML}K6$fEaLsg1MWxj7M!igJK7O9hF;{$nVp4`n zF$rL>>>Ky(WG;Lv`plYi-Ii46f99a6CH0h&u}awedp?2#rek+uoPQGsBd{-qyf5c_ z?u^Q9p$WQ3P2xc9;7HN%tB-886@~XGND9Nn!6KyGPv<+i_cqG9QI6_|$(_tf;lGM0 zO%Hdnc`dFbHPef9&zpxc)dh$dpPjle8^G9$Y!BY6r-GVF8QRjpH(O3c(NX!R^zi5T z>2(69&di^*AVSiWU0q2THt87@;tDsSyt8k`cv#{YPq|A4*bJw&P-o8SU?#qMlF^Lbpy!Pvkh%9t*&lTTUs(SMFRNEn zmglsv=qFCYSjFOui4@hAFV%?~b@$QQ=t+7be6cSOR|bQ&Nf>Pd6JglNB~onvmg!i0 z|D5}Y;H~)V$Zc0+$B~XD>#_I&3yDn~e|@z!jIDDT5a3zHNKd}7&z*nyPW@NNO0(EPK!Tj7@ zvT#K{+24)Nmyhw^6InP?Q7MgHObCD@70$~2C#j57Qu-5vD=Yn15e2#b#*k-)|I9@q zm1Y0Ll;!`uhw*=?Qih*t{8J1<=}$~f_N@OACnqZ}ubg%yX3Qv$lxKv&w6E(h{tFRK BgLwb| delta 4222 zcmZWqXE+>Ow@o1!L^pcxgc&oU6Nw&em=Pt~h!PScTGY{lsH2N9YDTZo+o;h(qDCi# zAc)@Mlkfhx_rCA>bM~{>+Uu-mKj-XIYA^(k(?WHd4wS9bQq64&AKXqpMJwM51Ktd} zai_{fT(*;RZF(|*O6-@&-nrn!`{vSkTRJkPQr9p9Z&4p|SMNBMx9oBO+g4_AMLS_} zRUINxa!TwN`RQy;Y*Uun>hjaVL`A;e5`gd$Jsu*uJMoGqfS4< zo68(yOzlr+C+nRXGCuj2VErvx{jY?4={Ux4$3NavJ`B-4I~sK!uuwDxDo#cK7C4V| ze$9MMr>2m}{=k@FQ{PyZRKQO_+NWakT?vko7%!-+r9QN0mQy6vJW|LGx1%WA)^7zkX{di&wwF)+N~c%v^wHiauI?b(x2U#d7{C`YD(M?YAvttss9?) zGARDL&jS5?uOW2Kj){-$Kw0~6cWzp9434|#=?X4zD1te^lyI0f+1<~K;D>dO<;+iM zju9u|PNm!b7jCx{{XYh%aKZ01yAleJ4zQgUybc8?FJWnKC=uEc8>OYH@-AnM1xIOxsoE&V@7LVfnbdN|1vU7ih&R`mZr`|D`-xJQ z0Yu_QeES}PDeyB)JAdu2$i{RobT71!pf${G@T67NL4e}-1e$yGIdeRAn`^H*qP}J& zyf9ite)&AyL4eBn9xswH#juBhEBU!qOHE?eb zmNtOIqeHimhZH-0boI#M0;~m)Ce2wWYP^1TpE|66(v1BE{Hkb|;8b_vSGo zqUk&(%=hZn&URbo%qvT;H=b9h*yH?4J#7aJXtfMDjc&?5DP#9r?swldQ|n)j!EbH1 zGDkVfDsAPFhB(RUZPHuoxGZ?Ge#3fj&`lp21EU%?Luyx?1X?%M`zzX`^#&K*m*32i z*H0OSISIJVkUtlkqx792KaoFLR+1dT&As8(Cv}~Sj%`TZHIo#-U@}9E`hST(%h&MG zOe*>~%<4`iMxOg_xX<>JN8YqM#8*`;w3^+Uh{@MKXt-^0{|*+rv_DWM=-YX5UZ?1P zR(33+nVKOE09Uf^xk~Zwt1sFclN@AWB zK)m?$X|y`iqXJc*mU?YX?Q!SSItQgB*U@@Pu*O6EcWTJ!vLXoY z@_wJNya-r&tp0f+#0VKc3bA3F@YKZJ2Q*}X8&4PJyM-YepfFAgucjXG{tW?cCP0USKY-4I=`T;(ACR`4w z^$n$o&B`~otit0}#}$+?ql^h81oNbA2o5^5olDU6DfK`IB2tL!x-rFZXXT{$LGNtZ$VC2_}-V*dPwz7KAOM zi~~ifmgNx{$#a?OHabD+KLD;%*0f3U*Y*u#I$7{blLmC!;?@A2A_9p4$L0YZ$Oc}+ zM%1=>ljpG@B@qOBS2!xr1_kvj_6Y2|b(||2j}t zv=5joM89>di!a$wG3@64)<8~_)3I|{wZudNm>pnP3&Q* z7f*l7xKpM7fz^U0V53BhRu7Egt8aQz5vy@KRUO0LHFQs9ilm=uDZG`3N5I(pQU8&a zrTpyVZ3%2v>0Oc@YRp*gJK4;2Uhj>_0m4=ETmkE2H&^g-4zG`lhra>tt0q#pn~sO| z$1O<{G5X|4`=H;C7^Cg{lPia=Q!qwshY4y0sx73m9hS;6~_i?3+^1A z0-rqPB(kvb8q(Z#MF1fl!wn;!p@U(!xT43ZTSDY->J#rInCev* z)k}Ze+QZu*1NUyK6+6an6y9blopqjVwWv%yj6kR1+T(4K>>OivEHg8rpGiC&RS1GH-9u z^Sd^8gd)Xxy5Y@~Oiz(ysgzC6JIFK$hWuK!A3)7R`FkkESjj~f4$5!Xg4zGD^Qq@) zHnjRqw!TB1f_<4|o%Qgs+JwG~+`d^^QC{?iw=?zMz}hs|t;mv?ioT)nCf*Dx`C$#I zl}T(REx3LwLM%gw%{d98mLA3qh3b~eJ3trzW%JJ_6AhgurH6|#!OSAc6sQaG~@^X!1Hkc{h5l|qJ z>xFx_oAQWJ-ldHC`FP*ya!tz7`;;TrFScA;Se=O-4YNGoCIw91f`Rl5EFH+tw>Cg( zn$wrtpSiYaaNwbwwBp`uO7Sr*;w-dS&0w`UVOI@aPYj2(Sti++sq)a}7H!MMn(_kf zB2b83{`!4;DC`RVwP=-Q*8L!IgN27ZDw~6{pZ& z@DmrlZ204hLTYyrs?z4TW8QuhkPu@>GowpIoFuJI32*T4!2v>Iih_f#sEQvHAe%I2 z`z!KFUtdOPfA^r(oE}<_LVTxhU`>{wN0q??oO$4dmWo!jzta^&5(VIe*K?Vs4DRf_ zSCCYAV2`N8z8Y2A67c+C9K!%F6s($$^Gk$$Y zz=RjM_PzHFj%d3YC+*)*d3JbHm3e@xdB-;lU-3Iz$&bWS;Sr_w)#9Vu_cKw8qISeW ztpmbsJ<3dFLX*Q%>%ZN{qy^5{EALv+%~bbd`+3l^h%)Re!_!fifTs)&Vt8uM?NMmZ z%(8HY5Y?j}Afw#$?5qpqDM@6*el4;vFU)wJk?b0{2N07tGr4m+v|#ES(YiMEP@Ofe82 znd;&*{?NIM&)*79BJ^p5_s2i_cwz*&5en*hyx1}f_vqfc4EIVDi>{DQTG>0HZ9-?7 z*ednUl9r#9hhn+H#XCl2=(jE<1k*#P*D*Y6%oVogL%Liyq_$Y=E_907U!f1b7}Rfi z>e!Wt2N=v76ST|gbYRNz8rHe zAY$&LC;aQR(D#C(M*RrHflFYnHY2WiCpZ-q*Sp+xPKC`o1NBH04|}?TmJ9D%F ziA$p>{_#V*7dTi<$0+~xjqqm5NBb~A{a|AQgnSX0FF`xEkNTbu*~&Sb1};pvD5Inw zNt{!}%{NgZpHnRix_r8Mo9rwzY;s9DM`*X@joZUNH&nPlW>{r{>n}kAAk;h{tfm

&0d@%ROuKX9>Na_#Ip(mfuJXPb-d7R1!orYER;-W>4Q|UBqy;0xxgB=5{%Kw)KZQ zF9}w2`B2QO=#*d!qno_%McNg9BrCE_5pWES^rn_%XZ3o^w}sOpq9dRdh4IYCayU2j zpZR%h6H(h&lZDT;Tq*t$SUmrMl0XW*2G}?~b9H-dWA*ZHd4EeVR +/** @file + @brief Contains a functions for parsing the command line, used by the applications. */ -using namespace NTL; +#include -/** @file - @brief Contains a ceiling function that was somehow missing from my NTL library. */ +using namespace std; + +//TODO: documentation + +char* get_cmd_option(char** begin, char** end, const string& option); -/** Funcion for ceiling NTL::RR values */ -void ceilPrec(RR& x, const RR& a, long p); +bool cmd_option_exists(char** begin, char** end, const string& option); #endif diff --git a/include/cleanbkz/version.hpp b/include/cleanbkz/version.hpp index dde0168..805dc13 100644 --- a/include/cleanbkz/version.hpp +++ b/include/cleanbkz/version.hpp @@ -2,11 +2,11 @@ #ifndef CBKZ_VERSION_H #define CBKZ_VERSION_H -#define CBKZ_VERSION "0.7.11" +#define CBKZ_VERSION "0.8.0" #define CBKZ_MAJOR_VERSION (0) -#define CBKZ_MINOR_VERSION (7) -#define CBKZ_REVISION (11) +#define CBKZ_MINOR_VERSION (8) +#define CBKZ_REVISION (0) #endif diff --git a/ratio.tmp b/ratio.tmp index ec739bd..39bc333 100644 --- a/ratio.tmp +++ b/ratio.tmp @@ -1 +1 @@ -0 0 1 1 3 6 11 18 29 39 64 87 143 182 284 361 543 682 1001 1247 1727 1867 2477 2449 2998 2880 3275 3049 3480 3033 3390 2982 3040 2470 2316 1697 1566 1128 915 629 535 384 297 190 145 93 71 38 23 15 11 6 3 1 0 0 0 0 0 \ No newline at end of file +1 1 1 1 1 1 2 4 5 3 6 9 14 17 21 22 40 58 101 143 160 161 167 153 155 132 99 68 53 38 40 32 29 23 28 25 31 24 37 51 64 57 75 80 72 55 45 33 21 15 8 5 4 3 1 1 1 1 1 1 \ No newline at end of file diff --git a/src/boundary.cpp b/src/boundary.cpp index 8c707cd..e48f7c3 100644 --- a/src/boundary.cpp +++ b/src/boundary.cpp @@ -25,7 +25,6 @@ #include #include #include -#include using namespace std; @@ -33,9 +32,9 @@ inline double fact(double x) { return (x < 2 ? 1 : x * fact(x - 1)); } -RR fact_RR(int x, int prec) { +RR fact_RR(int x) { RR ret; - ret.SetPrecision(prec); + ret.SetPrecision(RR_PRECISION); ret= 1; for(int i= 2; i <= x; i++) @@ -44,37 +43,21 @@ RR fact_RR(int x, int prec) { return ret; } -RR polytope_volume_RR(RR vols[], RR bounds[], int dim, int prec) { - RR tmp, pwr, ret; - - tmp.SetPrecision(prec); - pwr.SetPrecision(prec); - ret.SetPrecision(prec); - - ret= 1; - if(dim==0) - return ret; - - ret= 0; - vols[dim-1]= polytope_volume_RR(vols, bounds, dim-1, prec); - - for(int i= 0; i < dim; i++) { - tmp= dim-i; - pow(pwr, bounds[i], tmp); - ret+= pwr*vols[i]/fact_RR(dim-i,prec)*((dim-i+1)%2==0?1:-1); - } - - return ret; +double ball_vol(int k, double r) { + if (k%2==0) + return pow(M_PI, k/2) / fact(k/2) * pow(r, k); + else + return 2 * pow(4*M_PI, k/2) * fact(k/2) / fact(k) * pow(r, k); } RR RR_PI; -RR ball_vol_RR(int k, RR r, int prec) { +RR ball_vol_RR(int k, RR r) { RR pwr, exp, pipow; - pwr.SetPrecision(prec); - exp.SetPrecision(prec); - pipow.SetPrecision(prec); + pwr.SetPrecision(RR_PRECISION); + exp.SetPrecision(RR_PRECISION); + pipow.SetPrecision(RR_PRECISION); exp= k; pow(pwr, r, exp); @@ -82,183 +65,124 @@ RR ball_vol_RR(int k, RR r, int prec) { if (k%2==0) { exp= k/2; pow(pipow, RR_PI, exp); - return pipow / fact_RR(k/2, prec) * pwr; + return pipow / fact_RR(k/2) * pwr; } else { exp= k/2; pow(pipow, 4*RR_PI, exp); - return 2 * pipow * fact_RR(k/2, prec) / fact_RR(k, prec) * pwr; + return 2 * pipow * fact_RR(k/2) / fact_RR(k) * pwr; } } -double ball_vol(int k, double r) { - - /*cout << "#" << endl << "# ball_vol: R= " << r << " dim= " << k << endl; - cout << "# Gamma: " << pow(M_PI, k/2.0) / tgamma(k/2.0+1) * pow(r, k) << endl; - if (k%2==0) - cout << "# Form: " << pow(M_PI, k/2) / fact(k/2) * pow(r, k) << endl; - else - cout << "# Form: " << 2 * pow(4*M_PI, k/2) * fact(k/2) / fact(k) * pow(r, k) << endl << "#" << endl;*/ - - if (k%2==0) - return pow(M_PI, k/2) / fact(k/2) * pow(r, k); - else - return 2 * pow(4*M_PI, k/2) * fact(k/2) / fact(k) * pow(r, k); -} - -double integral_odd(int ltilde, int l, double tvec[], double vvec[]) { - double ret= 0; - - if(ltilde==0) - return 1; - - vvec[ltilde-1]= integral_odd(ltilde-1, l, tvec, vvec); +RR integral_even_RR(int h, int l, RR tvec[], RR vvec[]) { + RR tmp, pwr, ret; - ret-= pow(1-tvec[l-ltilde], (2*ltilde+1)/2.0)*pow(2,2*ltilde+1)*fact(ltilde+1)/fact(2*ltilde+2); - for(int i= 1; i < ltilde; i++) - ret+= pow(tvec[l-ltilde],ltilde-i)*vvec[i]/fact(ltilde-i)*((ltilde-i-1)%2==0?1:-1); + tmp.SetPrecision(RR_PRECISION); + pwr.SetPrecision(RR_PRECISION); + ret.SetPrecision(RR_PRECISION); - if(ltilde!=l) + ret= 1; + if(h==0) return ret; - return ret + pow(2,2*ltilde+1)*fact(ltilde+1)/fact(2*ltilde+2); -} - -double integral_even(int ltilde, int l, double tvec[], double vvec[]) { - double ret= 0; - - if(ltilde==0) - return 1; + ret= 0; - vvec[ltilde-1]= integral_even(ltilde-1, l, tvec, vvec); + vvec[h-1]= integral_even_RR(h-1, l, tvec, vvec); - for(int i= 0; i < ltilde; i++) - ret+= pow(tvec[l-ltilde],ltilde-i)*vvec[i]/fact(ltilde-i)*((ltilde-i-1)%2==0?1:-1); + for(int i= 0; i < h; i++) { + tmp= h-i; + pow(pwr, tvec[l-h], tmp); + ret+= pwr*vvec[i]/fact_RR(h-i)*((h-i-1)%2==0?1:-1); + } return ret; } -double polytope_volume(double vols[], double bounds[], int dim) { + +double integral_even(int h, int l, double tvec[], double vvec[]) { double ret= 0; - if(dim==0) + if(h==0) return 1; - vols[dim-1]= polytope_volume(vols, bounds, dim-1); + vvec[h-1]= integral_even(h-1, l, tvec, vvec); - for(int i= 0; i < dim; i++) - ret+= pow(bounds[i], dim-i)*vols[i]/fact(dim-i)*((dim-i+1)%2==0?1:-1); + for(int i= 0; i < h; i++) + ret+= pow(tvec[l-h],h-i)*vvec[i]/fact(h-i)*((h-i-1)%2==0?1:-1); return ret; } -RR t_full_RR(RR R, RR b_star_norm[], double t_node, double t_reduc, int n, int prec) { - RR N, V_act, denom; - N.SetPrecision(prec); - V_act.SetPrecision(prec); - denom.SetPrecision(prec); - denom= 1; - for(int i= 1; i<= n; i++) { - V_act= ball_vol_RR(i, R, prec); - - denom*= b_star_norm[n-i]; - N+= V_act/denom; - } - - N/= 2; - return t_reduc+t_node*N; - } +RR integral_odd_RR(int h, int l, RR tvec[], RR vvec[]) { + RR tmp, pwr, ret, pow2, two; -RR t_extreme_RR(RR Rvec[], RR b_star_norm[], double t_node, double t_reduc, int n, long prec) { - int pdim= (n%2==0)?n/2:n/2+1; - RR p_succ; - p_succ.SetPrecision(prec); - - RR* polytope_vols= new RR[pdim+1]; - RR* bounds= new RR[pdim]; + tmp.SetPrecision(RR_PRECISION); + pwr.SetPrecision(RR_PRECISION); + pow2.SetPrecision(RR_PRECISION); + ret.SetPrecision(RR_PRECISION); + two.SetPrecision(RR_PRECISION); - for(int i= 0; i< n; i+=2) { - bounds[i/2].SetPrecision(prec); - bounds[i/2]= Rvec[i]*Rvec[i]/(Rvec[n-1]*Rvec[n-1]); - } + two= 2; + ret= 1; + if(h==0) + return ret; - RR N, V_act, denom; - N.SetPrecision(prec); - V_act.SetPrecision(prec); - denom.SetPrecision(prec); - denom= 1; - polytope_vols[pdim]= polytope_volume_RR(polytope_vols, bounds, pdim, prec); - for(int i= 1; i<= n; i++) { - if(i%2==0) - V_act= ball_vol_RR(i, Rvec[i-1], prec) * polytope_vols[i/2] * fact_RR(i/2, prec); - else - V_act= ball_vol_RR(i, Rvec[i-1], prec) * (polytope_vols[i/2] * fact_RR(i/2,prec) + polytope_vols[i/2+1]*fact_RR(i/2+1,prec)) / 2; + ret= 0; - //cout << "# denom_" << i-1 << " = " << denom << endl; + vvec[h-1]= integral_odd_RR(h-1, l, tvec, vvec); - denom*= b_star_norm[n-i]; + tmp= (2*h+1)/2.0; + pow(pwr, 1-tvec[l-h], tmp); + tmp= 2*h+1; + pow(pow2, two, tmp); - /* cout << "# denom_" << i << " = " << denom << endl; - cout << "# bstarnorm_" << n-i << " = " << b_star_norm[n-i] << endl; - cout << "# V_" << i << " = " << V_act/denom << endl; - cout << "# Vball_" << i << "(" << Rvec[i-1] << ") = " << ball_vol_RR(i, Rvec[i-1], prec) << endl << endl;*/ + ret-= (pwr*pow2/fact_RR(2*h+2))*fact_RR(h+1); + for(int i= 1; i < h; i++) { + tmp= h-i; + pow(pwr, tvec[l-h], tmp); + ret+= pwr*vvec[i]/fact_RR(h-i)*((h-i-1)%2==0?1:-1); + } + + if(h!=l) + return ret; + return ret + pow2/fact_RR(2*h+2)*fact_RR(h+1); +} - N+= V_act/denom; - } - p_succ= polytope_vols[pdim-1]*fact_RR(pdim-1, prec); +double integral_odd(int h, int l, double tvec[], double vvec[]) { + double ret= 0; - delete [] bounds; - delete [] polytope_vols; + if(h==0) + return 1; - N/= 2; - return (t_reduc+t_node*N)/p_succ; -} + vvec[h-1]= integral_odd(h-1, l, tvec, vvec); -RR t_extreme(RR Rvec[], RR b_star_norm[], double t_node, double t_reduc, int n, long prec){ - return t_extreme_RR(Rvec, b_star_norm, t_node, t_reduc, n, prec); - } + ret-= pow(1-tvec[l-h], (2*h+1)/2.0)*pow(2,2*h+1)*fact(h+1)/fact(2*h+2); + for(int i= 1; i < h; i++) + ret+= pow(tvec[l-h],h-i)*vvec[i]/fact(h-i)*((h-i-1)%2==0?1:-1); -RR t_extreme_reference_RR(RR Rvec[], RR b_star_norm[], double t_node, double t_reduc, int n, int prec) { - int pdim= (n%2==0)?n/2:n/2+1; - RR p_succ; - p_succ.SetPrecision(prec); - - RR* polytope_vols= new RR[pdim+1]; - RR* bounds= new RR[pdim]; + if(h!=l) + return ret; - for(int i= 0; i< n; i+=2) { - bounds[i/2].SetPrecision(prec); - bounds[i/2]= Rvec[i]*Rvec[i]/(Rvec[n-1]*Rvec[n-1]); - } + return ret + pow(2,2*h+1)*fact(h+1)/fact(2*h+2); +} +RR t_full_RR(RR R, RR b_star_norm[], double t_node, double t_reduc, int n, int prec) { RR N, V_act, denom; N.SetPrecision(prec); V_act.SetPrecision(prec); denom.SetPrecision(prec); - polytope_vols[pdim]= polytope_volume_RR(polytope_vols, bounds, pdim, prec); + denom= 1; for(int i= 1; i<= n; i++) { - if(i%2==0) - V_act= ball_vol_RR(i, Rvec[i-1], prec) * polytope_vols[i/2] * fact_RR(i/2, prec); - else - V_act= ball_vol_RR(i, Rvec[i-1], prec) * (polytope_vols[i/2] * fact_RR(i/2,prec) + polytope_vols[i/2+1]*fact_RR(i/2+1,prec)) / 2; - - denom= 1; - for(int j= n-i; j < n; j++) - denom*= b_star_norm[j]; + V_act= ball_vol_RR(i, R); + denom*= b_star_norm[n-i]; N+= V_act/denom; } - p_succ= polytope_vols[pdim-1]*fact_RR(pdim-1, prec); - - delete [] bounds; - delete [] polytope_vols; - N/= 2; - return (t_reduc+t_node*N)/p_succ; -} - + return t_reduc+t_node*N; + } double n_full(double R, double b_star_norm[], int n) { double N, V_act, denom; @@ -279,43 +203,33 @@ double n_full(double R, double b_star_norm[], int n) { return N; } -inline double get_gsa(double scale, int dim, int bsize, int i) { - cout << "# get " << i << " exponent= " << (scale + gsa_slope[bsize]*(i-(dim+1)/2.0) + gsa_convexity[bsize]*(i*i-(dim+1)*i+(dim+1)*(dim+2)/6.0)) << endl; - return exp( (scale + gsa_slope[bsize]*(i-(dim+1)/2.0) + gsa_convexity[bsize]*(i*i-(dim+1)*i+(dim+1)*(dim+2)/6.0)) ); -} - -double n_full_gsa(double R, double scale, int bsize, int n) { - double N, V_act, denom; - - cout << "# " << "scale: " << scale << endl; - cout << "# " << "slope: " << gsa_slope[bsize] << endl; - cout << "# " << "convexity: " << gsa_convexity[bsize] << endl; - cout << "\"GSA\"" << endl << endl; +RR ci_prob_RR(RR Rvec[], int k) { + RR ret; + int l= k/2; + RR* bounds= new RR[l]; + RR* vvec= new RR[l]; - N= 0; - denom= 1; - for(int i= 1; i<= n; i++) { - V_act= ball_vol(i, R); + for(int i= 0; i< l; i++) + bounds[i]= Rvec[2*i]*Rvec[2*i]/(Rvec[k-1]*Rvec[k-1]); - denom*= sqrt(get_gsa(scale, n, bsize, n-i+1)); - N+= V_act/denom; + if(k%2==0) + ret= integral_even_RR(l, l, bounds, vvec)*fact_RR(l); + else + ret= integral_odd_RR(l, l, bounds, vvec)*fact_RR(2*l+2)/(power2_RR(2*l+1)*fact_RR(l+1)); - cout << "# " << n-i << " gsa: " << sqrt(get_gsa(scale, n, bsize, n-i+1)) << endl; - cout << i-1 << " " << V_act/denom << endl; - } - cout << endl << "# " << N/2 << endl; + delete [] bounds; + delete [] vvec; - N/= 2; - return N; + return ret; } double ci_prob(double Rvec[], int k) { + //TODO: ez a jó, a másikat is igy kijavitani double ret; int l= k/2; double* bounds= new double[l]; double* vvec= new double[l]; - /*cout << "# khalf: " << l << endl; cout << "# Rvec: [ "; for(int i= 0; i< 2*l; i++) @@ -341,12 +255,73 @@ double ci_prob(double Rvec[], int k) { return ret; } -// The case where n is even -double t_extreme_reference(double Rvec[], double b_star_norm[], double t_node, double t_reduc, int n) { +RR p_succ(RR Rvec[], int n){ + RR ret; + int l= (n%2==1)?n/2:n/2-1; + RR* bounds= new RR[l]; + RR* vvec= new RR[l]; + + for(int i= 0; i< l; i++) + bounds[i]= Rvec[2*i]*Rvec[2*i]/(Rvec[n-1]*Rvec[n-1]); + + if(n%2==0) + ret= integral_even_RR(l, l, bounds, vvec)*fact_RR(l); + else + ret= integral_odd_RR(l, l, bounds, vvec)*fact_RR(2*l+2)/(power2_RR(2*l+1)*fact_RR(l+1)); + + delete [] bounds; + delete [] vvec; + + return ret; +} + +RR t_extreme_RR(RR Rvec[], RR b_star_norm[], double t_node, double t_reduc, int n) { + RR N; + RR V_act; + RR denom; + + N= 0; + for(int k= 1; k<= n; k++) { + denom= 1; + + for(int i= n-k; i < n; i++) + denom*= b_star_norm[i]; + + V_act= ball_vol_RR(k, Rvec[k-1]) * ci_prob_RR(Rvec, k); + + N+= V_act/denom; + } + + + N/= 2; + return (t_reduc+t_node*N)/p_succ(Rvec, n); +} + +/*double t_extreme(double Rvec[], double b_star_norm[], double t_node, double t_reduc, int n) { + double N= 0; + double V_act; + double denom; + + for(int k= 1; k<= n; k++) { + denom= 1; + + for(int i= n-k; i < n; i++) + denom*= b_star_norm[i]; + + V_act= ball_vol(k, Rvec[k-1]) * ci_prob(Rvec, k); + + N+= V_act/denom; + } + + + N/= 2; + return (t_reduc+t_node*N)/ci_prob(Rvec, n/2); +}*/ + +void predict_nodes(double Rvec[], double b_star_norm[], int n) { double N= 0; double V_act; double denom; - //double N_next, N_prev; cout << "# T_extreme: " << endl << "# GS-lengths: [ "; for(int i= 0; i < n; i++) @@ -362,35 +337,8 @@ double t_extreme_reference(double Rvec[], double b_star_norm[], double t_node, d unsigned long nodes; unsigned long sum= 0; -// N_next= 0; -// N_prev= 1/b_star_norm[n-1]; for(int k= 1; k<= n; k++) { - /*if(k%2==0) { - //V_act= ball_vol(k, Rvec[k-1]) * pvol_and_scale(Rvec, k/2) * fact(k/2); - N_prev= N_next; - N+= N_prev; - cout << k-1 << " " << N_prev << endl; - } else { - denom= 1; - for(int i= n-k-1; i < n; i++) - denom*= b_star_norm[i]; - - N_next= ball_vol(k+1, Rvec[k]) * pvol_and_scale(Rvec, k/2+1) * fact(k/2+1) / denom; - N+= (N_prev+N_next)/2; - cout << k-1 << " " << (N_prev+N_next)/2 << endl; - - //V_act= ball_vol(k-1, Rvec[k-2]) * pvol_and_scale(Rvec, k/2)*fact(k/2) * 2 * Rvec[k]; // upper bound: V_{k-1}*2*R_k - //V_act= ball_vol(k, Rvec[k-1]) * ( pvol_and_scale(Rvec, k/2)*fact(k/2) + pvol_and_scale(Rvec, k/2+1)*fact(k/2+1)) / 2; - - cout << "# Simvol_" << k/2+1 << " = " << fact(k/2+1) << endl; - cout << "# Polvol_" << k/2+1 << " = " << pvol_and_scale(Rvec, k/2+1) << endl; - cout << "# V_ball(" << k+1 << ", " << Rvec[k] << ") = " << ball_vol(k+1, Rvec[k]) << endl; - cout << "# |b^*_" << n-k-1 << "| = " << b_star_norm[n-k-1] << endl; - cout << "# V_" << k+1 << " = " << ball_vol(k+1, Rvec[k]) * pvol_and_scale(Rvec, k/2+1) * fact(k/2+1) << endl; - cout << "# denom_" << k << " = " << denom << endl; -// }*/ denom= 1; - //for(int i= 0; i < k; i++) for(int i= n-k; i < n; i++) denom*= b_star_norm[i]; @@ -398,103 +346,74 @@ double t_extreme_reference(double Rvec[], double b_star_norm[], double t_node, d N+= V_act/denom/2; - //cout << "# Simvol_" << psd << " = " << fact(psd) << endl; - //cout << "# Polvol_" << psd << " = " << pvol_and_scale(Rvec, k) << endl; - //cout << "# Polprob_" << k/2 << " = " << polprob_even(Rvec, k) << endl; - cout << "# V_ball(" << k << ", " << Rvec[k-1] << ") = " << ball_vol(k, Rvec[k-1]) << endl; - //cout << "# |b^*_" << n-k << "| = " << b_star_norm[n-k] << endl; - //cout << "# V_" << k << " = " << V_act << endl; - cout << "# denom_" << k << " = " << denom << endl; - cout << "# ratio_" << k << " = " << pow(ball_vol(k, Rvec[k-1])/denom,1.0/k) << endl; - //cout << "# N_" << k << " = " << V_act/denom/2 << endl; + cout << "# N_" << k << " = " << V_act/denom/2 << endl; myfile >> nodes; sum+= nodes; cout << "# Measured nodes: " << nodes << endl; - //cout << "# Difference: " << nodes - V_act/denom/2 << endl; - //cout << "# Ratio: " << nodes/(V_act/denom/2) << endl; - //cout << "# GH_" << k << " = " << pow(denom/ball_vol(k, 1),1.0/k) << endl; - cout << "# R_" << k << " = " << Rvec[k-1] << endl; - if( Rvec[k-1] < sqrt(k)) - cout << "# WARNING! Gaussian Heuristics violation: " << k << "^(1/2) = " << sqrt(k) << endl; - //cout << "# R/GH in dim " << k << " = " << Rvec[k-1]/ pow(denom/ball_vol(k, 1),1.0/k) << endl; + cout << "# Difference: " << nodes - V_act/denom/2 << endl; + cout << "# Ratio: " << nodes/(V_act/denom/2) << endl; + cout << "# GH_" << k << " = " << pow(denom/ball_vol(k, 1),1.0/k) << endl; + cout << "# R/GH in dim " << k << " = " << Rvec[k-1]/ pow(denom/ball_vol(k, 1),1.0/k) << endl; - //cout << k-1 << " " << nodes/(V_act/denom/2) << endl; - //cout << k-1 << " " << nodes - V_act/denom/2 << endl; cout << k-1 << " " << V_act/denom/2 << endl; } myfile.close(); - //N/= 2; cout << "# Predicted nodes: " << N << endl; cout << "# Measured nodes: " << sum << endl; - return (t_reduc+t_node*N)/ci_prob(Rvec, n/2); } -double t_extreme(double Rvec[], double b_star_norm[], double t_node, double t_reduc, int n) { - int pdim= (n%2==0)?n/2:n/2+1; - - double* polytope_vols= new double[pdim+1]; - double* bounds= new double[pdim]; - - for(int i= 0; i< n; i+=2) - bounds[i/2]= Rvec[i]*Rvec[i]/(Rvec[n-1]*Rvec[n-1]); +void predict_nodes_RR(RR Rvec[], double b_star_norm[], int n) { + RR N; + N= 0; + RR V_act; + RR denom; - double N= 0; - double V_act, denom= 1; - polytope_vols[pdim]= polytope_volume(polytope_vols, bounds, pdim); - for(int i= 1; i<= n; i++) { - if(polytope_vols[i/2] * fact(i/2) > 1) - cout << " # " << polytope_vols[i/2] * fact(i/2) << endl; + RR_PI.SetPrecision(RR_PRECISION); + RR_PI= ComputePi_RR(); - if(i%2==0) - V_act= ball_vol(i, Rvec[i-1]) * polytope_vols[i/2] * fact(i/2); - else - V_act= ball_vol(i, Rvec[i-1]) * (polytope_vols[i/2] + polytope_vols[i/2+1])*fact(i/2+1) / 2; - //V_act= ball_vol(i, Rvec[i-1]) * (polytope_vols[i/2] + polytope_vols[i/2+1])*fact(i/2+1) / 2; + cout << "# T_extreme: " << endl << "# GS-lengths: [ "; + for(int i= 0; i < n; i++) + cout << b_star_norm[i] << " "; + cout << "]" << endl << "# Bounds: [ "; + for(int i= 0; i < n; i++) + cout << Rvec[i] << " "; + cout << "]" << endl << "#" << endl; - denom*= b_star_norm[n-i]; - N+= V_act/denom; - /*cout << "Simvol_" << i/2 << " = " << fact(i/2) << endl; - cout << "Polvol_" << i/2 << " = " << polytope_vols[i/2] << endl; - cout << "V_ball(" << i << ", " << Rvec[i-1] << ") = " << ball_vol(i, Rvec[i-1])<< endl; - cout << "|b^*_" << n-i << "| = " << b_star_norm[n-i] << endl; - cout << "V_" << i << " = " << V_act << endl; - cout << "denom_" << i << " = " << denom << endl; - cout << "N_" << i << " = " << V_act/denom << endl << endl;*/ - } - - N/= 2; + cout << "\"Predicted\"" << endl << endl; - cout << "Predicted nodes: " << N << endl; - return (t_reduc+t_node*N)/(polytope_vols[pdim-1]*fact(pdim-1)); -} + ifstream myfile ("ratio.tmp"); + unsigned long nodes; + unsigned long sum= 0; -RR compute_p_succ(RR Rvec[], int n, int prec){ - int pdim= (n%2==0)?n/2:n/2+1; - RR p_succ; - p_succ.SetPrecision(prec); - - RR* polytope_vols= new RR[pdim+1]; - RR* bounds= new RR[pdim]; + for(int k= 1; k<= n; k++) { + denom= 1; + for(int i= n-k; i < n; i++) + denom*= b_star_norm[i]; - for(int i= 0; i< n; i+=2) { - bounds[i/2].SetPrecision(prec); - bounds[i/2]= Rvec[i]*Rvec[i]/(Rvec[n-1]*Rvec[n-1]); - } + V_act= ball_vol_RR(k, Rvec[k-1]) * ci_prob_RR(Rvec, k); - polytope_vols[pdim]= polytope_volume_RR(polytope_vols, bounds, pdim, prec); + N+= V_act/denom/2; - p_succ= polytope_vols[pdim-1]*fact_RR(pdim-1, prec); + cout << "# Ball_" << k << " = " << ball_vol_RR(k, Rvec[k-1]) << endl; + cout << "# P_" << k << " = " << ci_prob_RR(Rvec, k) << endl; + cout << "# N_" << k << " = " << V_act/denom/2 << endl; + myfile >> nodes; sum+= nodes; + cout << "# Measured nodes: " << nodes << endl; + cout << "# Difference: " << nodes - V_act/denom/2 << endl; + //cout << "# Ratio: " << nodes/(V_act/denom/2) << endl; + //cout << "# GH_" << k << " = " << pow(denom/ball_vol(k, 1),1.0/k) << endl; + //cout << "# R/GH in dim " << k << " = " << Rvec[k-1]/ pow(denom/ball_vol(k, 1),1.0/k) << endl; - delete [] bounds; - delete [] polytope_vols; + cout << k-1 << " " << V_act/denom/2 << endl; + } + + myfile.close(); - return p_succ; + cout << "# Predicted nodes: " << N << endl; + cout << "# Measured nodes: " << sum << endl; } -RR p_succ(RR Rvec[], int n, int prec){ - return compute_p_succ(Rvec, n, prec); - } static void generate_boundary_step(RR b_star_norm[], double t_node, double t_reduc, int n, RR* act, double delta, unsigned long iterations, RR& t_enum, int& changes) { @@ -540,7 +459,7 @@ void generate_boundary_step(RR b_star_norm[], double t_node, double t_reduc, int break; - time= t_extreme_RR(mod, b_star_norm, t_node, t_reduc, n, RR_PRECISION); + time= t_extreme_RR(mod, b_star_norm, t_node, t_reduc, n); if(time < t_enum) { t_enum= time; @@ -558,7 +477,7 @@ void generate_boundary_step(RR b_star_norm[], double t_node, double t_reduc, int delete [] mod; } -void generate_boundary(RR b_star_norm[], double t_node, double t_reduc, int n, double Rvec[], double R, double delta, unsigned long iterations, double& p_succ, double& t_enum_d, bool quiet) { +void generate_boundary(RR b_star_norm[], double t_node, double t_reduc, int n, double Rvec[], double R, double delta, unsigned long iterations, double& p_succ_v, double& t_enum_d, bool quiet) { int changes; RR* act= new RR[n]; @@ -577,7 +496,7 @@ void generate_boundary(RR b_star_norm[], double t_node, double t_reduc, int n, d } act[n-1].SetPrecision(RR_PRECISION); act[n-1]= R; - t_enum= t_extreme_RR(act, b_star_norm, t_node, t_reduc, n, RR_PRECISION); + t_enum= t_extreme_RR(act, b_star_norm, t_node, t_reduc, n); if(!quiet) { RR tmp; @@ -596,7 +515,7 @@ void generate_boundary(RR b_star_norm[], double t_node, double t_reduc, int n, d if(!quiet) cout << "# Changes: " << changes << endl; - conv(p_succ, compute_p_succ(act, n, RR_PRECISION)); + conv(p_succ_v, p_succ(act, n)); conv(t_enum_d, t_enum); delete [] act; diff --git a/src/cjloss.cpp b/src/cjloss.cpp index 66c748f..b594a19 100644 --- a/src/cjloss.cpp +++ b/src/cjloss.cpp @@ -22,38 +22,60 @@ #include #include #include +#include using namespace std; cjloss::cjloss(long dimension, double density, long seed) { - basis.SetDims(dimension, dimension); - values.SetLength(dimension-1); - solution.SetLength(dimension-1); + values.SetLength(dimension); + solution.SetLength(dimension); // computing the maximal value to reach the prescribed density ZZ max; - conv(max,floor(pow(2, (dimension-1)/density))); + conv(max,floor(pow(2, dimension/density))); ZZ s; conv(s,seed); SetSeed(s); - randomize(dimension-1,max); + randomize(dimension,max); while(!check()) - randomize(dimension-1,max); + randomize(dimension,max); - for(int i= 0; i #include #include +#include using namespace std; @@ -161,7 +162,7 @@ void enumerate_ntl(mat_ZZ& basis, int beta, double* prune, vec_RR& result) { enumerate_ntl(mu, c, prune, 0, mu1.NumRows()-1, mu1.NumRows(), result); } -void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result, unsigned long &termination, double &time) { +void profile_enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result) { cout << "# Enumerate: " << endl << "# GS lengths: [ "; for(int i= 0; i < n; i++) @@ -212,14 +213,13 @@ void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result, unsigned long nodes= 0; unsigned long* prof= new unsigned long[n]; for(int i= 0; i < n; i++) - prof[i]= 0; + prof[i]= 1; int k= 0; - clock_t end,begin= clock(); cout << "\"Measured\"" << endl << endl; - while((termination==0)||(nodes < termination)) { + while(true) { rhovec[k]= rhovec[k+1]+(vvec[k]-cvec[k])*(vvec[k]-cvec[k])*b[k]; //cout << nodes << "\tk: " << k << "\trhovec: " << rhovec[k] << "\tRvec (n-k-1): " << Rvec[n-k-1] << endl; @@ -276,11 +276,6 @@ void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result, cout << "# Nodes: " << nodes << endl; - end= clock(); - time= (double) (end-begin) / CLOCKS_PER_SEC / nodes; - //if(termination!=0) - termination-= nodes; - ofstream myfile; myfile.open ("ratio.tmp",ios::trunc); for(int i= 0; i < n; i++) @@ -290,7 +285,7 @@ void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result, for(int i= 0; i < n; i++) cout << i << " " << prof[i] << endl; - + cout << endl << endl; if(k==0) { result.SetLength(n); @@ -311,6 +306,124 @@ void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result, delete [] rhovec; } +void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result, unsigned long &termination, double &time) { + bool pruning= true; + if(Rvec==NULL) { + pruning= false; + Rvec= new double[n]; + for(int i= 0; i < n; i++) + Rvec[i]= b[0]; + } + + double** sigmamat= new double*[n+1]; + for(int i= 0; i < n+1; i++) + sigmamat[i]= new double[n]; + for(int i= 0; i < n+1; i++) + for(int j= 0; j < n; j++) + sigmamat[i][j]= 0; + + int* rvec= new int[n+1]; + for(int i= 0; i < n+1; i++) + rvec[i]= i-1; + + double* rhovec= new double[n+1]; + for(int i= 0; i < n+1; i++) + rhovec[i]= 0; + + int* vvec= new int[n]; + vvec[0]= 1; + for(int i= 1; i < n; i++) + vvec[i]= 0; + + double* cvec= new double[n]; + for(int i= 0; i < n; i++) + cvec[i]= 0; + + int* wvec= new int[n]; + for(int i= 0; i < n; i++) + wvec[i]= 0; + + int last_nonzero= 0; + + unsigned long nodes= 0; + unsigned long* prof= new unsigned long[n]; + for(int i= 0; i < n; i++) + prof[i]= 0; + + int k= 0; + clock_t end,begin= clock(); + + // C++ uses processor (float) precision even when comparing doubles + double epsilon= std::numeric_limits::epsilon(); + + while((termination==0)||(nodes < termination)) { + rhovec[k]= rhovec[k+1]+(vvec[k]-cvec[k])*(vvec[k]-cvec[k])*b[k]; + + if(rhovec[k] < Rvec[n-k-1]+epsilon) { + nodes++; + + if(k==0) + break; + else { + k--; + rvec[k]= rvec[k]>rvec[k+1]?rvec[k]:rvec[k+1]; + for(int i= rvec[k+1]; i>k; i--) + sigmamat[i][k]= sigmamat[i+1][k] + vvec[i]*mu[i][k]; + cvec[k]= -sigmamat[k+1][k]; + + vvec[k]= lround(cvec[k]); + wvec[k]= 1; + + } + } else { + k++; + if(k==n) + break; + rvec[k]= k; + + if(k>=last_nonzero){ + last_nonzero= k; + vvec[k]++; + } else { + if(vvec[k]>cvec[k]) + vvec[k]-= wvec[k]; + else + vvec[k]+= wvec[k]; + wvec[k]++; + } + } + } + + end= clock(); + time= (double) (end-begin) / CLOCKS_PER_SEC / nodes; + + termination-= nodes; + + if(k==0) { + result.SetLength(n); + for(int i= 0; i < result.length(); i++) + conv(result[i],vvec[i]); + } else + result.SetLength(0); + + if(!pruning) + delete [] Rvec; + for(int i= 0; i < n+1; i++) + delete [] sigmamat[i]; + delete [] sigmamat; + delete [] rvec; + delete [] vvec; + delete [] wvec; + delete [] cvec; + delete [] rhovec; +} + +void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result) { + double time; + unsigned long termination= 0; + return enumerate_epr(mu, b, Rvec, n, result, termination, time); +} + void enumerate_epr(mat_ZZ& basis, int beta, double* prune, vec_RR& result) { double time; @@ -347,5 +460,5 @@ void enumerate_epr(mat_ZZ& basis, int beta, double* prune, vec_RR& result) { unsigned long termination= 0; enumerate_epr(mu, c, prune, mu1.NumRows(), result, termination, time); - cout << "# Nodes processed: " << -termination << endl; +// cout << "# Nodes processed: " << -termination << endl; } diff --git a/src/tools.cpp b/src/tools.cpp index aeb53ac..ae67e52 100644 --- a/src/tools.cpp +++ b/src/tools.cpp @@ -38,23 +38,3 @@ bool cmd_option_exists(char** begin, char** end, const string& option) { return find(begin, end, option) != end; } -void ceilPrec(RR& x, const RR& a, long p) { - if (p < 1 || NTL_OVERFLOW(p, 1, 0)) - Error("CeilPrec: bad precision"); - - ConvPrec(x,a,p); - - cout << "x = " << x << endl; - - RR epsilon; - epsilon= 10; - RR exp; - exp= (x.exponent()+15); - cout << "exponent = " << x.exponent() << endl; - cout << "precision = " << x.precision() << endl; - epsilon= pow(epsilon, exp); - cout << "epsilon = " << epsilon << endl; - - x+= epsilon; -} - diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 0951312..37c7ae5 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,7 +1,7 @@ SET(TARGET_NAME1 unit_tests) SET(TARGET_NAME2 act_test) -SET(TARGET_NAME3 randtest) -SET(TARGET_NAME4 detprofile) +#SET(TARGET_NAME3 randtest) +#SET(TARGET_NAME4 detprofile) INCLUDE_DIRECTORIES(${PROJECT_SOURCE_DIR}/include) LINK_DIRECTORIES(${LIBRARY_OUTPUT_PATH}) @@ -11,8 +11,8 @@ TARGET_LINK_LIBRARIES(${TARGET_NAME1} cleanbkz ntl) ADD_EXECUTABLE(${TARGET_NAME2} act_test.cpp) TARGET_LINK_LIBRARIES(${TARGET_NAME2} cleanbkz ntl) -ADD_EXECUTABLE(${TARGET_NAME3} randtest.cpp) -TARGET_LINK_LIBRARIES(${TARGET_NAME3} cleanbkz ntl) +#ADD_EXECUTABLE(${TARGET_NAME3} randtest.cpp) +#TARGET_LINK_LIBRARIES(${TARGET_NAME3} cleanbkz ntl) -ADD_EXECUTABLE(${TARGET_NAME4} detprofile.cpp) -TARGET_LINK_LIBRARIES(${TARGET_NAME4} cleanbkz ntl) +#ADD_EXECUTABLE(${TARGET_NAME4} detprofile.cpp) +#TARGET_LINK_LIBRARIES(${TARGET_NAME4} cleanbkz ntl) diff --git a/test/act_test.cpp b/test/act_test.cpp index 82d241a..070156c 100644 --- a/test/act_test.cpp +++ b/test/act_test.cpp @@ -31,51 +31,49 @@ using namespace std; -double t_extreme_reference(double Rvec[], double b_star_norm[], double t_node, double t_reduc, int n); -double t_extreme(double Rvec[], double b_star_norm[], double t_node, double t_reduc, int n); -double n_full(double R, double b_star_norm[], int n); -double n_full_gsa(double R, double scale, int bsize, int n); -double ball_vol(int k, double r); +extern void predict_nodes_RR(RR Rvec[], double b_star_norm[], int n); + +extern void predict_nodes(double Rvec[], double b_star_norm[], int n); +extern double ball_vol(int k, double r); +extern void profile_enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result); void gen_randlat(mat_ZZ& basis, ZZ determinant, int dim); -ZZ autoscale(mat_ZZ& basis, double Rvec[]); int main(int argc, char** argv) { int dim, bsize= 20; stringstream ss; + bool l_cjloss= true; + mat_ZZ basis; ss << argv[1]; ss >> dim; ss.clear(); - -// cjloss l(dim, 0.94, 10); - cjloss l(dim, 0.94, time(NULL)); - BKZ_QP1(l.basis, 0.99, bsize); - - l.basis.SetDims(dim-1,dim); - -// cout << l.basis << endl; -// return 0; - - dim--; - -// ZZ determinant; -// GenPrime(determinant,dim); -// GenPrime(determinant,floor(2*dim)); - mat_ZZ basis; -// gen_randlat(basis,determinant,dim); - - basis= l.basis; + if (l_cjloss) { + cjloss l(dim, 0.94, time(NULL)); + basis= l.get_basis(bsize); + } else { + ZZ determinant; + GenPrime(determinant,dim); + gen_randlat(basis,determinant,dim); + BKZ_QP1(basis, 0.99, bsize); + } mat_RR mu1; vec_RR c1; - BKZ_QP1(basis, 0.99, bsize); ComputeGS(basis,mu1,c1); double* c= new double[mu1.NumRows()]; for(int i= 0; i < mu1.NumRows(); i++) - conv(c[i], sqrt(c1[i])); + conv(c[i], c1[i]); + + + double** mu= new double*[mu1.NumRows()]; + for(int i= 0; i < mu1.NumRows(); i++) + mu[i]= new double[mu1.NumCols()]; + for(int i= 0; i < mu1.NumRows(); i++) + for(int j= 0; j < mu1.NumCols(); j++) + conv(mu[i][j], mu1[i][j]); // Gaussian heuristic double gh= 1; @@ -85,13 +83,16 @@ int main(int argc, char** argv) { gsghs[i]= pow(gh/ball_vol(i+1, 1),1.0/(i+1)); } gh= pow(gh/ball_vol(dim, 1),1.0/dim); - //cout << "# Shortest vector length: " << sqrt(dim-1) << endl; - //cout << "# Gaussian heuristic: " << gh << endl; - //cout << "# GH/lambda Ratio: " << gh/sqrt(dim-1) << endl; + if(l_cjloss) + cout << "# Shortest vector length: " << sqrt(dim) << endl; + cout << "# Gaussian heuristic: " << gh << endl; // Bounding functions - //double R= dim-1; - double R= gh*gh; + double R; + if(l_cjloss) + R= dim; + else + R= gh*gh; double* full= new double[dim]; for(int i= 0; i < dim; i++) @@ -115,7 +116,7 @@ int main(int argc, char** argv) { double* step= new double[dim]; for(int i= 0; i < dim/2; i++) - if(iscale) - conv(scale, sqrt((i+1)/Rvec[i])); - - // conversion computes floor and ceil is needed - scale+= 1; - if(scale == 1) - return scale;*/ - - power2(scale, basis.NumRows()); - - for(int i= 0; i < basis.NumRows(); i++) - for(int j= 0; j < basis.NumCols(); j++) - basis[i][j]*= scale; - - double rscale; - conv(rscale, scale); - rscale*= rscale; - for(int i= 0; i < basis.NumRows(); i++) - Rvec[i]*= rscale; - - return scale; } void gen_randlat(mat_ZZ& basis, ZZ determinant, int dim) { diff --git a/test/unit_tests.cpp b/test/unit_tests.cpp index 018349e..c8284ea 100644 --- a/test/unit_tests.cpp +++ b/test/unit_tests.cpp @@ -39,32 +39,31 @@ inline double fact(double x) { extern RR polytope_volume_RR(RR vols[], RR bounds[], int dim, int prec); -extern RR fact_RR(int x, int prec); +extern RR fact_RR(int x); -extern RR ball_vol_RR(int k, RR r, int prec); +extern RR ball_vol_RR(int k, RR r); extern RR RR_PI; extern double ball_vol(int k, double r); extern double integral_even(int ltilde, int l, double tvec[], double vvec[]); +extern RR integral_even_RR(int h, int l, RR tvec[], RR vvec[]); + +extern double integral_odd(int ltilde, int l, double tvec[], double vvec[]); +extern RR integral_odd_RR(int h, int l, RR tvec[], RR vvec[]); void genbounds(int l, double* tvec); int main(int argc, char** argv) { - double* vols_d= new double[40]; + /* Unit tests for the polytope volume computation. Reference values are simplices or computed with vinci. */ + + double* vols_d= new double[70]; bool vinci_test= true; double vinci_vol_d; double* bounds_d= new double[40]; - /*srand ( 0 ); - for(int i= 0; i< 1000; i++) { - genbounds( 30, bounds_d); - if (abs(integral_even(30,30,bounds_d,vols_d)- polytope_volume(vols_d, bounds_d, 30)) > 1e-100) - cout << integral_even(30,30,bounds_d,vols_d)- polytope_volume(vols_d, bounds_d, 30) << endl; - }*/ - bounds_d[0]= 0.1; bounds_d[1]= 0.2; bounds_d[2]= 0.3; bounds_d[3]= 0.4; bounds_d[4]= 0.5; vinci_vol_d= 1.08e-04; if (abs(integral_even(5,5,bounds_d,vols_d)-vinci_vol_d) > 1e-19) { @@ -103,54 +102,13 @@ int main(int argc, char** argv) { cout << vinci_vol_d << endl; } + if(vinci_test) + cout << "\tVinci test PASSED." << endl; + RR* vols= new RR[71]; int prec= RR_PRECISION; - /* Unit tests for the polytope volume computation. Reference values are simplices or computed with vinci. */ - - cout << "Testing polytope volume computation:" << endl; - - vinci_test= true; - RR vinci_vol; - vinci_vol.SetPrecision(prec); - RR* bounds= new RR[15]; - - bounds[0]= 0.1; bounds[1]= 0.2; bounds[2]= 0.3; bounds[3]= 0.4; bounds[4]= 0.5; - vinci_vol= 1.08e-04; - if (abs(polytope_volume_RR(vols, bounds, 5, prec)-vinci_vol) > 1e-19) { - vinci_test= false; - cout << "\tVinci test FAILED in dimension 5" << endl; - } - - bounds[0]= 0.1; bounds[1]= 0.2; bounds[2]= 0.3; bounds[3]= 0.4; bounds[4]= 0.5; bounds[5]= 0.6; - vinci_vol= 2.334305555556e-05; - if (abs(polytope_volume_RR(vols, bounds, 6, prec)-vinci_vol) > 1e-17) { - vinci_test= false; - cout << "\tVinci test FAILED in dimension 6" << endl; - } - - bounds[0]= 0.1; bounds[1]= 0.2; bounds[2]= 0.3; bounds[3]= 0.4; bounds[4]= 0.5; bounds[5]= 0.6; - bounds[6]= 0.7; bounds[7]= 0.8; bounds[8]= 0.9; - vinci_vol= 2.755731922399e-07; - if (abs(polytope_volume_RR(vols, bounds, 9, prec)-vinci_vol) > 1e-19) { - vinci_test= false; - cout << "\tVinci test FAILED in dimension 9" << endl; - } - - bounds[0]= 0.1; bounds[1]= 0.2; bounds[2]= 0.3; bounds[3]= 0.4; bounds[4]= 0.5; bounds[5]= 0.6; - bounds[6]= 0.7; bounds[7]= 0.8; bounds[8]= 0.9; bounds[9]= 0.95; - vinci_vol= 5.120005762235e-08; - if (abs(polytope_volume_RR(vols, bounds, 10, prec)-vinci_vol) > 1e-20) { - vinci_test= false; - cout << "\tVinci test FAILED in dimension 10" << endl; - } - - //TODO: dimension 15 multiple - - if(vinci_test) - cout << "\tVinci test PASSED." << endl; - RR epsilon, x, y; epsilon.SetPrecision(prec); x= 2; @@ -166,8 +124,8 @@ int main(int argc, char** argv) { RR pv, sv; bool simplex_test= true; for(int i= 1; i < 71; i++) { - pv= polytope_volume_RR(vols, simplex, i, prec); - sv= 1/fact_RR(i, prec); + pv= integral_even_RR(i, i, simplex, vols); + sv= 1/fact_RR(i); y= pv.exponent()+95; pow(epsilon, x, y); // Checking if relative error is greater than prec-95 binary digits @@ -185,6 +143,55 @@ int main(int argc, char** argv) { if(simplex_test) cout << "\tSimplex test PASSED." << endl; + RR* linear= new RR[70]; + double* linear_d= new double[70]; + for(int i= 0; i < 70; i++) + linear[i]= linear_d[i]= (i+1)/70.0; + + RR rv, dv; + bool double_test= true; + for(int i= 1; i < 40; i++) { + rv= integral_even_RR(i, i, linear, vols); + dv= integral_even(i, i, linear_d, vols_d); + y= rv.exponent()+175; + pow(epsilon, x, y); + // Checking if relative error is greater than prec-175 binary digits + if( abs(rv - dv) > epsilon) { + double_test= false; + cout << "\tDouble test (even) failed in dimension " << i << endl; + } + /*cout << "Dimension " << i << endl; + cout << "\t RR: " << rv << endl; + cout << "\t double: " << dv << endl; + cout << "\t diff: " << abs(rv-dv) << endl; + cout << "\t epsilon: " << epsilon << endl << endl;*/ + } + + if(double_test) + cout << "\tDouble test (even) PASSED." << endl; + + double_test= true; + for(int i= 1; i < 40; i++) { + rv= integral_odd_RR(i, i, linear, vols); + dv= integral_odd(i, i, linear_d, vols_d); + y= rv.exponent()+155; + pow(epsilon, x, y); + // Checking if relative error is greater than prec-155 binary digits + if( abs(rv - dv) > epsilon) { + double_test= false; + cout << "\tDouble test (odd) failed in dimension " << i << endl; + } + /*cout << "Dimension " << i << endl; + cout << "\t RR: " << rv << endl; + cout << "\t double: " << dv << endl; + cout << "\t diff: " << abs(rv-dv) << endl; + cout << "\t epsilon: " << epsilon << endl << endl;*/ + } + + if(double_test) + cout << "\tDouble test (odd) PASSED." << endl; + + // Unit tests for the n dimensional ball volume computation cout << "Testing ball volume computation:" << endl; @@ -198,7 +205,7 @@ int main(int argc, char** argv) { y= -50; pow(epsilon, x, y); for(int i= 1; i < 71; i++) { - if( abs(ball_vol_RR(i, r, prec)-ball_vol(i, r_d)) > epsilon) { + if( abs(ball_vol_RR(i, r)-ball_vol(i, r_d)) > epsilon) { ball_test= false; cout << "\tBall test failed in dimension " << i << endl; } @@ -207,71 +214,6 @@ int main(int argc, char** argv) { if(ball_test) cout << "\tBall test PASSED." << endl; - // Unit tests for node number estimation - - /* Unit tests for the enumeration running time estimator function - //TODO: cjloss destruktort megirni és itt meghivni - //TODO: setupot és teardownt irni a teszthez - - - cjloss l(80, 0.94, 0); - BKZ_QP1(l.basis, 0.99, 2); - - mat_RR mu1; - vec_RR c1; - ComputeGS(l.basis,mu1,c1); - - double* c_d= new double[mu1.NumRows()]; - for(int i= 0; i < mu1.NumRows(); i++) - conv(c_d[i], SqrRoot(c1[i])); - - double* boundary_d= new double[mu1.NumRows()]; - boundary_d[0]= boundary_d[1]= 2; - for(int i= 2; i < mu1.NumRows()-1; i+=2) - boundary_d[i]= boundary_d[i+1]= boundary_d[i-1]+boundary_d[0]; - boundary_d[mu1.NumRows()-1]= mu1.NumRows(); - - RR* c= new RR[mu1.NumRows()]; - for(int i= 0; i < mu1.NumRows(); i++) { - c[i].SetPrecision(prec); - c[i]= SqrRoot(c1[i]); - } - - RR* boundary= new RR[mu1.NumRows()]; - boundary[0].SetPrecision(prec); - boundary[1].SetPrecision(prec); - boundary[0]= boundary[1]= 2; - for(int i= 2; i < mu1.NumRows()-1; i+=2) { - boundary[i].SetPrecision(prec); - boundary[i]= boundary[i+1]= boundary[i-1]+boundary[0]; - } - boundary[mu1.NumRows()-1]= mu1.NumRows(); - - double t_node= 3.47193e-08; - double t_reduc= 0.101471; - if(t_extreme(boundary, c, t_node, t_reduc, mu1.NumRows())-t_extreme_reference(boundary, c, t_node, t_reduc, mu1.NumRows()) < 1e-10) - cout << "PASSED" << endl; - else cout << "FAILED" << endl; - - //cout << "Production (double): " << t_extreme(boundary_d, c_d, t_node, t_reduc, mu1.NumRows()) << endl; - cout << "Reference: " << t_extreme_reference_RR(boundary, c, t_node, t_reduc, mu1.NumRows(), prec) << endl; - cout << "Reference (double): " << t_extreme_reference(boundary_d, c_d, t_node, t_reduc, mu1.NumRows()) << endl; - cout << "Production: " << t_extreme(boundary, c, t_node, t_reduc, mu1.NumRows(), prec) << endl; - cout << "Diff: " << t_extreme(boundary, c, t_node, t_reduc, mu1.NumRows(), prec) - t_extreme_reference_RR(boundary, c, t_node, t_reduc, mu1.NumRows(), prec) << endl; */ - return 0; } -void genbounds(int l, double* tvec){ - for(int i= 0; i < l; i++) - tvec[i]= 1.0/(rand()%1000); - - double swap; - for(int i= 0; i < l-1; i++) - for(int j= 0; j < l-i-1; j++) - if (tvec[j] > tvec[j+1]) { - swap= tvec[j]; - tvec[j]= tvec[j+1]; - tvec[j+1]= swap; - } -} diff --git a/xaa b/xaa new file mode 100644 index 0000000..99dca22 --- /dev/null +++ b/xaa @@ -0,0 +1,460 @@ +/* + Copyright (C) 2007 Victor Shoup. + Copyright (C) 2014 Janos Follath. + + This file is part of cleanbkz. + + cleanbkz is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + cleanbkz is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with cleanbkz. If not, see . + +*/ + +#include +#include +#include +#include +#include + +using namespace std; + +static void enumerate_ntl(double** mu, double *c, double* prune, int jj, int kk, int m, vec_RR& result) { + int s, t; + //double cbar, eta; + double t1; + + double *ctilda; + ctilda = NTL_NEW_OP double[m+2]; + if (!ctilda) Error("ENUMERATE: out of memory"); + + double *vvec; + vvec = NTL_NEW_OP double[m+2]; + if (!vvec) Error("ENUMERATE: out of memory"); + + double *yvec; + yvec = NTL_NEW_OP double[m+2]; + if (!yvec) Error("ENUMERATE: out of memory"); + + double *uvec; + uvec = NTL_NEW_OP double[m+2]; + if (!uvec) Error("ENUMERATE: out of memory"); + + double *utildavec; + utildavec = NTL_NEW_OP double[m+2]; + if (!utildavec) Error("ENUMERATE: out of memory"); + + long *Deltavec; + Deltavec = NTL_NEW_OP long[m+2]; + if (!Deltavec) Error("ENUMERATE: out of memory"); + + long *deltavec; + deltavec = NTL_NEW_OP long[m+2]; + if (!deltavec) Error("ENUMERATE: out of memory"); + + //cbar = c[jj]; + utildavec[jj] = uvec[jj] = 1; + + yvec[jj] = vvec[jj] = 0; + Deltavec[jj] = 0; + + s = t = jj; + deltavec[jj] = 1; + + for (int i = jj+1; i <= kk+1; i++) { + ctilda[i] = uvec[i] = utildavec[i] = yvec[i] = 0; + Deltavec[i] = 0; + vvec[i] = 0; + deltavec[i] = 1; + } + + int nodes= 0; + while (t <= kk) { + + nodes++; + ctilda[t] = ctilda[t+1] + + (yvec[t]+utildavec[t])*(yvec[t]+utildavec[t])*c[t]; + + //cout << nodes << "\tt: " << t << "\tctilda: " << ctilda[t] << "\tprune (kk-t-1): " << prune[kk-t-2] << endl; + + ForceToMem(&ctilda[t]); // prevents an infinite loop + + if ((t>jj) && (ctilda[t] <= prune[kk-t])) { + // if (ctilda[t] <= prune[kk-t]) { + //if (t > jj) { + t--; + t1 = 0; + for (int i = t+1; i <= s; i++) + t1 += utildavec[i]*mu[i][t]; + yvec[t] = t1; + t1 = -t1; + if (t1 >= 0) + t1 = ceil(t1-0.5); + else + t1 = floor(t1+0.5); + utildavec[t] = vvec[t] = t1; + Deltavec[t] = 0; + if (utildavec[t] > -yvec[t]) + deltavec[t] = -1; + else + deltavec[t] = 1; + /*} + else { + //cbar = ctilda[jj]; + break; + for (int i = jj; i <= kk; i++) { + uvec[i] = utildavec[i]; + } + }*/ + } + else { + t++; + s = max(s, t); + if (t < s) Deltavec[t] = -Deltavec[t]; + if (Deltavec[t]*deltavec[t] >= 0) Deltavec[t] += deltavec[t]; + utildavec[t] = vvec[t] + Deltavec[t]; + } + } + + result.SetLength(m); + for(int i= 0; i < result.length(); i++) + conv(result[i],uvec[i]); + + cout << "# NTL nodes: " << nodes << endl; + + delete [] ctilda; + delete [] vvec; + delete [] yvec; + delete [] uvec; + delete [] utildavec; + delete [] Deltavec; + delete [] deltavec; +} + +void enumerate_ntl(mat_ZZ& basis, int beta, double* prune, vec_RR& result) { + + BKZ_QP1(basis, 0.99, beta); + + mat_RR mu1; + vec_RR c1; + ComputeGS(basis,mu1,c1); + + double** mu= new double*[mu1.NumRows()]; + for(int i= 0; i < mu1.NumRows(); i++) + mu[i]= new double[mu1.NumCols()]; + for(int i= 0; i < mu1.NumRows(); i++) + for(int j= 0; j < mu1.NumCols(); j++) + conv(mu[i][j], mu1[i][j]); + + double* c= new double[mu1.NumRows()]; + for(int i= 0; i < mu1.NumRows(); i++) + conv(c[i], c1[i]); + + enumerate_ntl(mu, c, prune, 0, mu1.NumRows()-1, mu1.NumRows(), result); +} + +void profile_enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result) { + + cout << "# Enumerate: " << endl << "# GS lengths: [ "; + for(int i= 0; i < n; i++) + cout << sqrt(b[i]) << " "; + cout << "]" << endl << "# Bounds: [ "; + for(int i= 0; i < n; i++) + cout << sqrt(Rvec[i]) << " "; + cout << "]" << endl << "#" << endl; + + bool pruning= true; + if(Rvec==NULL) { + pruning= false; + Rvec= new double[n]; + for(int i= 0; i < n; i++) + Rvec[i]= b[0]; + } + + double** sigmamat= new double*[n+1]; + for(int i= 0; i < n+1; i++) + sigmamat[i]= new double[n]; + for(int i= 0; i < n+1; i++) + for(int j= 0; j < n; j++) + sigmamat[i][j]= 0; + + int* rvec= new int[n+1]; + for(int i= 0; i < n+1; i++) + rvec[i]= i-1; + + double* rhovec= new double[n+1]; + for(int i= 0; i < n+1; i++) + rhovec[i]= 0; + + int* vvec= new int[n]; + vvec[0]= 1; + for(int i= 1; i < n; i++) + vvec[i]= 0; + + double* cvec= new double[n]; + for(int i= 0; i < n; i++) + cvec[i]= 0; + + int* wvec= new int[n]; + for(int i= 0; i < n; i++) + wvec[i]= 0; + + int last_nonzero= 0; + + unsigned long nodes= 0; + unsigned long* prof= new unsigned long[n]; + for(int i= 0; i < n; i++) + prof[i]= 0; + + int k= 0; + + cout << "\"Measured\"" << endl << endl; + + while(true) { + rhovec[k]= rhovec[k+1]+(vvec[k]-cvec[k])*(vvec[k]-cvec[k])*b[k]; + + //cout << nodes << "\tk: " << k << "\trhovec: " << rhovec[k] << "\tRvec (n-k-1): " << Rvec[n-k-1] << endl; + /*cout << "vvec: " << endl; + for(int i= 0; i < n; i++) + cout << vvec[i] << " "; + cout << endl;*/ + + //nodes++; + //prof[n-k-1]++; + + //if(rhovec[k] <= Rvec[n-k-1]) { + if((k!=0) && (rhovec[k] <= Rvec[n-k-1])) { + + nodes++; + prof[n-k-1]++; + + /*if(k==0) + break; + else {*/ + k--; + rvec[k]= rvec[k]>rvec[k+1]?rvec[k]:rvec[k+1]; + for(int i= rvec[k+1]; i>k; i--) + sigmamat[i][k]= sigmamat[i+1][k] + vvec[i]*mu[i][k]; + cvec[k]= -sigmamat[k+1][k]; + + vvec[k]= lround(cvec[k]); + wvec[k]= 1; + + //} + } else { + if((k==0) && (rhovec[0] <= Rvec[n-1])) { + nodes++; + prof[n-k-1]++; + } + + k++; + if(k==n) + break; + rvec[k]= k; + + if(k>=last_nonzero){ + last_nonzero= k; + vvec[k]++; + } else { + if(vvec[k]>cvec[k]) + vvec[k]-= wvec[k]; + else + vvec[k]+= wvec[k]; + wvec[k]++; + } + } + } + + cout << "# Nodes: " << nodes << endl; + + ofstream myfile; + myfile.open ("ratio.tmp",ios::trunc); + for(int i= 0; i < n; i++) + myfile << prof[i] << " "; + myfile.close(); + + + for(int i= 0; i < n; i++) + cout << i << " " << prof[i] << endl; + + + if(k==0) { + result.SetLength(n); + for(int i= 0; i < result.length(); i++) + conv(result[i],vvec[i]); + } else + result.SetLength(0); + + if(!pruning) + delete [] Rvec; + for(int i= 0; i < n+1; i++) + delete [] sigmamat[i]; + delete [] sigmamat; + delete [] rvec; + delete [] vvec; + delete [] wvec; + delete [] cvec; + delete [] rhovec; +} + +void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result, unsigned long &termination, double &time) { + bool pruning= true; + if(Rvec==NULL) { + pruning= false; + Rvec= new double[n]; + for(int i= 0; i < n; i++) + Rvec[i]= b[0]; + } + + double** sigmamat= new double*[n+1]; + for(int i= 0; i < n+1; i++) + sigmamat[i]= new double[n]; + for(int i= 0; i < n+1; i++) + for(int j= 0; j < n; j++) + sigmamat[i][j]= 0; + + int* rvec= new int[n+1]; + for(int i= 0; i < n+1; i++) + rvec[i]= i-1; + + double* rhovec= new double[n+1]; + for(int i= 0; i < n+1; i++) + rhovec[i]= 0; + + int* vvec= new int[n]; + vvec[0]= 1; + for(int i= 1; i < n; i++) + vvec[i]= 0; + + double* cvec= new double[n]; + for(int i= 0; i < n; i++) + cvec[i]= 0; + + int* wvec= new int[n]; + for(int i= 0; i < n; i++) + wvec[i]= 0; + + int last_nonzero= 0; + + unsigned long nodes= 0; + unsigned long* prof= new unsigned long[n]; + for(int i= 0; i < n; i++) + prof[i]= 0; + + int k= 0; + clock_t end,begin= clock(); + + while((termination==0)||(nodes < termination)) { + rhovec[k]= rhovec[k+1]+(vvec[k]-cvec[k])*(vvec[k]-cvec[k])*b[k]; + + if(rhovec[k] <= Rvec[n-k-1]) { + nodes++; + + if(k==0) + break; + else { + k--; + rvec[k]= rvec[k]>rvec[k+1]?rvec[k]:rvec[k+1]; + for(int i= rvec[k+1]; i>k; i--) + sigmamat[i][k]= sigmamat[i+1][k] + vvec[i]*mu[i][k]; + cvec[k]= -sigmamat[k+1][k]; + + vvec[k]= lround(cvec[k]); + wvec[k]= 1; + + } + } else { + k++; + if(k==n) + break; + rvec[k]= k; + + if(k>=last_nonzero){ + last_nonzero= k; + vvec[k]++; + } else { + if(vvec[k]>cvec[k]) + vvec[k]-= wvec[k]; + else + vvec[k]+= wvec[k]; + wvec[k]++; + } + } + } + + end= clock(); + time= (double) (end-begin) / CLOCKS_PER_SEC / nodes; + + termination-= nodes; + + if(k==0) { + result.SetLength(n); + for(int i= 0; i < result.length(); i++) + conv(result[i],vvec[i]); + } else + result.SetLength(0); + + if(!pruning) + delete [] Rvec; + for(int i= 0; i < n+1; i++) + delete [] sigmamat[i]; + delete [] sigmamat; + delete [] rvec; + delete [] vvec; + delete [] wvec; + delete [] cvec; + delete [] rhovec; +} + +void enumerate_epr(double** mu, double *b, double* Rvec, int n, vec_RR& result) { + double time; + unsigned long termination= 0; + return enumerate_epr(mu, b, Rvec, n, result, termination, time); +} + +void enumerate_epr(mat_ZZ& basis, int beta, double* prune, vec_RR& result) { + double time; + + +// cout << "Original: " << endl << basis << endl; + + BKZ_QP1(basis, 0.99, beta); + +// cout << "Rank: " << rank << endl; +// cout << "Reduced: " << endl << basis << endl; + + //cout << "Basis: " << endl << basis; + + mat_RR mu1; + vec_RR c1; + ComputeGS(basis,mu1,c1); + +// cout << "Reduced: " << endl << mu1 << endl; + + double** mu= new double*[mu1.NumRows()]; + for(int i= 0; i < mu1.NumRows(); i++) + mu[i]= new double[mu1.NumCols()]; + for(int i= 0; i < mu1.NumRows(); i++) + for(int j= 0; j < mu1.NumCols(); j++) + if(i==j) + mu[i][j]= 1; + else + conv(mu[i][j], mu1[i][j]); + + double* c= new double[mu1.NumRows()]; + for(int i= 0; i < mu1.NumRows(); i++) + conv(c[i], c1[i]); + + unsigned long termination= 0; + enumerate_epr(mu, c, prune, mu1.NumRows(), result, termination, time); + +// cout << "# Nodes processed: " << -termination << endl; +}