1.1 --- a/src/Tools/isac/Knowledge/GCD_Poly.thy Fri May 24 16:50:31 2013 +0200
1.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly.thy Fri May 31 17:49:34 2013 +0200
1.3 @@ -355,7 +355,18 @@
1.4
1.5 argumentns "a b d M P g p" named according to [1] p.93 *)
1.6 fun try_new_prime_up a b d M P g p =
1.7 - if P > M then g else
1.8 +(writeln ("try_new_prime_up: a = " ^ PolyML.makestring a ^
1.9 +", b = " ^ PolyML.makestring b ^
1.10 +", d = " ^ PolyML.makestring d ^
1.11 +", M = " ^ PolyML.makestring M ^
1.12 +", P = " ^ PolyML.makestring P ^
1.13 +", g = " ^ PolyML.makestring g ^
1.14 +", p = " ^ PolyML.makestring p);
1.15 + if P > M then
1.16 +(writeln ("try_new_prime_up 1 -----> " ^ PolyML.makestring g);
1.17 + g
1.18 +)
1.19 + else
1.20 let
1.21 val p = p next_prime_not_dvd d
1.22 val g_p = centr_up ( ( (normalise (mod_up_gcd a b p) p)
1.23 @@ -365,7 +376,11 @@
1.24 in
1.25 if deg_up g_p < deg_up g
1.26 then
1.27 - if (deg_up g_p) = 0 then [1]: unipoly else try_new_prime_up a b d M p g_p p
1.28 + if (deg_up g_p) = 0 then
1.29 +(writeln ("try_new_prime_up 2 -----> " ^ PolyML.makestring [1]);
1.30 + [1]: unipoly
1.31 +)
1.32 + else try_new_prime_up a b d M p g_p p
1.33 else
1.34 if deg_up g_p <> deg_up g then try_new_prime_up a b d M P g p else
1.35 let
1.36 @@ -373,7 +388,7 @@
1.37 val g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up P) P
1.38 in try_new_prime_up a b d M P g p end
1.39 end
1.40 -
1.41 +)
1.42 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly
1.43 HENSEL_lifting_up p1 p2 d M p = gcd
1.44 assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
1.45 @@ -382,19 +397,33 @@
1.46 yields gcd | p1 \<and> gcd | p2 \<and> gcd is primitive
1.47
1.48 argumentns "a b d M p" named according to [1] p.93 *)
1.49 - fun HENSEL_lifting_up a b d M p =
1.50 + fun HENSEL_lifting_up a b d M p =
1.51 +(writeln ("HENSEL_lifting_up: a = " ^ PolyML.makestring a ^
1.52 +", b = " ^ PolyML.makestring b ^
1.53 +", d = " ^ PolyML.makestring d ^
1.54 +", M = " ^ PolyML.makestring M ^
1.55 +", p = " ^ PolyML.makestring p);
1.56 let
1.57 val p = p next_prime_not_dvd d
1.58 - val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p (*see above*)
1.59 + val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p
1.60 + (*.. nesting of functions see above*)
1.61 in
1.62 - if deg_up g_p = 0 then [1] else
1.63 + if deg_up g_p = 0 then
1.64 +(writeln ("HENSEL_lifting_up 1 ----> " ^ PolyML.makestring [1]);
1.65 + [1]
1.66 +)
1.67 + else
1.68 let
1.69 val g = primitive_up (try_new_prime_up a b d M p g_p p)
1.70 in
1.71 - if (g %|% a) andalso (g %|% b) then g:unipoly else HENSEL_lifting_up a b d M p
1.72 + if (g %|% a) andalso (g %|% b) then
1.73 +(writeln ("HENSEL_lifting_up 2 -----> " ^ PolyML.makestring g);
1.74 + g:unipoly
1.75 +)
1.76 + else HENSEL_lifting_up a b d M p
1.77 end
1.78 end
1.79 -
1.80 +)
1.81 (* gcd_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> unipoly
1.82 gcd_up a b = c
1.83 assumes not (a = [] \<or> a = [0]) \<and> not (b = []\<or> b = [0]) \<and>
2.1 --- a/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Fri May 24 16:50:31 2013 +0200
2.2 +++ b/src/Tools/isac/Knowledge/GCD_Poly_FP.thy Fri May 31 17:49:34 2013 +0200
2.3 @@ -259,10 +259,9 @@
2.4 from primes_upto_1 have "n1 < last (primes_upto (Suc n1))" by auto
2.5 from this have "(int n1) - (int (last (primes_upto (Suc n1)) * q))
2.6 < (int (last (primes_upto (Suc n1)))) - (int (last (primes_upto (Suc n1)) * q))" by auto
2.7 - from this show ?thesis sorry
2.8 + from this show ?thesis sorry (*by (rule neg_less_iff_less)*)
2.9 qed (*?FLORIAN : lifting nat -> int ?*)
2.10
2.11 -
2.12 find_theorems "_ - _ < _ - _ = _" thm Nat.less_diff_iff
2.13 find_theorems "_ * _ < _ * _ = _" thm Nat.Suc_mult_less_cancel1
2.14 find_theorems "-_ < -_ = _" thm Groups.ordered_ab_group_add_class.neg_less_iff_less
2.15 @@ -439,9 +438,18 @@
2.16 value "[5, 4, -7, 8, -1] %* 5 = [25, 20, -35, 40, -5]"
2.17
2.18 (* scalar divison *)
2.19 +(*
2.20 definition swapf :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'c" where "swapf f a b = f b a"
2.21 +CODEGEN CAUSES ERROR:
2.22 +ML error (line 913 of "/home/neuper/devel/isac/codegen/gcd_univariate.sml"):
2.23 +Type error in function application. Function: div2 : inta -> inta -> inta
2.24 + Argument: (swapf, m) : (('a -> 'b -> 'c) -> 'b -> 'a -> 'c) * 'd
2.25 + Reason: Can't unify inta to (('a -> 'b -> 'c) -> 'b -> 'a -> 'c) * 'd (Incompatible types)
2.26 +THUS TYPE CONSTRAINED...
2.27 +*)
2.28 +definition swapf1 :: "(int \<Rightarrow> int \<Rightarrow> int) \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int" where "swapf1 f a b = f b a"
2.29 definition div_ups :: "unipoly \<Rightarrow> int \<Rightarrow> unipoly" (infixl "div'_ups" (* %/ error FLORIAN*) 70) where
2.30 -"p div_ups m = map (swapf op div2 m) p"
2.31 +"p div_ups m = map (swapf1 op div2 m) p"
2.32
2.33 value "[4, 3, 2, 5, 6] div_ups 3 = [1, 1, 0, 1, 2]"
2.34 value "[4, 3, 2, 0] div_ups 3 = [1, 1, 0, 0]"
2.35 @@ -561,12 +569,25 @@
2.36 value "centr_up [1, 2, 3, 4, 5] 10 = [1, 2, 3, 4, 5]"
2.37 value "centr_up [1, 2, 3, 4, 5] 11 = [1, 2, 3, 4, 5]"
2.38
2.39 +(*
2.40 +definition swapf :: "('a \<Rightarrow> 'b \<Rightarrow> 'c) \<Rightarrow> 'b \<Rightarrow> 'a \<Rightarrow> 'c" where "swapf f a b = f b a"
2.41 +CODEGEN CAUSES ERROR:
2.42 +ML error (line 913 of "/home/neuper/devel/isac/codegen/gcd_univariate.sml"):
2.43 +Type error in function application. Function: div2 : inta -> inta -> inta
2.44 + Argument: (swapf, m) : (('a -> 'b -> 'c) -> 'b -> 'a -> 'c) * 'd
2.45 + Reason: Can't unify inta to (('a -> 'b -> 'c) -> 'b -> 'a -> 'c) * 'd (Incompatible types)
2.46 +THUS TYPE CONSTRAINED...
2.47 +*)
2.48 +definition swapf2 :: "(int \<Rightarrow> nat \<Rightarrow> int) \<Rightarrow> nat \<Rightarrow> int \<Rightarrow> int" where "swapf2 f a b = f b a"
2.49 +
2.50 (* sum_lmb :: centr_up \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int
2.51 sum_lmb [p_0, .., p_k] e = s
2.52 assumes exp >= 0
2.53 yields. p_0^e + p_1^e + ... + p_k^e *)
2.54 definition sum_lmb :: "unipoly \<Rightarrow> nat \<Rightarrow> int" where
2.55 -"sum_lmb p exp = List.fold ((op +) o (swapf power exp)) p 0"
2.56 +"sum_lmb p exp = List.fold ((op +) o (swapf2 power exp)) p 0"
2.57 +
2.58 +value "power"
2.59
2.60 value "sum_lmb [-1, 2, -3, 4, -5] 1 = -3"
2.61 value "sum_lmb [-1, 2, -3, 4, -5] 2 = 55"
2.62 @@ -721,10 +742,10 @@
2.63
2.64 subsection {* gcd_up, code from [1] p.93 *}
2.65 (* try_new_prime_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> unipoly
2.66 - try_new_prime_up p1 p2 d M P g p = new_g
2.67 - assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
2.68 + try_new_prime_up a b d M P g p = new_g
2.69 + assumes d = gcd (lcoeff_up a, lcoeff_up b) \<and>
2.70 M = LANDAU_MIGNOTTE_bound \<and> p = prime \<and> p ~| d \<and> P \<ge> p \<and>
2.71 - p1 is primitive \<and> p2 is primitive
2.72 + a is primitive \<and> b is primitive
2.73 yields new_g = [1] \<or> (new_g \<ge> g \<and> P > M)
2.74
2.75 argumentns "a b d M P g p" named according to [1] p.93: "p" is "prime", not "poly" ! *)
2.76 @@ -760,6 +781,16 @@
2.77 *)
2.78 sorry
2.79
2.80 +value "try_new_prime_up
2.81 + [-18, -15, -20, 12, 20, -13, 2]
2.82 + [8, 28, 22, -11, -14, 1, 2]
2.83 + (nat \<bar>gcd (lcoeff_up [-18, -15, -20, 12, 20, -13, 2]) (lcoeff_up [8, 28, 22, -11, -14, 1, 2])\<bar>)
2.84 + ((2 *
2.85 + (nat \<bar>gcd (lcoeff_up [-18, -15, -20, 12, 20, -13, 2]) (lcoeff_up [8, 28, 22, -11, -14, 1, 2])\<bar>)
2.86 + * LANDAU_MIGNOTTE_bound [-18, -15, -20, 12, 20, -13, 2] [8, 28, 22, -11, -14, 1, 2]))
2.87 + (*p next_prime_not_dvd d ...thus this function cannot be tested immediately! FLORIAN!*)
2.88 +"
2.89 +
2.90 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly
2.91 HENSEL_lifting_up p1 p2 d M p = gcd
2.92 assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and>
3.1 --- a/test/Tools/isac/Knowledge/gcd_poly.sml Fri May 24 16:50:31 2013 +0200
3.2 +++ b/test/Tools/isac/Knowledge/gcd_poly.sml Fri May 31 17:49:34 2013 +0200
3.3 @@ -377,12 +377,40 @@
3.4 mod_up p')
3.5 p'
3.6 ;
3.7 +val (a, b, d, M, P, g, p) = (0, 0, 0, 0, 0, 0, 0) (*isolate test below*)
3.8 +
3.9 +val (a, b) = ([~18, ~15, ~20, 12, 20, ~13, 2], [8, 28, 22, ~11, ~14, 1, 2])
3.10 +val d = abs (Integer.gcd (lcoeff_up a) (lcoeff_up b))
3.11 +val M = 2 * d * LANDAU_MIGNOTTE_bound a b
3.12 +val p = 1 (*-> p in head of HENSEL_lifting_up*) next_prime_not_dvd d
3.13 +val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p (*see above*)
3.14 +val (a, b, d, M, P, g, p) = (a, b, d, M, p, g_p, p);
3.15 +
3.16 +if try_new_prime_up a b d M P g p = [~1, 1, ~1] then () else error "try_new_prime_up changed";
3.17
3.18 "----------- fun gcd_up ---------------------------------";
3.19 "----------- fun gcd_up ---------------------------------";
3.20 "----------- fun gcd_up ---------------------------------";
3.21 -if gcd_up [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] = [~2, ~1, 1] then ()
3.22 -else error "gcd_up [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] = [~2, ~1, 1] changed";
3.23 +(* output from changeset
3.24 +try_new_prime_up: a = [~18, ~15, ~20, 12, 20, ~13, 2], b = [8, 28, 22, ~11, ~14, 1, 2],
3.25 +d = 2, M = 10240, P = 12597, g = [~1, 1, ~1], p = 19
3.26 +try_new_prime_up 1 -----> [~1, 1, ~1] *)
3.27 +val (a, b, d, M, P, g, p) =
3.28 + ([~18, ~15, ~20, 12, 20, ~13, 2], [8, 28, 22, ~11, ~14, 1, 2],
3.29 + 2, 10240, 3, [~1, 1, ~1], 3);
3.30 +"-----";
3.31 +if try_new_prime_up a b d M P g p = [~1, 1, ~1]
3.32 + then () else error "try_new_prime_up changed";
3.33 +
3.34 +(* output from changeset
3.35 +try_new_prime_up: a = [~18, ~15, ~20, 12, 20, ~13, 2], b = [8, 28, 22, ~11, ~14, 1, 2],
3.36 +d = 2, M = 10240, P = 96577, g = [~4, ~2, 2], p = 23
3.37 +try_new_prime_up 1 -----> [~4, ~2, 2] *)
3.38 +val (a, b, d, M, P, g, p) =
3.39 + ([~18, ~15, ~20, 12, 20, ~13, 2], [8, 28, 22, ~11, ~14, 1, 2],
3.40 + 2, 10240, 5, [2, 2, 2, 2], 5);
3.41 +if try_new_prime_up a b d M P g p = [~4, ~2, 2]
3.42 + then () else error "try_new_prime_up changed";
3.43
3.44 "=========== Multivariate Case ==========================";
3.45 "=========== Multivariate Case ==========================";
4.1 --- a/test/Tools/isac/Test_Some2.thy Fri May 24 16:50:31 2013 +0200
4.2 +++ b/test/Tools/isac/Test_Some2.thy Fri May 31 17:49:34 2013 +0200
4.3 @@ -1,16 +1,42 @@
4.4
4.5 theory Test_Some2
4.6 imports "/usr/local/isabisac/src/Tools/isac/Knowledge/GCD_Poly"
4.7 + (*required for export_code ...*)
4.8 + "/usr/local/isabisac/src/Tools/isac/Knowledge/GCD_Poly_FP"
4.9 begin
4.10 -
4.11 ML_file "/usr/local/isabisac/test/Tools/isac/Knowledge/gcd_poly.sml"
4.12
4.13 -ML {*
4.14 -*}
4.15 -ML {*
4.16 -*}
4.17 -ML {*
4.18 -*}
4.19 +(*test codegen*)
4.20 +(*
4.21 +export_code "try_new_prime_up" in SML
4.22 +module_name GCD_Univariate file "/home/neuper/devel/isac/codegen/gcd_un_Test_Some2.ML"
4.23 +*)
4.24 +ML_file "/home/neuper/devel/isac/codegen/try_new_prime_up.ML" (*automat. generated *)
4.25 +ML_file "/home/neuper/devel/isac/codegen/try_new_prime_up.sml" (*included test:
4.26 +ML error (line 769 of "/home/neuper/devel/isac/codegen/try_new_prime_up.sml"):
4.27 +Type error in function application.
4.28 + Function: GCD_Univariate.try_new_prime_up :
4.29 + GCD_Univariate.inta list ->
4.30 + GCD_Univariate.inta list -> GCD_Univariate.nat -> GCD_Univariate.nat -> GCD_Univariate.nat ->
4.31 + GCD_Univariate.inta list -> GCD_Univariate.nat -> GCD_Univariate.inta list
4.32 + Argument: a : int list Reason: Can't unify GCD_Univariate.inta with int (*In Basis*)
4.33 + (Different type constructors)
4.34 +*)
4.35 +
4.36 +(*ML_file "/home/neuper/devel/isac/codegen/gcd_univariate-error1-swapf.sml" errors*)
4.37 +(*ML_file "/home/neuper/devel/isac/codegen/gcd_univariate-error2-div2.sml" errors*)
4.38 +(*
4.39 +export_code "gcd_up" in SML
4.40 +module_name GCD_Univariate file "/home/neuper/devel/isac/codegen/gcd_un_Test_Some2.ML"
4.41 +*)
4.42 +export_code "gcd_up" in SML (*WITHOUT IMPORT: Not a constant: gcd_up*)
4.43 +module_name GCD_Univariate file "/home/neuper/devel/isac/codegen/gcd_un_Test_Some2.ML"
4.44 +
4.45 +ML_file "/home/neuper/devel/isac/codegen/gcd_un_Test_Some2.ML" (*automat. generated *)
4.46 +(*ML error (line 915 of "/home/neuper/devel/isac/codegen/gcd_un_Test_Some2.ML"):
4.47 +Type error in function application. Function: div2 : inta -> inta -> inta
4.48 + Argument: (swapf1, m) : (('a -> 'b -> 'c) -> 'b -> 'a -> 'c) * 'd
4.49 + Reason: Can't unify inta to (('a -> 'b -> 'c) -> 'b -> 'a -> 'c) * 'd (Incompatible types)*)
4.50 ML {*
4.51 *}
4.52 ML {*
4.53 @@ -20,13 +46,35 @@
4.54 ML {*
4.55 *}
4.56 ML {*
4.57 +if gcd_up [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] = [~2, ~1, 1] then ()
4.58 +else error "gcd_up [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] = [~2, ~1, 1] changed";
4.59 +*}
4.60 +ML {*
4.61 +
4.62 +(* gcd (1 + x^2) (x + x^2) = (1 + x) *)
4.63 +if gcd_up [~1, 0 ,1] [0, 1, 1] = [1, 1]
4.64 + then () else error "gcd_up [~1, 0 ,1] ... changed";
4.65 *}
4.66 ML {*
4.67 +val (a, b, d, M, P, g, p) =
4.68 + ([~18, ~15, ~20, 12, 20, ~13, 2], [8, 28, 22, ~11, ~14, 1, 2],
4.69 + 2, 10240, 3, [~1, 1, ~1], 3);
4.70 +if try_new_prime_up a b d M P g p = [~1, 1, ~1]
4.71 + then () else error "try_new_prime_up changed";
4.72 +*}
4.73 +ML {*
4.74 +val (a, b, d, M, P, g, p) =
4.75 + ([~18, ~15, ~20, 12, 20, ~13, 2], [8, 28, 22, ~11, ~14, 1, 2],
4.76 + 2, 10240, 5, [2, 2, 2, 2], 5);
4.77 +if try_new_prime_up a b d M P g p = [~4, ~2, 2]
4.78 + then () else error "try_new_prime_up changed";
4.79 +*}
4.80 +ML {*
4.81 *}
4.82 ML {*
4.83 *}
4.84 ML {*
4.85 -*}
4.86 +*}
4.87 (*------------------------------------------------------------------------------------------
4.88 "~~~~~ fun , args:"; val () = ();
4.89 "~~~~~ to return val:"; val () = ();