(*
  Title:    Fast Fourier Transform
  Id:       $Id: Blatt7.thy,v 1.4 2005/06/28 12:33:27 obua Exp $
  Author:   Clemens Ballarin, started 12 April 2005
  Copyright 2005  Clemens Ballarin
*)

(*<*)theory Blatt7
imports Complex_Main
begin
<<<<<<< Blatt7.thy
(*>*)
=======
(*>*)

>>>>>>> 1.4
section {* Zeitplan *}

text {*
  \begin{itemize}
  \item[$\blacktriangleright$] Die komplexen Einheitswurzeln: Abgabe
  am 29.\ Juni.
  \item[$\blacktriangleright$] Diskrete Fourier-Transformation: Abgabe
  am 6.\ Juli.
  \item[$\blacktriangleright$] Schnelle Fourier-Transformation und
  Präsentation am 13.\ Juli.
  \end{itemize}
*}
  
section {* Vorbereitung *}

text {* In dieser Aufgabe werden Aussagen über die schnelle
  Fourier-Transformation über den komplexen Zahlen hergeleitet.  Viele
  der Lemmas über diesem Zahlenbereich, die Sie für diese Aufgabe
  benötigen, sind in Isabelle generisch formuliert.  Es ist deshalb
  zunächst erforderlich, dass Sie sich mit den Theorien für Körper
  \texttt{Ring\_and\_Field.thy} und Exponentiation \texttt{Power.thy}
  vertraut machen.  Ebenso mit dem Summationsoperator @{term setsum}
  aus der Theorie \texttt{Finite\_Set.thy}.

  Die Theorie der komplexen Zahlen befindet sich in einer speziellen
  Logik.  Sie müssen Isabelle mit der Option \texttt{-l HOL-Complex}
  starten. Benutzen Sie Ihre eigene Installation zuhause, so müssen Sie zuvor
  \texttt{HOL-Complex} compilieren. Wechseln Sie dazu in das Verzeichnis
  \texttt{isabelle/HOL}, wobei \texttt{isabelle} das Verzeichnis Ihrer Isabelle-Installation ist,
  und führen Sie dann den folgenden Befehl aus:
  \begin{center}
  \texttt{isatool make HOL-Complex}
  \end{center}
 *}

subsection {* Ringe und Körper *}

text {* Der Typ der komplexen Zahlen @{typ complex} ist Instanz der
  Typklassen @{text field}, @{text recpower} und @{text
  division_by_zero}.  Eine Beschreibung von Typklassen finden Sie im
  Isabelle/HOL-Tutorial, Kapitel~8.4.  Während @{text field} die gewöhnlichen
  Körper\-eigen\-schaften enthält, beschreibt @{text recpower} den
  Exponentiationsoperator @{text [quotes] "^"} und @{text
  division_by_zero} erleichtert den Umgang mit der Division.  Da in
  dieser Aufgabe letztlich nur Aussagen über die komplexen Zahlen
  bewiesen werden, können Sie Ihre Lemmas alternativ auch für den Typ
  @{text complex} einschränken. *}

text {*
  \begin{itemize}
  \item[$\blacktriangleright$] Zeigen Sie folgendes Lemma.

  Hinweis: Induktion mit dem Theorem @{thm [source] diff_induct}.
  \end{itemize}
 *}

lemma power_diff:
  assumes nz: "a ~= 0"
  shows "n <= m ==> (a::'a::{recpower, field}) ^ (m-n) = (a^m) / (a^n)"
oops

lemma nat_dvd_not_less:
  "0 < m ==> m < n ==> \<not> n dvd (m::nat)"
  by (unfold dvd_def) auto

lemma of_nat_cplx:
  "of_nat n = Complex (of_nat n) 0"
  by (induct n) (simp_all add: complex_one_def)


subsection {* Summen *}

text {* @{term "setsum f A"} ist für endlichen Mengen definiert als
  die Summe aller @{term "f i"} für @{term "i \<in> A"}.  Je nachdem, in
  welcher syntaktischen Form Funktion und Menge vorliegen, wählt
  Isabelle anschaulichere Syntax.  So erhält man @{term "\<Sum>i\<in>A. f i + g
  i"}, wenn die Funktion der $\lambda$-Ausdruck @{term "\<lambda>i. f i + g
  i"} ist.  Ist darüber hinaus
  die Menge das Interval @{text "{0..<n}"}, so verbessert sich die
  Syntax zu @{text "\<Sum>i = 0..<n. f i + g i"}.  Intervalle sind in
  \texttt{SetInterval.thy} definiert. *}

text {*
  \begin{itemize}
  \item[$\blacktriangleright$]
   Zeigen Sie folgende Distributivgesetze für Summen. 

  Hinweis: Fallunterscheidung, ob @{term A} endlich, dann
  Induktion über endlichen Mengen. 
  \end{itemize}
  *}

lemma setsum_left_distrib:
  "setsum f A * (r::'a::semiring_0_cancel) = (\<Sum>n\<in>A. f n * r)"
  oops

lemma setsum_divide_distrib:
  "setsum f A / (r::'a::field) = (\<Sum>n\<in>A. f n / r)"
  oops

text {*
  \begin{itemize}
  \item[$\blacktriangleright$] Zeigen Sie die Summationsformel für geometrische Reihen. 
  \end{itemize}
  *}

lemma geometric_sum:
  "x ~= 1 ==> (\<Sum>i=0..<n. x ^ i) =
  (x ^ n - 1) / (x - 1::'a::{field, recpower, division_by_zero})"
  oops

text {*
  \begin{itemize}
  \item[$\blacktriangleright$] Zeigen Sie, dass der erste Summand
  einer Summe abgespalten werden kann. 

  Hinweis: Ein ähnlicher Satz wurde mit @{thm [source]
  setsum_add_nat_ivl} bereits bewiesen. 
  \end{itemize}
  *}

lemma setsum_add_nat_ivl_singleton:
  assumes less: "m < (n::nat)"
  shows "f m + setsum f {m<..<n} = setsum f {m..<n}"
  oops


text {* Die folgenden Lemmas erleichtern das Experimentieren mit den
  Operatoren der Fouriertransformation, die später definiert
  werden. *}

lemma Ivl4:
  "{0..<4::nat} = {0, 1, 2, 3}"
proof -
  have "{0..<4::nat} = {0..<Suc (Suc (Suc (Suc 0)))}" by (simp add: nat_number)
  also have "... = {0, 1, 2, 3}"
    by (simp add: atLeastLessThanSuc nat_number insert_commute)
  finally show ?thesis .
qed

lemma Sum4:
  "(\<Sum>i=0..<4::nat. x i) = x 0 + x 1 + x 2 + x 3"
  by (simp add: Ivl4 nat_number)


subsection {* Summation und die Kongruenzregel *}

text {* Unter Umständen werden zum Vereinfachen des Terms unter einer
  Summe die Summationsgrenzen benötigt.  Der Simplifier kann diese
  Information nur nutzen, wenn mittels @{text "cong: setsum_cong"} die
  Kongruenzregel @{thm [display] setsum_cong} angegeben wird. *}

text {* Hierzu muss der Simplifier manchmal in den Power-Modus geschaltet
  werden.  Diese macht man mit dem folgenden Kommando.  *}

ML_setup {*
  simpset_ref() := simpset() setsubgoaler asm_full_simp_tac;
*}

text {*
  \begin{itemize}
  \item[$\blacktriangleright$] Zeigen Sie die folgenden beiden Lemmas. 
  \end{itemize}
  *}

lemma setsum_add_split_nat_ivl_singleton:
  assumes less: "m < (n::nat)"
    and g: "!!i. [| m < i; i < n |] ==> g i = f i"
  shows "f m + setsum g {m<..<n} = setsum f {m..<n}"
  oops

lemma setsum_add_split_nat_ivl:
  assumes le: "m <= (k::nat)" "k <= n"
    and g: "!!i. [| m <= i; i < k |] ==> g i = f i"
    and h: "!!i. [| k <= i; i < n |] ==> h i = f i"
  shows "setsum g {m..<k} + setsum h {k..<n} = setsum f {m..<n}"
  oops

text {* Da es den Simplifier erheblich verlangsamt, schalten wir den
  Power-Modus wieder aus. *}

ML_setup {*
  simpset_ref() := simpset() setsubgoaler asm_simp_tac;
*}

text {*
\textbf{Achtung Änderung:} Das Lemma \texttt{setsum-splice} samt Hilfslemmata muss nicht mehr bewiesen werden,
sondern ist vorgegeben.
*}

lemma ivl_splice_Un:
  "{0..<2*n::nat} = (op * 2 ` {0..<n}) \<union> ((%i. Suc (2*i)) ` {0..<n})"
  apply (unfold image_def Bex_def)
  apply auto
  apply arith
  done

lemma ivl_splice_Int:
  "(op * 2 ` {0..<n}) \<inter> ((%i. Suc (2*i)) ` {0..<n}) = {}"
  by auto arith

lemma double_inj_on:
  "inj_on (%i. 2*i::nat) A"
  by (simp add: inj_onI)

lemma Suc_double_inj_on:
  "inj_on (%i. Suc (2*i)) A"
  by (rule inj_onI) simp

lemma setsum_splice:
  "(\<Sum>i::nat = 0..<2*n. f i) = (\<Sum>i = 0..<n. f (2*i)) + (\<Sum>i = 0..<n. f (2*i+1))"
proof -
  have "(\<Sum>i::nat = 0..<2*n. f i) =
    setsum f (op * 2 ` {0..<n}) + setsum f ((%i. 2*i+1) ` {0..<n})"
    by (simp add: ivl_splice_Un ivl_splice_Int setsum_Un_disjoint)
  also have "... = (\<Sum>i = 0..<n. f (2*i)) + (\<Sum>i = 0..<n. f (2*i+1))"
    by (simp add: setsum_reindex [OF double_inj_on]
      setsum_reindex [OF Suc_double_inj_on] cong: setsum_cong)
  finally show ?thesis .
qed

text {* Lemma über die Vertauschung von innerer und äusserer Summe. *}

lemma swap_inj_on:
  "inj_on (%(i, j). (j, i)) (A \<times> B)"
  by (unfold inj_on_def) fast

lemma swap_product:
  "(%(i, j). (j, i)) ` (A \<times> B) = B \<times> A"
  by (simp add: split_def image_def) blast

lemma setsum_commute:
  "(\<Sum>i\<in>A. \<Sum>j\<in>B. f i j) = (\<Sum>j\<in>B. \<Sum>i\<in>A. f i j)"
proof (simp add: setsum_cartesian_product)
  have "(\<Sum>z\<in>A \<times> B. f (fst z) (snd z)) =
    (\<Sum>z\<in>(%(i, j). (j, i)) ` (A \<times> B). f (snd z) (fst z))"
    (is "?s = _")
    apply (simp add: setsum_reindex [where f = "%(i, j). (j, i)"] swap_inj_on
      cong: setsum_cong)
    apply (simp add: split_def)
    done
  also have "... = (\<Sum>z\<in>B \<times> A. f (snd z) (fst z))"
    (is "_ = ?t")
    apply (simp add: swap_product)
    done
  finally show "?s = ?t" .
qed


section {* Die komplexen Einheitswurzeln *}

text {* Mit diesem Abschnitt beginnt die Formalisierung der
  Fourier-Transformation.  Lesen Sie zunächst eine Einführung.
  Empfehlenswert ist das Kapitel in Cormen \textit{et al.},
  \emph{Introduction to Algorithms}, MIT Press, 1990. *}

text {* Einheitswurzeln werden definiert mit Hilfe der Funktion @{term
  cis} aus \texttt{Complex.thy}, die einem Winkel den zugehörigen
  Punkt auf dem Einheitskreis zuordnet.  Neben der Definition
  \begin{center} @{thm [source] cis_def}: @{thm cis_def} \end{center}
  ist
  \begin{center} @{thm [source] DeMoivre}: @{thm DeMoivre}
  \end{center}
  die wichtigste Eigenschaft. *}

constdefs
  root :: "nat => complex"
  "root n == cis (2*pi/(real (n::nat)))"


subsection {* Grundlagen *}

text {* Hinweis: Ein Reihe von Aussagen darüber, wann @{term sin} und
  @{term cos} die Werte @{text 0} oder @{text 1} annehmen, sind bei
  den folgenden Beweisen neben der De Moirve'schen Formel
  hilfreich. *}

text {* @{thm [source] sin_zero_abs_cos_one}: @{thm
  sin_zero_abs_cos_one} *}

lemma sin_periodic_pi_diff [simp]: "sin (x - pi) = - sin x"
  by (simp add: sin_diff)

lemma sin_cos_between_zero_two_pi:
  assumes 0: "0 < x" and pi: "x < 2 * pi"
  shows "sin x \<noteq> 0 \<or> cos x \<noteq> 1"
proof -
  { assume "0 < x" and "x < pi"
    then have "sin x \<noteq> 0" by (auto dest: sin_gt_zero_pi) }
  moreover
  { assume "x = pi"
    then have "cos x \<noteq> 1" by simp }
  moreover
  { assume pi1: "pi < x" and pi2: "x < 2 * pi"
    then have "0 < x - pi" and "x - pi < pi" by arith+
    then have "sin (x - pi) \<noteq> 0" by (auto dest: sin_gt_zero_pi)
    with pi1 pi2 have "sin x \<noteq> 0" by simp }
  ultimately show ?thesis using 0 pi by arith
qed

text {*
  \begin{itemize}
  \item[$\blacktriangleright$] Zeigen Sie die grundlegenden Eigenschaften der
  Einheitswurzeln.  Da unser Exponentiationsoperator nur für
  natürliche, nicht ganze Zahlen definiert ist, müssen manche Lemmas,
  insbesondere das Summationslemma, auch für inverse Einheitswurzeln
  formuliert werden. 
  \end{itemize}
  *}

lemma root_nonzero:
  "root n ~= 0"
  oops

lemma root_unity:
  "root n ^ n = 1"
  oops

lemma root_cancel:
  "0 < d ==> root (d * n) ^ (d * k) = root n ^ k"
  oops

lemma root_summation:
  assumes k: "0 < k" "k < n"
  shows "(\<Sum>i=0..<n. (root n ^ k) ^ i) = 0"
  oops

lemma root_summation_inv:
  assumes k: "0 < k" "k < n"
  shows "(\<Sum>i=0..<n. ((1 / root n) ^ k) ^ i) = 0"
  oops

lemma root0 [simp]:
  "root 0 = 1"
  oops

lemma root1 [simp]:
  "root 1 = 1"
  oops

lemma root2 [simp]:
  "root 2 = Complex -1 0"
  oops

lemma root4 [simp]:
  "root 4 = ii"
  oops


section {* Diskrete Fourier-Transformation *}

text {*
  \begin{itemize}
  \item[$\blacktriangleright$] Definieren Sie Operatoren @{text DFT}
  und @{text IDFT} für die
  diskrete Fourier-Trans\-for\-mation und ihr Inverses.  Stellen Sie
  Vektoren als Funktionen vom Typ @{text "nat => complex"} dar. 

  Hinweis: Berücksichtigen Sie, dass die Transformation für
  Vektoren bestimmter aber beliebiger Größe definiert werden soll. 
  \end{itemize}
  *}

text {*
  \begin{itemize}
  \item[$\blacktriangleright$]
  Als Vorbereitung auf den Korrektheitsbeweis der schnellen
  Fourier-Transformation, formulieren Sie Lemmas, mit denen die
  Transformation eines Vektors der Länge @{text "2 * 2 ^ k"} auf die
  Transformation von Vektoren der Länge @{text "2 ^ k"} zurückgeführt
  wird.
  \end{itemize}  *}

text {*
  \textbf{Achtung Änderung: Dies ist eine Zusatzaufgabe, die Sie überspringen können.}
  \begin{itemize}
  \item[$\blacktriangleright$] Zeigen Sie, dass @{text DFT} und @{text
  IDFT} zueinander
  invers sind.  D.h.\ genauer, dass die Anwendung von zunächst @{text
  DFT} und dann @{text IDFT} nur zu einer Skalierung um die Vektorgröße
  führt.
  \end{itemize}  *}


section {* Diskrete, schnelle Fourier-Transformation *}

text {*
  \begin{itemize}
  \item[$\blacktriangleright$] Definieren Sie mit Hilfe primitiver Rekursion Programme für
  die schnelle Fourier-Transformation @{text FFT} und das Inverse
  @{text IFFT}.  
  \end{itemize}
  *}

text {*
  \begin{itemize}
  \item[$\blacktriangleright$] Zeigen Sie, für Vektoren der Größe @{text "2 ^ k"}, die
  Äquivalenz von @{text DFT} und @{text FFT} und von @{text IDFT} und
  @{text IFFT}. 
  \end{itemize}
  *}

(*<*)
end
(*>*)