wneuper/isa: comparison src/Tools/isac/Knowledge/GCD_Poly

equal deleted inserted replaced

-:84da1e3e4600
+:679c95f19808
 value "gcd (15::int) (6::int) = 3"
 value "gcd (10::int) (3::int) = 1"
 value "gcd (6::int) (24::int) = 6"
-(* why has this been removed since 2002 ? *)
-definition last_elem :: "'a list \<Rightarrow> 'a" where
-"last_elem ls = (hd o rev) ls"
 (* drop tail elements equal 0 *)
 primrec drop_hd_zeros :: "int list \<Rightarrow> int list" where
 "drop_hd_zeros (p # ps) = (if p = 0 then drop_hd_zeros ps else (p # ps))"
 definition drop_tl_zeros :: "int list \<Rightarrow> int list" where
 "drop_tl_zeros p = (rev o drop_hd_zeros o rev) p"
 subsection {* modulo calculations for integers *}
 (* modi is just local for mod_inv *)
 function modi :: "int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat" where
 "modi r rold m rinv =
 (if m \<le> rinv
 then 0 else
 if r mod (int m) = 1
 then rinv
 else modi (rold * ((int rinv) + 1)) rold m (rinv + 1))"
 by auto
 termination by (relation "measure (\<lambda>(r, rold, m, rinv). m - rinv)") auto
 (* mod_inv :: int \<Rightarrow> nat \<Rightarrow> nat
 mod_inv r m = s
 assumes True
 yields r * s mod m = 1 *)
-definition mod_inv :: "int \<Rightarrow> nat \<Rightarrow> nat" where
+definition mod_inv :: "int \<Rightarrow> nat \<Rightarrow> nat" where "mod_inv r m = modi r r m 1"
-"mod_inv r m = modi r r m 1"
 value "modi 5 5 7 1   = 3"
 value "modi 3 3 7 1   = 5"
 value "modi 4 4 339 1 = 85"
 (* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat
 mod_div a b m = c
 assumes True
 yields a * b ^(-1) mod m = c   <\<Longrightarrow>  a mod m = *)
 definition mod_div :: "int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat" where
 "mod_div a b m = (nat a) * (mod_inv b m) mod m"
 (* root1 is just local to approx_root *)
 function root1 :: "int \<Rightarrow> nat \<Rightarrow> nat" where
 "root1 a n = (if (int (n * n)) < a then root1 a (n + 1) else n)"
 by auto
 termination by (relation "measure (\<lambda>(a, n). nat (a - (int (n * n))))") auto
 (* required just for one approximation:
 approx_root :: nat \<Rightarrow> nat
 approx_root a = r
 assumes True
 yields r * r \<ge> a *)
-definition approx_root :: "int \<Rightarrow> nat" where
+definition approx_root :: "int \<Rightarrow> nat" where "approx_root a = root1 a 1"
-"approx_root a = root1 a 1"
 (* chinese_remainder :: int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int
 chinese_remainder (r1, m1) (r2, m2) = r
 assumes True
 yields r = r1 mod m1 \<and> r = r2 mod m2 *)
 definition chinese_remainder :: "int \<Rightarrow> nat \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat" where
 "chinese_remainder r1 m1 r2 m2 =
 ((nat (r1 mod (int m1))) +
 (nat (((r2 - (r1 mod (int m1))) * (int (mod_inv (int m1) m2))) mod (int m2))) * m1)"
 value "chinese_remainder 17 9 3 4  = 35"
 value "chinese_remainder  7 2 6 11 = 17"
 subsection {* creation of lists of primes for efficiency *}
 assumes max ps < n  \<and>  n \<le> (max ps)^2  \<and>
 (*     FIXME: really ^^^^^^^^^^^^^^^? *)
 (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and> (\<forall>p. p < n \<and> Primes.prime p \<longrightarrow> List.member ps p)
 yields Primes.prime n *)
 fun is_prime :: "nat list \<Rightarrow> nat \<Rightarrow> bool" where
 "is_prime ps n =
 (if List.length ps > 0
 then
 if (n mod (List.hd ps)) = 0
 then False
 else is_prime (List.tl ps) n
 else True)"
 declare is_prime.simps [simp del] -- "make_primes, next_prime_not_dvd"
 value "is_prime [2, 3] 2  = False"    --"... precondition!"
 value "is_prime [2, 3] 3  = False"    --"... precondition!"
 value "is_prime [2, 3] 4  = False"
 assumes last_p = maxs ps  \<and>  (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p)  \<and>
 (\<forall>p. p < last_p \<and> Primes.prime p \<longrightarrow> List.member ps p)
 yields n \<le> maxs pps  \<and>  (\<forall>p. List.member pps p \<longrightarrow> Primes.prime p)  \<and>
 (\<forall>p. (p < maxs pps \<and> Primes.prime p) \<longrightarrow> List.member pps p)*)
 function make_primes :: "nat list \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat list" where
 "make_primes ps last_p n =
