Botan 3.9.0
Crypto and TLS for C&
pcurves_algos.h
Go to the documentation of this file.
1/*
2* (C) 2024,2025 Jack Lloyd
3*
4* Botan is released under the Simplified BSD License (see license.txt)
5*/
6
7#ifndef BOTAN_PCURVES_ALGOS_H_
8#define BOTAN_PCURVES_ALGOS_H_
9
10#include <botan/types.h>
11#include <botan/internal/ct_utils.h>
12#include <span>
13#include <vector>
14
15namespace Botan {
16
17/**
18* Field inversion concept
19*
20* This concept checks if the curve class supports fe_invert2
21*/
22template <typename C>
23concept curve_supports_fe_invert2 = requires(const typename C::FieldElement& fe) {
24 { C::fe_invert2(fe) } -> std::same_as<typename C::FieldElement>;
25};
26
27/**
28* Field inversion
29*
30* Uses the specialized fe_invert2 if available, or otherwise the standard
31* (FLT-based) field inversion.
32*/
33template <typename C>
34inline constexpr auto invert_field_element(const typename C::FieldElement& fe) {
35 if constexpr(curve_supports_fe_invert2<C>) {
36 return C::fe_invert2(fe) * fe;
37 } else {
38 return fe.invert();
39 }
40}
41
42/**
43* Field square root
44*
45* This concept checks if the curve class supports fe_sqrt
46*/
47template <typename C>
48concept curve_supports_fe_sqrt = requires(const typename C::FieldElement& fe) {
49 { C::fe_sqrt(fe) } -> std::same_as<typename C::FieldElement>;
50};
51
52/**
53* Field square root
54*
55* Uses the specialized fe_sqrt if available, or otherwise the standard
56* square root
57*/
58template <typename C>
59inline constexpr CT::Option<typename C::FieldElement> sqrt_field_element(const typename C::FieldElement& fe) {
60 if constexpr(curve_supports_fe_sqrt<C>) {
61 auto z = C::fe_sqrt(fe);
62 // Zero out the return value if it would otherwise be incorrect
63 const CT::Choice correct = (z.square() == fe);
64 z.conditional_assign(!correct, C::FieldElement::zero());
65 return CT::Option(z, correct);
66 } else {
67 return fe.sqrt();
68 }
69}
70
71/**
72* Convert a projective point into affine
73*/
74template <typename C>
75inline constexpr auto to_affine(const typename C::ProjectivePoint& pt) {
76 // Not strictly required right? - default should work as long
77 // as (0,0) is identity and invert returns 0 on 0
78
79 if constexpr(curve_supports_fe_invert2<C>) {
80 const auto z2_inv = C::fe_invert2(pt.z());
81 const auto z3_inv = z2_inv.square() * pt.z();
82 return typename C::AffinePoint(pt.x() * z2_inv, pt.y() * z3_inv);
83 } else {
84 const auto z_inv = invert_field_element<C>(pt.z());
85 const auto z2_inv = z_inv.square();
86 const auto z3_inv = z_inv * z2_inv;
87 return typename C::AffinePoint(pt.x() * z2_inv, pt.y() * z3_inv);
88 }
89}
90
91/**
92* Convert a projective point into affine and return x coordinate only
93*/
94template <typename C>
95auto to_affine_x(const typename C::ProjectivePoint& pt) {
96 if constexpr(curve_supports_fe_invert2<C>) {
97 return pt.x() * C::fe_invert2(pt.z());
98 } else {
99 const auto z_inv = invert_field_element<C>(pt.z());
100 const auto z2_inv = z_inv.square();
101 return pt.x() * z2_inv;
102 }
103}
104
105template <typename C>
106auto to_affine_batch(std::span<const typename C::ProjectivePoint> projective) {
107 using AffinePoint = typename C::AffinePoint;
108
109 const size_t N = projective.size();
110 std::vector<AffinePoint> affine;
111 affine.reserve(N);
112
113 CT::Choice any_identity = CT::Choice::no();
114
115 for(const auto& pt : projective) {
116 any_identity = any_identity || pt.is_identity();
117 }
118
119 // Conditional acceptable: N is public. State of points is not necessarily
120 // public, but we don't leak which point was the identity. In practice with
121 // the algorithms currently in use, the only time an identity can occur is
122 // during mul2 where the two points g/h have a small relation (ie h = g*k for
123 // some k < 16)
124
125 if(N <= 2 || any_identity.as_bool()) {
126 // If there are identity elements, using the batch inversion gets
127 // tricky. It can be done, but this should be a rare situation so
128 // just punt to the serial conversion if it occurs
129 for(size_t i = 0; i != N; ++i) {
130 affine.push_back(to_affine<C>(projective[i]));
131 }
132 } else {
133 std::vector<typename C::FieldElement> c;
134 c.reserve(N);
135
136 /*
137 Batch projective->affine using Montgomery's trick
138
139 See Algorithm 2.26 in "Guide to Elliptic Curve Cryptography"
140 (Hankerson, Menezes, Vanstone)
141 */
142
143 c.push_back(projective[0].z());
144 for(size_t i = 1; i != N; ++i) {
145 c.push_back(c[i - 1] * projective[i].z());
146 }
147
148 auto s_inv = invert_field_element<C>(c[N - 1]);
149
150 for(size_t i = N - 1; i > 0; --i) {
151 const auto& p = projective[i];
152
153 const auto z_inv = s_inv * c[i - 1];
154 const auto z2_inv = z_inv.square();
155 const auto z3_inv = z_inv * z2_inv;
156
157 s_inv = s_inv * p.z();
158
159 affine.push_back(AffinePoint(p.x() * z2_inv, p.y() * z3_inv));
160 }
161
162 const auto z2_inv = s_inv.square();
163 const auto z3_inv = s_inv * z2_inv;
164 affine.push_back(AffinePoint(projective[0].x() * z2_inv, projective[0].y() * z3_inv));
165 std::reverse(affine.begin(), affine.end());
166 return affine;
167 }
168
169 return affine;
170}
171
172/*
173Projective point addition
174
175https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-1998-cmo-2
176
177Cost: 12M + 4S + 6add + 1*2
178*/
179template <typename ProjectivePoint, typename FieldElement>
180inline constexpr ProjectivePoint point_add(const ProjectivePoint& a, const ProjectivePoint& b) {
181 const auto a_is_identity = a.is_identity();
182 const auto b_is_identity = b.is_identity();
183
184 const auto Z1Z1 = a.z().square();
185 const auto Z2Z2 = b.z().square();
186 const auto U1 = a.x() * Z2Z2;
187 const auto U2 = b.x() * Z1Z1;
188 const auto S1 = a.y() * b.z() * Z2Z2;
189 const auto S2 = b.y() * a.z() * Z1Z1;
190 const auto H = U2 - U1;
191 const auto r = S2 - S1;
192
193 /* Risky conditional
194 *
195 * This implementation uses projective coordinates, which do not have an efficient complete
196 * addition formula. We rely on the design of the multiplication algorithms to avoid doublings.
197 *
198 * This conditional only comes into play for the actual doubling case, not x + (-x) which
199 * is another exceptional case in some circumstances. Here if a == -b then H == 0 && r != 0,
200 * in which case at the end we'll set z to a.z * b.z * H = 0, resulting in the correct
201 * output (the identity element)
202 */
203 if((r.is_zero() && H.is_zero() && !(a_is_identity && b_is_identity)).as_bool()) {
204 return a.dbl();
205 }
206
207 const auto HH = H.square();
208 const auto HHH = H * HH;
209 const auto V = U1 * HH;
210 const auto t2 = r.square();
211 const auto t3 = V + V;
212 const auto t4 = t2 - HHH;
213 auto X3 = t4 - t3;
214 const auto t5 = V - X3;
215 const auto t6 = S1 * HHH;
216 const auto t7 = r * t5;
217 auto Y3 = t7 - t6;
218 const auto t8 = b.z() * H;
219 auto Z3 = a.z() * t8;
220
221 // if a is identity then return b
222 FieldElement::conditional_assign(X3, Y3, Z3, a_is_identity, b.x(), b.y(), b.z());
223
224 // if b is identity then return a
225 FieldElement::conditional_assign(X3, Y3, Z3, b_is_identity, a.x(), a.y(), a.z());
226
227 return ProjectivePoint(X3, Y3, Z3);
228}
229
230template <typename ProjectivePoint, typename AffinePoint, typename FieldElement>
231inline constexpr ProjectivePoint point_add_mixed(const ProjectivePoint& a,
232 const AffinePoint& b,
233 const FieldElement& one) {
234 const auto a_is_identity = a.is_identity();
235 const auto b_is_identity = b.is_identity();
236
237 /*
238 https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-1998-cmo-2
239
240 Cost: 8M + 3S + 6add + 1*2
241 */
242
243 const auto Z1Z1 = a.z().square();
244 const auto U2 = b.x() * Z1Z1;
245 const auto S2 = b.y() * a.z() * Z1Z1;
246 const auto H = U2 - a.x();
247 const auto r = S2 - a.y();
248
249 /* Risky conditional
250 *
251 * This implementation uses projective coordinates, which do not have an efficient complete
252 * addition formula. We rely on the design of the multiplication algorithms to avoid doublings.
253 *
254 * This conditional only comes into play for the actual doubling case, not x + (-x) which
255 * is another exceptional case in some circumstances. Here if a == -b then H == 0 && r != 0,
256 * in which case at the end we'll set z to a.z * H = 0, resulting in the correct output
257 * (the identity element)
258 */
259 if((r.is_zero() && H.is_zero() && !(a_is_identity && b_is_identity)).as_bool()) {
260 return a.dbl();
261 }
262
263 const auto HH = H.square();
264 const auto HHH = H * HH;
265 const auto V = a.x() * HH;
266 const auto t2 = r.square();
267 const auto t3 = V + V;
268 const auto t4 = t2 - HHH;
269 auto X3 = t4 - t3;
270 const auto t5 = V - X3;
271 const auto t6 = a.y() * HHH;
272 const auto t7 = r * t5;
273 auto Y3 = t7 - t6;
274 auto Z3 = a.z() * H;
275
276 // if a is identity then return b
277 FieldElement::conditional_assign(X3, Y3, Z3, a_is_identity, b.x(), b.y(), one);
278
279 // if b is identity then return a
280 FieldElement::conditional_assign(X3, Y3, Z3, b_is_identity, a.x(), a.y(), a.z());
281
282 return ProjectivePoint(X3, Y3, Z3);
283}
284
285template <typename ProjectivePoint, typename AffinePoint, typename FieldElement>
286inline constexpr ProjectivePoint point_add_or_sub_mixed(const ProjectivePoint& a,
287 const AffinePoint& b,
288 CT::Choice sub,
289 const FieldElement& one) {
290 const auto a_is_identity = a.is_identity();
291 const auto b_is_identity = b.is_identity();
292
293 /*
294 https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-1998-cmo-2
295
296 Cost: 8M + 3S + 6add + 1*2
297 */
298
299 auto by = b.y();
300 by.conditional_assign(sub, by.negate());
301
302 const auto Z1Z1 = a.z().square();
303 const auto U2 = b.x() * Z1Z1;
304 const auto S2 = by * a.z() * Z1Z1;
305 const auto H = U2 - a.x();
306 const auto r = S2 - a.y();
307
308 /* Risky conditional
309 *
310 * This implementation uses projective coordinates, which do not have an efficient complete
311 * addition formula. We rely on the design of the multiplication algorithms to avoid doublings.
312 *
313 * This conditional only comes into play for the actual doubling case, not x + (-x) which
314 * is another exceptional case in some circumstances. Here if a == -b then H == 0 && r != 0,
315 * in which case at the end we'll set z to a.z * H = 0, resulting in the correct output
316 * (the identity element)
317 */
318 if((r.is_zero() && H.is_zero() && !(a_is_identity && b_is_identity)).as_bool()) {
319 return a.dbl();
320 }
321
322 const auto HH = H.square();
323 const auto HHH = H * HH;
324 const auto V = a.x() * HH;
325 const auto t2 = r.square();
326 const auto t3 = V + V;
327 const auto t4 = t2 - HHH;
328 auto X3 = t4 - t3;
329 const auto t5 = V - X3;
330 const auto t6 = a.y() * HHH;
331 const auto t7 = r * t5;
332 auto Y3 = t7 - t6;
333 auto Z3 = a.z() * H;
334
335 // if a is identity then return b
336 FieldElement::conditional_assign(X3, Y3, Z3, a_is_identity, b.x(), by, one);
337
338 // if b is identity then return a
339 FieldElement::conditional_assign(X3, Y3, Z3, b_is_identity, a.x(), a.y(), a.z());
340
341 return ProjectivePoint(X3, Y3, Z3);
342}
343
344/*
345Point doubling
346
347Using https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-1998-cmo-2
348
349Cost (generic A): 4M + 6S + 4A + 2*2 + 1*3 + 1*4 + 1*8
350Cost (A == -3): 4M + 4S + 5A + 2*2 + 1*3 + 1*4 + 1*8
351Cost (A == 0): 3M + 4S + 3A + 2*2 + 1*3 + 1*4 + 1*8
352*/
353
354template <typename ProjectivePoint>
355inline constexpr ProjectivePoint dbl_a_minus_3(const ProjectivePoint& pt) {
356 /*
357 if a == -3 then
358 3*x^2 + a*z^4 == 3*x^2 - 3*z^4 == 3*(x^2-z^4) == 3*(x-z^2)*(x+z^2)
359 */
360 const auto z2 = pt.z().square();
361 const auto m = (pt.x() - z2).mul3() * (pt.x() + z2);
362
363 // Remaining cost: 3M + 3S + 3A + 2*2 + 1*4 + 1*8
364 const auto y2 = pt.y().square();
365 const auto s = pt.x().mul4() * y2;
366 const auto nx = m.square() - s.mul2();
367 const auto ny = m * (s - nx) - y2.square().mul8();
368 const auto nz = pt.y().mul2() * pt.z();
369
370 return ProjectivePoint(nx, ny, nz);
371}
372
373template <typename ProjectivePoint>
374inline constexpr ProjectivePoint dbl_a_zero(const ProjectivePoint& pt) {
375 // If a == 0 then 3*x^2 + a*z^4 == 3*x^2
376 // Cost: 1S + 1*3
377 const auto m = pt.x().square().mul3();
378
379 // Remaining cost: 3M + 3S + 3A + 2*2 + 1*4 + 1*8
380 const auto y2 = pt.y().square();
381 const auto s = pt.x().mul4() * y2;
382 const auto nx = m.square() - s.mul2();
383 const auto ny = m * (s - nx) - y2.square().mul8();
384 const auto nz = pt.y().mul2() * pt.z();
385
386 return ProjectivePoint(nx, ny, nz);
387}
388
389template <typename ProjectivePoint, typename FieldElement>
390inline constexpr ProjectivePoint dbl_generic(const ProjectivePoint& pt, const FieldElement& A) {
391 // Cost: 1M + 3S + 1A + 1*3
392 const auto z2 = pt.z().square();
393 const auto m = pt.x().square().mul3() + A * z2.square();
394
395 // Remaining cost: 3M + 3S + 3A + 2*2 + 1*4 + 1*8
396 const auto y2 = pt.y().square();
397 const auto s = pt.x().mul4() * y2;
398 const auto nx = m.square() - s.mul2();
399 const auto ny = m * (s - nx) - y2.square().mul8();
400 const auto nz = pt.y().mul2() * pt.z();
401
402 return ProjectivePoint(nx, ny, nz);
403}
404
405/*
406Repeated doubling using an adaptation of Algorithm 3.23 in
407"Guide To Elliptic Curve Cryptography" (Hankerson, Menezes, Vanstone)
408
409Curiously the book gives the algorithm only for A == -3, but
410the largest gains come from applying it to the generic A case,
411where it saves 2 squarings per iteration.
412
413For A == 0
414Pay 1*2 + 1half to save n*(1*4 + 1*8)
415
416For A == -3:
417Pay 2S + 1*2 + 1half to save n*(1A + 1*4 + 1*8) + 1M
418
419For generic A:
420Pay 2S + 1*2 + 1half to save n*(2S + 1*4 + 1*8)
421
422The value of n is assumed to be public and should be a constant
423*/
424template <typename ProjectivePoint>
425inline constexpr ProjectivePoint dbl_n_a_minus_3(const ProjectivePoint& pt, size_t n) {
426 auto nx = pt.x();
427 auto ny = pt.y().mul2();
428 auto nz = pt.z();
429 auto w = nz.square().square();
430
431 // Conditional ok: loop iteration count is public
432 while(n > 0) {
433 const auto ny2 = ny.square();
434 const auto ny4 = ny2.square();
435 const auto t1 = (nx.square() - w).mul3();
436 const auto t2 = nx * ny2;
437 nx = t1.square() - t2.mul2();
438 nz *= ny;
439 ny = t1 * (t2 - nx).mul2() - ny4;
440 n--;
441 // Conditional ok: loop iteration count is public
442 if(n > 0) {
443 w *= ny4;
444 }
445 }
446 return ProjectivePoint(nx, ny.div2(), nz);
447}
448
449template <typename ProjectivePoint>
450inline constexpr ProjectivePoint dbl_n_a_zero(const ProjectivePoint& pt, size_t n) {
451 auto nx = pt.x();
452 auto ny = pt.y().mul2();
453 auto nz = pt.z();
454
455 // Conditional ok: loop iteration count is public
456 while(n > 0) {
457 const auto ny2 = ny.square();
458 const auto ny4 = ny2.square();
459 const auto t1 = nx.square().mul3();
460 const auto t2 = nx * ny2;
461 nx = t1.square() - t2.mul2();
462 nz *= ny;
463 ny = t1 * (t2 - nx).mul2() - ny4;
464 n--;
465 }
466 return ProjectivePoint(nx, ny.div2(), nz);
467}
468
469template <typename ProjectivePoint, typename FieldElement>
470inline constexpr ProjectivePoint dbl_n_generic(const ProjectivePoint& pt, const FieldElement& A, size_t n) {
471 auto nx = pt.x();
472 auto ny = pt.y().mul2();
473 auto nz = pt.z();
474 auto w = nz.square().square() * A;
475
476 // Conditional ok: loop iteration count is public
477 while(n > 0) {
478 const auto ny2 = ny.square();
479 const auto ny4 = ny2.square();
480 const auto t1 = nx.square().mul3() + w;
481 const auto t2 = nx * ny2;
482 nx = t1.square() - t2.mul2();
483 nz *= ny;
484 ny = t1 * (t2 - nx).mul2() - ny4;
485 n--;
486 // Conditional ok: loop iteration count is public
487 if(n > 0) {
488 w *= ny4;
489 }
490 }
491 return ProjectivePoint(nx, ny.div2(), nz);
492}
493
494} // namespace Botan
495
496#endif
static constexpr Choice no()
Definition ct_utils.h:335
constexpr bool as_bool() const
Definition ct_utils.h:357
auto to_affine_batch(std::span< const typename C::ProjectivePoint > projective)
auto to_affine_x(const typename C::ProjectivePoint &pt)
constexpr CT::Option< typename C::FieldElement > sqrt_field_element(const typename C::FieldElement &fe)
constexpr ProjectivePoint dbl_n_generic(const ProjectivePoint &pt, const FieldElement &A, size_t n)
constexpr ProjectivePoint dbl_a_minus_3(const ProjectivePoint &pt)
constexpr ProjectivePoint dbl_n_a_zero(const ProjectivePoint &pt, size_t n)
constexpr ProjectivePoint dbl_a_zero(const ProjectivePoint &pt)
constexpr auto to_affine(const typename C::ProjectivePoint &pt)
constexpr ProjectivePoint point_add_mixed(const ProjectivePoint &a, const AffinePoint &b, const FieldElement &one)
constexpr ProjectivePoint point_add_or_sub_mixed(const ProjectivePoint &a, const AffinePoint &b, CT::Choice sub, const FieldElement &one)
constexpr auto invert_field_element(const typename C::FieldElement &fe)
constexpr ProjectivePoint dbl_n_a_minus_3(const ProjectivePoint &pt, size_t n)
constexpr ProjectivePoint dbl_generic(const ProjectivePoint &pt, const FieldElement &A)
constexpr ProjectivePoint point_add(const ProjectivePoint &a, const ProjectivePoint &b)