Lucas–Lehmer primality test by Coq

LucasLehmer.v

From mathcomp
Require Import all_ssreflect.
From mathcomp
Require Import ssralg ring_quotient finalg poly polydiv zmodp.
From mathcomp
Require Import all_fingroup cyclic.

Set Implicit Arguments.
Unset Strict Implicit.
Unset Printing Implicit Defensive.

Fixpoint lucas n :=
  if n is n'.+1 then (lucas n')^2 - 2 else 4.

Lemma lucas_leq : forall n, lucas n >= 4.
Proof.
  elim=> //=.
  move=> n H.
  rewrite ltn_subRL.
  apply: (@leq_trans (4^2)); first done.
  by rewrite leq_exp2r.
Qed.

Import GRing.Theory.
Open Scope ring_scope.

Section Sqrt3.
  Variable R : finComUnitRingType.
  Variable sqrt3 : R.
  Hypothesis sqrt3E: sqrt3 ^+ 2 = 3%:R.
  
  Definition alpha := 2%:R + sqrt3.
  Definition beta  := 2%:R - sqrt3.
  
  Lemma alpha_beta: alpha * beta = 1.
  Proof.
    rewrite /alpha /beta mulrC -subr_sqr.
    rewrite expr2 sqrt3E.
    rewrite mulrnAl mul1r.
    by rewrite -mulrnA -mulrnBr.
  Qed.
  
  Lemma alpha_beta_n n: alpha^+n * beta^+ n = 1.
  Proof.
    by rewrite -exprMn alpha_beta expr1n.
  Qed.

  Lemma lucasE n: (lucas n)%:R = alpha^+(2^n) + beta^+(2^n).
  Proof.
    elim: n.
    * rewrite !expr1 /lucas /alpha /beta.
      by rewrite addrCA addrK -mulrnDr.
    * move=> n /=.
      rewrite mulrnBr; last first.
      * apply: (@leq_trans (4^2)) => //.
        by rewrite leq_exp2r // lucas_leq.
      * rewrite expnS expn1 natrM.
        move=> ->.
        rewrite mulrDl !mulrDr.
        rewrite -!exprD expnS mulSn mul1n.
        rewrite -!addrA.
        congr (_ + _).
        rewrite alpha_beta_n.
        rewrite -exprMn mulrC exprMn.
        rewrite alpha_beta_n.
        rewrite addrCA addrCA.
        by rewrite [1 + _]addrA -mulr2n subrr addr0.
  Qed.
  
  Variable r : nat.
  Hypotheses lucas0: (lucas r)%:R = 0 :> R.
  
  Lemma alpha_r: alpha^+(2^r.+1) = -1.
  Proof.
    move: (lucasE r) => /eqP.
    rewrite lucas0 eq_sym addr_eq0.
    rewrite expnS mulSn mul1n exprD.
    move/eqP => {2}->.
    by rewrite mulrN alpha_beta_n.
  Qed.

  Lemma alpha_r_unit: alpha^+(2^r.+1) \is a GRing.unit.
  Proof.
    rewrite alpha_r.
    by apply: unitrN1.
  Qed.

  Lemma alpha_unit: alpha \is a GRing.unit.
  Proof.
    move: alpha_r_unit.
    rewrite unitrX_pos //.
    by rewrite expn_gt0.
  Qed.

  Definition ualpha := FinRing.unit R alpha_unit.

  Local Open Scope group_scope.
  
  Lemma ualpha_r: ualpha^+(2^r.+1) = FinRing.unit R (GRing.unitrN1 R).
  Proof.
    apply: val_inj => /=.
    rewrite FinRing.val_unitX /=.
    by rewrite alpha_r.
  Qed.
  
  Lemma ualpha_r_dvd: (#[ualpha] %| 2^r.+2)%N.
  Proof.
    rewrite order_dvdn.
    rewrite expnS mulSn mul1n.
    rewrite expgD ualpha_r.
    apply/eqP.
    apply: val_inj => /=.
    by rewrite mulrNN mul1r.
  Qed.
  
  Hypothesis nz2: 2%:R != 0 :> R.
  
  Lemma ualpha_r_ndvd: ~~(#[ualpha] %| 2^r.+1)%N.
  Proof.
    move: nz2.
    apply: contra.
    rewrite order_dvdn.
    rewrite ualpha_r.
    rewrite mulr2n.
    case/eqP => {2}<-.
    by rewrite subrr.
  Qed.
  
  Lemma ualpha_order: #[ualpha] >= 2^r.+2.
  Proof.
    have Hp: prime 2 by [].
    move: ualpha_r_ndvd => /dvdn_pfactor.
    move/(_ Hp) => Hn.
    move: ualpha_r_dvd.
    move/dvdn_pfactor.
    move/(_ Hp) => [] n H1 H2.
    rewrite leqNgt.
    apply/negP => H3.
    apply: Hn.
    exists n => //.
    move: H2 H3 => ->.
    by rewrite ltn_exp2l.
  Qed.
  
  Definition units_R := [set: {unit R}].

  Lemma Runit_card: #|units_R| >= 2^r.+2.
  Proof.
    move: ualpha_order.
    move/leq_trans.
    apply.
    rewrite orderE.
    apply: subset_leq_card.
    rewrite cycle_subG.
    by apply: in_setT.
  Qed.
  
  Lemma Rcard: 2^r.+2 < #|R|.
  Proof.
    apply: leq_ltn_trans.
    * by apply: Runit_card.
    * rewrite cardsT.
      rewrite card_sub.
      apply: proper_card.
      apply/properP.
      rewrite subset_predT.
      split=> //.
      exists 0.
      * by rewrite inE.
      * rewrite inE.
        by rewrite unitr0.
  Qed.
End Sqrt3.

Section AddSqrt3.
  Variable R : finComUnitRingType.
  
  Definition PR := {poly R}.
  
  Definition p : PR := 'X ^+ 2 - 3%:R%:P.
  
  Lemma pE: p = Poly [:: -3%:R; 0; 1].
  Proof.
    apply/polyP => i.
    rewrite /p coefD coefN coefC expr2 coefMX coefX coef_Poly.
    case: i => [|[|[|i]]] /=.
    * by rewrite add0r.
    * by rewrite add0r oppr0.
    * by rewrite subr0.
    * by rewrite subr0 nth_nil.
  Qed.
  
  Lemma pseqE: p = [:: -3%:R; 0; 1] :> seq R.
  Proof.
    rewrite pE /=.
    rewrite !polyseq_cons !nil_poly eqxx /=.
    by rewrite nil_poly polyseqC !oner_neq0.
  Qed.
  
  Lemma size_p: size p = 3.
  Proof.
    by rewrite pseqE.
  Qed.

  Lemma pne0: p != 0.
  Proof.
    by rewrite -size_poly_gt0 size_p.
  Qed.
  
  Lemma p_is_monic: p \is monic.
  Proof.
    by rewrite monicE /lead_coef size_p pseqE.
  Qed.

  Import Pdiv.Ring.
  Import Pdiv.RingMonic.
  
  Definition pI_pred : pred_class := rdvdp p.
  Fact pI_key : pred_key pI_pred. Proof. by []. Qed.
  Canonical pI := KeyedPred pI_key.
  
  Lemma pI_zmod_closed: zmod_closed pI.
  Proof.
    rewrite /pI.
    split.
    * rewrite unfold_in.
      apply/eqP.
      by exact: rmod0p.
    * move=> a b.
      rewrite !unfold_in.
      rewrite (rmodp_add p_is_monic).
      rewrite -mulN1r.
      rewrite -(rmodp_mulmr p_is_monic).
      move=> /eqP -> /eqP ->.
      by rewrite add0r mulr0 rmod0p.
  Qed.
  
  Canonical pI_opprPred := OpprPred pI_zmod_closed.
  Canonical pI_addrPred := AddrPred pI_zmod_closed.
  Canonical pI_zmodPred := ZmodPred pI_zmod_closed.
  
  Lemma pI_proper_ideal: proper_ideal pI.
  Proof.
    split.
    * rewrite unfold_in.
      by rewrite rmodp_small ?size_polyC ?size_p oner_neq0.
    * move=> a b.
      rewrite !unfold_in.
      rewrite -(rmodp_mulmr p_is_monic).
      move/eqP => ->.
      by rewrite mulr0 rmod0p.
  Qed.
  
  Canonical pI_idealr := MkIdeal pI_zmodPred pI_proper_ideal.
  
  Local Open Scope quotient_scope.
  
  Definition QpI : predArgType := {ideal_quot pI}.
  
  Definition pair_of_QpI (x : QpI) : (R * R) :=
    let: m := rmodp (generic_quotient.repr x) p in
      (m`_0, m`_1).

  Definition QpI_of_pair (x : R * R) : QpI :=
    let: (a, b) := x in
      let: y := Poly [:: a; b] in
        \pi y.
  
  Lemma pair_of_QpI_cancel: cancel pair_of_QpI QpI_of_pair.
  Proof.
    move=> x.
    apply/eqP.
    rewrite -[x]reprK !piE.
    rewrite /pair_of_QpI reprK.
    move: (generic_quotient.repr x) => {x} x.
    rewrite Quotient.equivE unfold_in.
    set z := rmodp x p.
    have ->: Poly [:: z`_0; z`_1] = z.
    * apply/polyP=> i.
      rewrite coef_Poly.
      case: i => [|[|i]] //.
      rewrite nth_default //.
      rewrite nth_default //.
      apply: (@leq_trans 2) => //.
      by rewrite -ltnS -size_p ltn_rmodp pne0.
    * rewrite (rmodp_add p_is_monic).
      rewrite /z.
      rewrite rmodp_small.
      * by rewrite -(rmodp_add p_is_monic) addrN rmod0p.
      * by rewrite ltn_rmodp pne0.
  Qed.
  
  Definition QpI_countMixin := CanCountMixin pair_of_QpI_cancel.
  Canonical QpI_countType := CountType QpI QpI_countMixin.
  Definition QpI_finMixin := CanFinMixin pair_of_QpI_cancel.
  Canonical QpI_finType := FinType QpI QpI_finMixin.

  Canonical Quotient.rquot_comRingType.
  
  Definition QpI_unitRingMixin :=
    Eval hnf in FinRing.Ring.UnitMixin [finRingType of QpI].
  Canonical QpI_unitRingType := UnitRingType QpI QpI_unitRingMixin.
  Canonical QpI_comUnitRingType := Eval hnf in [comUnitRingType of QpI].

  Canonical QpI_finComUnitRingType := Eval hnf in [finComUnitRingType of QpI].
  
  Definition RaS3 := QpI_finComUnitRingType.
  
  Definition sqrt3 : RaS3 := \pi 'X.
  
  Lemma sqrt3E: sqrt3 ^+ 2 = 3%:R.
  Proof.
    apply/eqP.
    rewrite expr2 /sqrt3.
    rewrite !piE Quotient.equivE.
    rewrite unfold_in.
    rewrite -(rmodpp p_is_monic).
    by rewrite -expr2 -polyC1 -!polyC_add -mulr2n -mulrS.
  Qed.
  
  Lemma card_RaS3: #|RaS3| <= #|R|^2.
  Proof.
    rewrite expnS expn1 -card_prod.
    rewrite -(@card_image _ _ pair_of_QpI).
    apply: subset_leq_card.
    apply: subset_predT.
    apply: can_inj.
    apply: pair_of_QpI_cancel.
  Qed.
  
  Lemma Ras3_0 n: (n%:R == 0 :> R) = (n%:R == 0 :> RaS3).
  Proof.
    rewrite -pi_oner /=.
    rewrite -raddfMn /=.
    rewrite piE Quotient.equivE.
    rewrite unfold_in subr0.
    rewrite -raddfMn /=.
    rewrite rmodp_small.
    * by rewrite polyC_eq0.
    * rewrite size_p size_polyC.
      by case: eqP.
  Qed.
End AddSqrt3.

Section Primality.
  Variable r : nat.

  Local Open Scope nat_scope.
  
  Definition Mp := 2^r.+2 - 1.
  Hypothesis LLT: Mp %| lucas r.

  Definition q := pdiv Mp.
  Definition Zq := [finComUnitRingType of 'Z_q].
  
  Definition ZqS3 := RaS3 Zq.
  Definition ZqS3_sqrt3 := sqrt3 Zq.
  
  Lemma Mp_gtn_1: Mp > 1.
  Proof.
    rewrite /Mp expnS.
    rewrite ltn_subRL /addn /=.
    apply: (@leq_ltn_trans (2 * 1)) => //.
    rewrite ltn_mul2l /=.
    apply: (@leq_ltn_trans (2^0)) => //.
    rewrite ltn_exp2l //.
  Qed.
  
  Lemma q_gtn_1: q > 1.
  Proof.
    rewrite /q.
    move: (@pdiv_prime Mp).
    move/(_ Mp_gtn_1).
    apply: prime_gt1.
  Qed.
  
  Lemma Mp_odd: odd Mp.
  Proof.
    rewrite /Mp.
    rewrite odd_sub /=.
    rewrite odd_exp //=.
    rewrite expn_gt0 //=.
  Qed.
  
  Lemma q_odd: odd q.
  Proof.
    rewrite /q.
    apply: dvdn_odd.
    apply: pdiv_dvd.
    apply: Mp_odd.
  Qed.
  
  Lemma Zq3_card: 2^r.+2 < #|ZqS3|.
  Proof.
    apply: (Rcard (sqrt3E _)).
    * apply/eqP.
      rewrite -Ras3_0.
      move: LLT.
      move/dvdnP => [k ->].
      move: (pdiv_dvd Mp).
      move/dvdnP => [t ->].
      rewrite !natrM.
      rewrite char_Zp ?q_gtn_1 //.
      by rewrite mulr0 mulr0.
    * rewrite -Ras3_0.
      rewrite eqE /=.
      rewrite Zp_cast ?q_gtn_1 //.
      rewrite [1 %% q]modn_small ?q_gtn_1 //.
      rewrite modn_small //.
      rewrite ltn_neqAle /addn /= q_gtn_1 andbT.
      by apply/eqP => H; move: H q_odd => <-.
  Qed.
  
  Lemma exp_r_ltn: 2^r.+2 < q^2.
  Proof.
    apply: leq_trans.
    apply: Zq3_card.
    apply: leq_trans.
    apply: card_RaS3.
    rewrite -(@card_Zp q).
    rewrite /Zq /Zp /=.
    rewrite q_gtn_1 cardsE cardT //=.
    apply: (@leq_trans 2) => //.
    apply: q_gtn_1.
  Qed.
  
  Lemma prime_Mp: prime Mp.
  Proof.
    apply: ltn_pdiv2_prime.
    apply: (@leq_trans 2) => //.
    apply: Mp_gtn_1.
    apply: (@leq_ltn_trans (2^r.+2)).
    rewrite /Mp.
    apply: leq_subr.
    apply: exp_r_ltn.
  Qed.
End Primality.

Local Open Scope nat_scope.

Theorem LucasLehmerPrimalityTest:
  forall p, Mp p %| lucas p -> prime (Mp p).
Proof.
  apply: prime_Mp.
Qed.