Before simd refactor

diff --git a/Makefile b/Makefile
index 17ab78a98bc54b6aa8953461ae926a9b6f3ffa36..9e03c09ddeb0c69056539868313743c786f74794 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,10 +1,27 @@
+CC ?= gcc
+CXX ?= g++
+
+CFLAGS ?= -O3 -lgmp -msse -msse2 -msse3 -msse4.2 -mavx -mavx2
+CPPFLAGS ?= -std=c++11
+
 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
+OBJECTS=slib/bitpacking.o slib/integratedbitpacking.o slib/simdbitpacking.o slib/usimdbitpacking.o slib/simdintegratedbitpacking.o slib/intersection.o slib/varintdecode.o slib/streamvbyte.o slib/simdpackedsearch.o slib/simdpackedselect.o slib/frameofreference.o slib/for.o
+
+HEADERS=$(shell ls slib/include/*h)
+
+.c.o:
+       $(CC) -c $(CFLAGS) $(CPPFLAGS) -o $@ $<
+.cpp.o:
+       $(CXX) -c $(CFLAGS) $(CPPFLAGS) -o $@ $<
+
+
+ng: ng.cpp $(OBJECTS) $(HEADERS)
+       g++ -O3 -std=c++11 -o ng ng.cpp $(OBJECTS) -lgmp -msse -msse2 -msse3 -msse4 -msse4.2 -mavx -mavx2 -ffast-math -Islib/include
 
 debug:
+       g++ -std=c++11 -g -o ng-dbg ng.cpp -lgmp
        g++ -g -o primoral-debug main.cpp -lgmp
 

diff --git a/ng.cpp b/ng.cpp
index d33adb69cf9c446d09f6837c8ecea1e1c2a26df2..de9da278854fb0a05134ac1822a1656d0a25c059 100644 (file)
--- a/ng.cpp
+++ b/ng.cpp
@@ -12,24 +12,65 @@
 #include <cstdlib>
 #include <cmath>
 #include <gmp.h>
+#include "codecfactory.h"
+#include "intersection.h"
+#include <ctime>
 
 using namespace std;
+using namespace SIMDCompressionLib;
 
 double target_mult = 1e300;
 double bestest = 1e300;
-double tt = 0.0;
-vector<set<pair<double, short> > > bsets;
+double total_target = 0.0;
+vector<set<int>> 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<pair<double, bitset<36> > > ncombs;
 vector<double> prods;
 vector<long long> primes;
 int n;
 
+void test_slib() {
+  int N = 1000;
+  intersectionfunction inter =
+    IntersectionFactory::getFromName("simd");
+  vector<uint32_t> mydata1(N);
+  vector<uint32_t> mydata2(N);
+  vector<uint32_t> mydata3(N);
+  for (uint32_t i = 0; i < N; ++i) {
+    mydata1[i] = 3 * i;
+    mydata2[i] = 2 * i;
+  }
+  clock_t start = clock();
+  for (int j = 0; j < 100000; ++j) {
+    size_t intersize = inter(mydata1.data(), mydata1.size(), mydata2.data(),
+                             mydata2.size(), mydata3.data());
+
+    mydata3.resize(intersize);
+    mydata3.shrink_to_fit();
+  }
+  clock_t end = clock();
+  cout << "TIME FOR SIMD: " << end - start << endl;
+
+  set<uint32_t> s1(mydata1.begin(), mydata1.end());
+  set<uint32_t> s2(mydata2.begin(), mydata2.end());
+  set<uint32_t> s3;
+  start = clock();
+  for (int j = 0; j < 100000; ++j) {
+    set_intersection(s1.begin(), s1.end(), s2.begin(), s2.end(), inserter(s3, s3.begin()));
+  }
+  end = clock();
+  cout << "TIME FOR SET INTERSECTION: " << end - start << endl;
+  cout << "Intersection size: " << mydata3.size() << "  integers. " << endl;
+
+  for (uint32_t i = 0; i < mydata3.size(); ++i) {
+    cout << mydata3[i] << ", ";
+  }
+  cout << endl;
+}
 
 vector<long long> get_primes(int n) {
        vector <long long> primes;
@@ -67,10 +108,10 @@ double get_target_mult(vector<long long>& primes, int n) {
 double upper = 100;
 double lower = 100;
 
-multimap<double, bitset<55> > combs;
+multimap<double, bitset<36> > combs;
 
 int siz = 0;
-int try_comb(int placed, int start, double prod, bitset<55> b) {
+int try_comb(int placed, int start, double prod, bitset<36> b) {
        //cout << "placed: " << placed << endl;
        if (placed == n-1) {
                if (prod > lower && prod < upper) {
@@ -94,7 +135,7 @@ int try_comb(int placed, int start, double prod, bitset<55> b) {
 }
 
 void smart_search() {
-       bitset<55> b;
+       bitset<36> b;
        try_comb(0, 0, 1.0, b);
 }
 
@@ -105,7 +146,7 @@ void exhaustive_search(vector<long long>& primes, int n, double target) {
        // print integers and permute bitmask
        do {
                double prod = 1.0;
-               bitset<55> b;
+               bitset<36> b;
                for (int i = 0; i < primes.size(); ++i) // [0..N-1] integers
                {
                        if (bitmask[i]) {
@@ -145,7 +186,7 @@ double score_arr(vector<int>& arr) {
 }
 
 
-void print_bitset(bitset<55> b) {
+void print_bitset(bitset<36> b) {
        double prod = 1.0;
        cout << "{";
        int printed = 0;
@@ -159,7 +200,7 @@ void print_bitset(bitset<55> b) {
        cout << "}";
 }
 
-double score_bitset(bitset<55> b) {
+double score_bitset(bitset<36> b) {
        double prod = 1.0;
        for (int j = 0; j < b.size(); ++j) {
                if (b[j]) {
@@ -193,112 +234,73 @@ void complete_rest(
        }
 }
 
-//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,
+               set<int>& allowed,
+               bitset<36> rest,
+               bitset<36> current,
+               bitset<36> not_allowed,
                int maxi,
                int mults) {
 
        if (mults >= n) {
                double ss = score_arr(arr);
-               if (ss - tt < bestest) {
+               if (ss - total_target < bestest) {
                        cout << "YEEEEEE: " << endl;
-                       bestest = ss - tt;
+                       bestest = ss - total_target;
                        print_arr(arr);
-                       cout << " points: " << ss - tt << endl;
+                       cout << " points: " << ss - total_target << endl;
                        return 1;
                } else {
                        cout << "At the end but not good enough..." << endl;
                        return 0;
                }
-       } else if (score - tt > bestest) {
+       } else if (score - total_target > 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) {
+               if (ss - total_target < bestest) {
                        cout << "YEEEEEE: " << endl;
-                       bestest = ss - tt;
+                       bestest = ss - total_target;
                        print_arr(arr);
                        cout << ", ";
                        print_bitset(rest);
                        cout << endl;
-                       cout << " points: " << ss - tt << endl;
+                       cout << " points: " << ss - total_target << endl;
                        return 1;
                } else {
-                       cout << "We tried but it wasn't good enough..." << endl;
+                       //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(),
+      if (*it <= maxi) continue;
+                       if ((not_allowed & ncombs[*it].second).count() > 0) continue;
+
+                       auto p = ncombs[*it];
+
+      set<int> new_allowed;
+
+                       set_intersection(allowed.begin(), allowed.end(),
+                                       bsets[*it].begin(), bsets[*it].end(),
                                        inserter(new_allowed, new_allowed.begin()));
 
-                       arr.push_back(it->second);
+                       arr.push_back(*it);
+
                        int ret = find_best(
-                                       score + prods[it->second],
+                                       score + prods[*it],
                                        arr,
                                        new_allowed,
                                        p.second ^ rest,
-                                       current | ncombs[it->second].second,
-                                       not_allowed | (ncombs[it->second].second & current),
-                                       it->second,
+                                       current | ncombs[*it].second,
+                                       not_allowed | (ncombs[*it].second & current),
+                                       *it,
                                        mults + 1);
-                       arr.pop_back();
-                       num_tried++;
 
-                       if (ret == 0 && mults > 1 && num_tried > beam_width) {
-                               return 0;
-                       }
+                       arr.pop_back();
                }
                return 0;
        }
@@ -309,12 +311,12 @@ void get_sets() {
        double sum = 0;
        int outer = 0;
        for (auto it = ncombs.begin(); it != ncombs.end(); ++it) {
-               set<pair<double, short> > b;
+               set<int> 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));
+                               b.insert(i);
                        }
                }
                sum += b.size();
@@ -342,8 +344,8 @@ void print_graph(multimap<double, bitset<BSIZE> >& combs) {
 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;
+                       if (*it < i)
+                               cout << i+1 << " " << *it+1 << endl;
                }
        }
 }
@@ -351,27 +353,28 @@ void print_all() {
 int main(int argc, char* argv[]) {
        n = atoi(argv[1]);
        primes = get_primes((n*(n-1))/2);
+  test_slib();
+  char c;
+  cin >> c;
        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;
+       total_target = target_mult;
        target_mult /= n;
-       cout << "TARGET ENERGY: " << tt << endl;
+       cout << "TARGET ENERGY: " << total_target << 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);
+       //auto combs = exhaustive_search<36>(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;
@@ -384,7 +387,7 @@ int main(int argc, char* argv[]) {
                prods.push_back(s);
                if (min_score > s) min_score = s;
        }
-       cout << "FOund " << ncombs.size() << " combs" << endl;
+       cout << "Found " << ncombs.size() << " combs" << endl;
 
        cout << "GETTING SETS" << endl;
        get_sets();
@@ -394,7 +397,7 @@ int main(int argc, char* argv[]) {
 
        vector<int> arr;
        i = 0;
-       bitset<55> em;
+       bitset<36> em;
        for (auto it = ncombs.begin(); it != ncombs.end(); ++it) {
                if (i % 10 == 0) cout << i << endl;
                arr.push_back(i);