<html><head><meta http-equiv="Content-Type" content="text/html charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class="">Hello Curves,<div class=""><br class=""></div><div class="">I was porting an embedded EC library to work with other curves, and I was wondering if anyone knew a compact, fast, deterministic algorithm for square roots mod NIST p224 (or other annoying primes)? I’m not worried about constant-time operation, since I probably have to blind anyway to reduce other side channels, but a deterministic algorithm is preferable because of possible real-time constraints.</div><div class=""><br class=""></div><div class="">OpenSSL’s implementation takes something like 4700 multiplications on average, which is more than a point*scalar multiply.</div><div class=""><br class=""></div><div class="">Dan’s paper [1] gives a variant of the Tonelli-Shanks algorithm which is fast and deterministic, requiring only 364 multiplies. But it’s very heavy, with 1024 precomputed elements. Reducing the table size makes the algorithm considerably slower.</div><div class=""><br class=""></div><div class="">I’ve tried to combine Dan’s Tonelli-Shanks variant with Sutherland’s algorithm, below. A variant with 4 precomputed points, a 64-element discrete log table (not as big as it sounds) and 10 or so field elements in memory takes 900 multiplies worst-case (~760 average-case).</div><div class=""><br class=""></div><div class="">Cipolla’s algorithm is nondeterministic.</div><div class=""><br class=""></div><div class="">Is Pocklington the way to go? It’s technically nondeterministic as well, but if I read [1] correctly the failure probability is 2^-96. A straightforward implementation probably takes more than 900 multiplies, but hopefully it’s simpler and uses less memory, right?</div><div class=""><br class=""></div><div class="">Is there some other algorithm I’m missing?</div><div class=""><br class=""></div><div class="">Thanks,</div><div class="">— Mike</div><div class=""><br class=""></div><div class="">[1] <a href="http://cr.yp.to/papers/sqroot.pdf" class="">http://cr.yp.to/papers/sqroot.pdf</a></div><div class=""><br class=""></div><div class="">####################################################</div><div class=""><br class=""></div><div class=""><div class="">p = 2^224 - 2^96 + 1</div><div class="">F = GF(p)</div><div class="">qnr = F(11)^(2^128-1)</div><div class=""><br class=""></div><div class=""># 3-element precomputed table with roots of unity</div><div class="">precmp = dict((</div><div class=""> (3*2^i, qnr^2^(96-3*2^i))</div><div class=""> for i in [2,3,4]</div><div class="">))</div><div class=""><br class=""></div><div class=""># 64-element discrete log table</div><div class="">tbits = 6</div><div class="">tlogger = qnr^2^(96-tbits)</div><div class="">dlog_table = dict((</div><div class=""> (tlogger^i,2^tbits-i) for i in xrange(2^tbits)</div><div class="">))</div><div class=""><br class=""></div><div class="">def isqrt(x):</div><div class=""> # inverse square root of x, assuming x is a QR</div><div class=""> # First, DJB’s addition chain</div><div class=""> s = x^3</div><div class=""> s = s^2 * x # 2^3 - 1</div><div class=""> u = s^2^3 * s # 2^6 - 1</div><div class=""> s = u^2^6 * u # 2^12 - 1</div><div class=""> t = s^2^12 * s # 2^24 - 1</div><div class=""> s = t^2^24 * t # 2^48 - 1</div><div class=""> s = s^2^48 * s # 2^96 - 1</div><div class=""> s = s^2^24 * t # 2^120 - 1</div><div class=""> s = s^2^6 * u # 2^126 - 1</div><div class=""> s = s^2 * x</div><div class=""> # here s == x^(2^127-1)</div><div class=""> </div><div class=""> ret = [s,s^2*x] # working powers of return value</div><div class=""> rou = [qnr,qnr^2] # working roots of unity</div> stack = [95]<div class=""> </div><div class=""> def leaf(n):</div><div class=""> dlog = dlog_table[ret.pop()]>>(tbits-n)</div><div class=""> for i in xrange(len(rou)):</div><div class=""> for j in xrange(n):</div><div class=""> if dlog & (1<<j): ret[i] = ret[i] * rou[i]</div><div class=""> rou[i] = rou[i]^2</div><div class=""> </div><div class=""> # max stack depth is 4, so 5+5 field elements in ret+rou</div><div class=""> while len(stack):</div><div class=""> n = stack.pop()</div><div class=""> ret.append(ret[-1]^2^floor(n/2))</div><div class=""> </div><div class=""> if n <= 2*tbits:</div><div class=""> leaf(ceil(n/2))</div><div class=""> rou.pop()</div><div class=""> leaf(floor(n/2))</div><div class=""> else:</div><div class=""> rou.append(precmp[ceil(n/2)])</div><div class=""> stack.append(floor(n/2))</div><div class=""> stack.append(ceil(n/2))</div><div class=""> </div><div class=""> return ret[0]</div></div></body></html>