-(if n <= last_elem ps then ps else
+(if n <= last ps then ps else
 if is_prime ps (last_p + 2)
 then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n
 else make_primes  ps                   (last_p + 2) n)"
 by pat_completeness auto --"simp del: is_prime.simps <-- declare"
 (*Calls:
 a) (ps, last_p, n) ~> (ps @ [last_p + 2], last_p + 2, n)
 b) (ps, last_p, n) ~> (ps,                last_p + 2, n)
 Measures:
 assumes True
 yields (\<forall>p. List.member ps p \<longrightarrow> Primes.prime p) \<and>
 n \<le> maxs ps \<and> maxs ps \<le> Fact.fact n + 1 \<and>
 (\<forall>p. p \<le> maxs ps \<and> Primes.prime p \<longrightarrow> List.member ps p) *)
 definition primes_upto :: "nat \<Rightarrow> nat list" where
 "primes_upto n = (if n < 3 then [2] else make_primes [2, 3] 3 n)"
 value "primes_upto  2 = [2]"
 value "primes_upto  3 = [2, 3]"
 value "primes_upto  4 = [2, 3, 5]"
 value "primes_upto  5 = [2, 3, 5]"
 value "primes_upto  9 = [2, 3, 5, 7, 11]"
 value "primes_upto 10 = [2, 3, 5, 7, 11]"
 value "primes_upto 11 = [2, 3, 5, 7, 11]"
 (* max's' is analogous to Integer.gcds *)
-definition maxs :: "nat list \<Rightarrow> nat" where
+definition maxs :: "nat list \<Rightarrow> nat" where "maxs ns = List.fold max ns (List.hd ns)"
-"maxs ns = List.fold max ns (List.hd ns)"
 value "maxs [5, 3, 7, 1, 2, 4] = 7"
 (* find the next prime greater p not dividing the number n:
 next_prime_not_dvd :: int \<Rightarrow> int \<Rightarrow> int (infix)
 n1 next_prime_not_dvd n2 = p
 assumes True
 yields p is_prime  \<and>  n1 < p  \<and>  \<not> p dvd n2  \<and> (* smallest with these properties... *)
 (\<forall> p'. (p' is_prime  \<and>  n1 < p'  \<and>  \<not> p' dvd n2)  \<longrightarrow>  p \<le> p') *)
 function next_prime_not_dvd :: "nat \<Rightarrow> nat \<Rightarrow> nat" (infixl "next'_prime'_not'_dvd" 70) where
 "n1 next_prime_not_dvd n2 =
 (let
 ps = primes_upto (n1 + 1);
 nxt = last ps
 in
 if n2 mod nxt \<noteq> 0
 then nxt
 else nxt next_prime_not_dvd n2)"
 by auto --"simp del: is_prime.simps, make_primes.simps, primes_upto.simps <-- declare"
 termination (*next_prime_not_dvd ..Inner syntax error*) (*not finished*)
 apply (relation "measure (\<lambda>(n1, n2). n2 - n1)")
 apply auto
 sorry (*by lexicographic_order unfinished subgoals*)
 subsection {* basic notions for univariate polynomials *}
 (* not in List.thy, copy from library.ML *)
 fun nth_drop :: "nat \<Rightarrow> 'a list \<Rightarrow> 'a list" where
 "nth_drop n xs = List.take n xs @ List.drop (n + 1) xs"
 value "nth_drop 0 []              = []"
 value "nth_drop 0 [1, 2, 3::int]  = [2, 3]"
 value "nth_drop 2 [1, 2, 3::int]  = [1, 2]"
 value "nth_drop 77 [1, 2, 3::int] = [1, 2, 3]"
 (* leading coefficient *)
-definition lcoeff_up :: "unipoly \<Rightarrow> int" where
+definition lcoeff_up :: "unipoly \<Rightarrow> int" where "lcoeff_up p = (last o drop_tl_zeros) p"
-"lcoeff_up p = (last_elem o drop_tl_zeros) p"
 value "lcoeff_up [3, 4, 5, 6]    = 6"
 value "lcoeff_up [3, 4, 5, 6, 0] = 6"
 (* drop leading coefficients equal 0 *)
 "drop_lc0_up p = drop_tl_zeros p"
 (**)fun drop_lc0_up :: "unipoly \<Rightarrow> unipoly" where
 "drop_lc0_up p = drop_tl_zeros p"
 *)
 fun drop_lc0_up :: "unipoly \<Rightarrow> unipoly" where
