From aee4ba14e6e3367807764b1597aa95662d7a0b6e Mon Sep 17 00:00:00 2001 From: =?utf8?q?Bj=C3=B8rn=20Rustad?= Date: Fri, 3 Nov 2017 19:57:34 +0100 Subject: [PATCH 1/1] Initial commit --- Makefile | 10 ++ main.cpp | 447 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ ng.cpp | 407 +++++++++++++++++++++++++++++++++++++++++++++++++ scores.sh | 11 ++ test.cpp | 273 +++++++++++++++++++++++++++++++++ 5 files changed, 1148 insertions(+) create mode 100644 Makefile create mode 100644 main.cpp create mode 100644 ng.cpp create mode 100755 scores.sh create mode 100644 test.cpp diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..17ab78a --- /dev/null +++ b/Makefile @@ -0,0 +1,10 @@ +all: + g++ -O3 -o primoral main.cpp -lgmp -msse -msse2 -msse3 + g++ -O3 -std=c++11 -o test test.cpp -lgmp -msse -msse2 -msse3 + +ng: ng.cpp + g++ -O3 -std=c++11 -o ng ng.cpp -lgmp -msse -msse2 -msse3 + +debug: + g++ -g -o primoral-debug main.cpp -lgmp + diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..66e24ea --- /dev/null +++ b/main.cpp @@ -0,0 +1,447 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class Graph; +std::ostream& operator<<(std::ostream& o, Graph const& g); + +class Coord { +public: + int row; + int col; + + Coord(int row, int col) : row(row), col(col) {} + Coord() {} +}; + +Coord rand_tuple(int n) { + int row = rand() % (n-1); + int col = rand() % (n-1-row) + 1 + row; + + return Coord(row, col); +} + +Coord rand_perm(Coord c, int n) { + int row_or_col = rand() % 2; + int row, col; + row = c.row; + col = c.col; + int i = 0; + do { + if (row_or_col) { + row = c.row + (rand() % 2)*2 - 1; + } else { + col = c.col + (rand() % 2)*2 - 1; + } + ++i; + } while (i < 2 && (col < row + 1 || row < 0 || col < 0 || row >= n || col >= n)); + + if (i == 2) return c; + + return Coord(row, col); +} + +std::ostream& operator<<(std::ostream& o, Coord& c) { + o << c.row << ", " << c.col; + return o; +} + +class Graph { +public: + long long **graph; + long long *primes; + long long n; + int n_primes; + double *row_pd; + mpz_t *row_prod; + mpz_t target_energy; + double best; + + Graph(long long n) : n(n) { + n_primes = n * (n-1) / 2; + graph = new long long*[n]; + for (int i = 0; i < n; ++i) { + graph[i] = new long long[n]; + } + + for (int i = 0; i < n; ++i) { + graph[i][i] = 1; + } + + row_prod = new mpz_t[n]; + row_pd = new double[n]; + for (int i = 0; i < n; ++i) { + mpz_init(row_prod[i]); + } + int i = 0; + primes = new long long[n*(n-1)/2]; + for (int p = 2; ; ++p) { + bool isprime = true; + for (int j = 0; j < i; ++j) { + if (p % primes[j] == 0) { + isprime = false; + break; + } + } + if (isprime) { + primes[i] = p; + ++i; + } + if (i >= n*(n-1)/2) break; + } + + long long index = 0; + for (int i = 0; i < n; ++i) { + for (int j = i+1; j < n; ++j) { + graph[i][j] = primes[index]; + graph[j][i] = primes[index]; + index++; + } + } + + calc_row_prod(); + calc_target_energy(); + best = read_best(); + shuffle(1000); + std::cout << "Energy to beat from file is " << best << std::endl; + gmp_printf("Target energy is %Zd\n", target_energy); + } + + double read_best() { + double best; + std::stringstream ss; + ss << "scores/" << n; + std::ifstream infile(ss.str().c_str()); + if (infile.good()) { + infile >> best; + infile.close(); + return best; + } + return std::numeric_limits::max(); + } + + void write() { + std::stringstream ss; + ss << "scores/" << n; + std::ofstream outfile(ss.str().c_str()); + if (outfile.good()) { + outfile << std::setprecision(90); + outfile << score() << std::endl; + outfile << *this << std::endl; + } else { + std::cout << "PROBLEM OPENING FILE" << std::endl; + } + outfile.close(); + } + + void shuffle(int k) { + for (int i = 0; i < k; ++i) { + Coord from = rand_tuple(n); + Coord to = rand_tuple(n); + if (from.row == to.row && from.col == to.col) continue; + jiggle(from, to); + } + } + + void calc_target_energy() { + mpz_init(target_energy); + mpz_set_si(target_energy, 1); + for (int i = 0; i < n_primes; ++i) { + mpz_mul_si(target_energy, target_energy, primes[i]); + } + mpz_pow_ui(target_energy, target_energy, 2); + mpz_root(target_energy, target_energy, n); + mpz_mul_si(target_energy, target_energy, n); + } + + void get_candidates(double s, std::vector& max_cand, std::vector& min_cand) { + double avg = s / n; + double min = std::numeric_limits::max(); + double max = 0; + int mini, maxi; + for (int i = 0; i < n; ++i) { + double d = mpz_get_d(row_prod[i]); + if (d < min) { + min = d; + mini = i; + } + if (d > max) { + max = d; + maxi = i; + } + if (d < avg * 0.3) { + for (int j = 0; j < n; ++j) { + min_cand.push_back(Coord(i, j)); + } + } + if (d > avg * 1.7) { + for (int j = 0; j < n; ++j) { + max_cand.push_back(Coord(i, j)); + } + } + } + for (int i = 0; i < n; ++i) { + double d = mpz_get_d(row_prod[i]); + if (maxi != i && d > avg * 2) { + max_cand.push_back(Coord(maxi, i)); + } + if (mini != i && d < avg * 0.5) { + min_cand.push_back(Coord(mini, i)); + } + } + } + + double get_target_energy() { + return mpz_get_d(target_energy); + } + + void calc_row_prod() { + for (int i = 0; i < n; ++i) { + mpz_t mult; + mpz_init(mult); + mpz_set_si(mult, 1); + for (int j = 0; j < n; ++j) { + mpz_mul_si(mult, mult, graph[i][j]); + } + mpz_set(row_prod[i], mult); + mpz_clear(mult); + } + } + + double score() { + mpz_t score; + mpz_init(score); + mpz_set_si(score, 0); + for (int i = 0; i < n; ++i) { + mpz_add(score, score, row_prod[i]); + } + mpz_sub(score, score, target_energy); + + double d = mpz_get_d(score); + mpz_clear(score); + return d; + } + + void ddiv(int i, long long div) { + mpz_divexact_ui(row_prod[i], row_prod[i], div); + } + void mmul(int i, long long mul) { + mpz_mul_si(row_prod[i], row_prod[i], mul); + } + + void calc_approximate_row_prods() { + for (int i = 0; i < n; ++i) { + row_pd[i] = mpz_get_d(row_prod[i]); + } + } + + double approximate_score_change(Coord from, Coord to) { + double change = 0; + long long fromn = graph[from.row][from.col]; + long long ton = graph[to.row][to.col]; + + change += row_pd[from.row] * (double(ton) / fromn - 1.0); + change += row_pd[from.col] * (double(ton) / fromn - 1.0); + change += row_pd[to.row] * (double(fromn) / ton - 1.0); + change += row_pd[to.col] * (double(fromn) / ton - 1.0); + + return change; + } + + void jiggle(Coord from, Coord to) { + long long fromn = graph[from.row][from.col]; + long long ton = graph[to.row][to.col]; + mmul(from.row, ton); + mmul(from.col, ton); + mmul(to.row, fromn); + mmul(to.col, fromn); + ddiv(from.row, fromn); + ddiv(from.col, fromn); + ddiv(to.row, ton); + ddiv(to.col, ton); + //row_pd[from.row] *= double(ton)/fromn; + //row_pd[from.col] *= double(ton)/fromn; + //row_pd[to.row] *= double(fromn)/ton; + //row_pd[to.col] *= double(fromn)/ton; + graph[to.row][to.col] = fromn; + graph[to.col][to.row] = fromn; + graph[from.row][from.col] = ton; + graph[from.col][from.row] = ton; + } + + void mac(mpz_t *c, int idx, long long mul, long long div) { + mpz_t change; + mpz_t res; + mpz_init(change); + mpz_init(res); + + mpz_mul_si(res, row_prod[idx], mul); + mpz_divexact_ui(res, res, div); + + mpz_set(change, res); + mpz_sub(change, change, row_prod[idx]); + + mpz_add(*c, *c, change); + mpz_clear(change); + mpz_clear(res); + } + + double get_change(Coord from, Coord to) { + mpz_t change; + mpz_init(change); + mpz_set_si(change, 0); + long long fromn = graph[from.row][from.col]; + long long ton = graph[to.row][to.col]; + + mac(&change, from.row, ton, fromn); + mac(&change, from.col, ton, fromn); + mac(&change, to.row, fromn, ton); + mac(&change, to.col, fromn, ton); + //change += row_pd[from.row] * (double(ton) / fromn - 1.0); + //change += row_pd[from.col] * (double(ton) / fromn - 1.0); + //change += row_pd[to.row] * (double(fromn) / ton - 1.0); + //change += row_pd[to.col] * (double(fromn) / ton - 1.0); + + double c = mpz_get_d(change); + mpz_clear(change); + return c; + } + + bool optimize() { + double before = score(); + for (int i = 0; i < n; i++) { + for (int j = i+1; j < n; j++) { + for (int u = 0; u < n; u++) { + for (int v = u+1; v < n; v++) { + if (i == u && j == v) + continue; + double s = score(); + Coord from(i,j); + Coord to(u,v); + jiggle(from, to); + double new_s = score(); + if (new_s > s) + jiggle(to, from); + } + } + } + } + double after = score(); + if (after < before) return true; + else return false; + } +}; + +std::ostream& operator<<(std::ostream& o, Graph const& g) { + for (int i = 0; i < g.n; ++i) { + o << "{"; + for (int j = 0; j < g.n; ++j) { + if (i == j) continue; + o << g.graph[i][j]; + if (j >= g.n-1) continue; + if (i == g.n-1 && j == g.n-2) break; + o << ", "; + } + o << "}"; + if (i != g.n-1) o << ", "; + } + + return o; +} + + +int main(int argc, char* argv[]) { + Graph g(atoi(argv[1])); + srand(time(NULL)); + std::cout << g << std::endl; + std::cout << "SCORE: " << g.score() << std::endl; + std::cout << std::setprecision(10); + double cur = g.score(); + double min = g.score(); + double T = atof(argv[2]); + int i = 0; + double last = 0; + int num_picks = 100;; + int moves = 0; + std::vector cmin, cmax; + while (true) { + Coord from; + Coord to; + if (double(rand())/RAND_MAX > 0.1) { + if (num_picks > 100) { + cmin.clear(); + cmax.clear(); + g.get_candidates(cur, cmax, cmin); + num_picks = 0; + } + //std::cout << "Got " << cand.size() << " candidates" << std::endl; + //for (int j = 0; j < cand.size(); ++j) { + // std::cout << cand[j] << ";"; + //} + if (cmin.size() <= 0) continue; + if (cmax.size() <= 0) continue; + int u = rand() % (cmax.size()); + int v = rand() % (cmin.size()); + from = cmax[u]; + to = cmin[v]; + } else { + from = rand_tuple(g.n); + to = rand_perm(from, g.n); + } + //std::cout << "FROM = " << from << " TO = " << to << std::endl; + if (from.row == to.row && from.col == to.col) continue; + g.jiggle(from, to); + double score = g.score(); + double change = score - cur; + //double getcc = g.get_change(from, to); + //if (abs(change - getcc) > 0.0001) { + // std::cout << "BIG DIFF: " << change << " - " << getcc << std::endl; + //} + double acc_prob = exp(-(change)/T); + if (change < 0) { + cur = score; + moves++; + } else if (acc_prob > double(rand()) / RAND_MAX) { + cur = score; + moves++; + } else if (double(rand()) / RAND_MAX > 0.99 || T > 1e240) { + cur = score; + } else { + g.jiggle(to, from); + } + if (cur < min) { + min = cur; + std::cout << "NEW MIN: " << min << " <> " << g.best << std::endl; + std::cout << "\ttemp = " << T << std::endl; + } + if (i % 100000 == 0) { + std::cout << "CURRENT: " << cur << " <> " << g.best << std::endl; + std::cout << "\ttemp = " << T << std::endl; + } + if (cur < g.best) { + g.best = cur; + std::cout << "NEW BEST MAN!!!!!!!!!!!!! " << g.best << std::endl; + g.write(); + } + T *= atof(argv[3]); + if (i % 10000 == 0) { + //if (fabs(last - cur) < 0.01) { + if (moves == 0) { + T *= pow(2.0 - atof(argv[3]), 1000000) * 3; + g.shuffle(g.n); + std::cout << "No moves lately " << T << std::endl; + } + last = cur; + moves = 0; + } + i++; + } + return 0; +} + diff --git a/ng.cpp b/ng.cpp new file mode 100644 index 0000000..d33adb6 --- /dev/null +++ b/ng.cpp @@ -0,0 +1,407 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +double target_mult = 1e300; +double bestest = 1e300; +double tt = 0.0; +vector > > bsets; +long long dbad = 0; +long long dbad_limit = 0; +long long without_success = 0; +long long success_limit = 0; +double min_score = 1e100; +int beam_width = 10; +vector > > ncombs; +vector prods; +vector primes; +int n; + + +vector get_primes(int n) { + vector primes; + for (int p = 2; ; ++p) { + bool isprime = true; + for (int j = 0; j < primes.size(); ++j) { + if (p % primes[j] == 0) { + isprime = false; + break; + } + } + if (isprime) { + primes.push_back(p); + } + if (primes.size() >= n) break; + } + return primes; +} + +double get_target_mult(vector& primes, int n) { + mpz_t target_energy; + mpz_init(target_energy); + mpz_set_si(target_energy, 1); + for (int i = 0; i < primes.size(); ++i) { + mpz_mul_si(target_energy, target_energy, primes[i]); + } + mpz_pow_ui(target_energy, target_energy, 2); + mpz_root(target_energy, target_energy, n); + mpz_mul_si(target_energy, target_energy, n); + double ret = mpz_get_d(target_energy); + mpz_clear(target_energy); + return ret; +} + +double upper = 100; +double lower = 100; + +multimap > combs; + +int siz = 0; +int try_comb(int placed, int start, double prod, bitset<55> b) { + //cout << "placed: " << placed << endl; + if (placed == n-1) { + if (prod > lower && prod < upper) { + combs.insert(make_pair(fabs(prod - target_mult), b)); + siz++; + if (siz % 100 == 0) cout << "Added 1, now: " << siz << endl; + return 0; + } + return 0; + } else if (prod > upper) { + return 1; + } + for (int i = start; i < primes.size(); ++i) { + b[i] = 1; + int ret = try_comb(placed + 1, i+1, prod * primes[i], b); + b[i] = 0; + if (ret == 1) // too big + return 0; + } + return 0; +} + +void smart_search() { + bitset<55> b; + try_comb(0, 0, 1.0, b); +} + +void exhaustive_search(vector& primes, int n, double target) { + string bitmask(n-1, 1); // K leading 1's + bitmask.resize(primes.size(), 0); // N-K trailing 0's + + // print integers and permute bitmask + do { + double prod = 1.0; + bitset<55> b; + for (int i = 0; i < primes.size(); ++i) // [0..N-1] integers + { + if (bitmask[i]) { + prod *= primes[i]; + b[i] = 1; + } + } + //combs.insert(make_pair(prod, b)); + if (prod > lower && prod < upper) + combs.insert(make_pair(fabs(prod - target), b)); + } while (prev_permutation(bitmask.begin(), bitmask.end())); +} + +void print_arr(vector& arr) { + + for (auto it = arr.begin(); it != arr.end(); ++it) { + cout << "{"; + int prints = 0; + for (int j = 0; j < ncombs[*it].second.size(); ++j) { + if (ncombs[*it].second[j]) { + cout << primes[j]; + prints++; + if (prints <= n-2) cout << ", "; + } + } + cout << "}"; + if (it + 1 != arr.end()) cout << ", "; + } +} + +double score_arr(vector& arr) { + double res = 0.0; + for (auto it = arr.begin(); it != arr.end(); ++it) { + res += prods[*it]; + } + return res; +} + + +void print_bitset(bitset<55> b) { + double prod = 1.0; + cout << "{"; + int printed = 0; + for (int j = 0; j < b.size(); ++j) { + if (b[j]) { + cout << primes[j]; + printed++; + if (printed < n-1) cout << ", "; + } + } + cout << "}"; +} + +double score_bitset(bitset<55> b) { + double prod = 1.0; + for (int j = 0; j < b.size(); ++j) { + if (b[j]) { + prod *= primes[j]; + } + } + return prod; +} + +void complete_rest( + vector >& restprimes, + double sum, + double prod, + set picked, + int cur, + int anode + ) { + for (int j = 0; j < restprimes[cur].size(); ++j) { + int prime = primes[restprimes[cur][j]]; + if (picked.find(prime) != picked.end()) continue; + + set newpicked = picked; + newpicked.insert(prime); + complete_rest( + restprimes, + sum, + prod, + newpicked, + cur, + anode + 1); + } +} + +//void complete(vector& arr, double score) { +// vector > restprimes; +// bitset<55> lastprimes; +// vector lprimes; +// lastprimes.set(); +// for (auto it = arr.begin(); it != arr.end(); ++it) { +// lastprimes &= ~(it->second); +// vector a; +// restprimes.push_back(a); +// bitset<55> bb; +// for (auto jt = arr.begin(); jt != arr.end(); ++jt) { +// if (*it == *jt) continue; +// bb |= ncombs[*jt].second; +// } +// auto restbits = ncombs[*it].second & (~bb); +// +// for (int j = 0; j < restbits.size(); ++j) { +// if (restbits[j]) +// restprimes.back().push_back(primes[j]); +// } +// } +// for (int j = 0; j < lastprimes.size(); ++j) { +// lprimes.push_back(primes[j]); +// } +// +// set empt; +// complete_rest(restprimes, score, 1.0, empt, 0); +//} + +int find_best( + double score, + vector& arr, + set >& allowed, + bitset<55> rest, + bitset<55> current, + bitset<55> not_allowed, + int maxi, + int mults) { + + if (mults >= n) { + double ss = score_arr(arr); + if (ss - tt < bestest) { + cout << "YEEEEEE: " << endl; + bestest = ss - tt; + print_arr(arr); + cout << " points: " << ss - tt << endl; + return 1; + } else { + cout << "At the end but not good enough..." << endl; + return 0; + } + } else if (score - tt > bestest) { + return 0; + } else if (mults == n - 1) { + //cout << "ONLY ONE LEFT: " << rest << endl; + double ss = score_arr(arr) + score_bitset(rest); + if (ss - tt < bestest) { + cout << "YEEEEEE: " << endl; + bestest = ss - tt; + print_arr(arr); + cout << ", "; + print_bitset(rest); + cout << endl; + cout << " points: " << ss - tt << endl; + return 1; + } else { + cout << "We tried but it wasn't good enough..." << endl; + return 0; + } + return 0; + } else if (score - tt > 8691987*1) { + return 0; + } else if (score + min_score * (n - mults) - tt > 8691987*2) { + return 0; + } else if (allowed.size() == 0) { + return 0; + } else { + int num_tried = 0; + set > new_allowed; + for (auto it = allowed.begin(); it != allowed.end(); ++it) { + if (it->second <= maxi) continue; + if (score + prods[it->second] - tt > 8691987*2) continue; + if ((not_allowed & ncombs[it->second].second).count() > 0) continue; + auto p = ncombs[it->second]; + //if ((p.second & cur).count() != mults) continue; + new_allowed.clear(); + set_intersection(it, allowed.end(), + bsets[it->second].begin(), bsets[it->second].end(), + inserter(new_allowed, new_allowed.begin())); + + arr.push_back(it->second); + int ret = find_best( + score + prods[it->second], + arr, + new_allowed, + p.second ^ rest, + current | ncombs[it->second].second, + not_allowed | (ncombs[it->second].second & current), + it->second, + mults + 1); + arr.pop_back(); + num_tried++; + + if (ret == 0 && mults > 1 && num_tried > beam_width) { + return 0; + } + } + return 0; + } + return 0; +} + +void get_sets() { + double sum = 0; + int outer = 0; + for (auto it = ncombs.begin(); it != ncombs.end(); ++it) { + set > b; + short i = -1; + for (auto jt = ncombs.begin(); jt != ncombs.end(); ++jt) { + ++i; + if ((it->second & jt->second).count() == 1) { + b.insert(make_pair(prods[i], i)); + } + } + sum += b.size(); + bsets.push_back(b); + outer++; + } + cout << "AVERAGE SET SIZE: " << sum / bsets.size() << endl; +} + +template +void print_graph(multimap >& combs) { + + for (auto it = combs.begin(); it != combs.end(); ++it) { + for (auto jt = combs.begin(); jt != combs.end(); ++jt) { + if ((it->second & jt->second).count() == 1) { + cout << "1 "; + } else { + cout << "0 "; + } + } + cout << endl; + } +} + +void print_all() { + for (int i = 0; i < ncombs.size(); ++i) { + for (auto it = bsets[i].begin(); it != bsets[i].end(); ++it) { + if (it->second < i) + cout << i+1 << " " << it->second+1 << endl; + } + } +} + +int main(int argc, char* argv[]) { + n = atoi(argv[1]); + primes = get_primes((n*(n-1))/2); + cout << "PRIMES:" << endl; + cout << setprecision(90); + for (auto it = primes.begin(); it != primes.end(); ++it) { + cout << *it << endl; + } + target_mult = get_target_mult(primes, n); + tt = target_mult; + target_mult /= n; + cout << "TARGET ENERGY: " << tt << endl; + double diff = target_mult * atof(argv[2]); //000001; + upper = target_mult + diff; + lower = target_mult - diff; + beam_width = atoi(argv[3]); + + cout << "Number of primes: " << primes.size() << endl; + cout << "Target node product: " << target_mult << endl; + //auto combs = exhaustive_search<55>(primes, n, target_mult); + smart_search(); + cout << "GOT COMBS: " << combs.size() << endl; + + + int i = 0; + for (auto it = combs.begin(); it != combs.end(); ++it) { + cout << it->first << " -> " << it->second << endl; + if (i >= 10) break; + ++i; + } + for (auto it = combs.begin(); it != combs.end(); ++it) { + ncombs.push_back(make_pair(it->first, it->second)); + double s = score_bitset(it->second); + prods.push_back(s); + if (min_score > s) min_score = s; + } + cout << "FOund " << ncombs.size() << " combs" << endl; + + cout << "GETTING SETS" << endl; + get_sets(); + //print_graph(combs); + cout << "GOTEM " << bsets.size() << endl; + //print_all(); + + vector arr; + i = 0; + bitset<55> em; + for (auto it = ncombs.begin(); it != ncombs.end(); ++it) { + if (i % 10 == 0) cout << i << endl; + arr.push_back(i); + int ret = find_best(prods[i], arr, bsets[i], it->second, it->second, em, i, 1); + arr.clear(); + ++i; + } + return 0; +} + diff --git a/scores.sh b/scores.sh new file mode 100755 index 0000000..214fabb --- /dev/null +++ b/scores.sh @@ -0,0 +1,11 @@ +#!/bin/bash +for i in `seq 4 28` +do + a=0 + if [ -e scores/$i ] + then + a=`head -n 1 scores/$i` + fi + echo "$i $a" +done + diff --git a/test.cpp b/test.cpp new file mode 100644 index 0000000..a02d5bf --- /dev/null +++ b/test.cpp @@ -0,0 +1,273 @@ +#include +#include +#include +#include +#include +#include +#include + +class Edge { +public: + int u; + int v; + + Edge(int u, int v) : u(u), v(v) {} + Edge() {} +}; + +std::ostream& operator<<(std::ostream& o, Edge& e) { + o << e.u << "-" << e.v; + return o; +} + +class Graph { +public: + long long *primes; + long long n; + double *row_pd; + mpz_t *row_prod; + Edge *edges; + + Graph(long long n) : n(n) { + row_prod = new mpz_t[n]; + row_pd = new double[n]; + for (int i = 0; i < n; ++i) { + mpz_init(row_prod[i]); + } + int i = 0; + primes = new long long[n*(n-1)/2]; + for (int p = 2; ; ++p) { + bool isprime = true; + for (int j = 0; j < i; ++j) { + if (p % primes[j] == 0) { + isprime = false; + break; + } + } + if (isprime) { + primes[i] = p; + ++i; + } + if (i >= n*(n-1)/2) break; + } + + edges = new Edge[n*(n-1)/2]; + int idx = 0; + for (int i = 0; i < n; i++) { + for (int j = i + 1; j < n; ++j) { + edges[idx] = Edge(i, j); + idx++; + } + } + + calc_row_prod(); + } + + void print() { + std::cout << "THE PRIMES" << std::endl; + for (int i = 0; i < n*(n-1)/2; ++i) { + std::cout << primes[i] << std::endl; + } + + std::cout << "THE EDGES" << std::endl; + for (int i = 0; i < n*(n-1)/2; ++i) { + std::cout << edges[i] << std::endl; + } + } + + void calc_row_prod() { + for (int i = 0; i < n; ++i) { + mpz_t mult; + mpz_init(mult); + mpz_set_si(mult, 1); + for (int j = 0; j < n*(n-1)/2; ++j) { + if (edges[j].u == i || edges[j].v == i) + mpz_mul_si(mult, mult, primes[j]); + } + mpz_set(row_prod[i], mult); + mpz_clear(mult); + } + } + + double score() { + mpz_t score; + mpz_init(score); + mpz_set_si(score, 0); + for (int i = 0; i < n; ++i) { + mpz_add(score, score, row_prod[i]); + } + + double d = mpz_get_d(score); + mpz_clear(score); + return d; + } + + void print_info() { + std::vector v; + for (int i = 0; i < n; ++i) { + double d = mpz_get_d(row_prod[i]); + v.push_back(d); + std::cout << d << ", "; + } + std::cout << std::endl; + double sum = std::accumulate(v.begin(), v.end(), 0.0); + double mean = sum / v.size(); + + double sq_sum = std::inner_product(v.begin(), v.end(), v.begin(), 0.0); + double stdev = std::sqrt(sq_sum / v.size() - mean * mean); + std::cout << "AVG: " << mean << std::endl; + std::cout << "STD: " << stdev << std::endl; + std::cout << "MIN: " << *std::min_element(v.begin(), v.end()) << std::endl; + std::cout << "MAX: " << *std::max_element(v.begin(), v.end()) << std::endl; + } + + void ddiv(int i, long long div) { + //std::cout << "dividing by " << div << std::endl; + mpz_divexact_ui(row_prod[i], row_prod[i], div); + } + void mmul(int i, long long mul) { + //std::cout << "mul by " << mul << std::endl; + mpz_mul_si(row_prod[i], row_prod[i], mul); + } + + void shuffle(int k) { + for (int i = 0; i < k; ++i) { + int from = rand() % (n * (n-1) / 2); + int to = rand() % (n * (n-1) / 2); + if (to == from) continue; + jiggle(from, to); + } + } + + void jiggle(int from, int to) { + Edge e = edges[from]; + mmul(edges[from].u, primes[to]); + mmul(edges[from].v, primes[to]); + mmul(edges[to].u, primes[from]); + mmul(edges[to].v, primes[from]); + ddiv(edges[from].u, primes[from]); + ddiv(edges[from].v, primes[from]); + ddiv(edges[to].u, primes[to]); + ddiv(edges[to].v, primes[to]); + edges[from] = edges[to]; + edges[to] = e; + } + +}; + +std::ostream& operator<<(std::ostream& o, const Graph& g) { + for (int i = 0; i < g.n; ++i) { + std::cout << "{"; + int pp = 0; + for (int j = 0; j < (g.n * (g.n-1)/2); ++j) { + if (g.edges[j].u == i || g.edges[j].v == i) { + std::cout << g.primes[j]; + pp++; + if (pp < g.n-1) { + std::cout << ", "; + } + } + } + std::cout << "}"; + if (i != g.n-1) std::cout << ", "; + } + + return o; +} + +double get_target_mult(long long *primes, int n) { + mpz_t target_energy; + mpz_init(target_energy); + mpz_set_si(target_energy, 1); + for (int i = 0; i < n*(n-1) / 2; ++i) { + mpz_mul_si(target_energy, target_energy, primes[i]); + } + mpz_pow_ui(target_energy, target_energy, 2); + mpz_root(target_energy, target_energy, n); + mpz_mul_si(target_energy, target_energy, n); + double ret = mpz_get_d(target_energy); + mpz_clear(target_energy); + return ret; +} + + + +int main(int argc, char* argv[]) { + Graph g(atoi(argv[1])); + srand(time(NULL)); + std::cout << g << std::endl; + std::cout << "SCORE: " << g.score() << std::endl; + std::cout << std::setprecision(10); + std::cout << "SHUFFLIGN" << std::endl; + g.shuffle(10000000); + std::cout << "SHUFFLIGN DONE" << std::endl; + double cur = g.score(); + double min = g.score(); + double T = atof(argv[2]); + double target = get_target_mult(g.primes, atoi(argv[1])); + std::cout << "TARGET IS: " << target << std::endl; + int i = 0; + int moves = 1; + g.print(); + int fromi = 0; + std::vector take_from; + for (int ii = 0; ii < g.n*(g.n-1)/2; ++ii) { + take_from.push_back(ii); + } + while (true) { + fromi++; + if (fromi >= g.n * (g.n-1) / 2) { + std::random_shuffle(take_from.begin(), take_from.end()); + fromi = 0; + } + int from = take_from[fromi]; + //int from = rand() % (g.n * (g.n-1) / 2); + //int to = rand() % (g.n * (g.n-1) / 2); + //int to2 = rand() % (g.n * (g.n-1) / 2); + int to = from + (rand() % 7); + int to2 = from - (rand() % 7); + if (to >= g.n * (g.n-1)/2) to = g.n * (g.n-1)/2 - 1; + if (to2 >= g.n * (g.n-1)/2) to2 = g.n * (g.n-1)/2 - 1; + if (to < 0) to = 0; + if (to2 < 0) to2 = 0; + if (to == to2) continue; + if (to == from) continue; + if (to2 == from) continue; + g.jiggle(from, to); + g.jiggle(from, to2); + double score = g.score(); + double change = log(score) - log(cur); + double acc_prob = exp(-(change)/T); + if (change < 0.0 || (acc_prob > double(rand()) / RAND_MAX)) { + if (fabs(change) > 0.0000000000001) + moves += 1; + cur = score; + } else { + g.jiggle(to2, from); + g.jiggle(to, from); + } + if (cur < min) { + min = cur; + std::cout << "NEW MIN: " << min - target << std::endl; + std::cout << g << std::endl; + std::cout << "temp = " << T << std::endl; + std::cout << "prob = " << acc_prob << std::endl; + //g.print(); + g.print_info(); + } + T *= atof(argv[3]); + if (i % 10000000 == 0) { + if (moves <= 0) { + T /= pow(atof(argv[3]), 100000000); + std::cout << "No moves lately " << T << ", score:" << cur - target << std::endl; + std::cout << "Shuffling" << std::endl; + g.shuffle(10); + cur = g.score(); + } + moves = 0; + } + i++; + } + return 0; +} + -- 2.47.3