Line data    Source code

       1             : #include "./fd_bn254.h"
       2             : 
       3             : /* G1 */
       4             : 
       5             : static inline int
       6     6500019 : fd_bn254_g1_is_zero( fd_bn254_g1_t const * p ) {
       7     6500019 :   return fd_bn254_fp_is_zero( &p->Z );
       8     6500019 : }
       9             : 
      10             : static inline fd_bn254_g1_t *
      11             : fd_bn254_g1_set( fd_bn254_g1_t *       r,
      12       60213 :                  fd_bn254_g1_t const * p ) {
      13       60213 :   fd_bn254_fp_set( &r->X, &p->X );
      14       60213 :   fd_bn254_fp_set( &r->Y, &p->Y );
      15       60213 :   fd_bn254_fp_set( &r->Z, &p->Z );
      16       60213 :   return r;
      17       60213 : }
      18             : 
      19             : static inline fd_bn254_g1_t *
      20          54 : fd_bn254_g1_set_zero( fd_bn254_g1_t * r ) {
      21             :   // fd_bn254_fp_set_zero( &r->X );
      22             :   // fd_bn254_fp_set_zero( &r->Y );
      23          54 :   fd_bn254_fp_set_zero( &r->Z );
      24          54 :   return r;
      25          54 : }
      26             : 
      27             : static inline fd_bn254_g1_t *
      28             : fd_bn254_g1_to_affine( fd_bn254_g1_t *       r,
      29       60165 :                        fd_bn254_g1_t const * p ) {
      30       60165 :   if( FD_UNLIKELY( fd_bn254_fp_is_zero( &p->Z ) || fd_bn254_fp_is_one( &p->Z ) ) ) {
      31       30060 :     return fd_bn254_g1_set( r, p );
      32       30060 :   }
      33             : 
      34       30105 :   fd_bn254_fp_t iz[1], iz2[1];
      35       30105 :   fd_bn254_fp_inv( iz, &p->Z );
      36       30105 :   fd_bn254_fp_sqr( iz2, iz );
      37             : 
      38             :   /* X / Z^2, Y / Z^3 */
      39       30105 :   fd_bn254_fp_mul( &r->X, &p->X, iz2 );
      40       30105 :   fd_bn254_fp_mul( &r->Y, &p->Y, iz2 );
      41       30105 :   fd_bn254_fp_mul( &r->Y, &r->Y, iz );
      42       30105 :   fd_bn254_fp_set_one( &r->Z );
      43       30105 :   return r;
      44       60165 : }
      45             : 
      46             : uchar *
      47             : fd_bn254_g1_tobytes( uchar                 out[64],
      48             :                      fd_bn254_g1_t const * p,
      49       60183 :                      int                   big_endian ) {
      50       60183 :   if( FD_UNLIKELY( fd_bn254_g1_is_zero( p ) ) ) {
      51          18 :     fd_memset( out, 0, 64UL );
      52             :     /* no flags */
      53          18 :     return out;
      54          18 :   }
      55             : 
      56       60165 :   fd_bn254_g1_t r[1];
      57       60165 :   fd_bn254_g1_to_affine( r, p );
      58             : 
      59       60165 :   fd_bn254_fp_from_mont( &r->X, &r->X );
      60       60165 :   fd_bn254_fp_from_mont( &r->Y, &r->Y );
      61             : 
      62       60165 :   fd_bn254_fp_tobytes_nm( &out[ 0], &r->X, big_endian );
      63       60165 :   fd_bn254_fp_tobytes_nm( &out[32], &r->Y, big_endian );
      64             :   /* no flags */
      65       60165 :   return out;
      66       60183 : }
      67             : 
      68             : /* fd_bn254_g1_affine_add computes r = p + q.
      69             :    Both p, q are affine, i.e. Z==1. */
      70             : fd_bn254_g1_t *
      71             : fd_bn254_g1_affine_add( fd_bn254_g1_t *       r,
      72             :                         fd_bn254_g1_t const * p,
      73       60180 :                         fd_bn254_g1_t const * q ) {
      74             :   /* p==0, return q */
      75       60180 :   if( FD_UNLIKELY( fd_bn254_g1_is_zero( p ) ) ) {
      76          21 :     return fd_bn254_g1_set( r, q );
      77          21 :   }
      78             :   /* q==0, return p */
      79       60159 :   if( FD_UNLIKELY( fd_bn254_g1_is_zero( q ) ) ) {
      80           9 :     return fd_bn254_g1_set( r, p );
      81           9 :   }
      82             : 
      83       60150 :   fd_bn254_fp_t lambda[1], x[1], y[1];
      84             : 
      85             :   /* same X, either the points are equal or opposite */
      86       60150 :   if( fd_bn254_fp_eq( &p->X, &q->X ) ) {
      87           6 :     if( fd_bn254_fp_eq( &p->Y, &q->Y ) ) {
      88             :       /* p==q => point double: lambda = 3 * x1^2 / (2 * y1) */
      89           6 :       fd_bn254_fp_sqr( x, &p->X ); /* x =   x1^2 */
      90           6 :       fd_bn254_fp_add( y, x, x );  /* y = 2 x1^2 */
      91           6 :       fd_bn254_fp_add( x, x, y );  /* x = 3 x1^2 */
      92           6 :       fd_bn254_fp_add( y, &p->Y, &p->Y );
      93           6 :       fd_bn254_fp_inv( lambda, y );
      94           6 :       fd_bn254_fp_mul( lambda, lambda, x );
      95           6 :     } else {
      96             :       /* p==-q => r=0 */
      97             :       /* COV: this may never happen with real data */
      98           0 :       return fd_bn254_g1_set_zero( r );
      99           0 :     }
     100       60144 :   } else {
     101             :     /* point add: lambda = (y1 - y2) / (x1 - x2) */
     102       60144 :     fd_bn254_fp_sub( x, &p->X, &q->X );
     103       60144 :     fd_bn254_fp_sub( y, &p->Y, &q->Y );
     104       60144 :     fd_bn254_fp_inv( lambda, x );
     105       60144 :     fd_bn254_fp_mul( lambda, lambda, y );
     106       60144 :   }
     107             : 
     108             :   /* x3 = lambda^2 - x1 - x2 */
     109       60150 :   fd_bn254_fp_sqr( x, lambda );
     110       60150 :   fd_bn254_fp_sub( x, x, &p->X );
     111       60150 :   fd_bn254_fp_sub( x, x, &q->X );
     112             : 
     113             :   /* y3 = lambda * (x1 - x3) - y1 */
     114       60150 :   fd_bn254_fp_sub( y, &p->X, x );
     115       60150 :   fd_bn254_fp_mul( y, y, lambda );
     116       60150 :   fd_bn254_fp_sub( y, y, &p->Y );
     117             : 
     118       60150 :   fd_bn254_fp_set( &r->X, x );
     119       60150 :   fd_bn254_fp_set( &r->Y, y );
     120       60150 :   fd_bn254_fp_set_one( &r->Z );
     121       60150 :   return r;
     122       60150 : }
     123             : 
     124             : /* fd_bn254_g1_dbl computes r = 2p.
     125             :    https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-dbl-2009-l */
     126             : fd_bn254_g1_t *
     127             : fd_bn254_g1_dbl( fd_bn254_g1_t *       r,
     128     3790419 :                  fd_bn254_g1_t const * p ) {
     129             :   /* p==0, return 0 */
     130     3790419 :   if( FD_UNLIKELY( fd_bn254_g1_is_zero( p ) ) ) {
     131           0 :     return fd_bn254_g1_set_zero( r );
     132           0 :   }
     133             : 
     134     3790419 :   fd_bn254_fp_t a[1], b[1], c[1];
     135     3790419 :   fd_bn254_fp_t d[1], e[1], f[1];
     136             : 
     137             :   /* A = X1^2 */
     138     3790419 :   fd_bn254_fp_sqr( a, &p->X );
     139             :   /* B = Y1^2 */
     140     3790419 :   fd_bn254_fp_sqr( b, &p->Y );
     141             :   /* C = B^2 */
     142     3790419 :   fd_bn254_fp_sqr( c, b );
     143             :   /* D = 2*((X1+B)^2-A-C)
     144             :      (X1+B)^2 = X1^2 + 2*X1*B + B^2
     145             :      D = 2*(X1^2 + 2*X1*B + B^2 - A    - C)
     146             :      D = 2*(X1^2 + 2*X1*B + B^2 - X1^2 - B^2)
     147             :             ^               ^     ^      ^
     148             :             |---------------|-----|      |
     149             :                             |------------|
     150             :      These terms cancel each other out, and we're left with:
     151             :      D = 2*(2*X1*B) */
     152     3790419 :   fd_bn254_fp_mul( d, &p->X, b );
     153     3790419 :   fd_bn254_fp_add( d, d, d );
     154     3790419 :   fd_bn254_fp_add( d, d, d );
     155             :   /* E = 3*A */
     156     3790419 :   fd_bn254_fp_add( e, a, a );
     157     3790419 :   fd_bn254_fp_add( e, a, e );
     158             :   /* F = E^2 */
     159     3790419 :   fd_bn254_fp_sqr( f, e );
     160             :   /* X3 = F-2*D */
     161     3790419 :   fd_bn254_fp_add( &r->X, d, d );
     162     3790419 :   fd_bn254_fp_sub( &r->X, f, &r->X );
     163             :   /* Z3 = (Y1+Z1)^2-YY-ZZ
     164             :      note: compute Z3 before Y3 because it depends on p->Y,
     165             :      that might be overwritten if r==p. */
     166             :   /* Z3 = 2*Y1*Z1 */
     167     3790419 :   fd_bn254_fp_mul( &r->Z, &p->Y, &p->Z );
     168     3790419 :   fd_bn254_fp_add( &r->Z, &r->Z, &r->Z );
     169             :   /* Y3 = E*(D-X3)-8*C */
     170     3790419 :   fd_bn254_fp_sub( &r->Y, d, &r->X );
     171     3790419 :   fd_bn254_fp_mul( &r->Y, e, &r->Y );
     172     3790419 :   fd_bn254_fp_add( c, c, c ); /* 2*c */
     173     3790419 :   fd_bn254_fp_add( c, c, c ); /* 4*y */
     174     3790419 :   fd_bn254_fp_add( c, c, c ); /* 8*y */
     175     3790419 :   fd_bn254_fp_sub( &r->Y, &r->Y, c );
     176     3790419 :   return r;
     177     3790419 : }
     178             : 
     179             : /* fd_bn254_g1_add_mixed computes r = p + q, when q->Z==1.
     180             :    http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-madd-2007-bl */
     181             : fd_bn254_g1_t *
     182             : fd_bn254_g1_add_mixed( fd_bn254_g1_t *       r,
     183             :                        fd_bn254_g1_t const * p,
     184     2376864 :                        fd_bn254_g1_t const * q ) {
     185             :   /* p==0, return q */
     186     2376864 :   if( FD_UNLIKELY( fd_bn254_g1_is_zero( p ) ) ) {
     187           0 :     return fd_bn254_g1_set( r, q );
     188           0 :   }
     189     2376864 :   fd_bn254_fp_t zz[1], u2[1], s2[1];
     190     2376864 :   fd_bn254_fp_t h[1], hh[1];
     191     2376864 :   fd_bn254_fp_t i[1], j[1];
     192     2376864 :   fd_bn254_fp_t rr[1], v[1];
     193             :   /* Z1Z1 = Z1^2 */
     194     2376864 :   fd_bn254_fp_sqr( zz, &p->Z );
     195             :   /* U2 = X2*Z1Z1 */
     196     2376864 :   fd_bn254_fp_mul( u2, &q->X, zz );
     197             :   /* S2 = Y2*Z1*Z1Z1 */
     198     2376864 :   fd_bn254_fp_mul( s2, &q->Y, &p->Z );
     199     2376864 :   fd_bn254_fp_mul( s2, s2, zz );
     200             : 
     201             :   /* if p==q, call fd_bn254_g1_dbl */
     202     2376864 :   if( FD_UNLIKELY( fd_bn254_fp_eq( u2, &p->X ) && fd_bn254_fp_eq( s2, &p->Y ) ) ) {
     203             :     /* COV: this may never happen with real data */
     204           0 :     return fd_bn254_g1_dbl( r, p );
     205           0 :   }
     206             : 
     207             :   /* H = U2-X1 */
     208     2376864 :   fd_bn254_fp_sub( h, u2, &p->X );
     209             :   /* HH = H^2 */
     210     2376864 :   fd_bn254_fp_sqr( hh, h );
     211             :   /* I = 4*HH */
     212     2376864 :   fd_bn254_fp_add( i, hh, hh );
     213     2376864 :   fd_bn254_fp_add( i, i, i );
     214             :   /* J = H*I */
     215     2376864 :   fd_bn254_fp_mul( j, h, i );
     216             :   /* r = 2*(S2-Y1) */
     217     2376864 :   fd_bn254_fp_sub( rr, s2, &p->Y );
     218     2376864 :   fd_bn254_fp_add( rr, rr, rr );
     219             :   /* V = X1*I */
     220     2376864 :   fd_bn254_fp_mul( v, &p->X, i );
     221             :   /* X3 = r^2-J-2*V */
     222     2376864 :   fd_bn254_fp_sqr( &r->X, rr );
     223     2376864 :   fd_bn254_fp_sub( &r->X, &r->X, j );
     224     2376864 :   fd_bn254_fp_sub( &r->X, &r->X, v );
     225     2376864 :   fd_bn254_fp_sub( &r->X, &r->X, v );
     226             :   /* Y3 = r*(V-X3)-2*Y1*J
     227             :      note: i no longer used */
     228     2376864 :   fd_bn254_fp_mul( i, &p->Y, j ); /* i =   Y1*J */
     229     2376864 :   fd_bn254_fp_add( i, i, i );     /* i = 2*Y1*J */
     230     2376864 :   fd_bn254_fp_sub( &r->Y, v, &r->X );
     231     2376864 :   fd_bn254_fp_mul( &r->Y, &r->Y, rr );
     232     2376864 :   fd_bn254_fp_sub( &r->Y, &r->Y, i );
     233             :   /* Z3 = (Z1+H)^2-Z1Z1-HH */
     234     2376864 :   fd_bn254_fp_add( &r->Z, &p->Z, h );
     235     2376864 :   fd_bn254_fp_sqr( &r->Z, &r->Z );
     236     2376864 :   fd_bn254_fp_sub( &r->Z, &r->Z, zz );
     237     2376864 :   fd_bn254_fp_sub( &r->Z, &r->Z, hh );
     238     2376864 :   return r;
     239     2376864 : }
     240             : 
     241             : /* fd_bn254_g1_scalar_mul computes r = [s]P.
     242             :    p must be in affine form (p->Z == 1).
     243             :    The result is in projective coordinates. */
     244             : fd_bn254_g1_t *
     245             : fd_bn254_g1_scalar_mul( fd_bn254_g1_t *           r,
     246             :                         fd_bn254_g1_t const *     p,
     247       30126 :                         fd_bn254_scalar_t const * s ) {
     248       30126 :   if( FD_UNLIKELY( fd_uint256_is_zero( s ) || fd_bn254_g1_is_zero( p ) ) ) {
     249           3 :     return fd_bn254_g1_set_zero( r );
     250           3 :   }
     251       30123 :   const ulong g1_const[ 3 ] = { 0x5398fd0300ff6565UL, 0x4ccef014a773d2d2UL, 0x0000000000000002UL };
     252       30123 :   ulong b1[ 3 ];
     253       30123 :   ulong b2[ 2 ];
     254       30123 :   fd_bn254_glv_sxg3( b1, s, g1_const );
     255       30123 :   fd_bn254_glv_sxg2( b2, s, g2_const );
     256             : 
     257             :   /* k1 = s - b1*N_A - b2*N_B (always non-negative for G1) */
     258       30123 :   fd_uint256_t k1[1];
     259       30123 :   {
     260       30123 :     ulong p11[ 4 ];
     261             :     /* b2*nb will produce at most 3 limbs, but we want the 4th zeroed for the addition. */
     262       30123 :     ulong p21[ 4 ] = {0};
     263       30123 :     ulong  t[ 4 ];
     264       30123 :     fd_bn254_glv_mul3x2( p11, b1, na );
     265       30123 :     fd_bn254_glv_mul2x1( p21, b2, nb );
     266       30123 :     fd_bn254_glv_add4( t, p11, p21 );
     267       30123 :     fd_bn254_glv_sub4( k1->limbs, s->limbs, t );
     268       30123 :   }
     269             : 
     270             :   /* k2 = b1*N_B - b2*N_C (may be negative) */
     271       30123 :   fd_uint256_t k2_abs[1];
     272       30123 :   int k2_neg = 0;
     273       30123 :   {
     274       30123 :     ulong pos[ 4 ], neg[ 4 ];
     275       30123 :     fd_bn254_glv_mul3x1( pos, b1, nb );
     276       30123 :     fd_bn254_glv_mul2x2( neg, b2, nc );
     277       30123 :     ulong borrow = fd_bn254_glv_sub4( k2_abs->limbs, pos, neg );
     278       30123 :     if( borrow ) {
     279           0 :       k2_neg = 1;
     280           0 :       fd_bn254_glv_negate4( k2_abs->limbs );
     281           0 :     }
     282       30123 :   }
     283             : 
     284             :   /* pt2 = phi(P) = (beta * P.x, P.y). If k2 < 0, negate pt2. */
     285       30123 :   fd_bn254_g1_t pt2[1];
     286       30123 :   fd_bn254_fp_mul    ( &pt2->X, &p->X, fd_bn254_const_beta_mont );
     287       30123 :   fd_bn254_fp_set    ( &pt2->Y, &p->Y );
     288       30123 :   fd_bn254_fp_set_one( &pt2->Z );
     289       30123 :   if( k2_neg ) {
     290           0 :     fd_bn254_fp_neg( &pt2->Y, &pt2->Y );
     291           0 :   }
     292             : 
     293       30123 :   fd_bn254_g1_t pt12[1];
     294       30123 :   fd_bn254_g1_affine_add( pt12, p, pt2 );
     295             : 
     296             :   /* Shamir's trick: simultaneous double-and-add on k1, k2. */
     297       30123 :   int i = 255;
     298     3921069 :   for( ; i>=0; i-- ) {
     299     3921069 :     int k1b = !!fd_uint256_bit( k1, i );
     300     3921069 :     int k2b = !!fd_uint256_bit( k2_abs, i );
     301     3921069 :     if( k1b || k2b ) {
     302       30123 :       fd_bn254_g1_set( r, ( k1b && k2b ) ? pt12 : ( k1b ? p : pt2 ) );
     303       30123 :       break;
     304       30123 :     }
     305     3921069 :   }
     306       30123 :   if( FD_UNLIKELY( i<0 ) ) {
     307           0 :     return fd_bn254_g1_set_zero( r );
     308           0 :   }
     309     3820542 :   for( i--; i >= 0; i-- ) {
     310     3790419 :     fd_bn254_g1_dbl( r, r );
     311     3790419 :     int k1b = !!fd_uint256_bit( k1, i );
     312     3790419 :     int k2b = !!fd_uint256_bit( k2_abs, i );
     313     3790419 :     if( k1b && k2b ) {
     314     1262370 :       fd_bn254_g1_add_mixed( r, r, pt12 );
     315     2528049 :     } else if( k1b ) {
     316      452739 :       fd_bn254_g1_add_mixed( r, r, p );
     317     2075310 :     } else if( k2b ) {
     318      661755 :       fd_bn254_g1_add_mixed( r, r, pt2 );
     319      661755 :     }
     320     3790419 :   }
     321             : 
     322       30123 :   return r;
     323       30123 : }
     324             : 
     325             : /* fd_bn254_g1_frombytes_internal extracts (x, y) and performs basic checks.
     326             :    This is used by fd_bn254_g1_compress() and fd_bn254_g1_frombytes_check_subgroup().
     327             :    https://github.com/arkworks-rs/algebra/blob/v0.4.2/ec/src/models/short_weierstrass/mod.rs#L173-L178 */
     328             : static inline fd_bn254_g1_t *
     329             : fd_bn254_g1_frombytes_internal( fd_bn254_g1_t * p,
     330             :                                 uchar const     in[64],
     331      121215 :                                 int             big_endian ) {
     332             :   /* Special case: all zeros => point at infinity */
     333      121215 :   const uchar zero[64] = { 0 };
     334      121215 :   if( FD_UNLIKELY( fd_memeq( in, zero, 64 ) ) ) {
     335          48 :     return fd_bn254_g1_set_zero( p );
     336          48 :   }
     337             : 
     338             :   /* Check x < p */
     339      121167 :   if( FD_UNLIKELY( !fd_bn254_fp_frombytes_nm( &p->X, &in[0], big_endian, NULL, NULL ) ) ) {
     340           0 :     return NULL;
     341           0 :   }
     342             : 
     343             :   /* Check flags and y < p */
     344      121167 :   int is_inf, is_neg;
     345      121167 :   if( FD_UNLIKELY( !fd_bn254_fp_frombytes_nm( &p->Y, &in[32], big_endian, &is_inf, &is_neg ) ) ) {
     346           0 :     return NULL;
     347           0 :   }
     348             : 
     349      121167 :   if( FD_UNLIKELY( is_inf ) ) {
     350           3 :     return fd_bn254_g1_set_zero( p );
     351           3 :   }
     352             : 
     353      121164 :   fd_bn254_fp_set_one( &p->Z );
     354      121164 :   return p;
     355      121167 : }
     356             : 
     357             : /* fd_bn254_g1_frombytes_check_subgroup performs frombytes AND checks subgroup membership. */
     358             : static inline fd_bn254_g1_t *
     359             : fd_bn254_g1_frombytes_check_subgroup( fd_bn254_g1_t * p,
     360             :                                       uchar const     in[64],
     361       91116 :                                       int             big_endian ) {
     362       91116 :   if( FD_UNLIKELY( !fd_bn254_g1_frombytes_internal( p, in, big_endian ) ) ) {
     363           0 :     return NULL;
     364           0 :   }
     365       91116 :   if( FD_UNLIKELY( fd_bn254_g1_is_zero( p ) ) ) {
     366          45 :     return p;
     367          45 :   }
     368             : 
     369       91071 :   fd_bn254_fp_to_mont( &p->X, &p->X );
     370       91071 :   fd_bn254_fp_to_mont( &p->Y, &p->Y );
     371       91071 :   fd_bn254_fp_set_one( &p->Z );
     372             : 
     373             :   /* Check that y^2 = x^3 + b */
     374       91071 :   fd_bn254_fp_t y2[1], x3b[1];
     375       91071 :   fd_bn254_fp_sqr( y2, &p->Y );
     376       91071 :   fd_bn254_fp_sqr( x3b, &p->X );
     377       91071 :   fd_bn254_fp_mul( x3b, x3b, &p->X );
     378       91071 :   fd_bn254_fp_add( x3b, x3b, fd_bn254_const_b_mont );
     379       91071 :   if( FD_UNLIKELY( !fd_bn254_fp_eq( y2, x3b ) ) ) {
     380           0 :     return NULL;
     381           0 :   }
     382             : 
     383             :   /* G1 has prime order, so we don't need to do any further checks. */
     384             : 
     385       91071 :   return p;
     386       91071 : }