-"drop_lc0_up [] = []"
+"drop_lc0_up [] = []" |
-| "drop_lc0_up p =
+"drop_lc0_up p =
 (let l = List.length p - 1
 in
 if p ! l = 0
 then drop_lc0_up (nth_drop l p)
 else p)"
 value "drop_lc0_up [0, 1, 2, 3, 4, 5, 0, 0] = [0, 1, 2, 3, 4, 5]"
 value "drop_lc0_up [0, 1, 2, 3, 4, 5]       = [0, 1, 2, 3, 4, 5]"
 (* degree *)
 (**)definition deg_up :: "unipoly \<Rightarrow> nat" where
 "deg_up p = ((op - 1) o length o drop_lc0_up) p"
 *)
 function deg_up :: "unipoly \<Rightarrow> nat" where
 "deg_up p =
 (let len = List.length p - 1
 in
 if p ! len \<noteq> 0 then len else deg_up (nth_drop len p))"
 by (*pat_completeness*) auto
 lemma termin_1:
 fixes p::"nat list"
 assumes "length p - Suc 0 = 0"
 value "deg_up [3, 4, 5, 6, 0] = 3"
 value "deg_up [1, 0, 3, 0, 0] = 2"
 (* norm is just local to normalise *)
 fun norm :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> unipoly" where
 "norm p nrm m lcp i =
 (if i = 0
 then [int (mod_div (p ! i) lcp m)] @ nrm
 else norm p ([int (mod_div (p ! i) lcp m)] @ nrm) m lcp (i - 1))"
 (* normalise a unipoly such that lcoeff_up mod m = 1.
 normalise :: unipoly \<Rightarrow> nat \<Rightarrow> unipoly
 normalise [p_0, .., p_k] m = [q_0, .., q_k]
 assumes
 yields \<exists> t. 0 \<le> i \<le> k \<Rightarrow> (p_i * t) mod m = q_i  \<and>  (p_k * t) mod m = 1 *)
 fun normalise :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where
 "normalise p m =
 (let
 p = drop_lc0_up p;
 lcp = lcoeff_up p
 in
 if List.length p = 0 then [] else norm p [] m lcp (List.length p - 1))"
 declare normalise.simps [simp del] --"HENSEL_lifting_up"
 value "normalise [-18, -15, -20, 12, 20, -13, 2] 5 = [1, 0, 0, 1, 0, 1, 1]"
 value "normalise [9, 6, 3] 10                      = [3, 2, 1]"
 subsection {* addition, multiplication, division *}
