improved math

7 files changed, 137 insertions, 68 deletions

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<bool> 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<ll> 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<double>; // Eigene Implementierung ist schneller.
+using cplx = complex<double>;
 
 void fft(vector<cplx>& 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);
-}}