1  # BN254 Miller loop implementation for the optimal ate pairing.
  2  #
  3  # Implements the Miller loop core of the BN254 pairing. The ate loop is shorter
  4  # than the Tate loop, improving performance.
  5  #
  6  # References:
  7  # - Aranha et al., "The Realm of the Pairings", 2013.
  8  # - LambdaClass, "How we implemented the BN254 Ate pairing in lambdaworks",
  9  #   2024.
 10  # - Arkworks `ark-bn254` for sparse multiplication helpers.
 11  
 12  # using GrothAlgebra
 13  # using StaticArrays
 14  
 15  # Import what we need
 16  import GrothCurves: conjugate, BN254_PRIME
 17  
 18  # BN254 parameter u = 4965661367192848881
 19  # For optimal ate pairing with Frobenius corrections, we iterate over u only
 20  # The factor of 6 and addition of 2 are handled by the two Frobenius correction lines
 21  const BN254_U = BigInt(4965661367192848881)
 22  const BN254_U_IS_NEGATIVE = false  # u > 0 for BN254
 23  # NAF representation of u for the ate loop (from arkworks)
 24  # This is the signed binary representation where digits are in {0, 1, -1}
 25  const ATE_LOOP_COUNT_NAF = Int8[
 26      0, 0, 0, 1, 0, 1, 0, -1, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1, 0, -1, 0, 0, 0, 1, 0, -1, 0, 0, 0,
 27      0, -1, 0, 0, 1, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 1, 0, -1, 0, 0, 0, -1, 0,
 28      -1, 0, 0, 0, 1, 0, 1, 1,
 29  ]
 30  
 31  # ============================================================================
 32  # p-power Endomorphism Constants for G2 Pairing
 33  # ============================================================================
 34  
 35  # These are the correct p-power endomorphism coefficients from arkworks
 36  # PSI_X = (u+9)^((p-1)/3) for x-coordinate multiplication after Frobenius
 37  const P_POWER_ENDOMORPHISM_COEFF_0 = Fp2Element(
 38      bn254_field(parse(BigInt, "21575463638280843010398324269430826099269044274347216827212613867836435027261")),
 39      bn254_field(parse(BigInt, "10307601595873709700152284273816112264069230130616436755625194854815875713954"))
 40  )
 41  
 42  # PSI_Y = (u+9)^((p-1)/2) for y-coordinate multiplication after Frobenius
 43  const P_POWER_ENDOMORPHISM_COEFF_1 = Fp2Element(
 44      bn254_field(parse(BigInt, "2821565182194536844548159561693502659359617185244120367078079554186484126554")),
 45      bn254_field(parse(BigInt, "3505843767911556378687030309984248845540243509899259641013678093033130930403"))
 46  )
 47  
 48  # Note: For π²(Q), we compose the endomorphism twice rather than using squared coefficients
 49  
 50  """
 51      LineCoeffs
 52  
 53  Raw coefficients for a line function in the Miller loop, before evaluation at P.
 54  For D-twist, the line has form: a*y + b*x + c where x,y are indeterminates.
 55  Note the order difference from M-twist!
 56  """
 57  struct LineCoeffs
 58      a::Fp2Element  # y coefficient (for D-twist)
 59      b::Fp2Element  # x coefficient (for D-twist)
 60      c::Fp2Element  # constant term
 61  end
 62  
 63  """
 64      doubling_step(T::G2Point) -> (G2Point, LineCoeffs)
 65  
 66  Compute the doubling step in the Miller loop.
 67  Returns 2T and the raw line coefficients for the tangent line at T.
 68  
 69  Uses D-twist formulas for BN254 with Jacobian coordinates.
 70  """
 71  function doubling_step(T::G2Point)
 72      # Handle point at infinity
 73      if iszero(T)
 74          return (T, LineCoeffs(zero(Fp2Element), zero(Fp2Element), one(Fp2Element)))
 75      end
 76  
 77      # Get coordinates of T in Jacobian form
 78      X, Y, Z = x_coord(T), y_coord(T), z_coord(T)
 79  
 80      # Precomputations
 81      X2 = X^2
 82      X3 = X2 * X
 83      Y2 = Y^2
 84      Y4 = Y2^2
 85      Z2 = Z^2
 86      Z3 = Z2 * Z
 87      Z6 = Z2^3  # More efficient than Z3 * Z3
 88  
 89      # Step 1: Compute line coefficients for D-twist
 90      # For point T = (X:Y:Z) where x = X/Z², y = Y/Z³
 91      # D-twist coefficients (note the different order from M-twist):
 92      a = -Fp2Element(2) * Y * Z3                          # y coefficient (D-twist)
 93      b = Fp2Element(3) * X2 * Z2                          # x coefficient (D-twist)
 94      c = -X3 + Fp2Element(2) * G2_B_TWIST * Z6            # constant term
 95  
 96      # Step 2: NOW compute 2T
 97      S = Fp2Element(4) * X * Y2
 98      M = Fp2Element(3) * X2
 99      T_val = M^2 - Fp2Element(2) * S