-(* scalar multiplication, %* in GCD_Poly.thy *)
+(* scalar multiplication *)
-definition  mult_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "mult'_ups" 70) where
+definition  mult_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "%*" 70) where
-"p mult_ups m = List.map (op * m) p"
+"p %* m = List.map (op * m) p"
-value "[5, 4, 7, 8, 1] mult_ups 5   = [25, 20, 35, 40, 5]"
+value "[5, 4, 7, 8, 1] %* 5   = [25, 20, 35, 40, 5]"
-value "[5, 4, -7, 8, -1] mult_ups 5 = [25, 20, -35, 40, -5]"
+value "[5, 4, -7, 8, -1] %* 5 = [25, 20, -35, 40, -5]"
-(* scalar divison, %/ in GCD_Poly.thy *)
+(* scalar divison *)
-definition swapf :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'c" where
+definition swapf :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'c" where "swapf f a b = f b a"
-"swapf f a b = f b a"
+definition div_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "div'_ups" (* %/ error FLORIAN*) 70) where
-definition div_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "div'_ups" 70) where
+"p div_ups m = map (swapf op div2 m) p"
-"p div_ups m = map (swapf op div2 m) p"
 value "[4, 3, 2, 5, 6] div_ups 3 = [1, 1, 0, 1, 2]"
 value "[4, 3, 2, 0] div_ups 3    = [1, 1, 0, 0]"
 (* not in List.thy, copy from library.ML *)
 fun map2 :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'a list \<Rightarrow> 'b list \<Rightarrow> 'c list" where
-"map2 _ [] [] = []"
+"map2 _ [] [] = []" |
-| "map2 f (x # xs) (y # ys) = f x y # map2 f xs ys"
+"map2 f (x # xs) (y # ys) = f x y # map2 f xs ys" |
-| "map2 _ _ _ = []" (*raise ListPair.UnequalLengths*)
+"map2 _ _ _ = []" (*raise ListPair.UnequalLengths*)
-(* minus of polys,  %-% in GCD_Poly.thy *)
+(* minus of polys *)
-definition minus_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" (infixl "minus'_up" 70) where
+definition minus_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" (infixl "%-%" 70) where
-"p1 minus_up p2 = map2 (op -) p1 p2"
+"p1 %-% p2 = map2 (op -) p1 p2"
-value "[8, -7, 0, 1] minus_up [-2, 2, 3, 0]     = [10, -9, -3, 1]"
+value "[8, -7, 0, 1] %-% [-2, 2, 3, 0]     = [10, -9, -3, 1]"
-value "[8, 7, 6, 5, 4] minus_up [2, 2, 3, 1, 1] = [6, 5, 3, 4, 3]"
+value "[8, 7, 6, 5, 4] %-% [2, 2, 3, 1, 1] = [6, 5, 3, 4, 3]"
 (* lemmata for pattern compatibility in dvd_up *)
