--- /dev/null
+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
+
--- /dev/null
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <sstream>
+#include <limits>
+#include <fstream>
+#include <cstdlib>
+#include <cmath>
+#include <gmp.h>
+
+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<double>::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<Coord>& max_cand, std::vector<Coord>& min_cand) {
+ double avg = s / n;
+ double min = std::numeric_limits<double>::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<Coord> 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;
+}
+
--- /dev/null
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <map>
+#include <set>
+#include <bitset>
+#include <algorithm>
+#include <iterator>
+#include <sstream>
+#include <limits>
+#include <fstream>
+#include <cstdlib>
+#include <cmath>
+#include <gmp.h>
+
+using namespace std;
+
+double target_mult = 1e300;
+double bestest = 1e300;
+double tt = 0.0;
+vector<set<pair<double, short> > > 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<pair<double, bitset<55> > > ncombs;
+vector<double> prods;
+vector<long long> primes;
+int n;
+
+
+vector<long long> get_primes(int n) {
+ vector <long long> 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<long long>& 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<double, bitset<55> > 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<long long>& 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<int>& 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<int>& 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<vector<int> >& restprimes,
+ double sum,
+ double prod,
+ set<int> 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<int> newpicked = picked;
+ newpicked.insert(prime);
+ complete_rest(
+ restprimes,
+ sum,
+ prod,
+ newpicked,
+ cur,
+ anode + 1);
+ }
+}
+
+//void complete(vector<int>& arr, double score) {
+// vector<vector<int> > restprimes;
+// bitset<55> lastprimes;
+// vector<int> lprimes;
+// lastprimes.set();
+// for (auto it = arr.begin(); it != arr.end(); ++it) {
+// lastprimes &= ~(it->second);
+// vector<int> 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<int> empt;
+// complete_rest(restprimes, score, 1.0, empt, 0);
+//}
+
+int find_best(
+ double score,
+ vector<int>& arr,
+ set<pair<double, short> >& 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<pair<double, short> > 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<pair<double, short> > 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 <size_t BSIZE>
+void print_graph(multimap<double, bitset<BSIZE> >& 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<int> 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;
+}
+
--- /dev/null
+#!/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
+
--- /dev/null
+#include <iostream>
+#include <iomanip>
+#include <algorithm>
+#include <cstdlib>
+#include <vector>
+#include <cmath>
+#include <gmp.h>
+
+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<double> 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<int> 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;
+}
+
projects / primorial-soup / commitdiff