fork(6) download
copy 
  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
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
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
https://ideone.com/SSB68P
language:
Octave (octave 4.4.1)
created:
8 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!