-lemma ex_fals0_1: "([d], [p]) = (v # vb # vc, ps) \<Longrightarrow> QUODLIBET"
+lemma ex_fals0_1: "([d], [p]) = (v # vb # vc, ps) \<Longrightarrow> QUODLIBET" by simp
-by simp
+lemma ex_fals0_2: "([], ps) = (v # vb # vc, psa) \<Longrightarrow> QUODLIBET" by simp
-lemma ex_fals0_2: "([], ps) = (v # vb # vc, psa) \<Longrightarrow> QUODLIBET"
+lemma ex_fals0_3: "(ds, []) = (dsa, v # vb # vc) \<Longrightarrow> QUODLIBET" by simp
-by simp
-lemma ex_fals0_3: "(ds, []) = (dsa, v # vb # vc) \<Longrightarrow> QUODLIBET"
+function (sequential) dvd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> bool" (infixl "%|%" 70) where
-by simp
+"[d] %|% [p] = ((\<bar>d\<bar> \<le> \<bar>p\<bar>) & (p mod d = 0))" |
+"ds %|% ps =
-function (sequential) dvd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> bool" (infixl "dvd'_up" 70) where
+(let
-"[d] dvd_up [p] = ((\<bar>d\<bar> \<le> \<bar>p\<bar>) & (p mod d = 0))"
+ds = drop_lc0_up ds; ps = drop_lc0_up ps;
-| "ds dvd_up ps =
+d000 = (List.replicate (List.length ps - List.length ds) 0) @ ds;
-(let
+quot = (lcoeff_up ps) div2 (lcoeff_up d000);
-ds = drop_lc0_up ds; ps = drop_lc0_up ps;
+rest = drop_lc0_up (ps %-% (d000 %* quot))
-d000 = (List.replicate (List.length ps - List.length ds) 0) @ ds;
+in
-quot = (lcoeff_up ps) div2 (lcoeff_up d000);
+if rest = [] then True else
-rest = drop_lc0_up (ps minus_up (d000 mult_ups quot))
+if quot \<noteq> 0 & List.length ds \<le> List.length rest then ds %|% rest else False)"
-in
-if rest = [] then True else
-if quot \<noteq> 0 & List.length ds \<le> List.length rest then ds dvd_up rest else False)"
 apply pat_completeness
 apply simp
 apply simp
 defer (*a "mixed" obligation*)
 apply simp
 done (* > 1 sec IMPROVED BY declare simp del:
 centr_up_def normalise.simps mod_up_gcd.simps lcoeff_up.simps*)
 termination (*dvd_up ..Inner syntax error*) sorry (*by lexicographic_order LOOPS*) (*not finished*)
 (*cutting definition to 'fun' loops, so no Calls, Measure, Matrix ---?*)
-value "[4] dvd_up [6]                 = False"
+value "[4] %|% [6]                 = False"
-value "[8] dvd_up [16, 0]             = True"
+value "[8] %|% [16, 0]             = True"
-value "[3, 2] dvd_up [0, 0, 9, 12, 4] = True"
+value "[3, 2] %|% [0, 0, 9, 12, 4] = True"
-value "[8, 0] dvd_up [16]             = True"
+value "[8, 0] %|% [16]             = True"
 subsection {* normalisation and Landau-Mignotte bound *}
 (* centr is just local to centr_up *)
 definition centr :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int" where
 "centr m mid p_i = (if mid < p_i then p_i - (int m) else p_i)"
 (* normalise :: centr_up \<Rightarrow> unipoly => int \<Rightarrow> unipoly
 normalise [p_0, .., p_k] m = [q_0, .., q_k]
 assumes
 yields 0 \<le> i \<le> k \<Rightarrow> |^ ~m/2 ^| <= q_i <=|^ m/2 ^|
 (where |^ x ^| means round x up to the next greater integer) *)
 definition centr_up :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" where
 "centr_up p m =
 (let
 mi = (int m) div2 2;
 mid = if (int m) mod mi = 0 then mi else mi + 1
 in map (centr m mid) p)"
 value "centr_up [7, 3, 5, 8, 1, 3] 10 = [-3, 3, 5, -2, 1, 3]"
 value "centr_up [1, 2, 3, 4, 5] 2     = [1, 0, 1, 2, 3]"
 value "centr_up [1, 2, 3, 4, 5] 3     = [1, -1, 0, 1, 2]"
 value "centr_up [1, 2, 3, 4, 5] 4     = [1, 2, -1, 0, 1]"
 value "centr_up [1, 2, 3, 4, 5] 8     = [1, 2, 3, 4, -3]"
 value "centr_up [1, 2, 3, 4, 5] 9     = [1, 2, 3, 4, 5]"
 value "centr_up [1, 2, 3, 4, 5] 10    = [1, 2, 3, 4, 5]"
 value "centr_up [1, 2, 3, 4, 5] 11    = [1, 2, 3, 4, 5]"
