ordp.sage 1.9 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394
  1. from tabulate import tabulate
  2. # This is a more usable version of valuate.sage, less instructional
  3. # Return components for basis
  4. def decomp(f, basis):
  5. b0, b1, b2 = basis
  6. a0, r = f.quo_rem(b0)
  7. a1, r = r.quo_rem(b1)
  8. a2, r = r.quo_rem(b2)
  9. assert r == 0
  10. return [a0, a1, a2]
  11. def comp(comps, basis):
  12. return sum(a*b for a, b in zip(comps, basis))
  13. def apply_reduction(a, g, Ef, Eg):
  14. assert a[2] == 0
  15. a[0] = a[0]*Eg + a[1]*Ef
  16. a[1] = 0
  17. g[0] *= Eg
  18. # EC_A, EC_B must be defined before calling this function
  19. def ordp(P, original_f, debug=False):
  20. EC = y^2 - x^3 - EC_A*x - EC_B
  21. Px, Py = P
  22. if Py != 0:
  23. b0, b1, b2 = [(x - Px), (y - Py), 1]
  24. # so we can replace (y - Py) with this
  25. Ef = b0^2 + binomial(3,2)*Px*b0^1 + (3*Px^2 + EC_A)
  26. Eg = (y + Py)
  27. assert EC == b1*Eg - b0*Ef
  28. else:
  29. b0, b1, b2 = [y, x, 1]
  30. # we can replace x with this
  31. Ef = b0
  32. Eg = x^2 + EC_A + EC_B
  33. basis = [b0, b1, b2]
  34. k = 0
  35. a = [original_f, 0, 0]
  36. g = [1]
  37. table = []
  38. table.append(("", "a", "g", "k"))
  39. def log(step_name, a, g, k):
  40. table.append((step_name, str(a), str(g), k))
  41. log("start", a, g, k)
  42. while True:
  43. f = a[0]
  44. a = decomp(f, basis)
  45. log("decomp", a, g, k)
  46. # Check remainder
  47. if a[2] != 0:
  48. break
  49. # We can apply a reduction
  50. k += 1
  51. apply_reduction(a, g, Ef, Eg)
  52. log("reduce", a, g, k)
  53. if debug:
  54. print(tabulate(table))
  55. u = b0
  56. f = comp(a, basis)
  57. g = g[0]
  58. assert u(Px, Py) == 0
  59. assert f(Px, Py) != 0
  60. assert g(Px, Py) != 0
  61. if debug:
  62. print(f"u = {u}")
  63. print(f"f = {f}")
  64. print(f"g = {g}")
  65. return k
  66. #K.<x, y> = GF(11)[]
  67. #EC_A = 4
  68. #EC_B = 0
  69. #P = (2, 4)
  70. #f = y - 2*x
  71. #k = ordp(P, f, debug=True)
  72. #print(f"k = {k}")