diff options
Diffstat (limited to 'content/math')
63 files changed, 2167 insertions, 0 deletions
diff --git a/content/math/berlekampMassey.cpp b/content/math/berlekampMassey.cpp new file mode 100644 index 0000000..29e084f --- /dev/null +++ b/content/math/berlekampMassey.cpp @@ -0,0 +1,31 @@ +constexpr ll mod = 1'000'000'007; +vector<ll> BerlekampMassey(const vector<ll>& s) { + int n = sz(s), L = 0, m = 0; + vector<ll> C(n), B(n), T; + C[0] = B[0] = 1; + + ll b = 1; + for (int i = 0; i < n; i++) { + m++; + ll d = s[i] % mod; + for (int j = 1; j <= L; j++) { + d = (d + C[j] * s[i - j]) % mod; + } + if (!d) continue; + T = C; + ll coef = d * powMod(b, mod-2, mod) % mod; + for (int j = m; j < n; j++) { + C[j] = (C[j] - coef * B[j - m]) % mod; + } + if (2 * L > i) continue; + L = i + 1 - L; + swap(B, T); + b = d; + m = 0; + } + + C.resize(L + 1); + C.erase(C.begin()); + for (auto& x : C) x = (mod - x) % mod; + return C; +} diff --git a/content/math/bigint.cpp b/content/math/bigint.cpp new file mode 100644 index 0000000..1b3b953 --- /dev/null +++ b/content/math/bigint.cpp @@ -0,0 +1,271 @@ +// base and base_digits must be consistent +constexpr ll base = 1'000'000; +constexpr ll base_digits = 6; +struct bigint { + using vll = vector<ll>; + vll a; ll sign; + + bigint() : sign(1) {} + + bigint(ll v) {*this = v;} + + bigint(const string &s) {read(s);} + + void operator=(ll v) { + sign = 1; + if (v < 0) sign = -1, v = -v; + a.clear(); + for (; v > 0; v = v / base) + a.push_back(v % base); + } + + bigint operator+(const bigint& v) const { + if (sign == v.sign) { + bigint res = v; + for (ll i = 0, carry = 0; i < max(sz(a), sz(v.a)) || carry; ++i) { + if (i == sz(res.a)) + res.a.push_back(0); + res.a[i] += carry + (i < sz(a) ? a[i] : 0); + carry = res.a[i] >= base; + if (carry) + res.a[i] -= base; + } + return res; + } + return *this - (-v); + } + + bigint operator-(const bigint& v) const { + if (sign == v.sign) { + if (abs() >= v.abs()) { + bigint res = *this; + for (ll i = 0, carry = 0; i < sz(v.a) || carry; ++i) { + res.a[i] -= carry + (i < sz(v.a) ? v.a[i] : 0); + carry = res.a[i] < 0; + if (carry) res.a[i] += base; + } + res.trim(); + return res; + } + return -(v - *this); + } + return *this + (-v); + } + + void operator*=(ll v) { + if (v < 0) sign = -sign, v = -v; + for (ll i = 0, carry = 0; i < sz(a) || carry; ++i) { + if (i == sz(a)) a.push_back(0); + ll cur = a[i] * v + carry; + carry = cur / base; + a[i] = cur % base; + } + trim(); + } + + bigint operator*(ll v) const { + bigint res = *this; + res *= v; + return res; + } + + friend pair<bigint, bigint> divmod(const bigint& a1, const bigint& b1) { + ll norm = base / (b1.a.back() + 1); + bigint a = a1.abs() * norm; + bigint b = b1.abs() * norm; + bigint q, r; + q.a.resize(sz(a.a)); + for (ll i = sz(a.a) - 1; i >= 0; i--) { + r *= base; + r += a.a[i]; + ll s1 = sz(r.a) <= sz(b.a) ? 0 : r.a[sz(b.a)]; + ll s2 = sz(r.a) <= sz(b.a) - 1 ? 0 : r.a[sz(b.a) - 1]; + ll d = (base * s1 + s2) / b.a.back(); + r -= b * d; + while (r < 0) r += b, --d; + q.a[i] = d; + } + q.sign = a1.sign * b1.sign; + r.sign = a1.sign; + q.trim(); + r.trim(); + return make_pair(q, r / norm); + } + + bigint operator/(const bigint& v) const { + return divmod(*this, v).first; + } + + bigint operator%(const bigint& v) const { + return divmod(*this, v).second; + } + + void operator/=(ll v) { + if (v < 0) sign = -sign, v = -v; + for (ll i = sz(a) - 1, rem = 0; i >= 0; --i) { + ll cur = a[i] + rem * base; + a[i] = cur / v; + rem = cur % v; + } + trim(); + } + + bigint operator/(ll v) const { + bigint res = *this; + res /= v; + return res; + } + + ll operator%(ll v) const { + if (v < 0) v = -v; + ll m = 0; + for (ll i = sz(a) - 1; i >= 0; --i) + m = (a[i] + m * base) % v; + return m * sign; + } + + void operator+=(const bigint& v) { + *this = *this + v; + } + void operator-=(const bigint& v) { + *this = *this - v; + } + void operator*=(const bigint& v) { + *this = *this * v; + } + void operator/=(const bigint& v) { + *this = *this / v; + } + + bool operator<(const bigint& v) const { + if (sign != v.sign) return sign < v.sign; + if (sz(a) != sz(v.a)) + return sz(a) * sign < sz(v.a) * v.sign; + for (ll i = sz(a) - 1; i >= 0; i--) + if (a[i] != v.a[i]) + return a[i] * sign < v.a[i] * sign; + return false; + } + + bool operator>(const bigint& v) const { + return v < *this; + } + bool operator<=(const bigint& v) const { + return !(v < *this); + } + bool operator>=(const bigint& v) const { + return !(*this < v); + } + bool operator==(const bigint& v) const { + return !(*this < v) && !(v < *this); + } + bool operator!=(const bigint& v) const { + return *this < v || v < *this; + } + + void trim() { + while (!a.empty() && !a.back()) a.pop_back(); + if (a.empty()) sign = 1; + } + + bool isZero() const { + return a.empty() || (sz(a) == 1 && a[0] == 0); + } + + bigint operator-() const { + bigint res = *this; + res.sign = -sign; + return res; + } + + bigint abs() const { + bigint res = *this; + res.sign *= res.sign; + return res; + } + + ll longValue() const { + ll res = 0; + for (ll i = sz(a) - 1; i >= 0; i--) + res = res * base + a[i]; + return res * sign; + } + + void read(const string& s) { + sign = 1; + a.clear(); + ll pos = 0; + while (pos < sz(s) && (s[pos] == '-' || s[pos] == '+')) { + if (s[pos] == '-') sign = -sign; + ++pos; + } + for (ll i = sz(s) - 1; i >= pos; i -= base_digits) { + ll x = 0; + for (ll j = max(pos, i - base_digits + 1); j <= i; j++) + x = x * 10 + s[j] - '0'; + a.push_back(x); + } + trim(); + } + + friend istream& operator>>(istream& stream, bigint& v) { + string s; + stream >> s; + v.read(s); + return stream; + } + + friend ostream& operator<<(ostream& stream, const bigint& v) { + if (v.sign == -1) stream << '-'; + stream << (v.a.empty() ? 0 : v.a.back()); + for (ll i = sz(v.a) - 2; i >= 0; --i) + stream << setw(base_digits) << setfill('0') << v.a[i]; + return stream; + } + + static vll karatsubaMultiply(const vll& a, const vll& b) { + ll n = sz(a); + vll res(n + n); + if (n <= 32) { + for (ll i = 0; i < n; i++) + for (ll j = 0; j < n; j++) + res[i + j] += a[i] * b[j]; + return res; + } + ll k = n >> 1; + vll a1(a.begin(), a.begin() + k); + vll a2(a.begin() + k, a.end()); + vll b1(b.begin(), b.begin() + k); + vll b2(b.begin() + k, b.end()); + vll a1b1 = karatsubaMultiply(a1, b1); + vll a2b2 = karatsubaMultiply(a2, b2); + for (ll i = 0; i < k; i++) a2[i] += a1[i]; + for (ll i = 0; i < k; i++) b2[i] += b1[i]; + vll r = karatsubaMultiply(a2, b2); + for (ll i = 0; i < sz(a1b1); i++) r[i] -= a1b1[i]; + for (ll i = 0; i < sz(a2b2); i++) r[i] -= a2b2[i]; + for (ll i = 0; i < sz(r); i++) res[i + k] += r[i]; + for (ll i = 0; i < sz(a1b1); i++) res[i] += a1b1[i]; + for (ll i = 0; i < sz(a2b2); i++) res[i + n] += a2b2[i]; + return res; + } + + bigint operator*(const bigint& v) const { + vll ta(a.begin(), a.end()); + vll va(v.a.begin(), v.a.end()); + while (sz(ta) < sz(va)) ta.push_back(0); + while (sz(va) < sz(ta)) va.push_back(0); + while (sz(ta) & (sz(ta) - 1)) + ta.push_back(0), va.push_back(0); + vll ra = karatsubaMultiply(ta, va); + bigint res; + res.sign = sign * v.sign; + for (ll i = 0, carry = 0; i < sz(ra); i++) { + ll cur = ra[i] + carry; + res.a.push_back(cur % base); + carry = cur / base; + } + res.trim(); + return res; + } +}; diff --git a/content/math/binomial0.cpp b/content/math/binomial0.cpp new file mode 100644 index 0000000..f37aea5 --- /dev/null +++ b/content/math/binomial0.cpp @@ -0,0 +1,14 @@ +constexpr ll lim = 10'000'000; +ll fac[lim], inv[lim]; + +void precalc() { + fac[0] = inv[0] = 1; + for (int i = 1; i < lim; i++) fac[i] = fac[i-1] * i % mod; + inv[lim - 1] = multInv(fac[lim - 1], mod); + for (int i = lim - 1; i > 0; i--) inv[i-1] = inv[i] * i % mod; +} + +ll calc_binom(ll n, ll k) { + if (n < 0 || n < k || k < 0) return 0; + return (fac[n] * inv[n-k] % mod) * inv[k] % mod; +} diff --git a/content/math/binomial1.cpp b/content/math/binomial1.cpp new file mode 100644 index 0000000..dab20b3 --- /dev/null +++ b/content/math/binomial1.cpp @@ -0,0 +1,8 @@ +ll calc_binom(ll n, ll k) { + if (k > n) return 0; + ll r = 1; + for (ll d = 1; d <= k; d++) {// Reihenfolge => Teilbarkeit + r *= n--, r /= d; + } + return r; +} diff --git a/content/math/binomial2.cpp b/content/math/binomial2.cpp new file mode 100644 index 0000000..4531505 --- /dev/null +++ b/content/math/binomial2.cpp @@ -0,0 +1,32 @@ +constexpr ll mod = 1'000'000'009; + +ll binomPPow(ll n, ll k, ll p) { + ll res = 1; + if (p > n) { + } else if (p > n - k || (p * p > n && n % p < k % p)) { + res *= p; + res %= mod; + } else if (p * p <= n) { + ll c = 0, tmpN = n, tmpK = k; + while (tmpN > 0) { + if (tmpN % p < tmpK % p + c) { + res *= p; + res %= mod; + c = 1; + } else c = 0; + tmpN /= p; + tmpK /= p; + }} + return res; +} + +ll calc_binom(ll n, ll k) { + if (k > n) return 0; + ll res = 1; + k = min(k, n - k); + for (ll i = 0; primes[i] <= n; i++) { + res *= binomPPow(n, k, primes[i]); + res %= mod; + } + return res; +} diff --git a/content/math/binomial3.cpp b/content/math/binomial3.cpp new file mode 100644 index 0000000..7a6ab4e --- /dev/null +++ b/content/math/binomial3.cpp @@ -0,0 +1,10 @@ +ll calc_binom(ll n, ll k, ll p) { + assert(n < p); //wichtig: sonst falsch! + if (k > n) return 0; + ll x = k % 2 != 0 ? p-1 : 1; + for (ll c = p-1; c > n; c--) { + x *= c - k; x %= p; + x *= multInv(c, p); x %= p; + } + return x; +} diff --git a/content/math/chineseRemainder.cpp b/content/math/chineseRemainder.cpp new file mode 100644 index 0000000..ccbc5dc --- /dev/null +++ b/content/math/chineseRemainder.cpp @@ -0,0 +1,14 @@ +struct CRT { + using lll = __int128; + lll M = 1, sol = 0; // Solution unique modulo M + bool hasSol = true; + + // Adds congruence x = a (mod m) + void add(ll a, ll m) { + auto [d, s, t] = extendedEuclid(M, m); + if((a - sol) % d != 0) hasSol = false; + lll z = M/d * s; + M *= m/d; + sol = (z % M * (a-sol) % M + sol + M) % M; + } +}; diff --git a/content/math/cycleDetection.cpp b/content/math/cycleDetection.cpp new file mode 100644 index 0000000..5e68c0c --- /dev/null +++ b/content/math/cycleDetection.cpp @@ -0,0 +1,18 @@ +pair<ll, ll> cycleDetection(ll x0, function<ll(ll)> f) { + ll a = x0, b = f(x0), length = 1; + for (ll power = 1; a != b; b = f(b), length++) { + if (power == length) { + power *= 2; + length = 0; + a = b; + }} + ll start = 0; + a = x0; b = x0; + for (ll i = 0; i < length; i++) b = f(b); + while (a != b) { + a = f(a); + b = f(b); + start++; + } + return {start, length}; +} diff --git a/content/math/discreteLogarithm.cpp b/content/math/discreteLogarithm.cpp new file mode 100644 index 0000000..68866e0 --- /dev/null +++ b/content/math/discreteLogarithm.cpp @@ -0,0 +1,17 @@ +ll dlog(ll a, ll b, ll m) { //a > 0! + ll bound = sqrtl(m) + 1; //memory usage bound < p + vector<pair<ll, ll>> vals(bound); + for (ll i = 0, e = 1; i < bound; i++, e = (e * a) % m) { + vals[i] = {e, i}; + } + vals.emplace_back(m, 0); + sort(all(vals)); + ll fact = powMod(a, m - bound - 1, m); + + for (ll i = 0; i < m; i += bound, b = (b * fact) % m) { + auto it = lower_bound(all(vals), pair<ll, ll>{b, 0}); + if (it->first == b) { + return (i + it->second) % m; + }} + return -1; +} diff --git a/content/math/discreteNthRoot.cpp b/content/math/discreteNthRoot.cpp new file mode 100644 index 0000000..403cb3b --- /dev/null +++ b/content/math/discreteNthRoot.cpp @@ -0,0 +1,5 @@ +ll root(ll a, ll b, ll m) { // a > 0! + ll g = findPrimitive(m); + ll c = dlog(powMod(g, a, m), b, m); + return c < 0 ? -1 : powMod(g, c, m); +} diff --git a/content/math/divisors.cpp b/content/math/divisors.cpp new file mode 100644 index 0000000..5afd4fb --- /dev/null +++ b/content/math/divisors.cpp @@ -0,0 +1,11 @@ +ll countDivisors(ll n) { + ll res = 1; + for (ll i = 2; i * i * i <= n; i++) { + ll c = 0; + while (n % i == 0) {n /= i; c++;} + res *= c + 1; + } + if (isPrime(n)) res *= 2; + else if (n > 1) res *= isSquare(n) ? 3 : 4; + return res; +} diff --git a/content/math/extendedEuclid.cpp b/content/math/extendedEuclid.cpp new file mode 100644 index 0000000..ecf4a16 --- /dev/null +++ b/content/math/extendedEuclid.cpp @@ -0,0 +1,6 @@ +// a*x + b*y = ggt(a, b) +array<ll, 3> extendedEuclid(ll a, ll b) { + if (a == 0) return {b, 0, 1}; + auto [d, x, y] = extendedEuclid(b % a, a); + return {d, y - (b / a) * x, x}; +} diff --git a/content/math/gauss.cpp b/content/math/gauss.cpp new file mode 100644 index 0000000..8129fd2 --- /dev/null +++ b/content/math/gauss.cpp @@ -0,0 +1,36 @@ +void normalLine(int line) { + double factor = mat[line][line]; + for (double& x : mat[line]) x /= factor; +} + +void takeAll(int n, int line) { + for (int i = 0; i < n; i++) { + if (i == line) continue; + double diff = mat[i][line]; + for (int j = 0; j < sz(mat[i]); j++) { + mat[i][j] -= diff * mat[line][j]; +}}} + +int gauss(int n) { + vector<bool> done(n, false); + for (int i = 0; i < n; i++) { + int swappee = i; // Sucht Pivotzeile für bessere Stabilität. + for (int j = 0; j < n; j++) { + if (done[j]) continue; + if (abs(mat[j][i]) > abs(mat[i][i])) swappee = j; + } + swap(mat[i], mat[swappee]); + if (abs(mat[i][i]) > EPS) { + normalLine(i); + takeAll(n, i); + done[i] = true; + }} + // Ab jetzt nur checks bzgl. Eindeutigkeit/Existenz der Lösung. + for (int i = 0; i < n; i++) { + bool allZero = true; + for (int j = i; j < n; j++) allZero &= abs(mat[i][j]) <= EPS; + if (allZero && abs(mat[i][n]) > EPS) return INCONSISTENT; + if (allZero && abs(mat[i][n]) <= EPS) return MULTIPLE; + } + return UNIQUE; +} diff --git a/content/math/gcd-lcm.cpp b/content/math/gcd-lcm.cpp new file mode 100644 index 0000000..a1c63c8 --- /dev/null +++ b/content/math/gcd-lcm.cpp @@ -0,0 +1,2 @@ +ll gcd(ll a, ll b) {return b == 0 ? a : gcd(b, a % b);} +ll lcm(ll a, ll b) {return a * (b / gcd(a, b));} diff --git a/content/math/goldenSectionSearch.cpp b/content/math/goldenSectionSearch.cpp new file mode 100644 index 0000000..28ee4c3 --- /dev/null +++ b/content/math/goldenSectionSearch.cpp @@ -0,0 +1,15 @@ +template<typename F> +ld gss(ld l, ld r, F&& f) { + ld inv = (sqrt(5.0l) - 1) / 2; + ld x1 = r - inv*(r-l), x2 = l + inv*(r-l); + ld f1 = f(x1), f2 = f(x2); + for (int i = 0; i < 200; i++) { + if (f1 < f2) { //change to > to find maximum + r = x2; x2 = x1; f2 = f1; + x1 = r - inv*(r-l); f1 = f(x1); + } else { + l = x1; x1 = x2; f1 = f2; + x2 = l + inv*(r-l); f2 = f(x2); + }} + return l; +} diff --git a/content/math/inversions.cpp b/content/math/inversions.cpp new file mode 100644 index 0000000..9e47f9b --- /dev/null +++ b/content/math/inversions.cpp @@ -0,0 +1,9 @@ +ll inversions(const vector<ll>& v) { + Tree<pair<ll, ll>> t; //ordered statistics tree @\sourceref{datastructures/pbds.cpp}@ + ll res = 0; + for (ll i = 0; i < sz(v); i++) { + res += i - t.order_of_key({v[i], i}); + t.insert({v[i], i}); + } + return res; +} diff --git a/content/math/inversionsMerge.cpp b/content/math/inversionsMerge.cpp new file mode 100644 index 0000000..8235b11 --- /dev/null +++ b/content/math/inversionsMerge.cpp @@ -0,0 +1,27 @@ +// Laufzeit: O(n*log(n)) +ll merge(vector<ll>& v, vector<ll>& left, vector<ll>& right) { + int a = 0, b = 0, i = 0; + ll inv = 0; + while (a < sz(left) && b < sz(right)) { + if (left[a] < right[b]) v[i++] = left[a++]; + else { + inv += sz(left) - a; + v[i++] = right[b++]; + } + } + while (a < sz(left)) v[i++] = left[a++]; + while (b < sz(right)) v[i++] = right[b++]; + return inv; +} + +ll mergeSort(vector<ll> &v) { // Sortiert v und gibt Inversionszahl zurück. + int n = sz(v); + vector<ll> left(n / 2), right((n + 1) / 2); + for (int i = 0; i < n / 2; i++) left[i] = v[i]; + for (int i = n / 2; i < n; i++) right[i - n / 2] = v[i]; + + ll result = 0; + if (sz(left) > 1) result += mergeSort(left); + if (sz(right) > 1) result += mergeSort(right); + return result + merge(v, left, right); +} diff --git a/content/math/kthperm.cpp b/content/math/kthperm.cpp new file mode 100644 index 0000000..504f09c --- /dev/null +++ b/content/math/kthperm.cpp @@ -0,0 +1,14 @@ +vector<ll> kthperm(ll n, ll k) { + Tree<ll> t; + vector<ll> res(n); + for (ll i = 1; i <= n; k /= i, i++) { + t.insert(i - 1); + res[n - i] = k % i; + } + for (ll& x : res) { + auto it = t.find_by_order(x); + x = *it; + t.erase(it); + } + return res; +} diff --git a/content/math/legendre.cpp b/content/math/legendre.cpp new file mode 100644 index 0000000..b85ea2a --- /dev/null +++ b/content/math/legendre.cpp @@ -0,0 +1,4 @@ +ll legendre(ll a, ll p) { // p prim >= 2 + ll s = powMod(a, p / 2, p); + return s < 2 ? s : -1ll; +} diff --git a/content/math/lgsFp.cpp b/content/math/lgsFp.cpp new file mode 100644 index 0000000..0241742 --- /dev/null +++ b/content/math/lgsFp.cpp @@ -0,0 +1,26 @@ +void normalLine(int line, ll p) { + ll factor = multInv(mat[line][line], p); + for (ll& x : mat[line]) x = (x * factor) % p; +} + +void takeAll(int n, int line, ll p) { + for (int i = 0; i < n; i++) { + if (i == line) continue; + ll diff = mat[i][line]; + for (int j = 0; j < sz(mat[i]); j++) { + mat[i][j] -= (diff * mat[line][j]) % p; + mat[i][j] = (mat[i][j] + p) % p; +}}} + +void gauss(int n, ll mod) { + vector<bool> done(n, false); + for (int i = 0; i < n; i++) { + int j = 0; + while (j < n && (done[j] || mat[j][i] == 0)) j++; + if (j == n) continue; + swap(mat[i], mat[j]); + normalLine(i, mod); + takeAll(n, i, mod); + done[i] = true; +}} +// für Eindeutigkeit, Existenz etc. siehe LGS über R @\sourceref{math/gauss.cpp}@ diff --git a/content/math/linearCongruence.cpp b/content/math/linearCongruence.cpp new file mode 100644 index 0000000..cdb5a37 --- /dev/null +++ b/content/math/linearCongruence.cpp @@ -0,0 +1,5 @@ +ll solveLinearCongruence(ll a, ll b, ll m) { + ll g = gcd(a, m); + if (b % g != 0) return -1; + return ((b / g) * multInv(a / g, m / g)) % (m / g); +} diff --git a/content/math/linearRecurrence.cpp b/content/math/linearRecurrence.cpp new file mode 100644 index 0000000..2501e64 --- /dev/null +++ b/content/math/linearRecurrence.cpp @@ -0,0 +1,33 @@ +constexpr ll mod = 1'000'000'007; +vector<ll> modMul(const vector<ll>& a, const vector<ll>& b, + const vector<ll>& c) { + ll n = sz(c); + vector<ll> res(n * 2 + 1); + for (int i = 0; i <= n; i++) { //a*b + for (int j = 0; j <= n; j++) { + res[i + j] += a[i] * b[j]; + res[i + j] %= mod; + }} + for (int i = 2 * n; i > n; i--) { //res%c + for (int j = 0; j < n; j++) { + res[i - 1 - j] += res[i] * c[j]; + res[i - 1 - j] %= mod; + }} + res.resize(n + 1); + return res; +} + +ll kthTerm(const vector<ll>& f, const vector<ll>& c, ll k) { + assert(sz(f) == sz(c)); + vector<ll> tmp(sz(c) + 1), a(sz(c) + 1); + tmp[0] = a[1] = 1; //tmp = (x^k) % c + + for (k++; k > 0; k /= 2) { + if (k & 1) tmp = modMul(tmp, a, c); + a = modMul(a, a, c); + } + + ll res = 0; + for (int i = 0; i < sz(c); i++) res += (tmp[i+1] * f[i]) % mod; + return res % mod; +} diff --git a/content/math/linearSieve.cpp b/content/math/linearSieve.cpp new file mode 100644 index 0000000..64440dd --- /dev/null +++ b/content/math/linearSieve.cpp @@ -0,0 +1,50 @@ +constexpr ll N = 10'000'000; +ll small[N], power[N], sieved[N]; +vector<ll> primes; + +//wird aufgerufen mit (p^k, p, k) für prime p und k > 0 +ll mu(ll pk, ll p, ll k) {return -(k == 1);} +ll phi(ll pk, ll p, ll k) {return pk - pk / p;} +ll div(ll pk, ll p, ll k) {return k+1;} +ll divSum(ll pk, ll p, ll k) {return (pk*p-1) / (p - 1);} +ll square(ll pk, ll p, ll k) {return k % 2 ? pk / p : pk;} +ll squareFree(ll pk, ll p, ll k) {return p;} + +void sieve() { // O(N) + small[1] = power[1] = sieved[1] = 1; + for (ll i = 2; i < N; i++) { + if (small[i] == 0) { + primes.push_back(i); + for (ll pk = i, k = 1; pk < N; pk *= i, k++) { + small[pk] = i; + power[pk] = pk; + sieved[pk] = mu(pk, i, k); // Aufruf ändern! + }} + for (ll j=0; i*primes[j] < N && primes[j] < small[i]; j++) { + ll k = i * primes[j]; + small[k] = power[k] = primes[j]; + sieved[k] = sieved[i] * sieved[primes[j]]; + } + if (i * small[i] < N && power[i] != i) { + ll k = i * small[i]; + small[k] = small[i]; + power[k] = power[i] * small[i]; + sieved[k] = sieved[power[k]] * sieved[k / power[k]]; +}}} + +ll naive(ll n) { // O(sqrt(n)) + ll res = 1; + for (ll p = 2; p * p <= n; p++) { + if (n % p == 0) { + ll pk = 1; + ll k = 0; + do { + n /= p; + pk *= p; + k++; + } while (n % p == 0); + res *= mu(pk, p, k); // Aufruf ändern! + }} + if (n > 1) res *= mu(n, n, 1); + return res; +} diff --git a/content/math/longestIncreasingSubsequence.cpp b/content/math/longestIncreasingSubsequence.cpp new file mode 100644 index 0000000..fcb63b4 --- /dev/null +++ b/content/math/longestIncreasingSubsequence.cpp @@ -0,0 +1,17 @@ +vector<int> lis(vector<ll>& a) { + int n = sz(a), len = 0; + vector<ll> dp(n, INF), dp_id(n), prev(n); + for (int i = 0; i < n; i++) { + int pos = lower_bound(all(dp), a[i]) - dp.begin(); + dp[pos] = a[i]; + dp_id[pos] = i; + prev[i] = pos ? dp_id[pos - 1] : -1; + len = max(len, pos + 1); + } + // reconstruction + vector<int> res(len); + for (int x = dp_id[len-1]; len--; x = prev[x]) { + res[len] = x; + } + return res; // indices of one LIS +} diff --git a/content/math/math.tex b/content/math/math.tex new file mode 100644 index 0000000..f670a70 --- /dev/null +++ b/content/math/math.tex @@ -0,0 +1,525 @@ +\section{Mathe}
+
+\begin{algorithm}{Longest Increasing Subsequence}
+ \begin{itemize}
+ \item \code{lower\_bound} $\Rightarrow$ streng monoton
+ \item \code{upper\_bound} $\Rightarrow$ monoton
+ \end{itemize}
+ \sourcecode{math/longestIncreasingSubsequence.cpp}
+\end{algorithm}
+\vfill\null\columnbreak
+
+\begin{algorithm}{Zykel Erkennung}
+ \begin{methods}
+ \method{cycleDetection}{findet Zyklus von $x_0$ und Länge in $f$}{b+l}
+ \end{methods}
+ \sourcecode{math/cycleDetection.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Permutationen}
+ \begin{methods}
+ \method{kthperm}{findet $k$-te Permutation \big($k \in [0, n!$)\big)}{n\*\log(n)}
+ \end{methods}
+ \sourcecode{math/kthperm.cpp}
+ \begin{methods}
+ \method{permIndex}{bestimmt Index der Permutation \big($\mathit{res} \in [0, n!$)\big)}{n\*\log(n)}
+ \end{methods}
+ \sourcecode{math/permIndex.cpp}
+\end{algorithm}
+\columnbreak
+
+\subsection{Mod-Exponent und Multiplikation über $\boldsymbol{\mathbb{F}_p}$}
+%\vspace{-1.25em}
+%\begin{multicols}{2}
+\method{mulMod}{berechnet $a \cdot b \bmod n$}{\log(b)}
+\sourcecode{math/modMulIterativ.cpp}
+% \vfill\null\columnbreak
+\method{powMod}{berechnet $a^b \bmod n$}{\log(b)}
+\sourcecode{math/modPowIterativ.cpp}
+%\end{multicols}
+%\vspace{-2.75em}
+\begin{itemize}
+ \item für $a > 10^9$ \code{__int128} oder \code{modMul} benutzten!
+\end{itemize}
+
+\begin{algorithm}{ggT, kgV, erweiterter euklidischer Algorithmus}
+ \runtime{\log(a) + \log(b)}
+ \sourcecode{math/extendedEuclid.cpp}
+\end{algorithm}
+
+\subsection{Multiplikatives Inverses von $\boldsymbol{x}$ in $\boldsymbol{\mathbb{Z}/m\mathbb{Z}}$}
+\textbf{Falls $\boldsymbol{m}$ prim:}\quad $x^{-1} \equiv x^{m-2} \bmod m$
+
+\textbf{Falls $\boldsymbol{\ggT(x, m) = 1}$:}
+\begin{itemize}
+ \item Erweiterter euklidischer Algorithmus liefert $\alpha$ und $\beta$ mit
+ $\alpha x + \beta m = 1$.
+ \item Nach Kongruenz gilt $\alpha x + \beta m \equiv \alpha x \equiv 1 \bmod m$.
+ \item $x^{-1} :\equiv \alpha \bmod m$
+\end{itemize}
+\textbf{Sonst $\boldsymbol{\ggT(x, m) > 1}$:}\quad Es existiert kein $x^{-1}$.
+% \sourcecode{math/multInv.cpp}
+\sourcecode{math/shortModInv.cpp}
+
+\paragraph{Lemma von \textsc{Bézout}}
+Sei $(x, y)$ eine Lösung der diophantischen Gleichung $ax + by = d$.
+Dann lassen sich wie folgt alle Lösungen berechnen:
+\[
+\left(x + k\frac{b}{\ggT(a, b)},~y - k\frac{a}{\ggT(a, b)}\right)
+\]
+
+\paragraph{\textsc{Pell}-Gleichungen}
+Sei $(\overline{x}, \overline{y})$ die Lösung von $x^2 - ny^2 = 1$, die $x>1$ minimiert.
+Sei $(\tilde{x}, \tilde{y})$ die Lösung von $x^2-ny^2 = c$, die $x>1$ minimiert. Dann lassen
+sich alle Lösungen von $x^2-ny^2=c$ berechnen durch:
+\begin{align*}
+ x_1&\coloneqq \tilde{x}, & y_1&\coloneqq\tilde{y}\\
+ x_{k+1}&\coloneqq \overline{x}x_k+n\overline{y}y_k, & y_{k+1}&\coloneqq\overline{x}y_k+\overline{y}x_k
+\end{align*}
+
+\begin{algorithm}{Lineare Kongruenz}
+ \begin{itemize}
+ \item Kleinste Lösung $x$ für $ax\equiv b\pmod{m}$.
+ \item Weitere Lösungen unterscheiden sich um \raisebox{2pt}{$\frac{m}{g}$}, es gibt
+ also $g$ Lösungen modulo $m$.
+ \end{itemize}
+ \sourcecode{math/linearCongruence.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Chinesischer Restsatz}
+ \begin{itemize}
+ \item Extrem anfällig gegen Overflows. Evtl. häufig 128-Bit Integer verwenden.
+ \item Direkte Formel für zwei Kongruenzen $x \equiv a \bmod n$, $x \equiv b \bmod m$:
+ \[
+ x \equiv a - y \cdot n \cdot \frac{a - b}{d} \bmod \frac{mn}{d}
+ \qquad \text{mit} \qquad
+ d := \ggT(n, m) = yn + zm
+ \]
+ Formel kann auch für nicht teilerfremde Moduli verwendet werden.
+ Sind die Moduli nicht teilerfremd, existiert genau dann eine Lösung,
+ wenn $a\equiv~b \bmod \ggT(m, n)$.
+ In diesem Fall sind keine Faktoren
+ auf der linken Seite erlaubt.
+ \end{itemize}
+ \sourcecode{math/chineseRemainder.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Primzahltest \& Faktorisierung}
+ \method{isPrime}{prüft ob Zahl prim ist}{\log(n)^2}
+ \sourcecode{math/millerRabin.cpp}
+ \method{rho}{findet zufälligen Teiler}{\sqrt[\leftroot{3}\uproot{2}4]{n}}
+ \sourcecode{math/rho.cpp}
+ %\method{squfof}{findet zufälligen Teiler}{\sqrt[\leftroot{4}\uproot{2}4]{n}}
+ %\sourcecode{math/squfof.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Teiler}
+ \begin{methods}
+ \method{countDivisors}{Zählt Teiler von $n$}{\sqrt[\leftroot{3}\uproot{2}3]{n}}
+ \end{methods}
+ \sourcecode{math/divisors.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Matrix-Exponent}
+ \begin{methods}
+ \method{precalc}{berechnet $m^{2^b}$ vor}{\log(b)\*n^3}
+ \method{calc}{berechnet $m^b\cdot$}{\log(b)\cdot n^2}
+ \end{methods}
+ \textbf{Tipp:} wenn \code{v[x]=1} und \code{0} sonst, dann ist \code{res[y]} = $m^b_{y,x}$.
+ \sourcecode{math/matrixPower.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Lineare Rekurrenz}
+ \begin{methods}
+ \method{BerlekampMassey}{Berechnet eine lineare Rekurrenz $n$-ter Ordnung}{n^2}
+ \method{}{aus den ersten $2n$ Werte}{}
+ \end{methods}
+ \sourcecode{math/berlekampMassey.cpp}
+ Sei $f(n)=c_{0}f(n-1)+c_{1}f(n-2)+\dots + c_{n-1}f(0)$ eine lineare Rekurrenz.
+
+ \begin{methods}
+ \method{kthTerm}{Berechnet $k$-ten Term einer Rekurrenz $n$-ter Ordnung}{\log(k)\cdot n^2}
+ \end{methods}
+ \sourcecode{math/linearRecurrence.cpp}
+ Alternativ kann der \mbox{$k$-te} Term in \runtime{n^3\log(k)} berechnet werden:
+ $$\renewcommand\arraystretch{1.5}
+ \setlength\arraycolsep{3pt}
+ \begin{pmatrix}
+ c_{0} & c_{1} & \smash{\cdots} & \smash{\cdots} & c_{n-1} \\
+ 1 & 0 & \smash{\cdots} & \smash{\cdots} & 0 \\
+ 0 & \smash{\ddots} & \smash{\ddots} & & \smash{\vdots} \\
+ \smash{\vdots} & \smash{\ddots} & \smash{\ddots} & \smash{\ddots} & \smash{\vdots} \\
+ 0 & \smash{\cdots} & 0 & 1 & 0 \\
+ \end{pmatrix}^k
+ \times~~
+ \begin{pmatrix}
+ f(n-1) \\
+ f(n-2) \\
+ \smash{\vdots} \\
+ \smash{\vdots} \\
+ f(0) \\
+ \end{pmatrix}
+ ~~=~~
+ \begin{pmatrix}
+ f(n-1+k) \\
+ f(n-2+k) \\
+ \smash{\vdots} \\
+ \smash{\vdots} \\
+ f(k) \makebox[0pt][l]{\hspace{15pt}$\vcenter{\hbox{\huge$\leftarrow$}}$}\\
+ \end{pmatrix}
+ $$
+\end{algorithm}
+
+\begin{algorithm}{Diskreter Logarithmus}
+ \begin{methods}
+ \method{solve}{bestimmt Lösung $x$ für $a^x=b \bmod m$}{\sqrt{m}\*\log(m)}
+ \end{methods}
+ \sourcecode{math/discreteLogarithm.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Diskrete Quadratwurzel}
+ \begin{methods}
+ \method{sqrtMod}{bestimmt Lösung $x$ für $x^2=a \bmod p$ }{\log(p)}
+ \end{methods}
+ \textbf{Wichtig:} $p$ muss prim sein!
+ \sourcecode{math/sqrtModCipolla.cpp}
+\end{algorithm}
+%\columnbreak
+
+\begin{algorithm}{Primitivwurzeln}
+ \begin{itemize}
+ \item Primitivwurzel modulo $n$ existiert $\Leftrightarrow$ $n \in \{2,\ 4,\ p^\alpha,\ 2\cdot p^\alpha \mid\ 2 < p \in \mathbb{P},\ \alpha \in \mathbb{N}\}$
+ \item es existiert entweder keine oder $\varphi(\varphi(n))$ inkongruente Primitivwurzeln
+ \item Sei $g$ Primitivwurzel modulo $n$.
+ Dann gilt:\newline
+ Das kleinste $k$, sodass $g^k \equiv 1 \bmod n$, ist $k = \varphi(n)$.
+ \end{itemize}
+ \begin{methods}
+ \method{isPrimitive}{prüft ob $g$ eine Primitivwurzel ist}{\log(\varphi(n))\*\log(n)}
+ \method{findPrimitive}{findet Primitivwurzel (oder -1)}{\abs{ans}\*\log(\varphi(n))\*\log(n)}
+ \end{methods}
+ \sourcecode{math/primitiveRoot.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Diskrete \textrm{\textit{n}}-te Wurzel}
+ \begin{methods}
+ \method{root}{bestimmt Lösung $x$ für $x^a=b \bmod m$ }{\sqrt{m}\*\log(m)}
+ \end{methods}
+ Alle Lösungen haben die Form $g^{c + \frac{i \cdot \phi(n)}{\gcd(a, \phi(n))}}$
+ \sourcecode{math/discreteNthRoot.cpp}
+\end{algorithm}
+
+\begin{algorithm}{\textsc{Legendre}-Symbol}
+ Sei $p \geq 3$ eine Primzahl, $a \in \mathbb{Z}$:
+ \vspace{-0.15cm}\begin{align*}
+ \hspace*{3cm}\legendre{a}{p} &=
+ \begin{cases*}
+ \hphantom{-}0 & falls $p~\vert~a$ \\[-1ex]
+ \hphantom{-}1 & falls $\exists x \in \mathbb{Z}\backslash p\mathbb{Z} : a \equiv x^2 \bmod p$ \\[-1ex]
+ -1 & sonst
+ \end{cases*} \\
+ \legendre{-1}{p} = (-1)^{\frac{p - 1}{2}} &=
+ \begin{cases*}
+ \hphantom{-}1 & falls $p \equiv 1 \bmod 4$ \\[-1ex]
+ -1 & falls $p \equiv 3 \bmod 4$
+ \end{cases*} \\
+ \legendre{2}{p} = (-1)^{\frac{p^2 - 1}{8}} &=
+ \begin{cases*}
+ \hphantom{-}1 & falls $p \equiv \pm 1 \bmod 8$ \\[-1ex]
+ -1 & falls $p \equiv \pm 3 \bmod 8$
+ \end{cases*}
+ \end{align*}
+ \begin{align*}
+ \legendre{p}{q} \cdot \legendre{q}{p} = (-1)^{\frac{p - 1}{2} \cdot \frac{q - 1}{2}} &&
+ \legendre{a}{p} \equiv a^{\frac{p-1}{2}}\bmod p
+ \end{align*}
+ \vspace{-0.05cm}
+ \sourcecode{math/legendre.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Lineares Sieb und multiplikative Funktionen}
+ Eine (zahlentheoretische) Funktion $f$ heißt multiplikativ wenn $f(1)=1$ und $f(a\cdot b)=f(a)\cdot f(b)$, falls $\ggT(a,b)=1$.
+
+ $\Rightarrow$ Es ist ausreichend $f(p^k)$ für alle primen $p$ und alle $k$ zu kennen.
+
+ \begin{methods}
+ \method{sieve}{berechnet Primzahlen und co.}{N}
+ \method{sieved}{Wert der entsprechenden multiplikativen Funktion}{1}
+
+ \method{naive}{Wert der entsprechenden multiplikativen Funktion}{\sqrt{n}}
+ \end{methods}
+ \textbf{Wichtig:} Sieb rechts ist schneller für \code{isPrime} oder \code{primes}!
+
+ \sourcecode{math/linearSieve.cpp}
+ \textbf{\textsc{Möbius} Funktion:}
+ \begin{itemize}
+ \item $\mu(n)=+1$, falls $n$ quadratfrei ist und gerade viele Primteiler hat
+ \item $\mu(n)=-1$, falls $n$ quadratfrei ist und ungerade viele Primteiler hat
+ \item $\mu(n)=0$, falls $n$ nicht quadratfrei ist
+ \end{itemize}
+
+ \textbf{\textsc{Euler}sche $\boldsymbol{\varphi}$-Funktion:}
+ \begin{itemize}
+ \item Zählt die relativ primen Zahlen $\leq n$.
+ \item $p$ prim, $k \in \mathbb{N}$:
+ $~\varphi(p^k) = p^k - p^{k - 1}$
+
+ \item \textbf{\textsc{Euler}'s Theorem:}
+ Für $b \geq \varphi(c)$ gilt: $a^b \equiv a^{b \bmod \varphi(c) + \varphi(c)} \pmod{c}$. Darüber hinaus gilt: $\gcd(a, c) = 1 \Leftrightarrow a^b \equiv a^{b \bmod \varphi(c)} \pmod{c}$.
+ Falls $m$ prim ist, liefert das den \textbf{kleinen Satz von \textsc{Fermat}}:
+ $a^{m} \equiv a \pmod{m}$
+ \end{itemize}
+\end{algorithm}
+
+\begin{algorithm}{Primzahlsieb von \textsc{Eratosthenes}}
+ \begin{itemize}
+ \item Bis $10^8$ in unter 64MB Speicher (lange Berechnung)
+ \end{itemize}
+ \begin{methods}
+ \method{primeSieve}{berechnet Primzahlen und Anzahl}{N\*\log(\log(N))}
+ \method{isPrime}{prüft ob Zahl prim ist}{1}
+ \end{methods}
+ \sourcecode{math/primeSieve.cpp}
+\end{algorithm}
+
+\begin{algorithm}{\textsc{Möbius}-Inversion}
+ \begin{itemize}
+ \item Seien $f,g : \mathbb{N} \to \mathbb{N}$ und $g(n) := \sum_{d \vert n}f(d)$.
+ Dann ist $f(n) = \sum_{d \vert n}g(d)\mu(\frac{n}{d})$.
+ \item $\sum\limits_{d \vert n}\mu(d) =
+ \begin{cases*}
+ 1 & falls $n = 1$\\
+ 0 & sonst
+ \end{cases*}$
+ \end{itemize}
+ \textbf{Beispiel Inklusion/Exklusion:}
+ Gegeben sein eine Sequenz $A={a_1,\ldots,a_n}$ von Zahlen, $1 \leq a_i \leq N$. Zähle die Anzahl der \emph{coprime subsequences}.\newline
+ \textbf{Lösung}:
+ Für jedes $x$, sei $cnt[x]$ die Anzahl der Vielfachen von $x$ in $A$.
+ Es gibt $2^{[x]}-1$ nicht leere Subsequences in $A$, die nur Vielfache von $x$ enthalten.
+ Die Anzahl der Subsequences mit $\ggT=1$ ist gegeben durch $\sum_{i = 1}^N \mu(i) \cdot (2^{cnt[i]} - 1)$.
+\end{algorithm}
+
+\subsection{LGS über $\boldsymbol{\mathbb{F}_p}$}
+\method{gauss}{löst LGS}{n^3}
+\sourcecode{math/lgsFp.cpp}
+
+\subsection{LGS über $\boldsymbol{\mathbb{R}}$}
+\method{gauss}{löst LGS}{n^3}
+\sourcecode{math/gauss.cpp}
+
+\vfill\null\columnbreak
+
+\begin{algorithm}{Numerisch Extremstelle bestimmen}
+ \sourcecode{math/goldenSectionSearch.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Numerisch Integrieren, Simpsonregel}
+ \sourcecode{math/simpson.cpp}
+\end{algorithm}
+
+
+\begin{algorithm}{Polynome, FFT, NTT \& andere Transformationen}
+ Multipliziert Polynome $A$ und $B$.
+ \begin{itemize}
+ \item $\deg(A \cdot B) = \deg(A) + \deg(B)$
+ \item Vektoren \code{a} und \code{b} müssen mindestens Größe
+ $\deg(A \cdot B) + 1$ haben.
+ Größe muss eine Zweierpotenz sein.
+ \item Für ganzzahlige Koeffizienten: \code{(ll)round(real(a[i]))}
+ \item \emph{xor}, \emph{or} und \emph{and} Transform funktioniert auch mit \code{double} oder modulo einer Primzahl $p$ falls $p \geq 2^{\texttt{bits}}$
+ \end{itemize}
+ %\sourcecode{math/fft.cpp}
+ %\sourcecode{math/ntt.cpp}
+ \sourcecode{math/transforms/fft.cpp}
+ \sourcecode{math/transforms/ntt.cpp}
+ \sourcecode{math/transforms/bitwiseTransforms.cpp}
+ Multiplikation mit 2 transforms statt 3: (nur benutzten wenn nötig!)
+ \sourcecode{math/transforms/fftMul.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Operations on Formal Power Series}
+ \sourcecode{math/transforms/seriesOperations.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Inversionszahl}
+ \sourcecode{math/inversions.cpp}
+\end{algorithm}
+
+\subsection{Satz von \textsc{Sprague-Grundy}}
+Weise jedem Zustand $X$ wie folgt eine \textsc{Grundy}-Zahl $g\left(X\right)$ zu:
+\[
+g\left(X\right) := \min\left\{
+\mathbb{Z}_0^+ \setminus
+\left\{g\left(Y\right) \mid Y \text{ von } X \text{ aus direkt erreichbar}\right\}
+\right\}
+\]
+$X$ ist genau dann gewonnen, wenn $g\left(X\right) > 0$ ist.\\
+Wenn man $k$ Spiele in den Zuständen $X_1, \ldots, X_k$ hat, dann ist die \textsc{Grundy}-Zahl des Gesamtzustandes $g\left(X_1\right) \oplus \ldots \oplus g\left(X_k\right)$.
+
+\subsection{Kombinatorik}
+
+\paragraph{\textsc{Wilson}'s Theorem}
+A number $n$ is prime if and only if
+$(n-1)!\equiv -1\bmod{n}$.\\
+($n$ is prime if and only if $(m-1)!\cdot(n-m)!\equiv(-1)^m\bmod{n}$ for all $m$ in $\{1,\dots,n\}$)
+\begin{align*}
+ (n-1)!\equiv\begin{cases}
+ -1\bmod{n},&\mathrm{falls}~n \in \mathbb{P}\\
+ \hphantom{-}2\bmod{n},&\mathrm{falls}~n = 4\\
+ \hphantom{-}0\bmod{n},&\mathrm{sonst}
+ \end{cases}
+\end{align*}
+
+\paragraph{\textsc{Zeckendorf}'s Theorem}
+Jede positive natürliche Zahl kann eindeutig als Summe einer oder mehrerer
+verschiedener \textsc{Fibonacci}-Zahlen geschrieben werden, sodass keine zwei
+aufeinanderfolgenden \textsc{Fibonacci}-Zahlen in der Summe vorkommen.\\
+\emph{Lösung:} Greedy, nimm immer die größte \textsc{Fibonacci}-Zahl, die noch
+hineinpasst.
+
+\paragraph{\textsc{Lucas}'s Theorem}
+Ist $p$ prim, $m=\sum_{i=0}^km_ip^i$, $n=\sum_{i=0}^kn_ip^i$ ($p$-adische Darstellung),
+so gilt
+\vspace{-0.75\baselineskip}
+\[
+ \binom{m}{n} \equiv \prod_{i=0}^k\binom{m_i}{n_i} \bmod{p}.
+\]
+
+%\begin{algorithm}{Binomialkoeffizienten}
+\paragraph{Binomialkoeffizienten}
+ Die Anzahl der \mbox{$k$-elementigen} Teilmengen einer \mbox{$n$-elementigen} Menge.
+
+ \begin{methods}
+ \method{precalc}{berechnet $n!$ und $n!^{-1}$ vor}{\mathit{lim}}
+ \method{calc\_binom}{berechnet Binomialkoeffizient}{1}
+ \end{methods}
+ \sourcecode{math/binomial0.cpp}
+ Falls $n >= p$ for $\mathit{mod}=p^k$ berechne \textit{fac} und \textit{inv} aber teile $p$ aus $i$ und berechne die häufigkeit von $p$ in $n!$ als $\sum\limits_{i=1}\big\lfloor\frac{n}{p^i}\big\rfloor$
+
+ \begin{methods}
+ \method{calc\_binom}{berechnet Binomialkoeffizient $(n \le 61)$}{k}
+ \end{methods}
+ \sourcecode{math/binomial1.cpp}
+
+ \begin{methods}
+ \method{calc\_binom}{berechnet Binomialkoeffizient modulo Primzahl $p$}{p-n}
+ \end{methods}
+ \sourcecode{math/binomial3.cpp}
+
+% \begin{methods}
+% \method{calc\_binom}{berechnet Primfaktoren vom Binomialkoeffizient}{n}
+% \end{methods}
+% \textbf{WICHTIG:} braucht alle Primzahlen $\leq n$
+% \sourcecode{math/binomial2.cpp}
+%\end{algorithm}
+
+\paragraph{\textsc{Catalan}-Zahlen}
+\begin{itemize}
+ \item Die \textsc{Catalan}-Zahl $C_n$ gibt an:
+ \begin{itemize}
+ \item Anzahl der Binärbäume mit $n$ nicht unterscheidbaren Knoten.
+ \item Anzahl der validen Klammerausdrücke mit $n$ Klammerpaaren.
+ \item Anzahl der korrekten Klammerungen von $n+1$ Faktoren.
+ \item Anzahl Möglichkeiten ein konvexes Polygon mit $n + 2$ Ecken zu triangulieren.
+ \item Anzahl der monotonen Pfade (zwischen gegenüberliegenden Ecken) in
+ einem $n \times n$-Gitter, die nicht die Diagonale kreuzen.
+ \end{itemize}
+\end{itemize}
+\[C_0 = 1\qquad C_n = \sum\limits_{k = 0}^{n - 1} C_kC_{n - 1 - k} =
+\frac{1}{n + 1}\binom{2n}{n} = \frac{4n - 2}{n+1} \cdot C_{n-1}\]
+\begin{itemize}
+ \item Formel $1$ erlaubt Berechnung ohne Division in \runtime{n^2}
+ \item Formel $2$ und $3$ erlauben Berechnung in \runtime{n}
+\end{itemize}
+
+\paragraph{\textsc{Catalan}-Convolution}
+\begin{itemize}
+ \item Anzahl an Klammerausdrücken mit $n+k$ Klammerpaaren, die mit $(^k$ beginnen.
+\end{itemize}
+\[C^k_0 = 1\qquad C^k_n = \sum\limits_{\mathclap{a_0+a_1+\dots+a_k=n}} C_{a_0}C_{a_1}\cdots C_{a_k} =
+\frac{k+1}{n+k+1}\binom{2n+k}{n} = \frac{(2n+k-1)\cdot(2n+k)}{n(n+k+1)} \cdot C_{n-1}\]
+
+\paragraph{\textsc{Euler}-Zahlen 1. Ordnung}
+Die Anzahl der Permutationen von $\{1, \ldots, n\}$ mit genau $k$ Anstiegen.
+Für die $n$-te Zahl gibt es $n$ mögliche Positionen zum Einfügen.
+Dabei wird entweder ein Anstieg in zwei gesplitted oder ein Anstieg um $n$ ergänzt.
+\[\eulerI{n}{0} = \eulerI{n}{n-1} = 1 \quad
+\eulerI{n}{k} = (k+1) \eulerI{n-1}{k} + (n-k) \eulerI{n-1}{k-1}=
+\sum_{i=0}^{k} (-1)^i\binom{n+1}{i}(k+1-i)^n\]
+\begin{itemize}
+ \item Formel $1$ erlaubt Berechnung ohne Division in \runtime{n^2}
+ \item Formel $2$ erlaubt Berechnung in \runtime{n\log(n)}
+\end{itemize}
+
+\paragraph{\textsc{Euler}-Zahlen 2. Ordnung}
+Die Anzahl der Permutationen von $\{1,1, \ldots, n,n\}$ mit genau $k$ Anstiegen.
+\[\eulerII{n}{0} = 1 \qquad\eulerII{n}{n} = 0 \qquad\eulerII{n}{k} = (k+1) \eulerII{n-1}{k} + (2n-k-1) \eulerII{n-1}{k-1}\]
+\begin{itemize}
+ \item Formel erlaubt Berechnung ohne Division in \runtime{n^2}
+\end{itemize}
+
+\paragraph{\textsc{Stirling}-Zahlen 1. Ordnung}
+Die Anzahl der Permutationen von $\{1, \ldots, n\}$ mit genau $k$ Zyklen.
+Es gibt zwei Möglichkeiten für die $n$-te Zahl. Entweder sie bildet einen eigene Zyklus, oder sie kann an jeder Position in jedem Zyklus einsortiert werden.
+\[\stirlingI{0}{0} = 1 \qquad
+\stirlingI{n}{0} = \stirlingI{0}{n} = 0 \qquad
+\stirlingI{n}{k} = \stirlingI{n-1}{k-1} + (n-1) \stirlingI{n-1}{k}\]
+\begin{itemize}
+ \item Formel erlaubt berechnung ohne Division in \runtime{n^2}
+\end{itemize}
+\[\sum_{k=0}^{n}\pm\stirlingI{n}{k}x^k=x(x-1)(x-2)\cdots(x-n+1)\]
+\begin{itemize}
+ \item Berechne Polynom mit FFT und benutzte betrag der Koeffizienten \runtime{n\log(n)^2} (nur ungefähr gleich große Polynome zusammen multiplizieren beginnend mit $x-k$)
+\end{itemize}
+
+\paragraph{\textsc{Stirling}-Zahlen 2. Ordnung}
+Die Anzahl der Möglichkeiten $n$ Elemente in $k$ nichtleere Teilmengen zu zerlegen.
+Es gibt $k$ Möglichkeiten die $n$ in eine $n-1$-Partition einzuordnen.
+Dazu kommt der Fall, dass die $n$ in ihrer eigenen Teilmenge (alleine) steht.
+\[\stirlingII{n}{1} = \stirlingII{n}{n} = 1 \qquad
+\stirlingII{n}{k} = k \stirlingII{n-1}{k} + \stirlingII{n-1}{k-1} =
+\frac{1}{k!} \sum\limits_{i=0}^{k} (-1)^{k-i}\binom{k}{i}i^n\]
+\begin{itemize}
+ \item Formel $1$ erlaubt Berechnung ohne Division in \runtime{n^2}
+ \item Formel $2$ erlaubt Berechnung in \runtime{n\log(n)}
+\end{itemize}
+
+\paragraph{\textsc{Bell}-Zahlen}
+Anzahl der Partitionen von $\{1, \ldots, n\}$.
+Wie \textsc{Stirling}-Zahlen 2. Ordnung ohne Limit durch $k$.
+\[B_1 = 1 \qquad
+B_n = \sum\limits_{k = 0}^{n - 1} B_k\binom{n-1}{k}
+= \sum\limits_{k = 0}^{n}\stirlingII{n}{k}\qquad\qquad B_{p^m+n}\equiv m\cdot B_n + B_{n+1} \bmod{p}\]
+
+\paragraph{Partitions}
+Die Anzahl der Partitionen von $n$ in genau $k$ positive Summanden.
+Die Anzahl der Partitionen von $n$ mit Elementen aus ${1,\dots,k}$.
+\begin{align*}
+ p_0(0)=1 \qquad p_k(n)&=0 \text{ für } k > n \text{ oder } n \leq 0 \text{ oder } k \leq 0\\
+ p_k(n)&= p_k(n-k) + p_{k-1}(n-1)\\[2pt]
+ p(n)&=\sum_{k=1}^{n} p_k(n)=p_n(2n)=\sum\limits_{k\neq0}^\infty(-1)^{k+1}p\bigg(n - \frac{k(3k-1)}{2}\bigg)
+\end{align*}
+\begin{itemize}
+ \item in Formel $3$ kann abgebrochen werden wenn $\frac{k(3k-1)}{2} > n$.
+ \item Die Anzahl der Partitionen von $n$ in bis zu $k$ positive Summanden ist $\sum\limits_{i=0}^{k}p_i(n)=p_k(n+k)$.
+\end{itemize}
+
+\subsection{The Twelvefold Way \textnormal{(verteile $n$ Bälle auf $k$ Boxen)}}
+\input{math/tables/twelvefold}
+
+\optional{
+\subsection{Primzahlzählfunktion $\boldsymbol{\pi}$}
+\begin{methods}
+ \method{init}{berechnet $\pi$ bis $N$}{N\*\log(\log(N))}
+ \method{phi}{zählt zu $p_i$ teilerfremde Zahlen $\leq n$ für alle $i \leq k$}{???}
+ \method{pi}{zählt Primzahlen $\leq n$ ($n < N^2$)}{n^{2/3}}
+\end{methods}
+\sourcecode{math/piLehmer.cpp}
+}
+
+%\input{math/tables/numbers}
+
+\begin{algorithm}[optional]{Big Integers}
+ \sourcecode{math/bigint.cpp}
+\end{algorithm}
diff --git a/content/math/matrixPower.cpp b/content/math/matrixPower.cpp new file mode 100644 index 0000000..d981e6e --- /dev/null +++ b/content/math/matrixPower.cpp @@ -0,0 +1,14 @@ +vector<mat> pows; + +void precalc(mat m) { + pows = {mat(sz(m.m), 1), m}; + for (int i = 1; i < 60; i++) pows.push_back(pows[i] * pows[i]); +} + +auto calc(ll b, vector<ll> v) { + for (ll i = 1; b > 0; i++) { + if (b & 1) v = pows[i] * v; + b /= 2; + } + return v; +} diff --git a/content/math/millerRabin.cpp b/content/math/millerRabin.cpp new file mode 100644 index 0000000..cb27d29 --- /dev/null +++ b/content/math/millerRabin.cpp @@ -0,0 +1,19 @@ +constexpr ll bases32[] = {2, 7, 61}; +constexpr ll bases64[] = {2, 325, 9375, 28178, 450775, + 9780504, 1795265022}; +bool isPrime(ll n) { + if (n < 2 || n % 2 == 0) return n == 2; + ll d = n - 1, j = 0; + while (d % 2 == 0) d /= 2, j++; + for (ll a : bases64) { + if (a % n == 0) continue; + ll v = powMod(a, d, n); //with mulmod or int128 + if (v == 1 || v == n - 1) continue; + for (ll i = 1; i <= j; i++) { + v = ((lll)v * v) % n; + if (v == n - 1 || v <= 1) break; + } + if (v != n - 1) return false; + } + return true; +} diff --git a/content/math/modExp.cpp b/content/math/modExp.cpp new file mode 100644 index 0000000..2329a94 --- /dev/null +++ b/content/math/modExp.cpp @@ -0,0 +1,6 @@ +ll powMod(ll a, ll b, ll n) { + if(b == 0) return 1; + if(b == 1) return a % n; + if(b & 1) return (powMod(a, b - 1, n) * a) % n; + else return powMod((a * a) % n, b / 2, n); +} diff --git a/content/math/modMulIterativ.cpp b/content/math/modMulIterativ.cpp new file mode 100644 index 0000000..611f09a --- /dev/null +++ b/content/math/modMulIterativ.cpp @@ -0,0 +1,9 @@ +ll mulMod(ll a, ll b, ll n) { + ll res = 0; + while (b > 0) { + if (b & 1) res = (a + res) % n; + a = (a * 2) % n; + b /= 2; + } + return res; +} diff --git a/content/math/modPowIterativ.cpp b/content/math/modPowIterativ.cpp new file mode 100644 index 0000000..0dc3fb1 --- /dev/null +++ b/content/math/modPowIterativ.cpp @@ -0,0 +1,9 @@ +ll powMod(ll a, ll b, ll n) { + ll res = 1; + while (b > 0) { + if (b & 1) res = (a * res) % n; + a = (a * a) % n; + b /= 2; + } + return res; +} diff --git a/content/math/multInv.cpp b/content/math/multInv.cpp new file mode 100644 index 0000000..647dc2d --- /dev/null +++ b/content/math/multInv.cpp @@ -0,0 +1,4 @@ +ll multInv(ll x, ll m) { + auto [d, a, b] = extendedEuclid(x, m); // Implementierung von oben. + return ((a % m) + m) % m; +} diff --git a/content/math/permIndex.cpp b/content/math/permIndex.cpp new file mode 100644 index 0000000..4cffc12 --- /dev/null +++ b/content/math/permIndex.cpp @@ -0,0 +1,13 @@ +ll permIndex(vector<ll> v) { + Tree<ll> t; + reverse(all(v)); + for (ll& x : v) { + t.insert(x); + x = t.order_of_key(x); + } + ll res = 0; + for (int i = sz(v); i > 0; i--) { + res = res * i + v[i - 1]; + } + return res; +} diff --git a/content/math/piLegendre.cpp b/content/math/piLegendre.cpp new file mode 100644 index 0000000..21b974b --- /dev/null +++ b/content/math/piLegendre.cpp @@ -0,0 +1,23 @@ +constexpr ll cache = 500; // requires O(cache^3)
+vector<vector<ll>> memo(cache * cache, vector<ll>(cache));
+
+ll pi(ll n);
+
+ll phi(ll n, ll k) {
+ if (n <= 1 || k < 0) return 0;
+ if (n <= primes[k]) return n - 1;
+ if (n < N && primes[k] * primes[k] > n) return n - pi(n) + k;
+ bool ok = n < cache * cache;
+ if (ok && memo[n][k] > 0) return memo[n][k];
+ ll res = n/primes[k] - phi(n/primes[k], k - 1) + phi(n, k - 1);
+ if (ok) memo[n][k] = res;
+ return res;
+}
+
+ll pi(ll n) {
+ if (n < N) { // implement this as O(1) lookup for speedup!
+ return distance(primes.begin(), upper_bound(all(primes), n));
+ } else {
+ ll k = pi(sqrtl(n) + 1);
+ return n - phi(n, k) + k;
+}}
diff --git a/content/math/piLehmer.cpp b/content/math/piLehmer.cpp new file mode 100644 index 0000000..17df85e --- /dev/null +++ b/content/math/piLehmer.cpp @@ -0,0 +1,52 @@ +constexpr ll cacheA = 2 * 3 * 5 * 7 * 11 * 13 * 17;
+constexpr ll cacheB = 7;
+ll memoA[cacheA + 1][cacheB + 1];
+ll memoB[cacheB + 1];
+ll memoC[N];
+
+void init() {
+ primeSieve(); // @\sourceref{math/primeSieve.cpp}@
+ for (ll i = 0; i < N; i++) {
+ memoC[i] = memoC[i - 1];
+ if (isPrime(i)) memoC[i]++;
+ }
+ memoB[0] = 1;
+ for(ll i = 0; i <= cacheA; i++) memoA[i][0] = i;
+ for(ll i = 1; i <= cacheB; i++) {
+ memoB[i] = primes[i - 1] * memoB[i - 1];
+ for(ll j = 1; j <= cacheA; j++) {
+ memoA[j][i] = memoA[j][i - 1] - memoA[j /
+ primes[i - 1]][i - 1];
+}}}
+
+ll phi(ll n, ll k) {
+ if(k == 0) return n;
+ if(k <= cacheB)
+ return memoA[n % memoB[k]][k] +
+ (n / memoB[k]) * memoA[memoB[k]][k];
+ if(n <= primes[k - 1]*primes[k - 1]) return memoC[n] - k + 1;
+ if(n <= primes[k - 1]*primes[k - 1]*primes[k - 1] && n < N) {
+ ll b = memoC[(ll)sqrtl(n)];
+ ll res = memoC[n] - (b + k - 2) * (b - k + 1) / 2;
+ for(ll i = k; i < b; i++) res += memoC[n / primes[i]];
+ return res;
+ }
+ return phi(n, k - 1) - phi(n / primes[k - 1], k - 1);
+}
+
+ll pi(ll n) {
+ if (n < N) return memoC[n];
+ ll a = pi(sqrtl(sqrtl(n)));
+ ll b = pi(sqrtl(n));
+ ll c = pi(cbrtl(n));
+ ll res = phi(n, a) + (b + a - 2) * (b - a + 1) / 2;
+ for (ll i = a; i < b; i++) {
+ ll w = n / primes[i];
+ res -= pi(w);
+ if (i > c) continue;
+ ll bi = pi(sqrtl(w));
+ for (ll j = i; j < bi; j++) {
+ res -= pi(w / primes[j]) - j;
+ }}
+ return res;
+}
diff --git a/content/math/polynomial.cpp b/content/math/polynomial.cpp new file mode 100644 index 0000000..44f6207 --- /dev/null +++ b/content/math/polynomial.cpp @@ -0,0 +1,65 @@ +struct poly { + vector<ll> data; + + poly(int deg = 0) : data(max(1, deg)) {} + poly(initializer_list<ll> _data) : data(_data) {} + + int size() const {return sz(data);} + + void trim() { + for (ll& x : data) x = (x % mod + mod) % mod; + while (size() > 1 && data.back() == 0) data.pop_back(); + } + + ll& operator[](int x) {return data[x];} + const ll& operator[](int x) const {return data[x];} + + ll operator()(int x) const { + ll res = 0; + for (int i = size() - 1; i >= 0; i--) + res = (res * x + data[i]) % mod; + return res % mod; + } + + poly& operator+=(const poly& o) { + if (size() < o.size()) data.resize(o.size()); + for (int i = 0; i < o.size(); i++) + data[i] = (data[i] + o[i]) % mod; + return *this; + } + + poly operator*(const poly& o) const { + poly res(size() + o.size() - 1); + for (int i = 0; i < size(); i++) { + for (int j = 0; j < o.size(); j++) { + res[i + j] += (data[i] * o[j]) % mod; + }} + res.trim(); + return res; + } + + //return p(x+a) + poly operator<<(ll a) const { + poly res(size()); + for (int i = size() - 1; i >= 0; i--) { + for (int j = size() - i - 1; j >= 1; j--) + res[j] = (res[j] * a + res[j - 1]) % mod; + res[0] = (res[0] * a + res[i]) % mod; + } + return res; + } + + pair<poly, poly> divmod(const poly& d) const { + int i = size() - d.size(); + poly s(i + 1), r = *this; + ll inv = multInv(d.data.back(), mod); + for (; i >= 0; i--) { + s[i] = (r.data.back() * inv) % mod; + r.data.pop_back(); + for (int j = 0; i + j < r.size(); j++) { + r[i + j] = (r.data[i + j] - s[i] * d[j]) % mod; + }} + s.trim(); r.trim(); + return {s, r}; + } +}; diff --git a/content/math/primeSieve.cpp b/content/math/primeSieve.cpp new file mode 100644 index 0000000..1b0f514 --- /dev/null +++ b/content/math/primeSieve.cpp @@ -0,0 +1,16 @@ +constexpr ll N = 100'000'000; +bitset<N / 2> isNotPrime; +vector<ll> primes = {2}; + +bool isPrime(ll x) { + if (x < 2 || x % 2 == 0) return x == 2; + else return !isNotPrime[x / 2]; +} + +void primeSieve() { + for (ll i = 3; i < N; i += 2) {// i * i < N reicht für isPrime + if (!isNotPrime[i / 2]) { + primes.push_back(i); // optional + for (ll j = i * i; j < N; j+= 2 * i) { + isNotPrime[j / 2] = 1; +}}}} diff --git a/content/math/primitiveRoot.cpp b/content/math/primitiveRoot.cpp new file mode 100644 index 0000000..39a0f64 --- /dev/null +++ b/content/math/primitiveRoot.cpp @@ -0,0 +1,23 @@ +bool isPrimitive(ll g, ll n, ll phi, map<ll, int>& phiFacts) { + if (g == 1) return n == 2; + if (gcd(g, n) > 1) return false; + for (auto [f, _] : phiFacts) + if (powMod(g, phi / f, n) == 1) return false; + return true; +} + +bool isPrimitive(ll g, ll n) { + ll phin = phi(n); //isPrime(n) => phi(n) = n - 1 + map<ll, int> phiFacts; + factor(phin, phiFacts); + return isPrimitive(g, n, phin, phiFacts); +} + +ll findPrimitive(ll n) { //test auf existens geht schneller + ll phin = phi(n); //isPrime(n) => phi(n) = n - 1 + map<ll, int> phiFacts; + factor(phin, phiFacts); + for (ll res = 1; res < n; res++) // oder zufällige Reihenfolge + if (isPrimitive(res, n, phin, phiFacts)) return res; + return -1; +} diff --git a/content/math/rho.cpp b/content/math/rho.cpp new file mode 100644 index 0000000..ad640cd --- /dev/null +++ b/content/math/rho.cpp @@ -0,0 +1,19 @@ +using lll = __int128; +ll rho(ll n) { // Findet Faktor < n, nicht unbedingt prim. + if (n % 2 == 0) return 2; + ll x = 0, y = 0, prd = 2, i = n/2 + 7; + auto f = [&](lll c){return (c * c + i) % n;}; + for (ll t = 30; t % 40 || gcd(prd, n) == 1; t++) { + if (x == y) x = ++i, y = f(x); + if (ll q = (lll)prd * abs(x-y) % n; q) prd = q; + x = f(x); y = f(f(y)); + } + return gcd(prd, n); +} + +void factor(ll n, map<ll, int>& facts) { + if (n == 1) return; + if (isPrime(n)) {facts[n]++; return;} + ll f = rho(n); + factor(n / f, facts); factor(f, facts); +} diff --git a/content/math/shortModInv.cpp b/content/math/shortModInv.cpp new file mode 100644 index 0000000..f696cce --- /dev/null +++ b/content/math/shortModInv.cpp @@ -0,0 +1,3 @@ +ll multInv(ll x, ll m) { // x^{-1} mod m + return 1 < x ? m - multInv(m % x, x) * m / x : 1; +} diff --git a/content/math/simpson.cpp b/content/math/simpson.cpp new file mode 100644 index 0000000..7f237a4 --- /dev/null +++ b/content/math/simpson.cpp @@ -0,0 +1,12 @@ +//double f(double x) {return x;} + +double simps(double a, double b) { + return (f(a) + 4.0 * f((a + b) / 2.0) + f(b)) * (b - a) / 6.0; +} + +double integrate(double a, double b) { + double m = (a + b) / 2.0; + double l = simps(a, m), r = simps(m, b), tot = simps(a, b); + if (abs(l + r - tot) < EPS) return tot; + return integrate(a, m) + integrate(m, b); +} diff --git a/content/math/sqrtModCipolla.cpp b/content/math/sqrtModCipolla.cpp new file mode 100644 index 0000000..1fac0c5 --- /dev/null +++ b/content/math/sqrtModCipolla.cpp @@ -0,0 +1,14 @@ +ll sqrtMod(ll a, ll p) {// teste mit legendre ob lösung existiert + if (a < 2) return a; + ll t = 0; + while (legendre((t*t-4*a) % p, p) >= 0) t = rng() % p; + ll b = -t, c = -t, d = 1, m = p; + for (m++; m /= 2; b = (a+a-b*b) % p, a = (a*a) % p) { + if (m % 2) { + d = (c-d*b) % p; + c = (c*a) % p; + } else { + c = (d*a - c*b) % p; + }} + return (d + p) % p; +} diff --git a/content/math/squfof.cpp b/content/math/squfof.cpp new file mode 100644 index 0000000..1cb97de --- /dev/null +++ b/content/math/squfof.cpp @@ -0,0 +1,89 @@ +using lll = __int128;
+
+constexpr lll multipliers[] = {1, 3, 5, 7,
+ 11, 3*5, 3*7, 3*11,
+ 5*7, 5*11, 7*11,
+ 3*5*7, 3*5*11, 3*7*11,
+ 5*7*11, 3*5*7*11};
+
+lll root(lll x) {
+ lll r = sqrtl(x);
+ while(r*r < x) r++;
+ while(r*r > x) r--;
+ return r;
+}
+
+lll croot(lll x) {
+ lll r = cbrtl(x);
+ while(r*r*r < x) r++;
+ while(r*r*r > x) r--;
+ return r;
+}
+
+lll squfof(lll N) {
+ lll s = croot(N);
+ if (s*s*s == N) return s;
+ s = root(N);
+ if (s*s == N) return s;
+ for (lll k : multipliers) {
+ lll D = k * N;
+ lll Po, P, Pprev, q, b, r, i;
+ Po = Pprev = P = root(D);
+ lll Qprev = 1;
+ lll Q = D - Po*Po;
+ lll L = 2 * root(2 * s);
+ lll B = 3 * L;
+ for (i = 2; i < B; i++) {
+ b = (Po + P) / Q;
+ P = b*Q - P;
+ q = Q;
+ Q = Qprev + b * (Pprev - P);
+ r = root(Q);
+ if (!(i & 1) && r*r == Q) break;
+ Qprev = q;
+ Pprev = P;
+ }
+ if (i >= B) continue;
+ b = (Po - P) / r;
+ Pprev = P = b*r + P;
+ Qprev = r;
+ Q = (D-Pprev*Pprev)/Qprev;
+ i = 0;
+ do {
+ b = (Po + P) / Q;
+ Pprev = P;
+ P = b*Q - P;
+ q = Q;
+ Q = Qprev + b * (Pprev - P);
+ Qprev = q;
+ i++;
+ } while(P != Pprev);
+ r = gcd(N, Qprev);
+ if (r != 1 && r != N) return r;
+ }
+ exit(1);//try fallback to pollard rho
+}
+
+constexpr lll trialLim = 5'000;
+
+void factor(lll n, map<lll, int>& facts) {
+ for (lll i = 2; i * i <= n && i <= trialLim; i++) {
+ while (n % i == 0) {
+ facts[i]++;
+ n /= i;
+ }}
+ if (n > 1 && n < trialLim * trialLim) {
+ facts[n]++;
+ } else {
+ vector<lll> todo = {n};
+ while (!todo.empty()) {
+ lll c = todo.back();
+ todo.pop_back();
+ if (c == 1) continue;
+ if (isPrime(c)) {
+ facts[c]++;
+ } else {
+ lll d = squfof(c);
+ todo.push_back(d);
+ todo.push_back(c / d);
+}}}}
diff --git a/content/math/tables.tex b/content/math/tables.tex new file mode 100644 index 0000000..c422d73 --- /dev/null +++ b/content/math/tables.tex @@ -0,0 +1,22 @@ +\enlargethispage{0.2cm} +\begin{multicols*}{2} + \refstepcounter{subsection} + \subsectionmark{Tables} + \addcontentsline{toc}{subsection}{\protect\numberline{\thesubsection}Tables} + + \input{math/tables/binom} + \vfill + \input{math/tables/prime-composite} + \vfill + \input{math/tables/platonic} + \vfill + \input{math/tables/series} + + \columnbreak + + \input{math/tables/probability} + \vfill + \input{math/tables/stuff} + \vfill + \input{math/tables/nim} +\end{multicols*} diff --git a/content/math/tables/binom.tex b/content/math/tables/binom.tex new file mode 100644 index 0000000..878a6b0 --- /dev/null +++ b/content/math/tables/binom.tex @@ -0,0 +1,28 @@ +\begin{tabularx}{\linewidth}{|XXXX|} + \hline + \multicolumn{4}{|c|}{Binomialkoeffizienten} \\ + \hline + \multicolumn{4}{|c|}{ + $\frac{n!}{k!(n - k)!} \hfill=\hfill + \binom{n}{k} \hfill=\hfill + \binom{n}{n - k} \hfill=\hfill + \frac{n}{k}\binom{n - 1}{k - 1} \hfill=\hfill + \frac{n-k+1}{k}\binom{n}{k - 1} \hfill=\hfill + \binom{n - 1}{k} + \binom{n - 1}{k - 1} \hfill=\hfill + (-1)^k \binom{k - n - 1}{k} \hfill\approx\hfill + 2^{n} \cdot \frac{2}{\sqrt{2\pi n}}\cdot\exp\left(-\frac{2(x - \frac{n}{2})^2}{n}\right)$ + } \\ + \grayhline + + $\sum\limits_{k = 0}^n \binom{n}{k} = 2^n$ & + $\sum\limits_{k = 0}^n \binom{k}{m} = \binom{n + 1}{m + 1}$ & + $\sum\limits_{i = 0}^n \binom{n}{i}^2 = \binom{2n}{n}$ & + $\sum\limits_{k = 0}^n\binom{r + k}{k} = \binom{r + n + 1}{n}$\\ + + $\binom{n}{m}\binom{m}{k} = \binom{n}{k}\binom{n - k}{m - k}$ & + $\sum\limits_{k = 0}^n \binom{r}{k}\binom{s}{n - k} = \binom{r + s}{n}$ & + \multicolumn{2}{l|}{ + $\sum\limits_{i = 1}^n \binom{n}{i} F_i = F_{2n} \quad F_n = n\text{-th Fib.}$ + }\\ + \hline +\end{tabularx} diff --git a/content/math/tables/nim.tex b/content/math/tables/nim.tex new file mode 100644 index 0000000..8490d42 --- /dev/null +++ b/content/math/tables/nim.tex @@ -0,0 +1,96 @@ +\begin{tabularx}{\linewidth}{|p{0.37\linewidth}|X|} + \hline + \multicolumn{2}{|c|}{Nim-Spiele (\ding{182} letzter gewinnt (normal), \ding{183} letzter verliert)} \\ + \hline + Beschreibung & + Strategie \\ + \hline + + $M = [\mathit{pile}_i]$\newline + $[x] := \{1, \ldots, x\}$& + $\mathit{SG} = \oplus_{i = 1}^n \mathit{pile}_i$\newline + \ding{182} Nimm von einem Stapel, sodass $\mathit{SG}$ $0$ wird.\newline + \ding{183} Genauso. + Außer: Bleiben nur noch Stapel der Größe $1$, erzeuge ungerade Anzahl solcher Stapel.\\ + \hline + + $M = \{a^m \mid m \geq 0\}$ & + $a$ ungerade: $\mathit{SG}_n = n \% 2$\newline + $a$ gerade:\newline + $\mathit{SG}_n = 2$, falls $n \equiv a \bmod (a + 1) $\newline + $\mathit{SG}_n = n \% (a + 1) \% 2$, sonst.\\ + \hline + + $M_{\text{\ding{172}}} = \left[\frac{\mathit{pile}_i}{2}\right]$\newline + $M_{\text{\ding{173}}} = + \left\{\left\lceil\frac{\mathit{pile}_i}{2}\right\rceil,~ + \mathit{pile}_i\right\}$ & + \ding{172} + $\mathit{SG}_{2n} = n$, + $\mathit{SG}_{2n+1} = \mathit{SG}_n$\newline + \ding{173} + $\mathit{SG}_0 = 0$, + $\mathit{SG}_n = [\log_2 n] + 1$ \\ + \hline + + $M_{\text{\ding{172}}} = \text{Teiler von $\mathit{pile}_i$}$\newline + $M_{\text{\ding{173}}} = \text{echte Teiler von $\mathit{pile}_i$}$ & + \ding{172} + $\mathit{SG}_0 = 0$, + $\mathit{SG}_n = \mathit{SG}_{\text{\ding{173},n}} + 1$\newline + \ding{173} + $\mathit{ST}_1 = 0$, + $\mathit{SG}_n = \text{\#Nullen am Ende von $n_{bin}$}$\\ + \hline + + $M_{\text{\ding{172}}} = [k]$\newline + $M_{\text{\ding{173}}} = S$, ($S$ endlich)\newline + $M_{\text{\ding{174}}} = S \cup \{\mathit{pile}_i\}$ & + $\mathit{SG}_{\text{\ding{172}}, n} = n \bmod (k + 1)$\newline + \ding{182} Niederlage bei $\mathit{SG} = 0$\newline + \ding{183} Niederlage bei $\mathit{SG} = 1$\newline + $\mathit{SG}_{\text{\ding{174}}, n} = \mathit{SG}_{\text{\ding{173}}, n} + 1$\\ + \hline + + \multicolumn{2}{|l|}{ + Für jedes endliche $M$ ist $\mathit{SG}$ eines Stapels irgendwann periodisch. + } \\ + \hline + + \textsc{Moore}'s Nim:\newline + Beliebige Zahl von maximal $k$ Stapeln. & + \ding{182} + Schreibe $\mathit{pile}_i$ binär. + Addiere ohne Übertrag zur Basis $k + 1$. + Niederlage, falls Ergebnis gleich 0.\newline + \ding{183} + Wenn alle Stapel $1$ sind: + Niederlage, wenn $n \equiv 1 \bmod (k + 1)$. + Sonst wie in \ding{182}.\\ + \hline + + Staircase Nim:\newline + $n$ Stapel in einer Reihe. + Beliebige Zahl von Stapel $i$ nach Stapel $i-1$. & + Niederlage, wenn Nim der ungeraden Spiele verloren ist:\newline + $\oplus_{i = 0}^{(n - 1) / 2} \mathit{pile}_{2i + 1} = 0$\\ + \hline + + \textsc{Lasker}'s Nim:\newline + Zwei mögliche Züge:\newline + 1) Nehme beliebige Zahl.\newline + 2) Teile Stapel in zwei Stapel (ohne Entnahme).& + $\mathit{SG}_n = n$, falls $n \equiv 1,2 \bmod 4$\newline + $\mathit{SG}_n = n + 1$, falls $n \equiv 3 \bmod 4$\newline + $\mathit{SG}_n = n - 1$, falls $n \equiv 0 \bmod 4$\\ + \hline + + \textsc{Kayles}' Nim:\newline + Zwei mögliche Züge:\newline + 1) Nehme beliebige Zahl.\newline + 2) Teile Stapel in zwei Stapel (mit Entnahme).& + Berechne $\mathit{SG}_n$ für kleine $n$ rekursiv.\newline + $n \in [72,83]: \quad 4, 1, 2, 8, 1, 4, 7, 2, 1, 8, 2, 7$\newline + Periode ab $n = 72$ der Länge $12$.\\ + \hline +\end{tabularx} diff --git a/content/math/tables/numbers.tex b/content/math/tables/numbers.tex new file mode 100644 index 0000000..1dc9f38 --- /dev/null +++ b/content/math/tables/numbers.tex @@ -0,0 +1,59 @@ +\begin{expandtable} +\begin{tabularx}{\linewidth}{|l|X|} + \hline + \multicolumn{2}{|c|}{Berühmte Zahlen} \\ + \hline + \textsc{Fibonacci} & + $f(0) = 0 \quad + f(1) = 1 \quad + f(n+2) = f(n+1) + f(n)$ \\ + \grayhline + + \textsc{Catalan} & + $C_0 = 1 \qquad + C_n = \sum\limits_{k = 0}^{n - 1} C_kC_{n - 1 - k} = + \frac{1}{n + 1}\binom{2n}{n} = \frac{2(2n - 1)}{n+1} \cdot C_{n-1}$ \\ + \grayhline + + \textsc{Euler} I & + $\eulerI{n}{0} = \eulerI{n}{n-1} = 1 \qquad + \eulerI{n}{k} = (k+1) \eulerI{n-1}{k} + (n-k) \eulerI{n-1}{k-1} $ \\ + \grayhline + + \textsc{Euler} II & + $\eulerII{n}{0} = 1 \quad + \eulerII{n}{n} = 0 \quad$\\ + & $\eulerII{n}{k} = (k+1) \eulerII{n-1}{k} + (2n-k-1) \eulerII{n-1}{k-1}$ \\ + \grayhline + + \textsc{Stirling} I & + $\stirlingI{0}{0} = 1 \qquad + \stirlingI{n}{0} = \stirlingI{0}{n} = 0 \qquad + \stirlingI{n}{k} = \stirlingI{n-1}{k-1} + (n-1) \stirlingI{n-1}{k}$ \\ + \grayhline + + \textsc{Stirling} II & + $\stirlingII{n}{1} = \stirlingII{n}{n} = 1 \qquad + \stirlingII{n}{k} = k \stirlingII{n-1}{k} + \stirlingII{n-1}{k-1} = + \frac{1}{k!} \sum\limits_{j=0}^{k} (-1)^{k-j}\binom{k}{j}j^n$\\ + \grayhline + + \textsc{Bell} & + $B_1 = 1 \qquad + B_n = \sum\limits_{k = 0}^{n - 1} B_k\binom{n-1}{k} + = \sum\limits_{k = 0}^{n}\stirlingII{n}{k}$\\ + \grayhline + + \textsc{Partitions} & + $p(0,0) = 1 \quad + p(n,k) = 0 \text{ für } k > n \text{ oder } n \leq 0 \text{ oder } k \leq 0$ \\ + & $p(n,k) = p(n-k,k) + p(n-1,k-1)$\\ + \grayhline + + \textsc{Partitions} & + $f(0) = 1 \quad f(n) = 0~(n < 0)$ \\ + & $f(n)=\sum\limits_{k=1}^\infty(-1)^{k-1}f(n - \frac{k(3k+1)}{2})+\sum\limits_{k=1}^\infty(-1)^{k-1}f(n - \frac{k(3k-1)}{2})$\\ + + \hline +\end{tabularx} +\end{expandtable} diff --git a/content/math/tables/platonic.tex b/content/math/tables/platonic.tex new file mode 100644 index 0000000..f4ee554 --- /dev/null +++ b/content/math/tables/platonic.tex @@ -0,0 +1,39 @@ +\begin{tabularx}{\linewidth}{|X|CCCX|} + \hline + \multicolumn{5}{|c|}{Platonische Körper} \\ + \hline + Übersicht & Seiten & Ecken & Kanten & dual zu \\ + \hline + Tetraeder & 4 & 4 & 6 & Tetraeder \\ + Würfel/Hexaeder & 6 & 8 & 12 & Oktaeder \\ + Oktaeder & 8 & 6 & 12 & Würfel/Hexaeder\\ + Dodekaeder & 12 & 20 & 30 & Ikosaeder \\ + Ikosaeder & 20 & 12 & 30 & Dodekaeder \\ + \hline + \multicolumn{5}{|c|}{Färbungen mit maximal $n$ Farben (bis auf Isomorphie)} \\ + \hline + \multicolumn{3}{|l}{Ecken vom Oktaeder/Seiten vom Würfel} & + \multicolumn{2}{l|}{$(n^6 + 3n^4 + 12n^3 + 8n^2)/24$} \\ + + \multicolumn{3}{|l}{Ecken vom Würfel/Seiten vom Oktaeder} & + \multicolumn{2}{l|}{$(n^8 + 17n^4 + 6n^2)/24$} \\ + + \multicolumn{3}{|l}{Kanten vom Würfel/Oktaeder} & + \multicolumn{2}{l|}{$(n^{12} + 6n^7 + 3n^6 + 8n^4 + 6n^3)/24$} \\ + + \multicolumn{3}{|l}{Ecken/Seiten vom Tetraeder} & + \multicolumn{2}{l|}{$(n^4 + 11n^2)/12$} \\ + + \multicolumn{3}{|l}{Kanten vom Tetraeder} & + \multicolumn{2}{l|}{$(n^6 + 3n^4 + 8n^2)/12$} \\ + + \multicolumn{3}{|l}{Ecken vom Ikosaeder/Seiten vom Dodekaeder} & + \multicolumn{2}{l|}{$(n^{12} + 15n^6 + 44n^4)/60$} \\ + + \multicolumn{3}{|l}{Ecken vom Dodekaeder/Seiten vom Ikosaeder} & + \multicolumn{2}{l|}{$(n^{20} + 15n^{10} + 20n^8 + 24n^4)/60$} \\ + + \multicolumn{3}{|l}{Kanten vom Dodekaeder/Ikosaeder (evtl. falsch)} & + \multicolumn{2}{l|}{$(n^{30} + 15n^{16} + 20n^{10} + 24n^6)/60$} \\ + \hline +\end{tabularx} diff --git a/content/math/tables/prime-composite.tex b/content/math/tables/prime-composite.tex new file mode 100644 index 0000000..99b3348 --- /dev/null +++ b/content/math/tables/prime-composite.tex @@ -0,0 +1,26 @@ +\begin{tabularx}{\linewidth}{|r|r|rIr|rIrIr|C|} + \hline + \multicolumn{8}{|c|}{Important Numbers} \\ + \hline + $10^x$ & Highly Composite & \# Divs & $<$ Prime & $>$ Prime & \# Primes & primorial & \\ + \hline + 1 & 6 & 4 & $-3$ & $+1$ & 4 & 2 & \\ + 2 & 60 & 12 & $-3$ & $+1$ & 25 & 3 & \\ + 3 & 840 & 32 & $-3$ & $+9$ & 168 & 4 & \\ + 4 & 7\,560 & 64 & $-27$ & $+7$ & 1\,229 & 5 & \\ + 5 & 83\,160 & 128 & $-9$ & $+3$ & 9\,592 & 6 & \\ + 6 & 720\,720 & 240 & $-17$ & $+3$ & 78\,498 & 7 & \\ + 7 & 8\,648\,640 & 448 & $-9$ & $+19$ & 664\,579 & 8 & \\ + 8 & 73\,513\,440 & 768 & $-11$ & $+7$ & 5\,761\,455 & 8 & \\ + 9 & 735\,134\,400 & 1\,344 & $-63$ & $+7$ & 50\,847\,534 & 9 & \\ + 10 & 6\,983\,776\,800 & 2\,304 & $-33$ & $+19$ & 455\,052\,511 & 10 & \\ + 11 & 97\,772\,875\,200 & 4\,032 & $-23$ & $+3$ & 4\,118\,054\,813 & 10 & \\ + 12 & 963\,761\,198\,400 & 6\,720 & $-11$ & $+39$ & 37\,607\,912\,018 & 11 & \\ + 13 & 9\,316\,358\,251\,200 & 10\,752 & $-29$ & $+37$ & 346\,065\,536\,839 & 12 & \\ + 14 & 97\,821\,761\,637\,600 & 17\,280 & $-27$ & $+31$ & 3\,204\,941\,750\,802 & 12 & \\ + 15 & 866\,421\,317\,361\,600 & 26\,880 & $-11$ & $+37$ & 29\,844\,570\,422\,669 & 13 & \\ + 16 & 8\,086\,598\,962\,041\,600 & 41\,472 & $-63$ & $+61$ & 279\,238\,341\,033\,925 & 13 & \\ + 17 & 74\,801\,040\,398\,884\,800 & 64\,512 & $-3$ & $+3$ & 2\,623\,557\,157\,654\,233 & 14 & \\ + 18 & 897\,612\,484\,786\,617\,600 & 103\,680 & $-11$ & $+3$ & 24\,739\,954\,287\,740\,860 & 16 & \\ + \hline +\end{tabularx} diff --git a/content/math/tables/probability.tex b/content/math/tables/probability.tex new file mode 100644 index 0000000..f265d10 --- /dev/null +++ b/content/math/tables/probability.tex @@ -0,0 +1,27 @@ +\begin{tabularx}{\linewidth}{|LICIR|} + \hline + \multicolumn{3}{|c|}{ + Wahrscheinlichkeitstheorie ($A,B$ Ereignisse und $X,Y$ Variablen) + } \\ + \hline + $\E(X + Y) = \E(X) + \E(Y)$ & + $\E(\alpha X) = \alpha \E(X)$ & + $X, Y$ unabh. $\Leftrightarrow \E(XY) = \E(X) \cdot \E(Y)$\\ + + $\Pr[A \vert B] = \frac{\Pr[A \land B]}{\Pr[B]}$ & + $A, B$ disj. $\Leftrightarrow \Pr[A \land B] = \Pr[A] \cdot \Pr[B]$ & + $\Pr[A \lor B] = \Pr[A] + \Pr[B] - \Pr[A \land B]$ \\ + \hline +\end{tabularx} +\vfill +\begin{tabularx}{\linewidth}{|Xlr|lrX|} + \hline + \multicolumn{6}{|c|}{\textsc{Bertrand}'s Ballot Theorem (Kandidaten $A$ und $B$, $k \in \mathbb{N}$)} \\ + \hline + & $\#A > k\#B$ & $Pr = \frac{a - kb}{a + b}$ & + $\#B - \#A \leq k$ & $Pr = 1 - \frac{a!b!}{(a + k + 1)!(b - k - 1)!}$ & \\ + + & $\#A \geq k\#B$ & $Pr = \frac{a + 1 - kb}{a + 1}$ & + $\#A \geq \#B + k$ & $Num = \frac{a - k + 1 - b}{a - k + 1} \binom{a + b - k}{b}$ & \\ + \hline +\end{tabularx} diff --git a/content/math/tables/series.tex b/content/math/tables/series.tex new file mode 100644 index 0000000..3042781 --- /dev/null +++ b/content/math/tables/series.tex @@ -0,0 +1,33 @@ +\begin{tabularx}{\linewidth}{|XIXIXIX|} + \hline + \multicolumn{4}{|c|}{Reihen} \\ + \hline + $\sum\limits_{i = 1}^n i = \frac{n(n+1)}{2}$ & + $\sum\limits_{i = 1}^n i^2 = \frac{n(n + 1)(2n + 1)}{6}$ & + $\sum\limits_{i = 1}^n i^3 = \frac{n^2 (n + 1)^2}{4}$ & + $H_n = \sum\limits_{i = 1}^n \frac{1}{i}$ \\ + \grayhline + + $\sum\limits_{i = 0}^n c^i = \frac{c^{n + 1} - 1}{c - 1} \quad c \neq 1$ & + $\sum\limits_{i = 0}^\infty c^i = \frac{1}{1 - c} \quad \vert c \vert < 1$ & + $\sum\limits_{i = 1}^\infty c^i = \frac{c}{1 - c} \quad \vert c \vert < 1$ & + $\sum\limits_{i = 0}^\infty ic^i = \frac{c}{(1 - c)^2} \quad \vert c \vert < 1$ \\ + \grayhline + + \multicolumn{2}{|lI}{ + $\sum\limits_{i = 0}^n ic^i = \frac{nc^{n + 2} - (n + 1)c^{n + 1} + c}{(c - 1)^2} \quad c \neq 1$ + } & + \multicolumn{2}{l|}{ + $\sum\limits_{i = 1}^n iH_i = \frac{n(n + 1)}{2}H_n - \frac{n(n - 1)}{4}$ + } \\ + \grayhline + + \multicolumn{2}{|lI}{ + $\sum\limits_{i = 1}^n H_i = (n + 1)H_n - n$ + } & + \multicolumn{2}{l|}{ + $\sum\limits_{i = 1}^n \binom{i}{m}H_i = + \binom{n + 1}{m + 1} \left(H_{n + 1} - \frac{1}{m + 1}\right)$ + } \\ + \hline +\end{tabularx} diff --git a/content/math/tables/stuff.tex b/content/math/tables/stuff.tex new file mode 100644 index 0000000..3cf8b4c --- /dev/null +++ b/content/math/tables/stuff.tex @@ -0,0 +1,32 @@ +\begin{tabularx}{\linewidth}{|ll|} + \hline + \multicolumn{2}{|C|}{Verschiedenes} \\ + \hline + Türme von Hanoi, minimale Schirttzahl: & + $T_n = 2^n - 1$ \\ + + \#Regionen zwischen $n$ Geraden & + $\frac{n\left(n + 1\right)}{2} + 1$ \\ + + \#abgeschlossene Regionen zwischen $n$ Geraden & + $\frac{n^2 - 3n + 2}{2}$ \\ + + \#markierte, gewurzelte Bäume & + $n^{n-1}$ \\ + + \#markierte, nicht gewurzelte Bäume & + $n^{n-2}$ \\ + + \#Wälder mit $k$ gewurzelten Bäumen & + $\frac{k}{n}\binom{n}{k}n^{n-k}$ \\ + + \#Wälder mit $k$ gewurzelten Bäumen mit vorgegebenen Wurzelknoten& + $\frac{k}{n}n^{n-k}$ \\ + + Derangements & + $!n = (n - 1)(!(n - 1) + !(n - 2)) = \left\lfloor\frac{n!}{e} + \frac{1}{2}\right\rfloor$ \\ + & + $\lim\limits_{n \to \infty} \frac{!n}{n!} = \frac{1}{e}$ \\ + \hline +\end{tabularx} + diff --git a/content/math/tables/twelvefold.tex b/content/math/tables/twelvefold.tex new file mode 100644 index 0000000..18d3955 --- /dev/null +++ b/content/math/tables/twelvefold.tex @@ -0,0 +1,32 @@ +\begin{expandtable} +\begin{tabularx}{\linewidth}{|C|CICICIC|} + \hline + Bälle & identisch & verschieden & identisch & verschieden \\ + Boxen & identisch & identisch & verschieden & verschieden \\ + \hline + -- & + $p_k(n + k)$ & + $\sum\limits_{i = 0}^k \stirlingII{n}{i}$ & + $\binom{n + k - 1}{k - 1}$ & + $k^n$ \\ + \grayhline + + \makecell{Bälle pro\\Box $\geq 1$} & + $p_k(n)$ & + $\stirlingII{n}{k}$ & + $\binom{n - 1}{k - 1}$ & + $k! \stirlingII{n}{k}$ \\ + \grayhline + + \makecell{Bälle pro\\Box $\leq 1$} & + $[n \leq k]$ & + $[n \leq k]$ & + $\binom{k}{n}$ & + $n! \binom{k}{n}$ \\ + \hline + \multicolumn{5}{|l|}{ + $[\text{Bedingung}]$: \code{return Bedingung ? 1 : 0;} + } \\ + \hline +\end{tabularx} +\end{expandtable} diff --git a/content/math/transforms/andTransform.cpp b/content/math/transforms/andTransform.cpp new file mode 100644 index 0000000..1fd9f5c --- /dev/null +++ b/content/math/transforms/andTransform.cpp @@ -0,0 +1,8 @@ +void fft(vector<ll>& a, bool inv = false) { + int n = sz(a); + for (int s = 1; s < n; s *= 2) { + for (int i = 0; i < n; i += 2 * s) { + for (int j = i; j < i + s; j++) { + ll& u = a[j], &v = a[j + s]; + tie(u, v) = inv ? pair(v - u, u) : pair(v, u + v); +}}}} diff --git a/content/math/transforms/bitwiseTransforms.cpp b/content/math/transforms/bitwiseTransforms.cpp new file mode 100644 index 0000000..28561da --- /dev/null +++ b/content/math/transforms/bitwiseTransforms.cpp @@ -0,0 +1,12 @@ +void bitwiseConv(vector<ll>& a, bool inv = false) { + int n = sz(a); + for (int s = 1; s < n; s *= 2) { + for (int i = 0; i < n; i += 2 * s) { + for (int j = i; j < i + s; j++) { + ll& u = a[j], &v = a[j + s]; + tie(u, v) = inv ? pair(v - u, u) : pair(v, u + v); // AND + //tie(u, v) = inv ? pair(v, u - v) : pair(u + v, u); //OR + //tie(u, v) = pair(u + v, u - v); // XOR + }}} + //if (inv) for (ll& x : a) x /= n; // XOR (careful with MOD) +} diff --git a/content/math/transforms/fft.cpp b/content/math/transforms/fft.cpp new file mode 100644 index 0000000..2bd95b2 --- /dev/null +++ b/content/math/transforms/fft.cpp @@ -0,0 +1,23 @@ +using cplx = complex<double>; + +void fft(vector<cplx>& a, bool inv = false) { + int n = sz(a); + for (int i = 0, j = 1; j < n - 1; ++j) { + for (int k = n >> 1; k > (i ^= k); k >>= 1); + if (j < i) swap(a[i], a[j]); + } + static vector<cplx> ws(2, 1); + for (static int k = 2; k < n; k *= 2) { + ws.resize(n); + cplx w = polar(1.0, acos(-1.0) / k); + for (int i=k; i<2*k; i++) ws[i] = ws[i/2] * (i % 2 ? w : 1); + } + for (int s = 1; s < n; s *= 2) { + for (int j = 0; j < n; j += 2 * s) { + for (int k = 0; k < s; k++) { + cplx u = a[j + k], t = a[j + s + k]; + t *= (inv ? conj(ws[s + k]) : ws[s + k]); + a[j + k] = u + t; + a[j + s + k] = u - t; + if (inv) a[j + k] /= 2, a[j + s + k] /= 2; +}}}} diff --git a/content/math/transforms/fftMul.cpp b/content/math/transforms/fftMul.cpp new file mode 100644 index 0000000..660ed79 --- /dev/null +++ b/content/math/transforms/fftMul.cpp @@ -0,0 +1,15 @@ +vector<cplx> mul(vector<ll>& a, vector<ll>& b) { + int n = 1 << (__lg(sz(a) + sz(b) - 1) + 1); + vector<cplx> c(all(a)), d(n); + c.resize(n); + for (int i = 0; i < sz(b); i++) c[i] = {real(c[i]), b[i]}; + fft(c); + for (int i = 0; i < n; i++) { + int j = (n - i) & (n - 1); + cplx x = (c[i] + conj(c[j])) / cplx{2, 0}; //fft(a)[i]; + cplx y = (c[i] - conj(c[j])) / cplx{0, 2}; //fft(b)[i]; + d[i] = x * y; + } + fft(d, true); + return d; +} diff --git a/content/math/transforms/multiplyBitwise.cpp b/content/math/transforms/multiplyBitwise.cpp new file mode 100644 index 0000000..f7cf169 --- /dev/null +++ b/content/math/transforms/multiplyBitwise.cpp @@ -0,0 +1,8 @@ +vector<ll> mul(vector<ll> a, vector<ll> b) { + int n = 1 << (__lg(2 * max(sz(a), sz(b)) - 1)); + a.resize(n), b.resize(n); + bitwiseConv(a), bitwiseConv(b); + for (int i=0; i<n; i++) a[i] *= b[i]; // MOD? + bitwiseConv(a, true); + return a; +} diff --git a/content/math/transforms/multiplyFFT.cpp b/content/math/transforms/multiplyFFT.cpp new file mode 100644 index 0000000..0022d1f --- /dev/null +++ b/content/math/transforms/multiplyFFT.cpp @@ -0,0 +1,12 @@ +vector<ll> mul(vector<ll>& a, vector<ll>& b) { + int n = 1 << (__lg(sz(a) + sz(b) - 1) + 1); + vector<cplx> a2(all(a)), b2(all(b)); + a2.resize(n), b2.resize(n); + fft(a2), fft(b2); + for (int i=0; i<n; i++) a2[i] *= b2[i]; + fft(a2, true); + + vector<ll> ans(n); + for (int i=0; i<n; i++) ans[i] = llround(a2[i].real()); + return ans; +} diff --git a/content/math/transforms/multiplyNTT.cpp b/content/math/transforms/multiplyNTT.cpp new file mode 100644 index 0000000..806d124 --- /dev/null +++ b/content/math/transforms/multiplyNTT.cpp @@ -0,0 +1,8 @@ +vector<ll> mul(vector<ll> a, vector<ll> b) { + int n = 1 << (__lg(sz(a) + sz(b) - 1) + 1); + a.resize(n), b.resize(n); + ntt(a), ntt(b); + for (int i=0; i<n; i++) a[i] = a[i] * b[i] % mod; + ntt(a, true); + return a; +} diff --git a/content/math/transforms/ntt.cpp b/content/math/transforms/ntt.cpp new file mode 100644 index 0000000..ca605d3 --- /dev/null +++ b/content/math/transforms/ntt.cpp @@ -0,0 +1,23 @@ +constexpr ll mod = 998244353, root = 3; + +void ntt(vector<ll>& a, bool inv = false) { + int n = sz(a); + auto b = a; + ll r = inv ? powMod(root, mod - 2, mod) : root; + + for (int s = n / 2; s > 0; s /= 2) { + ll ws = powMod(r, (mod - 1) / (n / s), mod), w = 1; + for (int j = 0; j < n / 2; j += s) { + for (int k = j; k < j + s; k++) { + ll u = a[j + k], t = a[j + s + k] * w % mod; + b[k] = (u + t) % mod; + b[n/2 + k] = (u - t + mod) % mod; + } + w = w * ws % mod; + } + swap(a, b); + } + if (inv) { + ll div = powMod(n, mod - 2, mod); + for (auto& x : a) x = x * div % mod; +}} diff --git a/content/math/transforms/orTransform.cpp b/content/math/transforms/orTransform.cpp new file mode 100644 index 0000000..eb1da44 --- /dev/null +++ b/content/math/transforms/orTransform.cpp @@ -0,0 +1,8 @@ +void fft(vector<ll>& a, bool inv = false) { + int n = sz(a); + for (int s = 1; s < n; s *= 2) { + for (int i = 0; i < n; i += 2 * s) { + for (int j = i; j < i + s; j++) { + ll& u = a[j], &v = a[j + s]; + tie(u, v) = inv ? pair(v, u - v) : pair(u + v, u); +}}}} diff --git a/content/math/transforms/seriesOperations.cpp b/content/math/transforms/seriesOperations.cpp new file mode 100644 index 0000000..4743674 --- /dev/null +++ b/content/math/transforms/seriesOperations.cpp @@ -0,0 +1,56 @@ +vector<ll> poly_inv(const vector<ll>& a, int n) { + vector<ll> q = {powMod(a[0], mod-2, mod)}; + for (int len = 1; len < n; len *= 2){ + vector<ll> a2 = a, q2 = q; + a2.resize(2*len), q2.resize(2*len); + ntt(q2); + for (int j : {0, 1}) { + ntt(a2); + for (int i = 0; i < 2*len; i++) a2[i] = a2[i]*q2[i] % mod; + ntt(a2, true); + for (int i = 0; i < len; i++) a2[i] = 0; + } + for (int i = len; i < min(n, 2*len); i++) { + q.push_back((mod - a2[i]) % mod); + }} + return q; +} + +vector<ll> poly_deriv(vector<ll> a) { + for (int i = 1; i < sz(a); i++) + a[i-1] = a[i] * i % mod; + a.pop_back(); + return a; +} + +vector<ll> poly_integr(vector<ll> a) { + if (a.empty()) return {0}; + a.push_back(a.back() * powMod(sz(a), mod-2, mod) % mod); + for (int i = sz(a)-2; i > 0; i--) + a[i] = a[i-1] * powMod(i, mod-2, mod) % mod; + a[0] = 0; + return a; +} + +vector<ll> poly_log(vector<ll> a, int n) { + a = mul(poly_deriv(a), poly_inv(a, n)); + a.resize(n-1); + a = poly_integr(a); + return a; +} + +vector<ll> poly_exp(vector<ll> a, int n) { + vector<ll> q = {1}; + for (int len = 1; len < n; len *= 2) { + vector<ll> p = poly_log(q, 2*len); + for (int i = 0; i < 2*len; i++) + p[i] = (mod - p[i] + (i < sz(a) ? a[i] : 0)) % mod; + vector<ll> q2 = q; + q2.resize(2*len); + ntt(p), ntt(q2); + for (int i = 0; i < 2*len; i++) p[i] = p[i] * q2[i] % mod; + ntt(p, true); + for (int i = len; i < min(n, 2*len); i++) q.push_back(p[i]); + } + return q; +} diff --git a/content/math/transforms/xorTransform.cpp b/content/math/transforms/xorTransform.cpp new file mode 100644 index 0000000..f9d1d82 --- /dev/null +++ b/content/math/transforms/xorTransform.cpp @@ -0,0 +1,10 @@ +void fft(vector<ll>& a, bool inv = false) { + int n = sz(a); + for (int s = 1; s < n; s *= 2) { + for (int i = 0; i < n; i += 2 * s) { + for (int j = i; j < i + s; j++) { + ll& u = a[j], &v = a[j + s]; + tie(u, v) = pair(u + v, u - v); + }}} + if (inv) for (ll& x : a) x /= n; +} |