-(*  sum_lmb :: centr_up \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int
+(* sum_lmb :: centr_up \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int
 sum_lmb [p_0, .., p_k] e = s
 assumes exp >= 0
 yields. p_0^e + p_1^e + ... + p_k^e *)
 definition sum_lmb :: "unipoly \<Rightarrow> nat \<Rightarrow> int" where
 "sum_lmb p exp = List.fold ((op +) o (swapf power exp)) p 0"
 value "sum_lmb [-1, 2, -3, 4, -5] 1 = -3"
 value "sum_lmb [-1, 2, -3, 4, -5] 2 = 55"
 value "sum_lmb [-1, 2, -3, 4, -5] 3 = -81"
 value "sum_lmb [-1, 2, -3, 4, -5] 4 = 979"
 LANDAU_MIGNOTTE_bound [a_0, ..., a_m] [b_0, .., b_n] = lmb
 assumes True
 yields lmb = 2^min(m,n) * gcd(a_m,b_n) *
 min( 1/|a_m| * root(sum_lmb [a_0,...a_m] 2 , 1/|b_n| * root(sum_lmb [b_0,...b_n] 2)*)
 definition LANDAU_MIGNOTTE_bound :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat" where
 "LANDAU_MIGNOTTE_bound p1 p2 =
 ((power 2 (min (deg_up p1) (deg_up p2))) * (nat (gcd (lcoeff_up p1) (lcoeff_up p2))) *
 (nat (min (abs ((int (approx_root (sum_lmb p1 2))) div2 -(lcoeff_up p1)))
 (abs ((int (approx_root (sum_lmb p2 2))) div2 -(lcoeff_up p2))))))"
-(*                                      |------------|        -|------------|
-|--------------------------|
-|--------------------------------|
-|-------------------------------------------------------|  *)
 value "LANDAU_MIGNOTTE_bound [1] [4, 5]          = 1"
 value "LANDAU_MIGNOTTE_bound [1, 2] [4, 5]       = 2"
 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4, 5]    = 2"
 value "LANDAU_MIGNOTTE_bound [1, 2, 3] [4]       = 1"
 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4]        = 1"
 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5]    = 2"
 value "LANDAU_MIGNOTTE_bound [-1, 2, -3] [4, -5, 6] = 12"
 subsection {* modulo calculations for polynomials *}
+value "List.compound"
 (* pair is just local to chinese_remainder_up, is "op ~~" in library.ML *)
 fun pair :: "unipoly \<Rightarrow> unipoly \<Rightarrow> ((int \<times> int) list)" (infix "pair" 4) where
-"([] pair []) = []"
+"([] pair []) = []" |
-| "([] pair ys) = []" (*raise ListPair.UnequalLengths*)
+"([] pair ys) = []" | (*raise ListPair.UnequalLengths*)
-| "(xs pair []) = []" (*raise ListPair.UnequalLengths*)
+"(xs pair []) = []" | (*raise ListPair.UnequalLengths*)
-| "((x#xs) pair (y#ys)) = (x, y) # (xs pair ys)"
+"((x#xs) pair (y#ys)) = (x, y) # (xs pair ys)"
 fun chinese_rem :: "nat \<times> nat \<Rightarrow> int \<times> int \<Rightarrow> int" where
 "chinese_rem (m1, m2) (p1, p2) = (int (chinese_remainder p1 m1 p2 m2))"
 (* chinese_remainder_up :: int * int \<Rightarrow> unipoly * unipoly \<Rightarrow> unipoly
 chinese_remainder_up (m1, m2) (p1, p2) = p
 assume m1, m2 relatively prime
 yields p1 = p mod m1 \<and> p2 = p mod m2 *)
 fun chinese_remainder_up :: "nat \<times> nat \<Rightarrow> unipoly \<times> unipoly \<Rightarrow> unipoly" where
 "chinese_remainder_up (m1, m2) (p1, p2) = map (chinese_rem (m1, m2)) (p1 pair p2)"
 value "chinese_remainder_up (5, 7) ([2, 2, 4, 3], [3, 2, 3, 5]) = [17, 2, 24, 33]"
 (* mod_up :: unipoly \<Rightarrow> int \<Rightarrow> unipoly
 mod_up [p1, p2, ..., pk] m = up
 assume m > 0
 yields up = [p1 mod m, p2 mod m, ..., pk mod m]*)
 definition mod' :: "nat \<Rightarrow> int \<Rightarrow> int" where "mod' m i = i mod (int m)"
 definition mod_up :: "unipoly \<Rightarrow> nat \<Rightarrow> unipoly" (infixl "mod'_up" 70) where
 "p mod_up m = map (mod' m) p"
 value "[5, 4, 7, 8, 1] mod_up 5 = [0, 4, 2, 3, 1]"
 value "[5, 4,-7, 8,-1] mod_up 5 = [0, 4, 3, 3, 4]"
