Skip to content

Instantly share code, notes, and snippets.

@natsugiri

natsugiri/SubsetSumNewton.cpp

Created August 16, 2018 17:00
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save natsugiri/d9515abee22527f72c3a5f87b1604adf to your computer and use it in GitHub Desktop.
Save natsugiri/d9515abee22527f72c3a5f87b1604adf to your computer and use it in GitHub Desktop.
Download ZIP
Raw
SubsetSumNewton.cpp
#include<bitset>
#include<chrono>
#include<cassert>
#include<stdio.h>
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<string.h>
using namespace std;
typedef long long LL;
typedef vector<int> VI;
#define REP(i,n) for(int i=0, i##_len=(n); i<i##_len; ++i)
#define EACH(i,c) for(__typeof((c).begin()) i=(c).begin(),i##_end=(c).end();i!=i##_end;++i)
template<unsigned MOD> struct ModInt {
static const unsigned static_MOD = MOD;
unsigned x;
void undef() { x = (unsigned)-1; }
bool isnan() const { return x == (unsigned)-1; }
inline int geti() const { return (int)x; }
ModInt() { x = 0; }
ModInt(const ModInt &y) { x = y.x; }
ModInt(int y) { if (y<0 || (int)MOD<=y) y %= (int)MOD; if (y<0) y += MOD; x=y; }
ModInt(unsigned y) { if (MOD<=y) x = y % MOD; else x = y; }
ModInt(long long y) { if (y<0 || MOD<=y) y %= MOD; if (y<0) y += MOD; x=y; }
ModInt(unsigned long long y) { if (MOD<=y) x = y % MOD; else x = y; }
ModInt &operator+=(const ModInt y) { if ((x += y.x) >= MOD) x -= MOD; return *this; }
ModInt &operator-=(const ModInt y) { if ((x -= y.x) & (1u<<31)) x += MOD; return *this; }
ModInt &operator*=(const ModInt y) { x = (unsigned long long)x * y.x % MOD; return *this; }
ModInt &operator/=(const ModInt y) { x = (unsigned long long)x * y.inv().x % MOD; return *this; }
ModInt operator-() const { return (x ? MOD-x: 0); }
ModInt inv() const {
unsigned a = MOD, b = x; int u = 0, v = 1;
while (b) {
int t = a / b;
a -= t * b; swap(a, b);
u -= t * v; swap(u, v);
}
if (u < 0) u += MOD;
return ModInt(u);
}
ModInt pow(long long y) const {
ModInt b = *this, r = 1;
if (y < 0) { b = b.inv(); y = -y; }
for (; y; y>>=1) {
if (y&1) r *= b;
b *= b;
}
return r;
}
friend ModInt operator+(ModInt x, const ModInt y) { return x += y; }
friend ModInt operator-(ModInt x, const ModInt y) { return x -= y; }
friend ModInt operator*(ModInt x, const ModInt y) { return x *= y; }
friend ModInt operator/(ModInt x, const ModInt y) { return x *= y.inv(); }
friend bool operator<(const ModInt x, const ModInt y) { return x.x < y.x; }
friend bool operator==(const ModInt x, const ModInt y) { return x.x == y.x; }
friend bool operator!=(const ModInt x, const ModInt y) { return x.x != y.x; }
};
template<unsigned MOD, unsigned ROOT> struct NTT {
typedef ModInt<MOD> Mint;
static void ntt(vector<Mint> &x, const int inverse=false) {
int n = x.size(); // n == 1<<k;
static vector<Mint> y, e, er;
if (n != (int)e.size()) {
y.resize(n); e.resize(n); er.resize(n);
e[0] = 1; e[1] = Mint(ROOT).pow((MOD-1)/n);
er[0] = 1; er[1] = e[1].inv();
for (int i=2; i<n; i++) e[i] = e[i-1] * e[1];
for (int i=2; i<n; i++) er[i] = er[i-1] * er[1];
}
if (inverse) swap(e, er);
int nn = n, s = 1, eo = 0;
for (; nn>1;) {
const int m = nn / 2;
for (int p=0; p<m; p++) {
Mint wp = e[p*(n/nn)];
for (int q=0; q<s; q++) {
Mint a = x[q+s*p];
Mint b = x[q+s*(p+m)];
y[q+s*(p+p)] = a + b;
y[q+s*(p+p+1)] = (a - b) * wp;
}
}
nn = m; s += s; eo ^= 1;
swap(x, y);
}
if (eo) {
for (int q=0; q<s; q++) y[q] = x[q];
swap(x, y);
}
if (inverse) {
swap(e, er);
Mint d = Mint(n).inv();
for (int i=0; i<n; i++) x[i] = x[i] * d;
}
}
static vector<unsigned> convolution(const vector<unsigned> &f, const vector<unsigned> &g) {
int sz = 1;
while (sz < (int)(f.size() + g.size())) sz *= 2;
static vector<Mint> f1, g1;
f1.resize(sz); fill(f1.begin(), f1.end(), 0);
g1.resize(sz); fill(g1.begin(), g1.end(), 0);
REP (i, f.size()) f1[i] = f[i];
REP (i, g.size()) g1[i] = g[i];
ntt(f1);
ntt(g1);
REP (i, sz) f1[i] *= g1[i];
ntt(f1, true);
vector<unsigned> ret(f.size() + g.size());
REP (i, ret.size()) ret[i] = f1[i].x;
return ret;
}
};
LL extgcd(LL x, LL y, LL&s, LL&t) {
for (LL u=t=1, v=s=0; x; ) {
LL q = y / x;
swap(s -= q * u, u);
swap(t -= q * v, v);
swap(y -= q * x, x);
}
return y;
}
LL inv_mod(LL a, LL mod) {
LL x, y;
if (extgcd(a, mod, x, y) == 1) return (x + mod) % mod;
return 0; // unsolvable
}
vector<unsigned> convolution(const vector<unsigned> &f, const vector<unsigned> &g, const unsigned mod) {
const unsigned MOD0 = 998244353, ROOT0 = 3;
const unsigned MOD1 = 2013265921, ROOT1 = 31;
const unsigned MOD2 = 2113929217, ROOT2 = 5;
if (MOD0 == mod) return NTT<MOD0, ROOT0>::convolution(f, g);
if (MOD1 == mod) return NTT<MOD1, ROOT1>::convolution(f, g);
if (MOD2 == mod) return NTT<MOD2, ROOT2>::convolution(f, g);
static vector<unsigned> vs[3];
vs[0] = NTT<MOD0, ROOT0>::convolution(f, g);
vs[1] = NTT<MOD1, ROOT1>::convolution(f, g);
vs[2] = NTT<MOD2, ROOT2>::convolution(f, g);
int sz = vs[0].size();
vector<unsigned> ret(sz);
const LL inv0 = inv_mod(MOD0, MOD1);
const LL inv1 = inv_mod((LL)MOD0 * MOD1, MOD2);
const LL modx = (LL)MOD0 * MOD1 % mod;
REP (i, sz) {
LL v1 = ((LL)vs[1][i] - vs[0][i]) * inv0 % MOD1;
if (v1 < 0) v1 += MOD1;
LL v2 = (vs[2][i] - (vs[0][i] + MOD0 * v1) % MOD2) * inv1 % MOD2;
if (v2 < 0) v2 += MOD2;
LL ans = (vs[0][i] + MOD0 * v1 + modx * v2) % mod;
if (ans < 0) ans += mod;
ret[i] = ans;
}
return ret;
}
const LL MOD = 1000000007;
typedef ModInt<MOD> Mint;
typedef vector<Mint> Poly;
namespace OPERATIONS {
Poly mul(const Poly &f, const Poly &g, int prec) {
int sf = min(prec, (int)f.size());
int sg = min(prec, (int)g.size());
vector<unsigned> uf(sf), ug(sg);
REP (i, sf) uf[i] = f[i].x;
REP (i, sg) ug[i] = g[i].x;
while (!uf.empty() && uf.back() == 0) uf.pop_back();
while (!ug.empty() && ug.back() == 0) ug.pop_back();
if (uf.empty() || ug.empty()) return Poly();
vector<unsigned> r = ::convolution(uf, ug, MOD);
Poly ret(prec);
for (int i=0; i<prec && i<(int)r.size(); i++) ret[i] = r[i];
return ret;
}
Poly inverse(const Poly &f, int prec) {
assert(f[0] == 1);
Poly g(prec, 0); g[0] = 1;
int lo = 1, hi = 2;
while (lo < prec) {
Poly e = mul(f, g, min(prec, hi));
Poly tmp = mul(g, e, min(prec, hi));
for (int i=lo; i<prec && i<(int)tmp.size(); i++) g[i] -= tmp[i];
lo = hi;
hi *= 2;
}
return g;
}
Poly log(const Poly &f, int prec) {
assert(f[0] == 1);
Poly f_inv = inverse(f, prec);
Poly df(f.size()-1u);
REP (i, df.size()) df[i] = (i+1) * f[i+1];
Poly dg = mul(df, f_inv, prec-1);
Poly g(prec);
REP (i, dg.size()) g[i+1] = dg[i] / (i+1);
return g;
}
Poly exp(const Poly &f, int prec) {
assert(f[0] == 0);
Poly g(prec); g[0] = 1;
int lo = 1, hi = 2;
while (lo < prec) {
Poly h = log(g, min(prec, hi));
REP (i, min(f.size(), h.size())) h[i] -= f[i];
h = mul(g, h, min(prec, hi));
for (int i=lo; i<prec && i<(int)h.size(); i++) g[i] -= h[i];
lo = hi;
hi *= 2;
}
return g;
}
}; // namespace OPERATIONS;
vector<Mint> subsetsum_fft(VI s, int t) {
sort(s.begin(), s.end());
Poly b(t+1);
for (int i=0, j; i<(int)s.size(); i=j) {
j = i+1;
while (j < (int)s.size() && s[i] == s[j]) j++;
for (int x=1; x<=t/s[i]; x++) {
b[x*s[i]] += (x&1? Mint(j-i): Mint(i-j)) / x;
}
}
Poly a = OPERATIONS::exp(b, t+1);
return a;
}
void MAIN() {
const int SIZE = 100000;
int n = SIZE;
int t = SIZE;
VI s(n);
REP (i, n) s[i] = i+1;
auto start = chrono::steady_clock::now();
vector<Mint> a = subsetsum_fft(s, t);
auto diff = chrono::steady_clock::now() - start;
printf("time\t%.6f\n", diff.count() * 1e-9);
printf("answer\t%d\n", a[t].geti());
}
int main() {
MAIN();
return 0;
}
Raw
SubsetSumRecursive.cpp
#include<bitset>
#include<chrono>
#include<cassert>
#include<stdio.h>
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<string.h>
using namespace std;
typedef long long LL;
typedef vector<int> VI;
#define REP(i,n) for(int i=0, i##_len=(n); i<i##_len; ++i)
#define EACH(i,c) for(__typeof((c).begin()) i=(c).begin(),i##_end=(c).end();i!=i##_end;++i)
template<unsigned MOD> struct ModInt {
static const unsigned static_MOD = MOD;
unsigned x;
void undef() { x = (unsigned)-1; }
bool isnan() const { return x == (unsigned)-1; }
inline int geti() const { return (int)x; }
ModInt() { x = 0; }
ModInt(const ModInt &y) { x = y.x; }
ModInt(int y) { if (y<0 || (int)MOD<=y) y %= (int)MOD; if (y<0) y += MOD; x=y; }
ModInt(unsigned y) { if (MOD<=y) x = y % MOD; else x = y; }
ModInt(long long y) { if (y<0 || MOD<=y) y %= MOD; if (y<0) y += MOD; x=y; }
ModInt(unsigned long long y) { if (MOD<=y) x = y % MOD; else x = y; }
ModInt &operator+=(const ModInt y) { if ((x += y.x) >= MOD) x -= MOD; return *this; }
ModInt &operator-=(const ModInt y) { if ((x -= y.x) & (1u<<31)) x += MOD; return *this; }
ModInt &operator*=(const ModInt y) { x = (unsigned long long)x * y.x % MOD; return *this; }
ModInt &operator/=(const ModInt y) { x = (unsigned long long)x * y.inv().x % MOD; return *this; }
ModInt operator-() const { return (x ? MOD-x: 0); }
ModInt inv() const {
unsigned a = MOD, b = x; int u = 0, v = 1;
while (b) {
int t = a / b;
a -= t * b; swap(a, b);
u -= t * v; swap(u, v);
}
if (u < 0) u += MOD;
return ModInt(u);
}
ModInt pow(long long y) const {
ModInt b = *this, r = 1;
if (y < 0) { b = b.inv(); y = -y; }
for (; y; y>>=1) {
if (y&1) r *= b;
b *= b;
}
return r;
}
friend ModInt operator+(ModInt x, const ModInt y) { return x += y; }
friend ModInt operator-(ModInt x, const ModInt y) { return x -= y; }
friend ModInt operator*(ModInt x, const ModInt y) { return x *= y; }
friend ModInt operator/(ModInt x, const ModInt y) { return x *= y.inv(); }
friend bool operator<(const ModInt x, const ModInt y) { return x.x < y.x; }
friend bool operator==(const ModInt x, const ModInt y) { return x.x == y.x; }
friend bool operator!=(const ModInt x, const ModInt y) { return x.x != y.x; }
};
template<unsigned MOD, unsigned ROOT> struct NTT {
typedef ModInt<MOD> Mint;
static void ntt(vector<Mint> &x, const int inverse=false) {
int n = x.size(); // n == 1<<k;
static vector<Mint> y, e, er;
if (n != (int)e.size()) {
y.resize(n); e.resize(n); er.resize(n);
e[0] = 1; e[1] = Mint(ROOT).pow((MOD-1)/n);
er[0] = 1; er[1] = e[1].inv();
for (int i=2; i<n; i++) e[i] = e[i-1] * e[1];
for (int i=2; i<n; i++) er[i] = er[i-1] * er[1];
}
if (inverse) swap(e, er);
int nn = n, s = 1, eo = 0;
for (; nn>1;) {
const int m = nn / 2;
for (int p=0; p<m; p++) {
Mint wp = e[p*(n/nn)];
for (int q=0; q<s; q++) {
Mint a = x[q+s*p];
Mint b = x[q+s*(p+m)];
y[q+s*(p+p)] = a + b;
y[q+s*(p+p+1)] = (a - b) * wp;
}
}
nn = m; s += s; eo ^= 1;
swap(x, y);
}
if (eo) {
for (int q=0; q<s; q++) y[q] = x[q];
swap(x, y);
}
if (inverse) {
swap(e, er);
Mint d = Mint(n).inv();
for (int i=0; i<n; i++) x[i] = x[i] * d;
}
}
static vector<unsigned> convolution(const vector<unsigned> &f, const vector<unsigned> &g) {
int sz = 1;
while (sz < (int)(f.size() + g.size())) sz *= 2;
static vector<Mint> f1, g1;
f1.resize(sz); fill(f1.begin(), f1.end(), 0);
g1.resize(sz); fill(g1.begin(), g1.end(), 0);
REP (i, f.size()) f1[i] = f[i];
REP (i, g.size()) g1[i] = g[i];
ntt(f1);
ntt(g1);
REP (i, sz) f1[i] *= g1[i];
ntt(f1, true);
vector<unsigned> ret(f.size() + g.size());
REP (i, ret.size()) ret[i] = f1[i].x;
return ret;
}
};
LL extgcd(LL x, LL y, LL&s, LL&t) {
for (LL u=t=1, v=s=0; x; ) {
LL q = y / x;
swap(s -= q * u, u);
swap(t -= q * v, v);
swap(y -= q * x, x);
}
return y;
}
LL inv_mod(LL a, LL mod) {
LL x, y;
if (extgcd(a, mod, x, y) == 1) return (x + mod) % mod;
return 0; // unsolvable
}
vector<unsigned> convolution(const vector<unsigned> &f, const vector<unsigned> &g, const unsigned mod) {
const unsigned MOD0 = 998244353, ROOT0 = 3;
const unsigned MOD1 = 2013265921, ROOT1 = 31;
const unsigned MOD2 = 2113929217, ROOT2 = 5;
if (MOD0 == mod) return NTT<MOD0, ROOT0>::convolution(f, g);
if (MOD1 == mod) return NTT<MOD1, ROOT1>::convolution(f, g);
if (MOD2 == mod) return NTT<MOD2, ROOT2>::convolution(f, g);
static vector<unsigned> vs[3];
vs[0] = NTT<MOD0, ROOT0>::convolution(f, g);
vs[1] = NTT<MOD1, ROOT1>::convolution(f, g);
vs[2] = NTT<MOD2, ROOT2>::convolution(f, g);
int sz = vs[0].size();
vector<unsigned> ret(sz);
const LL inv0 = inv_mod(MOD0, MOD1);
const LL inv1 = inv_mod((LL)MOD0 * MOD1, MOD2);
const LL modx = (LL)MOD0 * MOD1 % mod;
REP (i, sz) {
LL v1 = ((LL)vs[1][i] - vs[0][i]) * inv0 % MOD1;
if (v1 < 0) v1 += MOD1;
LL v2 = (vs[2][i] - (vs[0][i] + MOD0 * v1) % MOD2) * inv1 % MOD2;
if (v2 < 0) v2 += MOD2;
LL ans = (vs[0][i] + MOD0 * v1 + modx * v2) % mod;
if (ans < 0) ans += mod;
ret[i] = ans;
}
return ret;
}
const LL MOD = 1000000007;
typedef ModInt<MOD> Mint;
typedef vector<Mint> Poly;
namespace OPERATIONS_2 {
typedef typename Poly::iterator Iter;
Poly convolution(Iter f_begin, Iter f_end, Iter g_begin, Iter g_end) {
vector<unsigned> f, g;
f.reserve(f_end - f_begin);
for (auto it=f_begin; it!=f_end; it++) f.push_back(it->x);
g.reserve(g_end - g_begin);
for (auto it=g_begin; it!=g_end; it++) g.push_back(it->x);
vector<unsigned> h = ::convolution(f, g, MOD);
return Poly(h.begin(), h.end());
}
void compute(Poly &f, Poly &g, int l, int r) {
if (l < r) {
int m = (l+r)/2;
compute(f, g, l, m);
// fft
// for (int i=m+1; i<=r; i++) {
// for (int j=l; j<=m; j++) {
// g[i] += f[i-j] * g[j] / i;
// }
// }
Poly h = convolution(f.begin(), f.begin()+r-l+1, g.begin()+l, g.begin()+m+1);
for (int i=m+1; i<=r; i++) {
g[i] += h[i-l]/i;
}
compute(f, g, m+1, r);
}
}
Poly exp(Poly f) {
int n = f.size();
REP (i, n) f[i] *= i;
Poly g(n);
g[0] = 1;
compute(f, g, 0, n-1);
return g;
}
}; // namespace OPERATIONS_2;
vector<Mint> subsetsum_fft(VI s, int t) {
sort(s.begin(), s.end());
Poly b(t+1);
for (int i=0, j; i<(int)s.size(); i=j) {
j = i+1;
while (j < (int)s.size() && s[i] == s[j]) j++;
for (int x=1; x<=t/s[i]; x++) {
b[x*s[i]] += (x&1? Mint(j-i): Mint(i-j)) / x;
}
}
Poly a = OPERATIONS_2::exp(b);
return a;
}
void MAIN() {
const int SIZE = 100000;
int n = SIZE;
int t = SIZE;
VI s(n);
REP (i, n) s[i] = i+1;
auto start = chrono::steady_clock::now();
vector<Mint> a = subsetsum_fft(s, t);
auto diff = chrono::steady_clock::now() - start;
printf("time\t%.6f\n", diff.count() * 1e-9);
printf("answer\t%d\n", a[t].geti());
}
int main() {
MAIN();
return 0;
}
@TheSQLGuru
Copy link

I am sorry to trouble you about this old code you wrote, but I was wondering if either or both methods could be adapted to actually output the solution(s) found. I have become fascinated by the SSP, especially programatic code that can be run on modern CPU/GPU to solve large problems in this space. Unfortunately for me, all of the maths I took for my physics degree were almost 4 decades ago now - and even in my prime I could not decipher the maths used in most of the best algorithms for the SSP problem. Any help you could provide would be greatly appreciated, and I would be happy to provide payment for your efforts!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment