Module strassen

Require Import Ncring Ncring_tac.
Require Import ssreflect ssrbool ssrfun eqtype ssrnat seq choice fintype.
Require Import div finfun bigop prime binomial ssralg finset fingroup finalg.
Require Import perm zmodp matrix seqmatrix cssralg.


Section Strassen_exp2.

Local Open Scope nat_scope.

Definition exp2 := fun i => iter i (fun x => x + x) 1.

Local Open Scope ring_scope.

Variable (R : ringType).


Section matrixRing.

Variable n : nat.

Global Instance Zops : @Ring_ops (matrix_zmodType R n n) 0%R
  (scalar_mx 1) (@addmx R _ _) mulmx (fun M N => M - N) (@oppmx R _ _) eq.

Global Instance Zr: (@Ring _ _ _ _ _ _ _ _ Zops).

End matrixRing.


Fixpoint Strassen_exp2 {k} :=
  match k return let M := 'M[R]_(exp2 k) in M -> M -> M with
  | 0%N => fun A B => A *m B
  | l.+1 => fun A B =>
    if l <= 5 then A *m B else
    let A11 := ulsubmx A in
    let A12 := ursubmx A in
    let A21 := dlsubmx A in
    let A22 := drsubmx A in
    let B11 := ulsubmx B in
    let B12 := ursubmx B in
    let B21 := dlsubmx B in
    let B22 := drsubmx B in
    let X := A11 - A21 in
    let Y := B22 - B12 in
    let C21 := Strassen_exp2 X Y in
    let X := A21 + A22 in
    let Y := B12 - B11 in
    let C22 := Strassen_exp2 X Y in
    let X := X - A11 in
    let Y := B22 - Y in
    let C12 := Strassen_exp2 X Y in
    let X := A12 - X in
    let C11 := Strassen_exp2 X B22 in
    let X := Strassen_exp2 A11 B11 in
    let C12 := X + C12 in
    let C21 := C12 + C21 in
    let C12 := C12 + C22 in
    let C22 := C21 + C22 in
    let C12 := C12 + C11 in
    let Y := Y - B21 in
    let C11 := Strassen_exp2 A22 Y in
    let C21 := C21 - C11 in
    let C11 := Strassen_exp2 A12 B21 in
    let C11 := X + C11 in
    block_mx C11 C12 C21 C22
  end.

Import GRing.Theory.

Lemma Strassen_exp2P : forall n (M N : 'M[R]_(exp2 n)), (Strassen_exp2 M N) = (M *m N).

Variable CR : cringType R.

Fixpoint Strassen_exp2_seqmx k :=
  match k return let M := seqmatrix CR in M -> M -> M with
  | 0%N => fun A B => mulseqmx (exp2 k) (exp2 k) A B
  | l.+1 => fun A B =>
    if l <= 5 then mulseqmx (exp2 k) (exp2 k) A B else
    let off := exp2 l in
    let A11 := ulsubseqmx off off A in
    let A12 := ursubseqmx off off A in
    let A21 := dlsubseqmx off off A in
    let A22 := drsubseqmx off off A in
    let B11 := ulsubseqmx off off B in
    let B12 := ursubseqmx off off B in
    let B21 := dlsubseqmx off off B in
    let B22 := drsubseqmx off off B in
    let X := subseqmx A11 A21 in
    let Y := subseqmx B22 B12 in
    let C21 := Strassen_exp2_seqmx l X Y in
    let X := addseqmx A21 A22 in
    let Y := subseqmx B12 B11 in
    let C22 := Strassen_exp2_seqmx l X Y in
    let X := subseqmx X A11 in
    let Y := subseqmx B22 Y in
    let C12 := Strassen_exp2_seqmx l X Y in
    let X := subseqmx A12 X in
    let C11 := Strassen_exp2_seqmx l X B22 in
    let X := Strassen_exp2_seqmx l A11 B11 in
    let C12 := addseqmx X C12 in
    let C21 := addseqmx C12 C21 in
    let C12 := addseqmx C12 C22 in
    let C22 := addseqmx C21 C22 in
    let C12 := addseqmx C12 C11 in
    let Y := subseqmx Y B21 in
    let C11 := Strassen_exp2_seqmx l A22 Y in
    let C21 := subseqmx C21 C11 in
    let C11 := Strassen_exp2_seqmx l A12 B21 in
    let C11 := addseqmx X C11 in
    block_seqmx C11 C12 C21 C22
  end.

Variable k : nat.

Lemma Strassen_exp2_seqmxP :
  {morph (@seqmx_of_mx _ CR (exp2 k) (exp2 k)) :
    M N / Strassen_exp2 M N >-> Strassen_exp2_seqmx k M N}.

End Strassen_exp2.


Section Strassen.

Local Open Scope ring_scope.

Variable R : ringType.
Let K := 32%positive.

Local Coercion nat_of_pos : positive >-> nat.

Lemma addpp p : xO p = (p + p)%N :> nat.

Lemma addSpp p : xI p = (p + p).+1%N :> nat.

Lemma addp1 p : xI p = (xO p + 1)%N :> nat.

Lemma addpp1 p : xI p = (p + p + 1)%N :> nat.

Lemma lt0p : forall p : positive, 0 < p.

Definition Strassen_step {p : positive} (A B : 'M[R]_(p + p)) f : 'M[R]_(p + p) :=
  let A11 := ulsubmx A in
  let A12 := ursubmx A in
  let A21 := dlsubmx A in
  let A22 := drsubmx A in
  let B11 := ulsubmx B in
  let B12 := ursubmx B in
  let B21 := dlsubmx B in
  let B22 := drsubmx B in
  let X := A11 - A21 in
  let Y := B22 - B12 in
  let C21 := f X Y in
  let X := A21 + A22 in
  let Y := B12 - B11 in
  let C22 := f X Y in
  let X := X - A11 in
  let Y := B22 - Y in
  let C12 := f X Y in
  let X := A12 - X in
  let C11 := f X B22 in
  let X := f A11 B11 in
  let C12 := X + C12 in
  let C21 := C12 + C21 in
  let C12 := C12 + C22 in
  let C22 := C21 + C22 in
  let C12 := C12 + C11 in
  let Y := Y - B21 in
  let C11 := f A22 Y in
  let C21 := C21 - C11 in
  let C11 := f A12 B21 in
  let C11 := X + C11 in
  block_mx C11 C12 C21 C22.

Lemma Strassen_stepP (p : positive) (A B : 'M[R]_(p + p)) f :
  f =2 mulmx -> Strassen_step A B f = A *m B.

Fixpoint Strassen {n : positive} {struct n} :=
  match n return let M := 'M[R]_n in M -> M -> M with
  | xH => fun M N => M *m N
  | xO p => fun A B =>
    if p <= K then A *m B else
    let A := castmx (addpp p,addpp p) A in
    let B := castmx (addpp p,addpp p) B in
    castmx (esym (addpp p),esym (addpp p)) (Strassen_step A B Strassen)
  | xI p => fun M N =>
    if p <= K then M *m N else
    let M := castmx (addpp1 p, addpp1 p) M in
    let N := castmx (addpp1 p, addpp1 p) N in
    let M11 := ulsubmx M in
    let M12 := ursubmx M in
    let M21 := dlsubmx M in
    let M22 := drsubmx M in
    let N11 := ulsubmx N in
    let N12 := ursubmx N in
    let N21 := dlsubmx N in
    let N22 := drsubmx N in
    let C := Strassen_step M11 N11 Strassen + M12 *m N21 in
    let R12 := M11 *m N12 + M12 *m N22 in
    let R21 := M21 *m N11 + M22 *m N21 in
    let R22 := M21 *m N12 + M22 *m N22 in
    castmx (esym (addpp1 p), esym (addpp1 p)) (block_mx C R12 R21 R22)
end.

Lemma mulmx_cast {R' : ringType} {m n p m' n' p'} {M:'M[R']_(m,p)} {N:'M_(p,n)}
  {eqm : m = m'} (eqp : p = p') {eqn : n = n'} :
  castmx (eqm,eqn) (M *m N) = castmx (eqm,eqp) M *m castmx (eqp,eqn) N.

Lemma StrassenP : forall (n : positive) (M N : 'M[R]_n), Strassen M N = M *m N.

Variable CR : cringType R.

Definition Strassen_step_seqmx (p : positive) (A B : seqmatrix CR) f : seqmatrix CR :=
  let A11 := ulsubseqmx p p A in
  let A12 := ursubseqmx p p A in
  let A21 := dlsubseqmx p p A in
  let A22 := drsubseqmx p p A in
  let B11 := ulsubseqmx p p B in
  let B12 := ursubseqmx p p B in
  let B21 := dlsubseqmx p p B in
  let B22 := drsubseqmx p p B in
  let X := subseqmx A11 A21 in
  let Y := subseqmx B22 B12 in
  let C21 := f p X Y in
  let X := addseqmx A21 A22 in
  let Y := subseqmx B12 B11 in
  let C22 := f p X Y in
  let X := subseqmx X A11 in
  let Y := subseqmx B22 Y in
  let C12 := f p X Y in
  let X := subseqmx A12 X in
  let C11 := f p X B22 in
  let X := f p A11 B11 in
  let C12 := addseqmx X C12 in
  let C21 := addseqmx C12 C21 in
  let C12 := addseqmx C12 C22 in
  let C22 := addseqmx C21 C22 in
  let C12 := addseqmx C12 C11 in
  let Y := subseqmx Y B21 in
  let C11 := f p A22 Y in
  let C21 := subseqmx C21 C11 in
  let C11 := f p A12 B21 in
  let C11 := addseqmx X C11 in
  block_seqmx C11 C12 C21 C22.

Lemma Strassen_step_seqmxP (p : positive) f fI :
  {morph (@seqmx_of_mx _ CR p p) : M N / f M N >-> fI p M N} ->
  {morph (@seqmx_of_mx _ CR (p + p) (p + p)) :
    M N / Strassen_step M N f >-> Strassen_step_seqmx p M N fI}.

Fixpoint Strassen_seqmx (n : positive) :=
  match n return let M := seqmatrix CR in M -> M -> M with
  | xH => fun A B => mulseqmx n n A B
  | xO p => fun A B =>
    if p <= K then mulseqmx n n A B else
    Strassen_step_seqmx p A B Strassen_seqmx
  | xI p => fun M N =>
    if p <= K then mulseqmx n n M N else
    let off := xO p in
    let M11 := ulsubseqmx off off M in
    let M12 := ursubseqmx off off M in
    let M21 := dlsubseqmx off off M in
    let M22 := drsubseqmx off off M in
    let N11 := ulsubseqmx off off N in
    let N12 := ursubseqmx off off N in
    let N21 := dlsubseqmx off off N in
    let N22 := drsubseqmx off off N in
    let R12 := addseqmx (mulseqmx off 1 M11 N12) (mulseqmx 1 1 M12 N22) in
    let R21 := addseqmx (mulseqmx off off M21 N11) (mulseqmx 1 off M22 N21) in
    let R22 := addseqmx (mulseqmx off 1 M21 N12) (mulseqmx 1 1 M22 N22) in
    let C := addseqmx (Strassen_step_seqmx p M11 N11 Strassen_seqmx) (mulseqmx 1 off M12 N21) in
    block_seqmx C R12 R21 R22
  end.

Lemma Strassen_seqmxP : forall (p : positive),
  {morph (@seqmx_of_mx _ CR p p) : M N / Strassen M N >-> Strassen_seqmx p M N}.

End Strassen.