fixed

Pure Haskell large fixed-width integers and Montgomery arithmetic (docs.ppad.tech/fixed).
git clone git://git.ppad.tech/fixed.git
Log | Files | Refs | README | LICENSE

generate_inv.sh (4983B)


      1 #!/usr/bin/env bash
      2 
      3 # generates a constant-time haskell function for performing modular
      4 # inversion with montgomery arithmetic on a secp256k1-derived field.
      5 #
      6 # fermat inversion is used: a^-1 = a ^ (m - 2) mod m, where m is the
      7 # secp256k1 field prime (curve) or the scalar group order (scalar).
      8 #
      9 # rather than a naive square-and-multiply over the exponent bits (one
     10 # montgomery-multiply per set bit), we use an addition chain that
     11 # coalesces runs of 1 bits. for a maximal run of L consecutive 1 bits
     12 # we precompute the block
     13 #
     14 #   x_L = a ^ (2^L - 1)
     15 #
     16 # and fold the whole run in with a single multiply, since
     17 #
     18 #   result^(2^L) * x_L
     19 #
     20 # appends L ones to the accumulated exponent. the x_L blocks are built
     21 # with a doubling ladder using the identities
     22 #
     23 #   x_(i+j) = (x_i)^(2^j) * x_j         (so x_2k = (x_k)^(2^k) * x_k)
     24 #
     25 # which costs O(log L) multiplies per block instead of O(L).
     26 #
     27 # the square-and-multiply schedule is still fixed and data-independent,
     28 # so given constant-time 'sqr#' and 'mul#', 'inv#' is constant-time by
     29 # construction.
     30 
     31 # secp256k1 field prime - 2
     32 exponent_curve="1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111011111111111111111111110000101101"
     33 
     34 # secp256k1 scalar group order - 2
     35 exponent_scalar="1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111010111010101011101101110011100110101011110100100010100000001110111011111111010010010111101000110011010000001101100100000100111111"
     36 
     37 # adjust me as desired; use 'curve' or 'scalar'
     38 target="curve"
     39 
     40 declare exponent
     41 if [[ $target == "curve" ]]; then
     42   exponent=$exponent_curve
     43 elif [[ $target == "scalar" ]]; then
     44   exponent=$exponent_scalar
     45 else
     46   echo "unknown target: $target" >&2
     47   exit 1
     48 fi
     49 
     50 # next free t-label, and a map from block length L to the variable
     51 # holding x_L = a ^ (2^L - 1). x_1 is the argument 'a' itself.
     52 idx=1
     53 declare -A blk
     54 blk[1]="a"
     55 
     56 # REPLY conventions: emit_* helpers return the resulting variable in
     57 # the global REPLY.
     58 
     59 emit_sqr_n() { # $1 = source var, $2 = count
     60   local src=$1 n=$2 i
     61   for (( i = 0; i < n; i++ )); do
     62     echo "      !t$idx = sqr# $src"
     63     src=t$idx
     64     idx=$((idx + 1))
     65   done
     66   REPLY=$src
     67 }
     68 
     69 emit_mul() { # $1, $2 = vars
     70   echo "      !t$idx = mul# $1 $2"
     71   REPLY=t$idx
     72   idx=$((idx + 1))
     73 }
     74 
     75 # ensure x_1, x_2, x_4, ..., x_hp (powers of two up to hp) are built
     76 build_powers() { # $1 = highest power of two needed
     77   local hp=$1 p=1 np acc
     78   while (( p < hp )); do
     79     np=$((p * 2))
     80     if [[ -z ${blk[$np]:-} ]]; then
     81       emit_sqr_n "${blk[$p]}" "$p"; acc=$REPLY # x_p ^ (2^p)
     82       emit_mul "$acc" "${blk[$p]}"; acc=$REPLY # * x_p  = x_2p
     83       blk[$np]=$acc
     84     fi
     85     p=$np
     86   done
     87 }
     88 
     89 # ensure x_L is built, composing it from power-of-two blocks
     90 build_block() { # $1 = run length L
     91   local L=$1 hp pw rem acc
     92   [[ -n ${blk[$L]:-} ]] && return
     93   hp=1
     94   while (( hp * 2 <= L )); do hp=$((hp * 2)); done
     95   build_powers "$hp"
     96   rem=$L; acc=""; pw=$hp
     97   while (( pw >= 1 )); do
     98     if (( rem >= pw )); then
     99       if [[ -z $acc ]]; then
    100         acc=${blk[$pw]}
    101       else
    102         emit_sqr_n "$acc" "$pw"; acc=$REPLY
    103         emit_mul "$acc" "${blk[$pw]}"; acc=$REPLY
    104       fi
    105       rem=$((rem - pw))
    106     fi
    107     pw=$((pw / 2))
    108   done
    109   blk[$L]=$acc
    110 }
    111 
    112 # parse the exponent into alternating runs ('R', ones) and gaps ('G')
    113 seg_type=(); seg_len=()
    114 i=0; n=${#exponent}
    115 while (( i < n )); do
    116   c=${exponent:i:1}
    117   j=$i
    118   while (( j < n )) && [[ ${exponent:j:1} == "$c" ]]; do j=$((j + 1)); done
    119   if [[ $c == 1 ]]; then seg_type+=("R"); else seg_type+=("G"); fi
    120   seg_len+=("$((j - i))")
    121   i=$j
    122 done
    123 
    124 echo "-- generated by etc/generate_inv.sh"
    125 echo "inv#"
    126 echo "  :: Limb4"
    127 echo "  -> Limb4"
    128 echo "inv# a ="
    129 echo "  let"
    130 
    131 # build blocks for the longest run first so the doubling ladder is
    132 # shared, then any shorter runs reuse the cached intermediates.
    133 maxrun=0
    134 for k in "${!seg_type[@]}"; do
    135   if [[ ${seg_type[k]} == R ]] && (( seg_len[k] > maxrun )); then
    136     maxrun=${seg_len[k]}
    137   fi
    138 done
    139 build_block "$maxrun"
    140 for k in "${!seg_type[@]}"; do
    141   if [[ ${seg_type[k]} == R ]]; then build_block "${seg_len[k]}"; fi
    142 done
    143 
    144 # fold the exponent: result starts at the leading run, then each
    145 # subsequent run is appended via result^(2^(gap+run)) * x_run.
    146 result=""; pending=0; first=1
    147 for k in "${!seg_type[@]}"; do
    148   if [[ ${seg_type[k]} == G ]]; then
    149     pending=$((pending + seg_len[k]))
    150   else
    151     L=${seg_len[k]}
    152     if (( first )); then
    153       result=${blk[$L]}; first=0
    154     else
    155       emit_sqr_n "$result" "$((pending + L))"; result=$REPLY
    156       emit_mul "$result" "${blk[$L]}"; result=$REPLY
    157     fi
    158     pending=0
    159   fi
    160 done
    161 # trailing zero bits, if any
    162 if (( pending > 0 )); then
    163   emit_sqr_n "$result" "$pending"; result=$REPLY
    164 fi
    165 
    166 echo "      !r = $result"
    167 echo "  in  r"
    168 echo '{-# INLINE inv# #-}'