100  
101      # New coordinates for 2T
102      X_2T = T_val
103      Y_2T = M * (S - T_val) - Fp2Element(8) * Y4
104      Z_2T = Fp2Element(2) * Y * Z
105  
106      result_point = G2Point(X_2T, Y_2T, Z_2T)
107      coeffs = LineCoeffs(a, b, c)
108  
109      return (result_point, coeffs)
110  end
111  
112  """
113      addition_step(T::G2Point, Q::G2Point) -> (G2Point, LineCoeffs)
114  
115  Compute the addition step in the Miller loop.
116  Returns T+Q and the raw line coefficients for the line through T and Q.
117  
118  This implements the first equation from page 14 of "The Realm of the Pairings".
119  Uses mixed addition with Q in affine coordinates.
120  """
121  function addition_step(T::G2Point, Q::G2Point)
122      # Handle special cases
123      if iszero(T)
124          return (Q, LineCoeffs(zero(Fp2Element), zero(Fp2Element), one(Fp2Element)))
125      end
126      if iszero(Q)
127          return (T, LineCoeffs(zero(Fp2Element), zero(Fp2Element), one(Fp2Element)))
128      end
129  
130      # If T == Q, use doubling instead
131      if T == Q
132          return doubling_step(T)
133      end
134  
135      # Get coordinates
136      X, Y, Z = x_coord(T), y_coord(T), z_coord(T)
137      xQ, yQ = to_affine(Q)  # Q in affine (common for pairing)
138  
139      # Mixed addition formulas (T projective, Q affine)
140      Z2 = Z^2
141      Z3 = Z2 * Z
142  
143      # Bring Q to same projective denominator as T
144      U2 = xQ * Z2
145      S2 = yQ * Z3
146  
147      # Differences
148      H = U2 - X  # x-coordinate difference
149      R = S2 - Y  # y-coordinate difference
150  
151      # Check if points are equal (shouldn't happen if we checked T == Q above)
152      if iszero(H) && iszero(R)
153          return doubling_step(T)
154      end
155  
156      # Step 1: Compute line coefficients for D-twist BEFORE mutating T
157      # For Jacobian coordinates where T = (X:Y:Z) with x = X/Z², y = Y/Z³
158      # and Q = (xQ, yQ) in affine
159      # D-twist coefficients (note the different order from M-twist):
160      a = -Z * H                       # coefficient for y (D-twist)
161      b = R                            # coefficient for x (D-twist)
162      c = Z * (H * yQ) - R * xQ        # constant term
163  
164      # Step 2: NOW compute T + Q (mutate the point)
165      H2 = H^2
166      H3 = H2 * H
167  
168      X3 = R^2 - H3 - Fp2Element(2) * X * H2
169      Y3 = R * (X * H2 - X3) - Y * H3
170      Z3 = Z * H
171  
172      result_point = G2Point(X3, Y3, Z3)
173      coeffs = LineCoeffs(a, b, c)
174  
175      return (result_point, coeffs)
176  end
177  
178  """
179      evaluate_line(coeffs::LineCoeffs, P::G1Point) -> Fp12Element
180  
181  Evaluate the line function at point P and embed into Fp12.
182  Uses sparse 0-1-4 embedding matching arkworks' mul_by_014.
183  
184  The line a*x + b*y + c is evaluated at P = (x_P, y_P) to give
185  a*x_P + b*y_P + c, which is then embedded into Fp12.
186  """
187  function evaluate_line(coeffs::LineCoeffs, P::G1Point)
188      # Get affine coordinates of P
189      P_x, P_y = if iszero(P)
190          (zero(BN254Field), zero(BN254Field))
191      else
192          to_affine(P)
193      end
194  
195      # Convert to Fp2 for computation
196      xP = Fp2Element(P_x, zero(BN254Field))
197      yP = Fp2Element(P_y, zero(BN254Field))
198  
199      # D-twist embedding - different from M-twist!
200      # For D-twist, the line evaluation gives a sparse Fp12 element
201      # with non-zero entries at positions (0, 3, 4) using mul_by_034 format
202      # Line ℓ(P) = a*yP + b*xP*w + c*w²
203      # This maps to Fp12 coordinates as:
204      #   c0 = (a*yP, 0, 0)
205      #   c1 = (b*xP, c, 0)
206  
207      # Build sparse Fp12 element with correct D-twist mapping
208      c0_fp6 = Fp6Element(coeffs.a * yP, zero(Fp2Element), zero(Fp2Element))
209      c1_fp6 = Fp6Element(coeffs.b * xP, coeffs.c, zero(Fp2Element))
210  
211      return Fp12Element(c0_fp6, c1_fp6)
212  end
213  
214  """
215      p_power_endomorphism(Q::G2Point) -> G2Point
216  
217  Apply the p-power endomorphism to a G2 point.
218  This is the untwist-Frobenius-twist endomorphism: ψ ∘ π_p ∘ ψ^(-1)
219  Maps (x,y) → (x^p * PSI_X, y^p * PSI_Y)
220  """
221  function p_power_endomorphism(Q::G2Point)
222      if iszero(Q)
223          return Q
224      end
225  
226      # Get affine coordinates
227      x, y = to_affine(Q)
228  
229      # Apply Frobenius (conjugation in Fp2)
230      x_frob = conjugate(x)
231      y_frob = conjugate(y)
232  
233      # Multiply by the p-power endomorphism coefficients
234      x_final = x_frob * P_POWER_ENDOMORPHISM_COEFF_0
235      y_final = y_frob * P_POWER_ENDOMORPHISM_COEFF_1
236  
237      return G2Point(x_final, y_final)
238  end
239  
240  """
241      frobenius_g2(Q::G2Point, power::Int) -> G2Point
242  
243  Apply powers of the p-power endomorphism to a G2 point.
244  For power=1: applies p_power_endomorphism once
245  For power=2: applies it twice (π²)
246  """
247  function frobenius_g2(Q::G2Point, power::Int)
248      if iszero(Q)
249          return Q
250      end
251  
252      if power == 1
253          # π(Q) = p-power endomorphism
254          return p_power_endomorphism(Q)
255  
256      elseif power == 2
257          # π²(Q) = apply p-power endomorphism twice
258          # We need to compose: π²(Q) = π(π(Q))
259          intermediate = p_power_endomorphism(Q)
260          return p_power_endomorphism(intermediate)
261  
262      else
263          # For other powers, compose
264          result = Q
265          for _ in 1:power
266              result = p_power_endomorphism(result)
267          end
268          return result
269      end
270  end
271  
272  """
273      miller_loop(::BN254Engine, P::G1Point, Q::G2Point) -> Fp12Element
274  
275  Compute the Miller loop for the optimal ate pairing using the BN254 backend.
276  This is the main computational component of the pairing.
277  
278  The loop iterates over u (not 6u+2), with two Frobenius correction lines
279  at the end to achieve the full optimal ate pairing.
280  """
281  function miller_loop(::BN254Engine, P::G1Point, Q::G2Point)
282      # Handle special cases
283      if iszero(P) || iszero(Q)
284          return one(Fp12Element)
285      end
286  
287      # Initialize
288      T = Q
289      f = one(Fp12Element)
290  
291      # Process NAF representation from MSB to LSB
292      # The loop matches arkworks: iterate from len-1 down to 1, indexing bits at i-1
293      for i in (length(ATE_LOOP_COUNT_NAF)):-1:2
294          # Square f (except on the first iteration)
295          if i != length(ATE_LOOP_COUNT_NAF)
296              f = f^2
297          end
298  
299          # Always perform doubling step
300          T_new, line_coeffs = doubling_step(T)
301          T = T_new
302          f = f * evaluate_line(line_coeffs, P)
303  
304          # Addition/subtraction step based on NAF digit at index i-1
305          bit = ATE_LOOP_COUNT_NAF[i-1]
306          if bit == 1
307              T_new, line_coeffs = addition_step(T, Q)
308              T = T_new
309              f = f * evaluate_line(line_coeffs, P)
310          elseif bit == -1
311              T_new, line_coeffs = addition_step(T, -Q)
312              T = T_new
313              f = f * evaluate_line(line_coeffs, P)
314          end
315          # If bit == 0, we do nothing extra
316      end
317  
318      # For BN254 ate pairing, we need two Frobenius-based correction steps
319      # These are essential for bilinearity to work correctly
320  
321      # Compute Frobenius images of Q
322      # π(Q) applies Frobenius to the x,y coordinates
323      Q_pi = frobenius_g2(Q, 1)   # π(Q)
324      Q_pi2 = frobenius_g2(Q, 2)  # π²(Q)
325  
326      # Correction step 1: Line through T and π(Q)
327      T1, line1 = addition_step(T, Q_pi)
328      f = f * evaluate_line(line1, P)
329  
330      # Correction step 2: Line through T1 and -π²(Q)
331      neg_Q_pi2 = -Q_pi2  # Use proper group negation
332      T2, line2 = addition_step(T1, neg_Q_pi2)
333      f = f * evaluate_line(line2, P)
334  
335      return f
336  end
337  
338  """
339      miller_loop(P::G1Point, Q::G2Point)
340  
341  Convenience overload that reuses the shared BN254 engine.
342  """
343  miller_loop(P::G1Point, Q::G2Point) = miller_loop(BN254_ENGINE, P, Q)
344  
345  # Export functions
346  export LineCoeffs, doubling_step, addition_step, evaluate_line, miller_loop
347  export frobenius_g2, p_power_endomorphism
348  export ATE_LOOP_COUNT_NAF