## Stream: Coq users

### Topic: Any idea how to prove this?

#### Lessness (Aug 08 2023 at 08:39):

I'm stuck with the proof of the last theorem here: https://raw.githubusercontent.com/LessnessRandomness/PE/main/temporary.v

Can someone help with some fresh ideas? If you have free time, of course.

Thank you in advance.

#### Paolo Giarrusso (Aug 08 2023 at 11:47):

IIUC, what you need to prove (restating things a bit) is that if `i + 1` is prime and divides `n`, then `i + 1` divides `n / i!` where `i!` is factorial, and `n / i!` is `repeated_repeated_div`? Sorry, I misunderstood, but maybe the sketch could be adapted.

#### Paolo Giarrusso (Aug 08 2023 at 11:48):

It _seems_ that _almost_ follows because `i + 1` is coprime with any `j <= i` (since it's prime) , so it's coprime with `i!`, and of the following fact:
If `a` and `b` are coprime, `a | n` and `b | n`, then `a | (n / b)`.

#### Paolo Giarrusso (Aug 08 2023 at 11:57):

in fact, it seems `repeated_repeated_div i n` divides `n` by all integers `j < i` _only as many times as the division is exact_. But then the ideas above should work even better.

#### Lessness (Aug 08 2023 at 14:47):

It's done! :partying_face:

Thank you very much for the hints! :)

#### Lessness (Aug 08 2023 at 23:20):

I will just put it here (to avoid multiplying topics unnecessarily).
https://github.com/LessnessRandomness/PE/blob/main/Verif_EulerProject3.v

I have smth wrong with loop invariant (at line 640), I believe, as I reach some unprovable state later in the proof. Will sleep on it and, hopefully, tomorrow fix it.

But if someone has some hints or corrections, I would be very grateful.
(Corresponding C source here: https://github.com/LessnessRandomness/PE/blob/main/EulerProject3.c )

Thank you!

(deleted)

#### Lessness (Aug 09 2023 at 12:29):

Made some changes. And got some more problems. :|

At line 762 I have two `forward_if`s that I should fix (I believe mistakes in them are cause for failing `forward`s after them).

If someone helps with them, I would be grateful.
Taking rest for now and not changing anything for some time.

#### Lessness (Aug 11 2023 at 17:16):

Still no idea about the loop invariant in the `find_proof` lemma of the `Verif_EulerProject3.v`file.
Cleaned up a bit, added some more theorems etc.

Is it possible to place bounty on the loop invariant? :D

#### Lessness (Aug 11 2023 at 19:48):

Some more mathy challenges dropped.

More specifically, the lemma

``````Theorem repeated_div_thm5 n (H: 1 <= n) a b (Ha: prime a) (Hb: prime b):
snd (repeated_div a (snd (repeated_div b n))) = snd (repeated_div b (snd (repeated_div a n))).
``````

where hypotheses `prime a` and `prime b` probably should be replaced with `rel_prime a b` and (?) `2 <= a`, `2 <= b`.

Again, if someone has some hints, I will be grateful! :)

#### Paolo Giarrusso (Aug 11 2023 at 21:31):

The lasu change makes sense if rel_prime means coprime

#### Lessness (Aug 11 2023 at 21:32):

Yes, it does mean coprime.

#### Lessness (Aug 11 2023 at 21:33):

``````rel_prime = fun a b : Z => Zis_gcd a b 1
``````

#### Lessness (Aug 13 2023 at 15:24):

Right now I'm just spamming theorems with hope that smth will stick. There is some progress. :sweat_smile:

And I have such seemingly easy lemma to prove:
``` Theorem aux6 n a b (H1: rel_prime a b) (H2: (a | n)) (H3: (b | n)): (a * b | n). ```

Thanks for any hints!

#### Paolo Giarrusso (Aug 13 2023 at 15:48):

Do you mean this one is hard? I wonder if you have available the theory of prime factorizations