-value "curry"
-value "curry"
-value "5 mod 3::int"
 (* euclidean algoritm in Z_p[x].
 mod_up_gcd :: unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly
 mod_up_gcd p1 p2 m = g
 assumes
 yields  gcd (p1 mod m) (p2 mod m) = g *)
 function mod_up_gcd :: "unipoly \<Rightarrow> unipoly  \<Rightarrow> nat \<Rightarrow> unipoly" where
 "mod_up_gcd p1 p2 m =
 (let
 p1m = p1 mod_up m;
 p2m = drop_lc0_up (p2 mod_up m);
 p2n = (replicate (List.length p1 - List.length p2m) 0) @ p2m;
 quot = mod_div (lcoeff_up p1m) (lcoeff_up p2n) m;
-rest = drop_lc0_up ((p1m minus_up (p2n mult_ups (int quot))) mod_up m)
+rest = drop_lc0_up ((p1m %-% (p2n %* (int quot))) mod_up m)
 in
 if rest = []
 then p2
 else
 if List.length rest < List.length p2
 then mod_up_gcd p2 rest m
 else mod_up_gcd rest p2 m)"
 by auto
 termination mod_up_gcd (*not finished*)
 apply (relation "measures [\<lambda>(p1, p2, m). length p1,
 \<lambda>(p1, p2, m). 00000000000000]")
 apply auto
 gcds :: int list \<Rightarrow> int
 gcds ns = d
 assumes True
 yields THE d. ((\<forall>n. List.member ns n \<longrightarrow> d dvd n) \<and>
 (\<forall>d'. (\<forall>n. List.member ns n \<and> d' dvd n) \<longrightarrow> d'modp \<le> d)) *)
-fun gcds :: "int list \<Rightarrow> int" where
+fun gcds :: "int list \<Rightarrow> int" where "gcds ns = List.fold gcd ns (List.hd ns)"
-"gcds ns = List.fold gcd ns (List.hd ns)"
 value "gcds [6, 9, 12] = 3"
 value "gcds [6, -9, 12] = 3"
 value "gcds [8, 12, 16] = 4"
 value "gcds [-8, 12, -16] = 4"
 (* prim_poly :: unipoly \<Rightarrow> unipoly
 prim_poly p = pp
 assumes True
 yields \<forall>p1, p2. (List.member pp p1 \<and> List.member pp p2 \<and> p1 \<noteq> p2) \<longrightarrow> gcd p1 p2 = 1 *)
 fun primitive_up :: "unipoly \<Rightarrow> unipoly" where
-"primitive_up [c] = (if c = 0 then [0] else [1])"
+"primitive_up [c] = (if c = 0 then [0] else [1])" |
-| "primitive_up p =
+"primitive_up p =
 (let d = gcds p
 in
 if d = 1 then p else p div_ups d)"
 value "primitive_up [12, 16, 32, 44] = [3, 4, 8, 11]"
 value "primitive_up [4, 5, 12] =  [4, 5, 12]"
 value "primitive_up [0] = [0]"
 value "primitive_up [6] = [1]"
 yields  new_g = [1] \<or> (new_g \<ge> g \<and> P > M)
 argumentns "a b d M P g p" named according to [1] p.93 *)
 (*declare [[show_types]]*)
 function try_new_prime_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> unipoly"
 where
