From b067d9880606143ac4b860beff4d5b02e7d349bd Mon Sep 17 00:00:00 2001 From: mzuenni Date: Mon, 27 Feb 2023 15:37:24 +0100 Subject: improved math --- math/gauss.cpp | 18 ++++----- math/lgsFp.cpp | 31 ++++++++------- math/linearSieve.cpp | 49 ++++++++++++++++++++++++ math/math.tex | 92 ++++++++++++++++++++++++++++++--------------- math/primeSieve.cpp | 5 +-- math/transforms/fft.cpp | 2 +- math/transforms/fftPerm.cpp | 8 ---- 7 files changed, 137 insertions(+), 68 deletions(-) create mode 100644 math/linearSieve.cpp delete mode 100644 math/transforms/fftPerm.cpp (limited to 'math') diff --git a/math/gauss.cpp b/math/gauss.cpp index 0c7ab03..3e3b7aa 100644 --- a/math/gauss.cpp +++ b/math/gauss.cpp @@ -1,8 +1,7 @@ -void normalLine(int n, int line) { +void normalLine(int line) { double factor = mat[line][line]; - for (int i = 0; i <= n; i++) { - mat[line][i] /= factor; -}} + for (double& x : mat[line]) x /= factor; +} void takeAll(int n, int line) { for (int i = 0; i < n; i++) { @@ -21,18 +20,17 @@ int gauss(int n) { if (abs(mat[j][i]) > abs(mat[i][i])) swappee = j; } swap(mat[i], mat[swappee]); - if (abs(mat[i][i]) > EPSILON) { - normalLine(n, i); + 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++) - if (abs(mat[i][j]) > EPSILON) allZero = false; - if (allZero && abs(mat[i][n]) > EPSILON) return INCONSISTENT; - if (allZero && abs(mat[i][n]) < EPSILON) return MULTIPLE; + 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/math/lgsFp.cpp b/math/lgsFp.cpp index 52442c6..13cdbd0 100644 --- a/math/lgsFp.cpp +++ b/math/lgsFp.cpp @@ -1,9 +1,7 @@ -void normalLine(int n, int line, ll p) { +void normalLine(int line, ll p) { ll factor = multInv(mat[line][line], p); - for (int i = 0; i <= n; i++) { - mat[line][i] *= factor; - mat[line][i] %= 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++) { @@ -11,17 +9,18 @@ void takeAll(int n, int line, ll p) { ll diff = mat[i][line]; for (int j = 0; j <= n; j++) { mat[i][j] -= (diff * mat[line][j]) % p; - mat[i][j] %= p; - if (mat[i][j] < 0) mat[i][j] += p; + mat[i][j] = (mat[i][j] + p) % p; }}} -void gauss(int n, ll p) { // nx(n+1)-Matrix, Körper F_p. - for (int line = 0; line < n; line++) { - int swappee = line; - while (swappee < n && mat[swappee][line] == 0) swappee++; - if (swappee == n) continue; - swap(mat[line], mat[swappee]); - normalLine(n, line, p); - takeAll(n, line, p); - // für Eindeutigkeit, Existenz etc. siehe LGS +void gauss(int n, ll mod) { // Nx(N+1)-Matrix, Körper F_p. + vector done(n, false); + for (int i = 0; i < n; i++) { + int j = i; + 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 diff --git a/math/linearSieve.cpp b/math/linearSieve.cpp new file mode 100644 index 0000000..e1fa7ff --- /dev/null +++ b/math/linearSieve.cpp @@ -0,0 +1,49 @@ +constexpr ll N = 10000000; +ll smallest[N], power[N], sieved[N]; +vector primes; + +//wird aufgerufen mit (p^k, p, k) für prime p +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 k % 2 ? pk : 1;} + +void sieve() { // O(N) + smallest[1] = power[1] = sieved[1] = 1; + for (ll i = 2; i < N; i++) { + if (smallest[i] == 0) { + primes.push_back(i); + for (ll pk = i, k = 1; pk < N; pk *= i, k++) { + smallest[pk] = i; + power[pk] = pk; + sieved[pk] = mu(pk, i, k); // Aufruf ändern! + }} + for (ll j = 0; i * primes[j] < N && primes[j] < smallest[i]; j++) { + ll k = i * primes[j]; + smallest[k] = power[k] = primes[j]; + sieved[k] = sieved[i] * sieved[primes[j]]; + } + if (i * smallest[i] < N && power[i] != i) { + ll k = i * smallest[i]; + smallest[k] = smallest[i]; + power[k] = power[i] * smallest[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! + }} + return res; +} diff --git a/math/math.tex b/math/math.tex index 3e5e3e0..e640814 100644 --- a/math/math.tex +++ b/math/math.tex @@ -149,23 +149,53 @@ sich alle Lösungen von $x^2-ny^2=c$ berechnen durch: \sourcecode{math/discreteNthRoot.cpp} \end{algorithm} +\begin{algorithm}{Linearessieb 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 endsprechenden Multiplikativen Funktion}{1} + + \method{naive}{Wert der endsprechenden 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} Funtkion:} + \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{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 Kann erweitert werden: Für jede Zahl den \{kleinsten, größten\} Primfaktor speichern, Liste von Primzahlen berechnen \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{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}-Funktion und \textsc{Möbius}-Inversion} +\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})$. @@ -181,38 +211,33 @@ sich alle Lösungen von $x^2-ny^2=c$ berechnen durch: 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)$. - \sourcecode{math/mobius.cpp} + %\sourcecode{math/mobius.cpp} \end{algorithm} -\columnbreak -\subsection{\textsc{Euler}sche $\boldsymbol{\varphi}$-Funktion} -\begin{itemize} - \item Zählt die relativ primen Zahlen $\leq n$. - - \item Multiplikativ: - $\gcd(a,b) = 1 \Longrightarrow \varphi(a) \cdot \varphi(b) = \varphi(ab)$ - - \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} -\sourcecode{math/phi.cpp} +%\columnbreak +%\subsection{\textsc{Euler}sche $\boldsymbol{\varphi}$-Funktion} +%\begin{itemize} +% \item Zählt die relativ primen Zahlen $\leq n$. +% +% \item Multiplikativ: +% $\gcd(a,b) = 1 \Longrightarrow \varphi(a) \cdot \varphi(b) = \varphi(ab)$ +% +% \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} +%\sourcecode{math/phi.cpp} -\subsection{LGS über $\boldsymbol{\mathbb{F}_p}$} -\method{gauss}{löst LGS}{n^3} -\sourcecode{math/lgsFp.cpp} - -\subsection{LGS über $\boldsymbol{\mathbb{R}}$} -\sourcecode{math/gauss.cpp} \begin{algorithm}{Numerisch Extremstelle bestimmen} \sourcecode{math/goldenSectionSearch.cpp} \end{algorithm} + \begin{algorithm}{Numerisch Integrieren, Simpsonregel} \sourcecode{math/simpson.cpp} \end{algorithm} @@ -231,11 +256,18 @@ sich alle Lösungen von $x^2-ny^2=c$ berechnen durch: %\lstinputlisting{math/ntt.cpp} %\textcolor{safeOrange}{$\blacksquare$} NTT code, %\textcolor{safeGreen}{$\blacksquare$} FFT code \sourcecode{math/transforms/all.cpp} - Für sehr viele transforms kann die Vertauschung vorberechnet werden: - \sourcecode{math/transforms/fftPerm.cpp} Multiplikation mit 2 transforms statt 3: (nur benutzten wenn nötig!) \sourcecode{math/transforms/fftMul.cpp} \end{algorithm} + +\subsection{LGS über $\boldsymbol{\mathbb{R}}$} +\method{gauss}{löst LGS}{n^3} +\sourcecode{math/gauss.cpp} + +\subsection{LGS über $\boldsymbol{\mathbb{F}_p}$} +\method{gauss}{löst LGS}{n^3} +\sourcecode{math/lgsFp.cpp} + \clearpage \subsection{Primzahlzählfunktion $\boldsymbol{\pi}$} diff --git a/math/primeSieve.cpp b/math/primeSieve.cpp index 68f7fcb..3898ab7 100644 --- a/math/primeSieve.cpp +++ b/math/primeSieve.cpp @@ -8,10 +8,9 @@ bool isPrime(ll x) { } void primeSieve() { - // i * i < N is enough for isPrime - for (ll i = 3; i < N; i += 2) { + for (ll i = 3; i < N; i += 2) {// i * i < N reicht für isPrime if (!isNotPrime[i / 2]) { - primes.push_back(i); + primes.push_back(i); // optional for (ll j = i * i; j < N; j+= 2 * i) { isNotPrime[j / 2] = 1; }}}} diff --git a/math/transforms/fft.cpp b/math/transforms/fft.cpp index 4540ed8..6c5f5a8 100644 --- a/math/transforms/fft.cpp +++ b/math/transforms/fft.cpp @@ -1,4 +1,4 @@ -using cplx = complex; // Eigene Implementierung ist schneller. +using cplx = complex; void fft(vector& a, bool inverse = 0) { int n = sz(a); diff --git a/math/transforms/fftPerm.cpp b/math/transforms/fftPerm.cpp deleted file mode 100644 index 2b6fb10..0000000 --- a/math/transforms/fftPerm.cpp +++ /dev/null @@ -1,8 +0,0 @@ -int perm[MAXN]; //perm[i] = j in Zeile 10 -void genPerm(int n){ - ull mask = ~0ull << (__lg(n) - 1); - for (int i = 0, j = 0; i < n; i++) { - perm[i] = j; //if (i < j) swap(a[i], a[j]); - ull y = mask >> __builtin_ctz(~i); - j ^= y & (n - 1); -}} -- cgit v1.2.3