#### Paolo Giarrusso (Aug 13 2023 at 15:49):

But generally, I feel this code this definitely calls for a pen-and-paper proof before formalization.

#### Paolo Giarrusso (Aug 13 2023 at 16:12):

For instance, `aux6` becomes easy given a theory of prime factorizations; I recall from discussions here taht math-comp should include some building blocks, the stdlib might include them too.

Say `primes i` is the enumeration of primes, and the factorization `fact n` of a number `n : nat` is pair `(max_i : n, f : nat -> nat)` such that `\prod_{i : nat := 0 .. max_i} (primes i)^(f i) = n`.

Using nontrivial properties of prime factorizations, `aux6` reduces to the following lemma on _results_ of prime factorizations.

``````From Coq Require Import Utf8 Lia.
Definition factor : nat -> nat * (nat -> nat).
Admitted.

Lemma aux6' n a b :
let (max_a, fact_a) := factor a in
let (max_b, fact_b) := factor b in
let (max_n, fact_n) := factor n in
(∀ i,
(fact_a i = 0 \/ fact_b i = 0) (* rel_prime a b *)
/\ fact_a i <= fact_n i (* (a | n) *)
/\ fact_b i <= fact_n i  (* (b | n) *)) ->
(∀ i, fact_a i + fact_b i <= fact_n i).
Proof.
destruct
(factor a) as [max_a fact_a],
(factor b) as [max_b fact_b],
(factor n) as [max_n fact_n].
intros Hyps i; specialize (Hyps i) as ([-> | ->] & AN & BN); lia.
Qed.
``````

#### Paolo Giarrusso (Aug 13 2023 at 16:14):

a real formalization might encode factorizations as finite lists — but lookups in those lists can use `nth 0 factors`, using 0 as default, to get similarly easy reasoning :-).

#### Lessness (Aug 13 2023 at 17:23):

I'm trying to make it work without theory of factorizations. :sweat_smile: Just because it's fun. :sweat_smile:

Some year ago (approximately) I had proved the FTA (the fundamental theorem or arithmetics) and was thinking about doing the proofs for the PE 3rd exercise using factorizations. If I reach dead end with my current approach, I will do as you say, because it is definitely more reasonable.

#### Lessness (Aug 13 2023 at 17:33):

Ow, aux6 was super easy actually.

#### Lessness (Aug 13 2023 at 17:42):

Now I have left only this:
`n = n / a ^ fst (repeated_div a n) / b ^ fst (repeated_div b n) * a ^ fst (repeated_div a n) * b ^ fst (repeated_div b n)`
It seems hard but maybe I will be lucky.

And then I can return to the proof of the `find` function's correctness (to finding correct loop invariant and fixing problems with `forward_if`).

#### Lessness (Aug 13 2023 at 18:46):

Ok, it is done too.

#### Lessness (Aug 13 2023 at 21:55):

Math part is done, now only "VST part" is left.

Loop invariant with its following proofs is one thing that's still left. And there are two `admit` at the end, where I can't get two `forward_if` to work properly.

Thanks in advance, if someone chooses to help.

#### Lessness (Aug 14 2023 at 10:24):

The C program is here: https://github.com/LessnessRandomness/PE/blob/main/EulerProject3.c
The loop in the question is inside `find` function.