-"try_new_prime_up           a          b          d      M      P     g          p   =
+"try_new_prime_up             a          b          d      M      P     g          p   =
 (if P > M then g else
 let p' = p next_prime_not_dvd d;
-g_p = centr_up (          (          (normalise (mod_up_gcd a b p') p')
+g_p = centr_up (          (      (normalise (mod_up_gcd a b p') p')
-mult_ups ((int d) mod (int p')))
+%* ((int d) mod (int p')))
 mod_up p')
 p'
 in
 if deg_up g_p < deg_up g
 then
 if (deg_up g_p) = 0 then [1] else try_new_prime_up a b d M p g_p p'
 else
 if deg_up g_p \<noteq> deg_up g then try_new_prime_up a b d M P g p' else
 let
 P = P * p';
 g = centr_up ((chinese_remainder_up (P, p') (g, g_p)) mod_up P) P
 in
-try_new_prime_up a b d M P g p'
+try_new_prime_up a b d M P g p')"
-)"
 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps"
 termination try_new_prime_up sorry (*by lexicographic_order LOOPS*) (*not finished*)
 (*cutting definition to 'fun' loops, so no Calls, Measure, Matrix ---?*)
 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly
 p1 is primitive  \<and>  p2 is primitive
 yields  gcd | p1  \<and>  gcd | p2  \<and>  gcd is primitive
 argumentns "a b d M p" named according to [1] p.93 *)
 function HENSEL_lifting_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> unipoly" where
 "HENSEL_lifting_up a b d M p =
 (let
 p = p next_prime_not_dvd d;
-g_p = centr_up (((normalise (mod_up_gcd a b p) p) mult_ups (int (d mod p))) mod_up p) p
+g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (int (d mod p))) mod_up p) p
 in
 if deg_up g_p = 0 then [1] else
 (let
 g = primitive_up (try_new_prime_up a b d M p g_p p)
 in
-if (g dvd_up a) \<and> (g dvd_up b) then g else HENSEL_lifting_up a b d M p))"
+if (g %|% a) \<and> (g %|% b) then g else HENSEL_lifting_up a b d M p))"
 by pat_completeness auto --"simp del: centr_up_def normalise.simps mod_up_gcd.simps"
 termination HENSEL_lifting_up sorry (*by lexicographic_order LOOPS*) (*not finished*)
 (*cutting definition to 'fun' loops, so no Calls, Measure, Matrix ---?*)
 (* gcd_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly
 gcd_up a b = c
 assumes not (a = [] \<or> a = [0]) \<and> not (b = []\<or> b = [0]) \<and>
 a is primitive \<and> b is primitive
 yields c dvd a \<and> c dvd b \<and> (\<forall>c'. (c' dvd a \<and> c' dvd b) \<longrightarrow> c' \<le> c) *)
 function gcd_up :: "unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly" where
 "gcd_up a b =
 (let d = \<bar>gcd (lcoeff_up a) (lcoeff_up b)\<bar>
 in
 if a = b then a else
 HENSEL_lifting_up a b (nat d) (2 * (nat d) * LANDAU_MIGNOTTE_bound a b) 1)"
 by pat_completeness auto --"simp del: lcoeff_up.simps ?+ others?"
 termination by lexicographic_order (*works*)
 value "gcd_up [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2] = [-2, -1, 1]"
 (*     gcd    (-1 + x^2) (x + x^2) = (1 + x) *)
 value "gcd_up [-1, 0 ,1] [0, 1, 1] = [1, 1]"
-primrec nth :: "'a list => nat => 'a" (infixl "!!" 100) where
-nth_Cons: "(x # xs) !! n = (case n of 0 \<Rightarrow> x | Suc k \<Rightarrow> xs !! k)"
-value "[1,2,3,4,5] ! 2"
-definition  mult_ups2 :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "%*" 70) where
-"p %* m = List.map (op * m) p"
 (*
 print_configs
 declare [[simp_trace_depth_limit = 99]]
 declare [[simp_trace = true]]

changeset 48864	679c95f19808
parent 48863	84da1e3e4600
child 48865	97408b42129d