EllCurveMult[P_,  Q_] := 
      If[P == 
          Q,  {(((((3  P[[1]]^2)/(2  P[[2]]))))^2) - 
            2  P[[1]],  -P[[2]] + (((3  P[[1]]^2)/(2  P[[2]]))) (P[[1]] - ((((((3  P[[1]]^2)/(2  P[[2]]))))^2) - 2  
                          P[[1]]))}, {(((((Q[[2]] - P[[2]])/(Q[[1]] - P[[1]]))))^2) - P[[1]] - 
            Q[[1]], -P[[2]] + (((Q[[2]] - P[[2]])/(Q[[1]] - P[[1]]))) (P[[1]] - ((((((Q[[2]] - P[[2]])/(Q[[1]] - P[[1]]))))^2) - 
                      P[[1]] - Q[[1]]))}]

Diffie-Hellman

Public info:

p = Prime[1000]

s = Random[Integer,  {1,  p - 1}]

s = 6875

Alice chooses a secret a:

a = Random[Integer, {1,  p - 1}]

a = 6055

Bob chooses a secret b:

b = Random[Integer, {1,  p - 1}]

b = 4504

Alice computes α ≡ (s^a) (mod p) using succesive squaring:

s

a

p

α = PowerMod[6875,  6055,  7919]

Bob computes β ≡ (s^b) (mod p) using succesive squaring:

b

β = PowerMod[6875,  4504,  7919]

Alice                                                             Bob
                               α = 4324
                                 ---->
                                 
                               β = 5586
                                 <----

Alice computes (β^a) (mod p):

PowerMod[5586,  6055,  7919]

Bob computes (α^b) (mod p):

PowerMod[4324,  4504,  7919]

Elliptic curves

<< Graphics`

ImplicitPlot[
        (y^2) == (x^3) + x,  {x,  -2,  2}];

ImplicitPlot[
        (y^2) == (x^3) - x,  {x,  -5,  2}];

ImplicitPlot[
        (y^2) == (x^3) + 17,  {x,  -5,  5}];

Show[
        Graphics[{PointSize[ .05], Point[{-2, 3}], Point[{2, 5}],  
            Point[{(1/4),  -(33/8)}],   Point[{(1/4),  (33/8)}],  
            Line[{{(1/4),  8},  {(1/4),  -8}}]}], 
        ImplicitPlot[
          (y^2) == (x^3) + 17,  {x,  -5,  
            4}], 
Plot[(1/2) x + 4,  {x,  -3,  
            4}], 
 Axes -> True,  
        AspectRatio -> Automatic];

Elliptic Curve Diffie-Hellman

Public info:

p = 7

S = {2,  5}

Alice chooses a secret a:

a = Random[Integer, {1,  p - 1}]

a = 3

Bob chooses a secret b:

b = Random[Integer, {1,  p - 1}]

b = 4

Alice computes A ≡ (S^a) (mod p):

S

a

p

EllCurveMult[S, S]

Show[
        Graphics[{PointSize[ .05], Point[{2, 5}],  
            Point[{-(64/25), (59/125)}],  
            Point[{-(64/25), -(59/125)}], 
            Line[{{-(64/25),  8},  {-(64/25),  -8}}]}], 
        ImplicitPlot[
          (y^2) == (x^3) + 17,  {x,  -5,  
            4}], 
Plot[(6/5) (x - 2) + 
            5,  {x,  -3,  4}], 
 
        Axes -> True,  AspectRatio -> Automatic];

EllCurveMult[S, EllCurveMult[S, S]]

A = {Mod[5023,  7]/Mod[3249,  7],  
        Mod[-842480,  7]/Mod[185193,  7]}

Check:

Mod[(5^2), 7] == Mod[(4^3) + 17,  7]

Bob computes B ≡ (S^b) (mod p) using succesive squaring:

b

EllCurveMult[EllCurveMult[S, S],  EllCurveMult[S, S]]

A = {Mod[38194304,  7]/Mod[87025,  7],  
        Mod[-236046706033,  7]/Mod[25672375,  7]}

Alice                                                             Bob
                               A = (4, 5)
                                 ---->
                                 
                               B = (1, 2)
                                 <----

Alice computes (B^a) (mod p):

EllCurveMult[EllCurveMult[{1, 2},  {1, 2}],  {1, 2}]

{Mod[1873,  7]/Mod[1521,  7],  
      Mod[-130870,  7]/Mod[59319,  7]}

Bob computes (A^b) (mod p):

EllCurveMult[EllCurveMult[{4, 5},  {4, 5}],  
      EllCurveMult[{4, 5},  {4, 5}]]

{Mod[5455043344,  7]/Mod[1313700025,  7], 
      Mod[-271860050010047,  7]/Mod[47615057406125,  7]}

PowerMod[6,  -1,  7]

Mod[{2, 5*6},  7]