And my current try to write loop invariant looks like this (it's commented out in the body of `find_proof`:

``````forward_loop (
EX (i: Z),
PROP (let W := repeated_repeated_div (6 * i + 3) n H1 in
0 <= i /\ if prime_dec W then (6 * i + 5) <= W else
if Z.eq_dec W 1 then True else (6 * i + 5) * (6 * i + 5) <= W)
LOCAL (temp _n (Vlong (Int64.repr (repeated_repeated_div (6 * i + 3) n H1)));
temp _i (Vlong (Int64.repr (6 * i + 5))); gvars gv)
SEP (data_at Ews tulong (Vlong (Int64.repr (value_of_highest (6 * i + 3) n H1))) (gv _highest))
).
``````

#### Lessness (Aug 15 2023 at 12:41):

The correct loop invariant is, probably, such:

``````forward_loop (
EX (i: Z),
PROP (let W := repeated_repeated_div (6 * i + 3) n H1 in
0 <= i /\ 6 * i + 5 <= brute_force W)
LOCAL (temp _n (Vlong (Int64.repr (repeated_repeated_div (6 * i + 3) n H1)));
temp _i (Vlong (Int64.repr (6 * i + 5))); gvars gv)
SEP (data_at Ews tulong (Vlong (Int64.repr (value_of_highest (6 * i + 3) n H1))) (gv _highest))
).
``````

But I have incorrect `forward_if` a little bit above/before it, I'm fairly sure. Now working on it (thinking and exploring behaviour of different values of `n`).

#### Lessness (Aug 18 2023 at 00:18):

Nah, my current loop invariant isn't correct.

Anyway, about that `forward_if`. Current version looks correct to me, even if a bit complex:

``````    set (repeated_repeated_div 3 n H1) as W.
set (W = 1 \/ (Z.divide ((brute_force W) ^ 2) W) \/ (prime (brute_force W - 2) /\ Z.divide (brute_force W - 2) W)) as P.
assert ({P} + {~ P}).
{ unfold P. destruct (Z.eq_dec W 1).
+ left. auto.
+ destruct (Zdivide_dec (brute_force W ^ 2) W).
- left. auto.
- destruct (prime_dec (brute_force W - 2)).
* destruct (Zdivide_dec (brute_force W - 2) W).
++ left. auto.
++ right. tauto.
* tauto. }
- forward_if (
if H2
then PROP ()
LOCAL (temp _n (Vlong (Int64.repr 1)); gvars gv)
SEP (data_at Ews tulong (Vlong (Int64.repr (brute_force W))) (gv _highest))
else PROP ()
LOCAL (temp _n (Vlong (Int64.repr (brute_force W))); gvars gv)
SEP (data_at Ews tulong (Vlong (Int64.repr (value_of_second_highest W (repeated_repeated_div_thm0 _ _ H1)))) (gv _highest))
).
``````

I tried to prove that it is correct, but it is hard (impossible?) to do without correct loop invariant. For now it's just an empirical observation. And getting to the correct loop invariant is hard because C has size limitations on the types and I have to prove that, dunno, `6 * i + 3 <= Int64.unsigned` very often. That means, loop invariant has to take that into account and that makes harder to find the correct loop invariant than for "pseudocode language" that have unlimited size integers.

#### Lessness (Aug 20 2023 at 15:24):

I'm trying to make an inductive type that kinda simulates the behavior of the `find` function (from the https://github.com/LessnessRandomness/PE/blob/main/EulerProject3.c), but without size limitations on the integers.

``````Inductive state N (H: 1 <= N): Z -> Z -> nat -> Prop :=
| Start: state N H (value_of_highest 3 N) (snd (repeated_div 3 (snd (repeated_div 2 N)))) 0
| Loop: forall h n i, (6 * Z.of_nat i + 5) * (6 * Z.of_nat i + 5) <= n ->
state N H h n i ->
state N H (value_of_highest (6 * Z.of_nat i + 7) N) (snd (repeated_div (6 * Z.of_nat i + 7) (snd (repeated_div (6 * Z.of_nat i + 5) n)))) (i + 1).
``````

Does this look good/correct enough?

`value_of_highest x n` is the highest prime divisor of `n` that's less or equal to `x`.
`snd (repeated_div x n)` is equal to `n` with `x` divided out of it.

The goal is to define the analogue of loop invariant for it and prove its correctness.

``````Definition is_loop_invariant (T: Z -> Z -> Prop) :=
forall N H n i h (s: state N H h n i), T n (Z.of_nat i).
``````

Last updated: Jun 23 2024 at 05:02 UTC