derive_frobenius.sage
1 #!/usr/bin/sage 2 # vim: syntax=python 3 # vim: set ts=2 sw=2 et: 4 5 # Constantine 6 # Copyright (c) 2018-2019 Status Research & Development GmbH 7 # Copyright (c) 2020-Present Mamy André-Ratsimbazafy 8 # Licensed and distributed under either of 9 # * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT). 10 # * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0). 11 # at your option. This file may not be copied, modified, or distributed except according to those terms. 12 13 # ############################################################ 14 # 15 # Frobenius constants 16 # 17 # ############################################################ 18 19 # Imports 20 # --------------------------------------------------------- 21 22 import os 23 import inspect, textwrap 24 25 # Working directory 26 # --------------------------------------------------------- 27 28 os.chdir(os.path.dirname(__file__)) 29 30 # Sage imports 31 # --------------------------------------------------------- 32 # Accelerate arithmetic by accepting probabilistic proofs 33 from sage.structure.proof.all import arithmetic 34 arithmetic(False) 35 36 load('curves.sage') 37 38 # Utilities 39 # --------------------------------------------------------- 40 41 def fp2_to_hex(a): 42 v = vector(a) 43 return '0x' + Integer(v[0]).hex() + ' + β * ' + '0x' + Integer(v[1]).hex() 44 45 def field_to_nim(value, field, curve, prefix = "", comment_above = "", comment_right = ""): 46 result = '# ' + comment_above + '\n' if comment_above else '' 47 comment_right = ' # ' + comment_right if comment_right else '' 48 49 if field == 'Fp2': 50 v = vector(value) 51 52 result += inspect.cleandoc(f""" 53 {prefix}Fp2[{curve}].fromHex( {comment_right} 54 "0x{Integer(v[0]).hex()}", 55 "0x{Integer(v[1]).hex()}" 56 )""") 57 elif field == 'Fp': 58 result += inspect.cleandoc(f""" 59 {prefix}Fp[{curve}].fromHex( {comment_right} 60 "0x{Integer(value).hex()}") 61 """) 62 else: 63 raise NotImplementedError() 64 65 return result 66 67 # Code generators 68 # --------------------------------------------------------- 69 70 def genFrobeniusMapConstants(curve_name, curve_config): 71 embdeg = curve_config[curve_name]['tower']['embedding_degree'] 72 twdeg = curve_config[curve_name]['tower']['twist_degree'] 73 g2field = f'Fp{embdeg//twdeg}' if (embdeg//twdeg) > 1 else 'Fp' 74 75 p = curve_config[curve_name]['field']['modulus'] 76 Fp = GF(p) 77 K.<u> = PolynomialRing(Fp) 78 if g2field == 'Fp2': 79 QNR_Fp = curve_config[curve_name]['tower']['QNR_Fp'] 80 Fp2.<beta> = Fp.extension(u^2 - QNR_Fp) 81 else: 82 SNR_Fp = curve_config[curve_name]['tower']['SNR_Fp'] 83 Fp2.<beta> = Fp.extension(u^2 - SNR_Fp) 84 85 if g2field == 'Fp2': 86 SNR = curve_config[curve_name]['tower']['SNR_Fp2'] 87 SNR = Fp2(SNR) 88 else: 89 # To build the Fp6 extension, since we use a SexticNonResidue 90 # to build Fp2, we can reuse it as a cubic non-residue 91 # It always has [0, 1] coordinates in Fp2 92 SNR = Fp2([0, 1]) 93 94 halfK = embdeg//2 95 96 print('\n----> Frobenius extension field constants <----\n') 97 buf = inspect.cleandoc(f""" 98 # Frobenius map - on extension fields 99 # ----------------------------------------------------------------- 100 101 # We start from base frobenius constant for a {embdeg} embedding degree. 102 # with 103 # - a sextic twist, SNR being the Sextic Non-Residue. 104 # - coef being the Frobenius coefficient "ID" 105 # c = SNR^((p-1)/{halfK})^coef 106 # 107 # On Fp2 frobenius(c) = conj(c) so we have 108 # For n=2, with n the number of Frobenius applications 109 # c2 = c * (c^p) = c * frobenius(c) = c * conj(c) 110 # c2 = (SNR * conj(SNR))^((p-1)/{halfK})^coef) 111 # c2 = (norm(SNR))^((p-1)/{halfK})^coef) 112 # For k=3 113 # c3 = c * c2^p = c * frobenius(c2) = c * conj(c2) 114 # with conj(norm(SNR)) = norm(SNR) as a norm is strictly on the base field. 115 # c3 = (SNR * norm(SNR))^((p-1)/{halfK})^coef) 116 # 117 # A more generic formula can be derived by observing that 118 # c3 = c * c2^p = c * (c * c^p)^p 119 # c3 = c * c^p * c^p² 120 # with 4, we have 121 # c4 = c * c3^p = c * (c * c^p * c^p²)^p 122 # c4 = c * c^p * c^p² * c^p³ 123 # with n we have 124 # cn = c * c^p * c^p² ... * c^p^(n-1) 125 # cn = c^(1+p+p² + ... + p^(n-1)) 126 # This is the sum of first n terms of a geometric series 127 # hence cn = c^((p^n-1)/(p-1)) 128 # We now expand c 129 # cn = SNR^((p-1)/{halfK})^coef^((p^n-1)/(p-1)) 130 # cn = SNR^((p^n-1)/{halfK})^coef 131 # cn = SNR^(coef * (p^n-1)/{halfK}) 132 133 const {curve_name}_FrobeniusMapCoefficients* = [ 134 """) 135 136 arr = "" 137 maxN = 3 # We only need up to f^(p^3) in final exponentiation 138 139 for n in range(1, maxN + 1): 140 for coef in range(halfK): 141 if coef == 0: 142 arr += f'\n# frobenius({n}) -----------------------\n' 143 arr += '[' 144 frobmapcoef = SNR^(coef*((p^n-1)/halfK)) 145 hatN = '^' + str(n) if n>1 else '' 146 arr += field_to_nim(frobmapcoef, 'Fp2', curve_name, comment_right = f'SNR^((p{hatN}-1)/{halfK})^{coef}') 147 if coef != halfK - 1: 148 arr += ',\n' 149 arr += '],\n' 150 151 buf += textwrap.indent(arr, ' ') 152 buf += ']' 153 return buf 154 155 def genFrobeniusPsiConstants(curve_name, curve_config): 156 embdeg = curve_config[curve_name]['tower']['embedding_degree'] 157 twdeg = curve_config[curve_name]['tower']['twist_degree'] 158 twkind = curve_config[curve_name]['tower']['twist'] 159 g2field = f'Fp{embdeg//twdeg}' if (embdeg//twdeg) > 1 else 'Fp' 160 161 p = curve_config[curve_name]['field']['modulus'] 162 Fp = GF(p) 163 K.<u> = PolynomialRing(Fp) 164 if g2field == 'Fp2': 165 QNR_Fp = curve_config[curve_name]['tower']['QNR_Fp'] 166 Fp2.<beta> = Fp.extension(u^2 - QNR_Fp) 167 168 if g2field == 'Fp2': 169 SNR = curve_config[curve_name]['tower']['SNR_Fp2'] 170 SNR = Fp2(SNR) 171 else: 172 SNR = curve_config[curve_name]['tower']['SNR_Fp'] 173 SNR = Fp(SNR) 174 175 print('\n----> ψ (Psi) - Untwist-Frobenius-Twist Endomorphism constants <----\n') 176 buf = inspect.cleandoc(f""" 177 # ψ (Psi) - Untwist-Frobenius-Twist Endomorphisms on twisted curves 178 # ----------------------------------------------------------------- 179 """) 180 buf += '\n' 181 if twkind == 'D_Twist': 182 buf += f'# {curve_name} is a D-Twist: psi1_coef1 = SNR^((p-1)/{twdeg})\n\n' 183 xi = SNR 184 snrUsed = 'SNR' 185 else: 186 buf += f'# {curve_name} is a M-Twist: psi1_coef1 = (1/SNR)^((p-1)/{twdeg})\n\n' 187 xi = 1/SNR 188 snrUsed = '(1/SNR)' 189 190 maxPsi = CyclotomicField(embdeg).degree() 191 192 for n in range(1, maxPsi+1): 193 for coef in range(2, 3+1): 194 # Same formula as FrobeniusMap constants 195 # except that 196 # - we only need 2 coefs for elliptic curve twists 197 # - xi = SNR or 1/SNR depending on D-Twist or M-Twist respectively 198 # - the divisor is the twist degree isntead of half the embedding degree 199 frobpsicoef = xi^(coef*(p^n - 1)/twdeg) 200 hatN = '^' + str(n) if n>1 else '' 201 buf += field_to_nim( 202 frobpsicoef, g2field, curve_name, 203 prefix = f'const {curve_name}_FrobeniusPsi_psi{n}_coef{coef}* = ', 204 comment_above = f'{snrUsed}^({coef}(p{hatN}-1)/{twdeg})' 205 ) + '\n' 206 207 buf += '\n' 208 209 buf += inspect.cleandoc(f""" 210 # For a sextic twist 211 # - p ≡ 1 (mod 2) 212 # - p ≡ 1 (mod 3) 213 # 214 # psi2_coef3 is always -1 (mod p^m) with m = embdeg/twdeg 215 # Recap, with ξ (xi) the sextic non-residue for D-Twist or 1/SNR for M-Twist 216 # psi_2 ≡ ξ^((p-1)/6)^2 ≡ ξ^((p-1)/3) 217 # psi_3 ≡ psi_2 * ξ^((p-1)/6) ≡ ξ^((p-1)/3) * ξ^((p-1)/6) ≡ ξ^((p-1)/2) 218 # 219 # In Fp² (i.e. embedding degree of 12, G2 on Fp2) 220 # - quadratic non-residues respect the equation a^((p²-1)/2) ≡ -1 (mod p²) by the Legendre symbol 221 # - sextic non-residues are also quadratic non-residues so ξ^((p²-1)/2) ≡ -1 (mod p²) 222 # - QRT(1/a) = QRT(a) with QRT the quadratic residuosity test 223 # 224 # We have psi2_3 ≡ psi_3 * psi_3^p ≡ psi_3^(p+1) 225 # ≡ (ξ^(p-1)/2)^(p+1) (mod p²) 226 # ≡ ξ^((p-1)(p+1)/2) (mod p²) 227 # ≡ ξ^((p²-1)/2) (mod p²) 228 # And ξ^((p²-1)/2) ≡ -1 (mod p²) since ξ is a quadratic non-residue 229 # So psi2_3 ≡ -1 (mod p²) 230 # 231 # 232 # In Fp (i.e. embedding degree of 6, G2 on Fp) 233 # - Fermat's Little Theorem gives us a^(p-1) ≡ 1 (mod p) 234 # 235 # psi2_3 ≡ ξ^((p-1)(p+1)/2) (mod p) 236 # ≡ ξ^((p+1)/2)^(p-1) (mod p) as we have 2|p+1 237 # ≡ 1 (mod p) by Fermat's Little Theorem 238 """) 239 240 return buf 241 242 # CLI 243 # --------------------------------------------------------- 244 245 if __name__ == "__main__": 246 # Usage 247 # BLS12-381 248 # sage sage/derive_frobenius.sage BLS12_381 249 250 from argparse import ArgumentParser 251 252 parser = ArgumentParser() 253 parser.add_argument("curve",nargs="+") 254 args = parser.parse_args() 255 256 curve = args.curve[0] 257 258 if curve not in Curves: 259 raise ValueError( 260 curve + 261 ' is not one of the available curves: ' + 262 str(Curves.keys()) 263 ) 264 else: 265 trace = Curves[curve]['field']['trace'] 266 print(f'trace of Frobenius ({int(trace).bit_length()}-bit): 0x{Integer(trace).hex()}') 267 268 FrobMap = genFrobeniusMapConstants(curve, Curves) 269 FrobPsi = genFrobeniusPsiConstants(curve, Curves) 270 271 with open(f'{curve.lower()}_frobenius.nim', 'w') as f: 272 f.write(copyright()) 273 f.write('\n\n') 274 275 embdeg = Curves[curve]['tower']['embedding_degree'] 276 twdeg = Curves[curve]['tower']['twist_degree'] 277 278 if embdeg//twdeg >= 2: 279 f.write(inspect.cleandoc(""" 280 import 281 ../config/curves, 282 ../extension_fields, 283 ../io/io_extfields 284 """)) 285 else: 286 f.write(inspect.cleandoc(""" 287 import 288 ../config/curves, 289 ../extension_fields, 290 ../io/[io_fields, io_extfields] 291 """)) 292 f.write('\n\n') 293 f.write(FrobMap) 294 f.write('\n\n') 295 f.write(FrobPsi) 296 297 print(f'Successfully created {curve}_frobenius.nim') 298 299 print(inspect.cleandoc("""\n 300 For testing you can verify the following invariants: 301 302 Galbraith-Lin-Scott, 2008, Theorem 1 303 Fuentes-Castaneda et al, 2011, Equation (2) 304 ψ(ψ(P)) - t*ψ(P) + p*P == Infinity 305 306 Galbraith-Scott, 2008, Lemma 1 307 The cyclotomic polynomial with GΦ(ψ(P)) == Infinity 308 Hence for embedding degree k=12 309 ψ⁴(P) - ψ²(P) + P == Infinity 310 for embedding degree k=6 311 ψ²(P) - ψ(P) + P == Infinity 312 """))