fork(4) download
  1. %%
  2. % This script is for adding two points on an elliptic curve over a finite
  3. % field of size of a prime p that is given in the short weierstrass form.
  4. % note that p should not be excessively large, as many computations here
  5. % rely on bruteforcing.
  6.  
  7.  
  8. clc;clear;
  9. % Field characteristic p ~= 2,3
  10. p = 17;
  11. % Parameters of the elliptic curve y^2 = x^3 + Ax + B
  12. A = 2;
  13. B = 2;
  14.  
  15. % Points you want to add
  16.  
  17. P = [5,1]; % two coordinates [x,y] or just Inf for the neutral element
  18. Q = [7,6]; % write P=Inf; Q=Inf; to find all points on the curve
  19.  
  20. % DO NOT CHANGE ANYTHING BELOW THIS LINE
  21. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  22.  
  23.  
  24. F=@(x,y) mod(x.^3 + A*x + B - y.^2,p);
  25.  
  26.  
  27. %% perform checks
  28. assert(mod(4*A^3 + 27*B^2, p) ~=0,'A,B do not define an elliptic curve');
  29. assert(p~=2 && p~=3 && isprime(p),'p=2,3 or p is not prime');
  30. assert(P(1) == Inf || F(P(1),P(2)) == 0 , 'P is not on the curve');
  31. assert(Q(1) == Inf || F(Q(1),Q(2)) == 0 , 'Q is not on the curve');
  32.  
  33. %% find all points
  34. all_points = [];
  35. for x=0:p-1;
  36. for y=0:p-1;
  37. if F(x,y)==0
  38. all_points = [all_points; x,y];
  39. end
  40. end
  41. end
  42. disp('all points on the curve are the following:');
  43. disp(' inf');
  44. disp(all_points)
  45.  
  46.  
  47. %% tree
  48.  
  49. if P == Inf % P is neutral element
  50. disp('Result is:')
  51. disp(Q)
  52. elseif Q == Inf % Q is neutral element
  53. disp('Result is:')
  54. disp(P)
  55. else
  56. disp('apparently none of the summands is the neutral element...');
  57. if P(1) == Q(1) && mod(P(2)+Q(2),p) == 0 % P+Q == 0
  58. disp('Result is:')
  59. disp(Inf)
  60. else % P ~= -Q -> we have a nonvertical line
  61. disp('P and Q are apparently not additive inverses to each other...');
  62. x0 = P(1);
  63. y0 = P(2);
  64. x1 = Q(1);
  65. y1 = Q(2);
  66. if x0 == x1 && y0 == y1; % P == Q -> calcualte tangent
  67. disp('we have to double P == Q ...');
  68. % m = (3* x0^2 + A)/ (2 * y0);
  69. % find modular inverse of (2*y0)
  70. k = find( mod( (1:p) * 2*y0, p) == 1);
  71. m = mod((3*x0^2 + A) * k, p);
  72. q = mod( (-x0^3+A*x0 + 2*B) * k, p);
  73.  
  74.  
  75. else % calculate line (nonvertical) throu those two points
  76. disp('we can calculate a nonvertical line throu P and Q...');
  77. % m = (y1-y0) / (x1-x0)
  78. % find modular inverse of (x1-x0)
  79. k = find( mod( (1:p) * (x1-x0), p) == 1);
  80. m = mod( (y1-y0) * k, p);
  81. % q = y0 - m*x0
  82. q = mod( y0 - m*x0, p);
  83. end
  84. % calculate the zeros of this polynomial:
  85. % P = x^3+Ax+B - (m*x+q)^2
  86. % = x^3+Ax+B - (m^2 * x^2 + 2*m*q*x+q^2)
  87. % = x^3 + (-m^2) * x^2 + (A-2*m*q) * x + (B-q^2)
  88. pl = mod([1, -m^2, A-2*m*q, B-q^2], p);
  89. pl_roots = mod(find(mod(polyval(pl,1:p),p)==0,3),p);
  90.  
  91. if numel(pl_roots) == 3
  92. disp('the polynomial has 3 distinct roots...');
  93. X = pl_roots( pl_roots ~= x0 & pl_roots ~= x1);
  94. Y = mod(m*X+q, p);
  95. disp('Result is:');
  96. disp(mod([X,-Y],p));
  97. elseif numel(pl_roots) == 2;
  98. disp('the polynomial has 2 distinct roots...');
  99. if all( mod(conv(conv([1,-x0],[1,-x1]),[1,-pl_roots(1)]),p) == pl)
  100. disp('Result is:')
  101. X = pl_roots(1);
  102. Y = mod(m*X+q, p);
  103. disp(mod([X,-Y],p));
  104. elseif all( mod(conv(conv([1,-x0],[1,-x1]),[1,-pl_roots(2)]),p) == pl)
  105. disp('Result is:');
  106. X = pl_roots(2);
  107. Y = mod(m*X+q, p);
  108. disp(mod([X,-Y],p));
  109. else
  110. disp('ERROR: something went horribly wrong: could not reconstruct polynomial');
  111. end
  112. elseif numel(pl_roots) == 1;
  113. disp('the polynomial has 1 distinct root...');
  114. if all( mod(conv(conv([1,-x0],[1,-x1]),[1,-pl_roots]),p) == pl)
  115. disp('Result is:')
  116. X = pl_roots;
  117. Y = mod(m*X+q, p);
  118. disp(mod([X,-Y],p))
  119. else
  120. disp('ERROR: somethign went horribly wrong: could not reconstruct polynomial');
  121. end
  122.  
  123. else
  124. disp('ERROR: something went horribly wrong (not 1 or 2 or 3 roots)');
  125. end
  126.  
  127. end
  128. end
  129.  
  130.  
  131.  
Success #stdin #stdout 0.5s 122240KB
stdin
Standard input is empty
stdout
all points on the curve are the following:
     inf
    0    6
    0   11
    3    1
    3   16
    5    1
    5   16
    6    3
    6   14
    7    6
    7   11
    9    1
    9   16
   10    6
   10   11
   13    7
   13   10
   16    4
   16   13
apparently none of the summands is the neutral element...
P and Q are apparently not additive inverses to each other...
we can calculate a nonvertical line throu P and Q...
the polynomial has 2 distinct roots...
Result is:
    7   11