112 if r mod m = 1 |
112 if r mod m = 1 |
113 then rinv |
113 then rinv |
114 else modi (rold * (rinv + 1), rold, m, rinv + 1) |
114 else modi (rold * (rinv + 1), rold, m, rinv + 1) |
115 in modi (r, r, m, 1) end |
115 in modi (r, r, m, 1) end |
116 |
116 |
117 (* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> nat |
117 (* mod_div :: int \<Rightarrow> int \<Rightarrow> nat \<Rightarrow> natO |
118 mod_div a b m = c |
118 mod_div a b m = c |
119 assumes True |
119 assumes True |
120 yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = *) |
120 yields a * b ^(-1) mod m = c <\<Longrightarrow> a mod m = (b * c) mod m*) |
121 fun mod_div a b m = a * (mod_inv b m) mod m |
121 fun mod_div a b m = a * (mod_inv b m) mod m |
122 |
122 |
123 (* required just for one approximation: |
123 (* required just for one approximation: |
124 approx_root :: nat \<Rightarrow> nat |
124 approx_root :: nat \<Rightarrow> nat |
125 approx_root a = r |
125 approx_root a = r |
354 |
354 |
355 argumentns "a b d M P g p" named according to [1] p.93 *) |
355 argumentns "a b d M P g p" named according to [1] p.93 *) |
356 fun try_new_prime_up a b d M P g p = |
356 fun try_new_prime_up a b d M P g p = |
357 if P > M then g else |
357 if P > M then g else |
358 let |
358 let |
359 val p' = p next_prime_not_dvd d |
359 val p = p next_prime_not_dvd d |
360 val g_p = centr_up (( (normalise (mod_up_gcd a b p') p') |
360 val g_p = centr_up ( ( (normalise (mod_up_gcd a b p) p) |
361 %* (d mod p')) |
361 %* (d mod p)) |
362 mod_up p') |
362 mod_up p) |
363 p' |
363 p |
364 in |
364 in |
365 if deg_up g_p < deg_up g |
365 if deg_up g_p < deg_up g |
366 then |
366 then |
367 if (deg_up g_p) = 0 then [1]: unipoly else try_new_prime_up a b d M p g_p p' |
367 if (deg_up g_p) = 0 then [1]: unipoly else try_new_prime_up a b d M p g_p p |
368 else |
368 else |
369 if deg_up g_p <> deg_up g then try_new_prime_up a b d M P g p' else |
369 if deg_up g_p <> deg_up g then try_new_prime_up a b d M P g p else |
370 let |
370 let |
371 val P = P * p' |
371 val P = P * p |
372 val g = centr_up ((chinese_remainder_up (P, p') (g, g_p)) mod_up P) P |
372 val g = centr_up ((chinese_remainder_up (P, p) (g, g_p)) mod_up P) P |
373 in try_new_prime_up a b d M P g p' end |
373 in try_new_prime_up a b d M P g p end |
374 end |
374 end |
375 |
375 |
376 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly |
376 (* HENSEL_lifting_up :: unipoly \<Rightarrow> unipoly \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> unipoly |
377 HENSEL_lifting_up p1 p2 d M p = gcd |
377 HENSEL_lifting_up p1 p2 d M p = gcd |
378 assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and> |
378 assumes d = gcd (lcoeff_up p1, lcoeff_up p2) \<and> |
382 |
382 |
383 argumentns "a b d M p" named according to [1] p.93 *) |
383 argumentns "a b d M p" named according to [1] p.93 *) |
384 fun HENSEL_lifting_up a b d M p = |
384 fun HENSEL_lifting_up a b d M p = |
385 let |
385 let |
386 val p = p next_prime_not_dvd d |
386 val p = p next_prime_not_dvd d |
387 val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p |
387 val g_p = centr_up (((normalise (mod_up_gcd a b p) p) %* (d mod p)) mod_up p) p (*see above*) |
388 in |
388 in |
389 if deg_up g_p = 0 then [1] else |
389 if deg_up g_p = 0 then [1] else |
390 let |
390 let |
391 val g = primitive_up (try_new_prime_up a b d M p g_p p) |
391 val g = primitive_up (try_new_prime_up a b d M p g_p p) |
392 in |
392 in |