From de377290ee94d7e5089c4d245b90a979cb402b81 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 4 Jan 2018 13:43:43 -0800 Subject: [PATCH 1/8] Implement batch inversion using a product tree. --- src/field.rs | 93 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) diff --git a/src/field.rs b/src/field.rs index 4532820..23362e8 100644 --- a/src/field.rs +++ b/src/field.rs @@ -178,6 +178,63 @@ impl FieldElement { (t19, t3) } + /// Given a slice of public `FieldElements`, replace each with its inverse. + /// + /// All input `FieldElements` **MUST** be nonzero. + /// + /// This function is most efficient when the batch size (slice + /// length) is a power of 2. + pub fn batch_invert(inputs: &mut [FieldElement]) { + // First, compute the product of all inputs using a product + // tree: + // + // Inputs: [x_0, x_1, x_2] + // + // Tree: + // + // x_0*x_1*x_2*1 tree[1] + // / \ + // x_0*x_1 x_2*1 tree[2,3] + // / \ / \ + // x_0 x_1 x_2 1 tree[4,5,6,7] + // + // The leaves of the tree are the inputs. We store the tree in + // an array of length 2*n, similar to a binary heap. + // + // To initialize the tree, set every node to 1, then fill in + // the leaf nodes with the input variables. Finally, set every + // non-leaf node to be the product of its children. + + let n = inputs.len().next_power_of_two(); + let mut tree = vec![FieldElement::one(); 2*n]; + tree[n..n+inputs.len()].copy_from_slice(inputs); + for i in (1..n).rev() { + tree[i] = &tree[2*i] * &tree[2*i+1]; + } + + // The root of the tree is the product of all inputs, and is + // stored at index 1. Compute its inverse. + let allinv = tree[1].invert(); + + // To compute y_i = 1/x_i, start at the i-th leaf node of the + // tree, and walk up to the root of the tree, multiplying + // `allinv` by each sibling. This computes + // + // y_i = y * (all x_j except x_i) + // + // using lg(n) multiplications for each y_i, taking n*lg(n) in + // total. + for i in 0..inputs.len() { + let mut inv = allinv; + let mut node = n + i; + while node > 1 { + inv *= &tree[node ^ 1]; + node = node >> 1; + } + inputs[i] = inv; + } + } + /// Given a nonzero field element, compute its inverse. /// /// The inverse is computed as self^(p-2), since @@ -375,6 +432,21 @@ mod test { assert_eq!(FieldElement::one(), &a * &should_be_inverse); } + #[test] + fn batch_invert_a_matches_nonbatched() { + let a = FieldElement::from_bytes(&A_BYTES); + let ap58 = FieldElement::from_bytes(&AP58_BYTES); + let asq = FieldElement::from_bytes(&ASQ_BYTES); + let ainv = FieldElement::from_bytes(&AINV_BYTES); + let a2 = &a + &a; + let a_list = vec![a, ap58, asq, ainv, a2]; + let mut ainv_list = a_list.clone(); + FieldElement::batch_invert(&mut ainv_list[..]); + for i in 0..5 { + assert_eq!(a_list[i].invert(), ainv_list[i]); + } + } + #[test] fn a_p58_vs_ap58_constant() { let a = FieldElement::from_bytes(&A_BYTES); @@ -470,4 +542,25 @@ mod bench { let a = FieldElement::from_bytes(&A_BYTES); b.iter(|| a.invert()); } + + #[bench] + fn batch_16_inv(b: &mut Bencher) { + let a = FieldElement::from_bytes(&A_BYTES); + let mut a_vec = vec![a; 16]; + b.iter(|| FieldElement::batch_invert(&mut a_vec)); + } + + #[bench] + fn batch_128_inv(b: &mut Bencher) { + let a = FieldElement::from_bytes(&A_BYTES); + let mut a_vec = vec![a; 128]; + b.iter(|| FieldElement::batch_invert(&mut a_vec)); + } + + #[bench] + fn batch_1024_inv(b: &mut Bencher) { + let a = FieldElement::from_bytes(&A_BYTES); + let mut a_vec = vec![a; 1024]; + b.iter(|| FieldElement::batch_invert(&mut a_vec)); + } } From 59837c6ecff02b77b9d5ff84dbc239d0cf33ef90 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Tue, 16 Jan 2018 14:30:15 -0800 Subject: [PATCH 2/8] Update ristretto.sage to match goldilocks b0af87 --- vendor/ristretto.sage | 390 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 318 insertions(+), 72 deletions(-) diff --git a/vendor/ristretto.sage b/vendor/ristretto.sage index f40bebd..04cf4f9 100644 --- a/vendor/ristretto.sage +++ b/vendor/ristretto.sage @@ -49,7 +49,11 @@ def isqrt(x,exn=InvalidEncodingException("Not on curve")): """Return 1/sqrt(x)""" if x==0: return 0 if not is_square(x): raise exn - return 1/sqrt(x) + s = sqrt(x) + #if negative(s): s=-s + return 1/s + +def inv0(x): return 1/x if x != 0 else 0 def isqrt_i(x): """Return 1/sqrt(x) or 1/sqrt(zeta * x)""" @@ -117,16 +121,20 @@ class QuotientEdwardsPoint(object): else: return self.__class__(-self.x, -self.y) + def doubleAndEncodeSpec(self): + return (self+self).encode() # Utility functions @classmethod - def bytesToGf(cls,bytes,mustBeProper=True,mustBePositive=False): + def bytesToGf(cls,bytes,mustBeProper=True,mustBePositive=False,maskHiBits=False): """Convert little-endian bytes to field element, sanity check length""" if len(bytes) != cls.encLen: raise InvalidEncodingException("wrong length %d" % len(bytes)) s = dec_le(bytes) - if mustBeProper and s >= cls.F.modulus(): + if mustBeProper and s >= cls.F.order(): raise InvalidEncodingException("%d out of range!" % s) + bitlen = int(ceil(log(cls.F.order())/log(2))) + if maskHiBits: s &= 2^bitlen-1 s = cls.F(s) if mustBePositive and negative(s): raise InvalidEncodingException("%d is negative!" % s) @@ -197,7 +205,42 @@ class RistrettoPoint(QuotientEdwardsPoint): if negative(isr^2*num*y*t): y = -y s = isr*y*(z-y) + return self.gfToBytes(s,mustBePositive=True) + + @optimized_version_of("doubleAndEncodeSpec") + def doubleAndEncode(self): + X,Y,Z,T = self.xyzt() + a,d,mneg = self.a,self.d,self.mneg + + if self.cofactor==8: + e = 2*X*Y + f = Z^2+d*T^2 + g = Y^2-a*X^2 + h = Z^2-d*T^2 + + inv1 = 1/(e*f*g*h) + z_inv = inv1*e*g # 1 / (f*h) + t_inv = inv1*f*h + if negative(e*g*z_inv): + if a==-1: sqrta = self.i + else: sqrta = -1 + e,f,g,h = g,h,-e,f*sqrta + factor = self.i + else: + factor = self.magic + + if negative(h*e*z_inv): g=-g + s = (h-g)*factor*g*t_inv + + else: + foo = Y^2+a*X^2 + bar = X*Y + den = 1/(foo*bar) + if negative(2*bar^2*den): tmp = a*X^2 + else: tmp = Y^2 + s = self.magic*(Z^2-tmp)*foo*den + return self.gfToBytes(s,mustBePositive=True) @classmethod @@ -238,8 +281,9 @@ class RistrettoPoint(QuotientEdwardsPoint): @classmethod def elligatorSpec(cls,r0): a,d = cls.a,cls.d - r = cls.qnr * cls.bytesToGf(r0)^2 + r = cls.qnr * cls.bytesToGf(r0,mustBeProper=False,maskHiBits=True)^2 den = (d*r-a)*(a*r-d) + if den == 0: return cls() n1 = cls.a*(r+1)*(a+d)*(d-a)/den n2 = r*n1 if is_square(n1): @@ -253,7 +297,7 @@ class RistrettoPoint(QuotientEdwardsPoint): @optimized_version_of("elligatorSpec") def elligator(cls,r0): a,d = cls.a,cls.d - r0 = cls.bytesToGf(r0) + r0 = cls.bytesToGf(r0,mustBeProper=False,maskHiBits=True) r = cls.qnr * r0^2 den = (d*r-a)*(a*r-d) num = cls.a*(r+1)*(a+d)*(d-a) @@ -278,15 +322,11 @@ class Decaf_1_1_Point(QuotientEdwardsPoint): if self.cofactor==8 and negative(x*y*self.isoMagic): x,y = self.torque() - - isr2 = isqrt(a*(y^2-1)) * sqrt(a*d-1) - + sr = xsqrt(1-a*x^2) - assert sr in [isr2*x*y,-isr2*x*y] - - altx = 1/isr2*self.isoMagic - if negative(altx): s = (1+x*y*isr2)/(a*x) - else: s = (1-x*y*isr2)/(a*x) + altx = x*y*self.isoMagic / sr + if negative(altx): s = (1+sr)/x + else: s = (1-sr)/x return self.gfToBytes(s,mustBePositive=True) @@ -297,52 +337,141 @@ class Decaf_1_1_Point(QuotientEdwardsPoint): s = cls.bytesToGf(s,mustBePositive=True) if s==0: return cls() - isr = isqrt(s^4 + 2*(a-2*d)*s^2 + 1) - altx = 2*s*isr*cls.isoMagic - if negative(altx): isr = -isr + t = xsqrt(s^4 + 2*(a-2*d)*s^2 + 1) + altx = 2*s*cls.isoMagic/t + if negative(altx): t = -t x = 2*s / (1+a*s^2) - y = (1-a*s^2) * isr + y = (1-a*s^2) / t if cls.cofactor==8 and (negative(x*y*cls.isoMagic) or y==0): raise InvalidEncodingException("x*y is invalid: %d, %d" % (x,y)) return cls(x,y) - @optimized_version_of("encodeSpec") - def encode(self): - """Encode, optimized version""" + def toJacobiQuartic(self,toggle_rotation=False,toggle_altx=False,toggle_s=False): + "Return s,t on jacobi curve" a,d = self.a,self.d x,y,z,t = self.xyzt() if self.cofactor == 8: # Cofactor 8 version + # Simulate IMAGINE_TWIST because that's how libdecaf does it + x = self.i*x + t = self.i*t + a = -a + d = -d + + # OK, the actual libdecaf code should be here num = (z+y)*(z-y) den = x*y - tmp = isqrt(num*(a-d)*den^2) - - if negative(tmp^2*den*num*(a-d)*t^2*self.isoMagic): - den,num = num,den - tmp *= sqrt(a-d) # witness that cofactor is 8 - yisr = x*sqrt(a) - toggle = (a==1) - else: - yisr = y*(a*d-1) - toggle = False + isr = isqrt(num*(a-d)*den^2) + + iden = isr * den * self.isoMagic # 1/sqrt((z+y)(z-y)) = 1/sqrt(1-Y^2) / z + inum = isr * num # sqrt(1-Y^2) * z / xysqrt(a-d) ~ 1/sqrt(1-ax^2)/z - tiisr = tmp*num - altx = tiisr*t*self.isoMagic - if negative(altx) != toggle: tiisr =- tiisr - s = tmp*den*yisr*(tiisr*z - 1) + if negative(iden*inum*self.i*t^2*(d-a)) != toggle_rotation: + iden,inum = inum,iden + fac = x*sqrt(a) + toggle=(a==-1) + else: + fac = y + toggle=False + + imi = self.isoMagic * self.i + altx = inum*t*imi + neg_altx = negative(altx) != toggle_altx + if neg_altx != toggle: inum =- inum + + tmp = fac*(inum*z + 1) + s = iden*tmp*imi + + negm1 = (negative(s) != toggle_s) != neg_altx + if negm1: m1 = a*fac + z + else: m1 = a*fac - z + + swap = toggle_s else: # Much simpler cofactor 4 version num = (x+t)*(x-t) isr = isqrt(num*(a-d)*x^2) - ratio = isr*num - if negative(ratio*self.isoMagic): ratio=-ratio - s = (a-d)*isr*x*(ratio*z - t) + ratio = isr*num + altx = ratio*self.isoMagic + + neg_altx = negative(altx) != toggle_altx + if neg_altx: ratio =- ratio + + tmp = ratio*z - t + s = (a-d)*isr*x*tmp + + negx = (negative(s) != toggle_s) != neg_altx + if negx: m1 = -a*t + x + else: m1 = -a*t - x + + swap = toggle_s + + if negative(s): s = -s - return self.gfToBytes(s,mustBePositive=True) + return s,m1,a*tmp,swap + + def invertElligator(self,toggle_r=False,*args,**kwargs): + "Produce preimage of self under elligator, or None" + a,d = self.a,self.d + + rets = [] + + tr = [False,True] if self.cofactor == 8 else [False] + for toggle_rotation in tr: + for toggle_altx in [False,True]: + for toggle_s in [False,True]: + for toggle_r in [False,True]: + s,m1,m12,swap = self.toJacobiQuartic(toggle_rotation,toggle_altx,toggle_s) + + #print + #print toggle_rotation,toggle_altx,toggle_s + #print m1 + #print m12 + + + if self == self.__class__(): + if self.cofactor == 4: + # Hacks for identity! + if toggle_altx: m12 = 1 + elif toggle_s: m1 = 1 + elif toggle_r: continue + ## BOTH??? + + else: + m12 = 1 + imi = self.isoMagic * self.i + if toggle_rotation: + if toggle_altx: m1 = -imi + else: m1 = +imi + else: + if toggle_altx: m1 = 0 + else: m1 = a-d + + rnum = (d*a*m12-m1) + rden = ((d*a-1)*m12+m1) + if swap: rnum,rden = rden,rnum + + ok,sr = isqrt_i(rnum*rden*self.qnr) + if not ok: continue + sr *= rnum + #print "Works! %d %x" % (swap,sr) + + if negative(sr) != toggle_r: sr = -sr + ret = self.gfToBytes(sr) + if self.elligator(ret) != self and self.elligator(ret) != -self: + print "WRONG!",[toggle_rotation,toggle_altx,toggle_s] + if self.elligator(ret) == -self and self != -self: print "Negated!",[toggle_rotation,toggle_altx,toggle_s] + rets.append(bytes(ret)) + return rets + + @optimized_version_of("encodeSpec") + def encode(self): + """Encode, optimized version""" + return self.gfToBytes(self.toJacobiQuartic()[0]) @classmethod @optimized_version_of("decodeSpec") @@ -351,7 +480,7 @@ class Decaf_1_1_Point(QuotientEdwardsPoint): a,d = cls.a,cls.d s = cls.bytesToGf(s,mustBePositive=True) - if s==0: return cls() + #if s==0: return cls() s2 = s^2 den = 1+a*s2 num = den^2 - 4*d*s2 @@ -374,13 +503,63 @@ class Decaf_1_1_Point(QuotientEdwardsPoint): x = 2*s / (1+a*s^2) y = (1-a*s^2) / t return cls(x,sgn*y) + + @optimized_version_of("doubleAndEncodeSpec") + def doubleAndEncode(self): + X,Y,Z,T = self.xyzt() + a,d = self.a,self.d + + if self.cofactor == 8: + # Cofactor 8 version + # Simulate IMAGINE_TWIST because that's how libdecaf does it + X = self.i*X + T = self.i*T + a = -a + d = -d + # TODO: This is only being called for a=-1, so could + # be wrong for a=1 + + e = 2*X*Y + f = Y^2+a*X^2 + g = Y^2-a*X^2 + h = Z^2-d*T^2 + + eim = e*self.isoMagic + inv = 1/(eim*g*f*h) + fh_inv = eim*g*inv*self.i + + if negative(eim*g*fh_inv): + idf = g*self.isoMagic*self.i + bar = f + foo = g + test = eim*f + else: + idf = eim + bar = h + foo = -eim + test = g*h + + if negative(test*fh_inv): bar =- bar + s = idf*(foo+bar)*inv*f*h + + else: + xy = X*Y + h = Z^2-d*T^2 + inv = 1/(xy*h) + if negative(inv*2*xy^2*self.isoMagic): tmp = Y + else: tmp = X + s = tmp^2*h*inv # = X/Y or Y/X, interestingly + + return self.gfToBytes(s,mustBePositive=True) @classmethod - def elligatorSpec(cls,r0): + def elligatorSpec(cls,r0,fromR=False): a,d = cls.a,cls.d - r = cls.qnr * cls.bytesToGf(r0)^2 + if fromR: r = r0 + else: r = cls.qnr * cls.bytesToGf(r0,mustBeProper=False,maskHiBits=True)^2 den = (d*r-(d-a))*((d-a)*r-d) + if den == 0: return cls() n1 = (r+1)*(a-2*d)/den n2 = r*n1 if is_square(n1): @@ -394,7 +573,7 @@ class Decaf_1_1_Point(QuotientEdwardsPoint): @optimized_version_of("elligatorSpec") def elligator(cls,r0): a,d = cls.a,cls.d - r0 = cls.bytesToGf(r0) + r0 = cls.bytesToGf(r0,mustBeProper=False,maskHiBits=True) r = cls.qnr * r0^2 den = (d*r-(d-a))*((d-a)*r-d) num = (r+1)*(a-2*d) @@ -408,6 +587,40 @@ class Decaf_1_1_Point(QuotientEdwardsPoint): if negative(s) == iss: s = -s return cls.fromJacobiQuartic(s,t) + def elligatorInverseBruteForce(self): + """Invert Elligator using SAGE's polynomial solver""" + a,d = self.a,self.d + R. = self.F[] + r = self.qnr * r0^2 + den = (d*r-(d-a))*((d-a)*r-d) + n1 = (r+1)*(a-2*d)/den + n2 = r*n1 + ret = set() + for s2,t in [(n1, -(r-1)*(a-2*d)^2 / den - 1), + (n2,r*(r-1)*(a-2*d)^2 / den - 1)]: + x2 = 4*s2/(1+a*s2)^2 + y = (1-a*s2) / t + + selfT = self + for i in xrange(self.cofactor/2): + xT,yT = selfT + polyX = xT^2-x2 + polyY = yT-y + sx = set(r for r,_ in polyX.numerator().roots()) + sy = set(r for r,_ in polyY.numerator().roots()) + ret = ret.union(sx.intersection(sy)) + + selfT = selfT.torque() + + ret = [self.gfToBytes(r) for r in ret] + + for r in ret: + assert self.elligator(r) in [self,-self] + + ret = [r for r in ret if self.elligator(r) == self] + + return ret + class Ed25519Point(RistrettoPoint): F = GF(2^255-19) d = F(-121665/121666) @@ -455,7 +668,7 @@ class IsoEd448Point(RistrettoPoint): @classmethod def base(cls): return cls( # RFC has it wrong - -345397493039729516374008604150537410266655260075183290216406970281645695073672344430481787759340633221708391583424041788924124567700732, + 345397493039729516374008604150537410266655260075183290216406970281645695073672344430481787759340633221708391583424041788924124567700732, -363419362147803445274661903944002267176820680343659030140745099590306164083365386343198191849338272965044442230921818680526749009182718 ) @@ -464,7 +677,6 @@ class TwistedEd448GoldilocksPoint(Decaf_1_1_Point): d = F(-39082) a = F(-1) qnr = -1 - magic = isqrt(a*d-1) cofactor = 4 encLen = 56 isoMagic = IsoEd448Point.magic @@ -478,14 +690,13 @@ class Ed448GoldilocksPoint(Decaf_1_1_Point): d = F(-39081) a = F(1) qnr = -1 - magic = isqrt(a*d-1) cofactor = 4 encLen = 56 isoMagic = IsoEd448Point.magic @classmethod def base(cls): - return -2*cls( # FIXME: make not negative + return 2*cls( 224580040295924300187604334099896036246789641632564134246125461686950415467406032909029192869357953282578032075146446173674602635247710, 298819210078481492676017930443930673437544040154080242095928241372331506189835876003536878655418784733982303233503462500531545062832660 ) @@ -532,19 +743,29 @@ def test(cls,n): P = cls.base() - print "base", list(P.encode()) - for i in xrange(16): - Q = P*i - print i, list(Q.encode()) - Q = cls() for i in xrange(n): - #print i, binascii.hexlify(Q.encode()) - QQ = cls.decode(Q.encode()) + #print binascii.hexlify(Q.encode()) + QE = Q.encode() + QQ = cls.decode(QE) if QQ != Q: raise TestFailedException("Round trip %s != %s" % (str(QQ),str(Q))) + + # Testing s -> 1/s: encodes -point on cofactor + s = cls.bytesToGf(QE) + if s != 0: + ss = cls.gfToBytes(1/s,mustBePositive=True) + try: + QN = cls.decode(ss) + if cls.cofactor == 8: + raise TestFailedException("1/s shouldnt work for cofactor 8") + if QN != -Q: + raise TestFailedException("s -> 1/s should negate point for cofactor 4") + except InvalidEncodingException as e: + # Should be raised iff cofactor==8 + if cls.cofactor == 4: + raise TestFailedException("s -> 1/s should work for cofactor 4") QT = Q - QE = Q.encode() for h in xrange(cls.cofactor): QT = QT.torque() if QT.encode() != QE: @@ -559,27 +780,26 @@ def test(cls,n): Q2 = Q0*(r+1) if Q1 + Q0 != Q2: raise TestFailedException("Scalarmul doesn't work") Q = Q1 - -test(Ed25519Point,100) -#test(NegEd25519Point,100) -#test(IsoEd25519Point,100) -#test(IsoEd448Point,100) -#test(TwistedEd448GoldilocksPoint,100) -#test(Ed448GoldilocksPoint,100) - def testElligator(cls,n): print "Testing elligator on %s" % cls.__name__ for i in xrange(n): r = randombytes(cls.encLen) - Q = cls.elligator(r) - print list(r), list(Q.encode()) - -testElligator(Ed25519Point,100) -#testElligator(NegEd25519Point,100) -#testElligator(IsoEd448Point,100) -#testElligator(Ed448GoldilocksPoint,100) -#testElligator(TwistedEd448GoldilocksPoint,100) + P = cls.elligator(r) + if hasattr(P,"invertElligator"): + iv = P.invertElligator() + modr = bytes(cls.gfToBytes(cls.bytesToGf(r,mustBeProper=False,maskHiBits=True))) + iv2 = P.torque().invertElligator() + if modr not in iv: print "Failed to invert Elligator!" + if len(iv) != len(set(iv)): + print "Elligator inverses not unique!", len(set(iv)), len(iv) + if iv != iv2: + print "Elligator is untorqueable!" + #print [binascii.hexlify(j) for j in iv] + #print [binascii.hexlify(j) for j in iv2] + #break + else: + pass # TODO def gangtest(classes,n): print "Gang test",[cls.__name__ for cls in classes] @@ -607,5 +827,31 @@ def gangtest(classes,n): for c,ret in zip(classes,rets): print c,binascii.hexlify(ret) print -gangtest([IsoEd448Point,TwistedEd448GoldilocksPoint,Ed448GoldilocksPoint],100) -gangtest([Ed25519Point,IsoEd25519Point],100) + +def testDoubleAndEncode(cls,n): + print "Testing doubleAndEncode on %s" % cls.__name__ + for i in xrange(n): + r1 = randombytes(cls.encLen) + r2 = randombytes(cls.encLen) + u = cls.elligator(r1) + cls.elligator(r2) + u.doubleAndEncode() + +testDoubleAndEncode(Ed25519Point,100) +testDoubleAndEncode(NegEd25519Point,100) +testDoubleAndEncode(IsoEd25519Point,100) +testDoubleAndEncode(IsoEd448Point,100) +testDoubleAndEncode(TwistedEd448GoldilocksPoint,100) +#test(Ed25519Point,100) +#test(NegEd25519Point,100) +#test(IsoEd25519Point,100) +#test(IsoEd448Point,100) +#test(TwistedEd448GoldilocksPoint,100) +#test(Ed448GoldilocksPoint,100) +#testElligator(Ed25519Point,100) +#testElligator(NegEd25519Point,100) +#testElligator(IsoEd25519Point,100) +#testElligator(IsoEd448Point,100) +#testElligator(Ed448GoldilocksPoint,100) +#testElligator(TwistedEd448GoldilocksPoint,100) +#gangtest([IsoEd448Point,TwistedEd448GoldilocksPoint,Ed448GoldilocksPoint],100) +#gangtest([Ed25519Point,IsoEd25519Point],100) From b20ccbd6854e5c487e7f06d4f3671b27ba6821b8 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Tue, 16 Jan 2018 15:41:05 -0800 Subject: [PATCH 3/8] Implement batched encoding for RistrettoPoints --- src/ristretto.rs | 118 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 118 insertions(+) diff --git a/src/ristretto.rs b/src/ristretto.rs index 61b451c..50f9b96 100644 --- a/src/ristretto.rs +++ b/src/ristretto.rs @@ -666,6 +666,86 @@ impl RistrettoPoint { CompressedRistretto(s.to_bytes()) } + /// Double-and-compress a batch of points. + pub fn double_and_compress_batch<'a, I>(points: I) -> Vec + where I: IntoIterator + { + #[derive(Copy, Clone, Debug)] + struct BatchCompressState { + e: FieldElement, + f: FieldElement, + g: FieldElement, + h: FieldElement, + eg: FieldElement, + fh: FieldElement, + } + + impl BatchCompressState { + fn efgh(&self) -> FieldElement { + &self.eg * &self.fh + } + } + + impl<'a> From<&'a RistrettoPoint> for BatchCompressState { + fn from(P: &'a RistrettoPoint) -> BatchCompressState { + let XX = P.0.X.square(); + let YY = P.0.Y.square(); + let ZZ = P.0.Z.square(); + let dTT = &P.0.T.square() * &constants::EDWARDS_D; + + let e = &P.0.X * &(&P.0.Y + &P.0.Y); // = 2*X*Y + let f = &ZZ + &dTT; // = Z^2 + d*T^2 + let g = &YY + &XX; // = Y^2 - a*X^2 + let h = &ZZ - &dTT; // = Z^2 - d*T^2 + + let eg = &e * &g; + let fh = &f * &h; + + BatchCompressState{ e: e, f: f, g: g, h: h, eg: eg, fh: fh } + } + } + + let states: Vec = points.into_iter().map(|P| BatchCompressState::from(P)).collect(); + + let mut invs: Vec = states.iter().map(|state| state.efgh()).collect(); + + FieldElement::batch_invert(&mut invs[..]); + + states.iter().zip(invs.iter()).map(|(state, inv): (&BatchCompressState, &FieldElement)| { + let Zinv = &state.eg * &inv; + let Tinv = &state.fh * &inv; + + let mut magic = constants::INVSQRT_A_MINUS_D; + + let negcheck1 = (&state.eg * &Zinv).is_negative(); + + let mut e = state.e; + let mut g = state.g; + let mut h = state.h; + + let minus_e = -&e; + let f_times_sqrta = &state.f * &constants::SQRT_M1; + + e.conditional_assign(&state.g, negcheck1); + g.conditional_assign(&minus_e, negcheck1); + h.conditional_assign(&f_times_sqrta, negcheck1); + + magic.conditional_assign(&constants::SQRT_M1, negcheck1); + + let negcheck2 = (&(&h * &e) * &Zinv).is_negative(); + + g.conditional_negate(negcheck2); + + let mut s = &(&h - &g) * &(&magic * &(&g * &Tinv)); + + let s_is_negative = s.is_negative(); + s.conditional_negate(s_is_negative); + + CompressedRistretto(s.to_bytes()) + }).collect() + } + + /// Return the coset self + E[4], for debugging. fn coset4(&self) -> [ExtendedPoint; 4] { [ self.0 @@ -1216,6 +1296,20 @@ mod test { } } + #[test] + fn double_and_compress_1024_random_points() { + let mut rng = OsRng::new().unwrap(); + + let points: Vec = + (0..1024).map(|_| RistrettoPoint::random(&mut rng)).collect(); + + let compressed = RistrettoPoint::double_and_compress_batch(&points); + + for (P, P2_compressed) in points.iter().zip(compressed.iter()) { + assert_eq!(*P2_compressed, (P + P).compress()); + } + } + #[test] fn random_is_valid() { let mut rng = OsRng::new().unwrap(); @@ -1254,4 +1348,28 @@ mod bench { let P = B * &Scalar::random(&mut rng); b.iter(|| P.compress()); } + + fn double_and_compress_n_random_points(n: usize, b: &mut Bencher) { + let mut rng = OsRng::new().unwrap(); + + let points: Vec = + (0..n).map(|_| RistrettoPoint::random(&mut rng)).collect(); + + b.iter(|| RistrettoPoint::double_and_compress_batch(&points) ); + } + + #[bench] + fn double_and_compress_16_random_points(b: &mut Bencher) { + double_and_compress_n_random_points(16, b); + } + + #[bench] + fn double_and_compress_128_random_points(b: &mut Bencher) { + double_and_compress_n_random_points(128, b); + } + + #[bench] + fn double_and_compress_1024_random_points(b: &mut Bencher) { + double_and_compress_n_random_points(1024, b); + } } From 1f821d34a54ee489a6347afe4578e6ce10d2a761 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Fri, 19 Jan 2018 17:08:12 -0800 Subject: [PATCH 4/8] Feature-gate batch inversion and compression on alloc --- src/field.rs | 1 + src/ristretto.rs | 1 + 2 files changed, 2 insertions(+) diff --git a/src/field.rs b/src/field.rs index 23362e8..e6f2587 100644 --- a/src/field.rs +++ b/src/field.rs @@ -184,6 +184,7 @@ impl FieldElement { /// /// This function is most efficient when the batch size (slice /// length) is a power of 2. + #[cfg(any(feature = "alloc", feature = "std"))] pub fn batch_invert(inputs: &mut [FieldElement]) { // First, compute the product of all inputs using a product // tree: diff --git a/src/ristretto.rs b/src/ristretto.rs index 50f9b96..ee6260e 100644 --- a/src/ristretto.rs +++ b/src/ristretto.rs @@ -667,6 +667,7 @@ impl RistrettoPoint { } /// Double-and-compress a batch of points. + #[cfg(any(feature = "alloc", feature = "std"))] pub fn double_and_compress_batch<'a, I>(points: I) -> Vec where I: IntoIterator { From d8d235fb485710886b8cbda5e6d85375e9933d29 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Fri, 19 Jan 2018 13:49:10 -0800 Subject: [PATCH 5/8] Move pow2k into the backends and use it to implement square() --- src/backend/u32/field.rs | 10 +++ src/backend/u64/field.rs | 153 ++++++++++++++++++++++++++------------- src/field.rs | 9 --- 3 files changed, 112 insertions(+), 60 deletions(-) diff --git a/src/backend/u32/field.rs b/src/backend/u32/field.rs index a7a6419..fef4615 100644 --- a/src/backend/u32/field.rs +++ b/src/backend/u32/field.rs @@ -264,6 +264,16 @@ impl FieldElement32 { ]) } + /// Given `k > 0`, return `self^(2^k)`. + pub fn pow2k(&self, k: u32) -> FieldElement32 { + debug_assert!( k > 0 ); + let mut z = self.square(); + for _ in 1..k { + z = z.square(); + } + z + } + /// Given unreduced coefficients `z[0], ..., z[9]` of any size, /// carry and reduce them mod p to obtain a `FieldElement32` /// whose coefficients have excess `b < 0.007`. diff --git a/src/backend/u64/field.rs b/src/backend/u64/field.rs index 70d7fa8..1239491 100644 --- a/src/backend/u64/field.rs +++ b/src/backend/u64/field.rs @@ -337,73 +337,124 @@ impl FieldElement64 { s } - #[inline(always)] - fn square_inner(&self) -> [u64; 5] { + /// Given `k > 0`, return `self^(2^k)`. + pub fn pow2k(&self, mut k: u32) -> FieldElement64 { + + debug_assert!( k > 0 ); + /// Multiply two 64-bit integers with 128 bits of output. #[inline(always)] fn m(x: u64, y: u64) -> u128 { (x as u128) * (y as u128) } - // Alias self, _rhs for more readable formulas - let a: &[u64; 5] = &self.0; + let mut a: [u64; 5] = self.0; - // Precomputation: 64-bit multiply by 19 - let a3_19 = 19 * a[3]; - let a4_19 = 19 * a[4]; + loop { + // Precondition: assume input limbs a[i] are bounded as + // + // a[i] < 2^(51 + b) + // + // where b is a real parameter measuring the "bit excess" of the limbs. - // Multiply to get 128-bit coefficients of output - let c0: u128 = m(a[0], a[0]) + 2*( m(a[1], a4_19) + m(a[2], a3_19) ); - let mut c1: u128 = m(a[3], a3_19) + 2*( m(a[0], a[1]) + m(a[2], a4_19) ); - let mut c2: u128 = m(a[1], a[1]) + 2*( m(a[0], a[2]) + m(a[4], a3_19) ); - let mut c3: u128 = m(a[4], a4_19) + 2*( m(a[0], a[3]) + m(a[1], a[2]) ); - let mut c4: u128 = m(a[2], a[2]) + 2*( m(a[0], a[4]) + m(a[1], a[3]) ); + // Precomputation: 64-bit multiply by 19. + // + // This fits into a u64 whenever 51 + b + lg(19) < 64. + // + // Since 51 + b + lg(19) < 51 + 4.25 + b + // = 55.25 + b, + // this fits if b < 8.75. + let a3_19 = 19 * a[3]; + let a4_19 = 19 * a[4]; - // Same bound as in multiply: - // c[i] < 2^2b * (1+i + (4-i)*19) < 2^(2b + lg(1+4*19)) < 2^(2b + 6.27) - // where b is the bitlength of the input limbs. - // - // The carry (c[i] >> 51) fits into a u64 iff 2b+6.27 < 64+51 iff b <= 54. - // After the first carry pass, all c[i] fit into u64. - debug_assert!(a[0] < (1 << 54)); - debug_assert!(a[1] < (1 << 54)); - debug_assert!(a[2] < (1 << 54)); - debug_assert!(a[3] < (1 << 54)); - debug_assert!(a[4] < (1 << 54)); + // Multiply to get 128-bit coefficients of output. + // + // The 128-bit multiplications by 2 turn into 1 slr + 1 slrd each, + // which doesn't seem any better or worse than doing them as precomputations + // on the 64-bit inputs. + let c0: u128 = m(a[0], a[0]) + 2*( m(a[1], a4_19) + m(a[2], a3_19) ); + let mut c1: u128 = m(a[3], a3_19) + 2*( m(a[0], a[1]) + m(a[2], a4_19) ); + let mut c2: u128 = m(a[1], a[1]) + 2*( m(a[0], a[2]) + m(a[4], a3_19) ); + let mut c3: u128 = m(a[4], a4_19) + 2*( m(a[0], a[3]) + m(a[1], a[2]) ); + let mut c4: u128 = m(a[2], a[2]) + 2*( m(a[0], a[4]) + m(a[1], a[3]) ); - // The 128-bit output limbs are stored in two 64-bit registers (low/high part). - // By rebinding the names after carrying, we free the upper registers for reuse. - let low_51_bit_mask = (1u64 << 51) - 1; - c1 += (c0 >> 51) as u128; - let mut c0: u64 = (c0 as u64) & low_51_bit_mask; - c2 += (c1 >> 51) as u128; - let c1: u64 = (c1 as u64) & low_51_bit_mask; - c3 += (c2 >> 51) as u128; - let c2: u64 = (c2 as u64) & low_51_bit_mask; - c4 += (c3 >> 51) as u128; - let c3: u64 = (c3 as u64) & low_51_bit_mask; - c0 += ((c4 >> 51) as u64) * 19; - let c4: u64 = (c4 as u64) & low_51_bit_mask; + // Same bound as in multiply: + // c[i] < 2^(102 + 2*b) * (1+i + (4-i)*19) + // < 2^(102 + lg(1 + 4*19) + 2*b) + // < 2^(108.27 + 2*b) + // + // The carry (c[i] >> 51) fits into a u64 when + // 108.27 + 2*b - 51 < 64 + // 2*b < 6.73 + // b < 3.365. + // + // So we require b < 3 to ensure this fits. + debug_assert!(a[0] < (1 << 54)); + debug_assert!(a[1] < (1 << 54)); + debug_assert!(a[2] < (1 << 54)); + debug_assert!(a[3] < (1 << 54)); + debug_assert!(a[4] < (1 << 54)); - // Now c_i all fit into u64, but are not yet bounded by 2^51. - [c0,c1,c2,c3,c4] + const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1; + + // Casting to u64 and back tells the compiler that the carry is bounded by 2^64, so + // that the addition is a u128 + u64 rather than u128 + u128. + c1 += ((c0 >> 51) as u64) as u128; + a[0] = (c0 as u64) & LOW_51_BIT_MASK; + + c2 += ((c1 >> 51) as u64) as u128; + a[1] = (c1 as u64) & LOW_51_BIT_MASK; + + c3 += ((c2 >> 51) as u64) as u128; + a[2] = (c2 as u64) & LOW_51_BIT_MASK; + + c4 += ((c3 >> 51) as u64) as u128; + a[3] = (c3 as u64) & LOW_51_BIT_MASK; + + let carry: u64 = (c4 >> 51) as u64; + a[4] = (c4 as u64) & LOW_51_BIT_MASK; + + // To see that this does not overflow, we need a[0] + carry * 19 < 2^64. + // + // c4 < a2^2 + 2*a0*a4 + 2*a1*a3 + (carry from c3) + // < 2^(102 + 2*b + lg(5)) + 2^64. + // + // When b < 3 we get + // + // c4 < 2^110.33 so that carry < 2^59.33 + // + // so that + // + // a[0] + carry * 19 < 2^51 + 19 * 2^59.33 < 2^63.58 + // + // and there is no overflow. + a[0] = a[0] + carry * 19; + + // Now a[1] < 2^51 + 2^(64 -51) = 2^51 + 2^13 < 2^(51 + eps). + a[1] += a[0] >> 51; + a[0] &= LOW_51_BIT_MASK; + + // Now all a[i] < 2^(51 + eps) and a = self^(2^k). + + k = k - 1; + if k == 0 { + break; + } + } + + FieldElement64(a) } /// Returns the square of this field element. pub fn square(&self) -> FieldElement64 { - FieldElement64::reduce(self.square_inner()) + self.pow2k(1) } /// Returns 2 times the square of this field element. pub fn square2(&self) -> FieldElement64 { - let mut limbs = self.square_inner(); - // For this to work, need to have 1 extra bit of headroom after carry - // --> max 53 bit inputs, not 54 - // - // XXX check that this is correct; I think it isn't -- hdevalence - limbs[0] *= 2; - limbs[1] *= 2; - limbs[2] *= 2; - limbs[3] *= 2; - limbs[4] *= 2; - FieldElement64::reduce(limbs) + let mut square = self.pow2k(1); + for i in 0..5 { + square.0[i] *= 2; + } + + square } } diff --git a/src/field.rs b/src/field.rs index e6f2587..756423e 100644 --- a/src/field.rs +++ b/src/field.rs @@ -124,15 +124,6 @@ impl FieldElement { byte_is_nonzero(x) } - #[inline] - #[allow(dead_code)] - /// Requires k > 0; raise self to the 2^(2^k)-th power. - fn pow2k(&self, k: u32) -> FieldElement { - let mut z = self.square(); - for _ in 1..k { z = z.square(); } - z - } - /// Compute (self^(2^250-1), self^11), used as a helper function /// within invert() and pow22523(). /// From ec50ec96f79335053716d123c7fc40d3ab969af9 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Fri, 19 Jan 2018 15:06:29 -0800 Subject: [PATCH 6/8] Use parallel carry-ins and carry-outs in FieldElement64::reduce() --- src/backend/u64/field.rs | 62 ++++++++++++++++++++++++++++------------ 1 file changed, 43 insertions(+), 19 deletions(-) diff --git a/src/backend/u64/field.rs b/src/backend/u64/field.rs index 1239491..398e48e 100644 --- a/src/backend/u64/field.rs +++ b/src/backend/u64/field.rs @@ -78,8 +78,7 @@ impl<'a, 'b> Sub<&'b FieldElement64> for &'a FieldElement64 { // just bigger than _rhs and avoid having to do a reduction. // // Since we don't yet have type-level integers to do this, we - // have to add an explicit reduction call here, which is a - // somewhat significant cost. + // have to add an explicit reduction call here. FieldElement64::reduce([ (self.0[0] + 36028797018963664u64) - _rhs.0[0], (self.0[1] + 36028797018963952u64) - _rhs.0[1], @@ -200,20 +199,39 @@ impl FieldElement64 { FieldElement64([2251799813685228, 2251799813685247, 2251799813685247, 2251799813685247, 2251799813685247]) } - /// Given 64-bit limbs, reduce to enforce the bound c_i < 2^51. + /// Given 64-bit input limbs, reduce to enforce the bound 2^(51 + eps). #[inline(always)] fn reduce(mut limbs: [u64; 5]) -> FieldElement64 { - let low_51_bit_mask = (1u64 << 51) - 1; - limbs[1] += limbs[0] >> 51; - limbs[0] = limbs[0] & low_51_bit_mask; - limbs[2] += limbs[1] >> 51; - limbs[1] = limbs[1] & low_51_bit_mask; - limbs[3] += limbs[2] >> 51; - limbs[2] = limbs[2] & low_51_bit_mask; - limbs[4] += limbs[3] >> 51; - limbs[3] = limbs[3] & low_51_bit_mask; - limbs[0] += (limbs[4] >> 51) * 19; - limbs[4] = limbs[4] & low_51_bit_mask; + const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1; + + // Since the input limbs are bounded by 2^64, the biggest + // carry-out is bounded by 2^13. + // + // The biggest carry-in is c4 * 19, resulting in + // + // 2^51 + 19*2^13 < 2^51.0000000001 + // + // Because we don't need to canonicalize, only to reduce the + // limb sizes, it's OK to do a "weak reduction", where we + // compute the carry-outs in parallel. + + let c0 = limbs[0] >> 51; + let c1 = limbs[1] >> 51; + let c2 = limbs[2] >> 51; + let c3 = limbs[3] >> 51; + let c4 = limbs[4] >> 51; + + limbs[0] &= LOW_51_BIT_MASK; + limbs[1] &= LOW_51_BIT_MASK; + limbs[2] &= LOW_51_BIT_MASK; + limbs[3] &= LOW_51_BIT_MASK; + limbs[4] &= LOW_51_BIT_MASK; + + limbs[0] += c4 * 19; + limbs[1] += c0; + limbs[2] += c1; + limbs[3] += c2; + limbs[4] += c3; FieldElement64(limbs) } @@ -260,18 +278,24 @@ impl FieldElement64 { /// Serialize this `FieldElement64` to a 32-byte array. The /// encoding is canonical. pub fn to_bytes(&self) -> [u8; 32] { - // This reduces to the range [0,2^255), but we need [0,2^255-19). - let mut limbs = FieldElement64::reduce(self.0).0; - // Let h = limbs[0] + limbs[1]*2^51 + ... + limbs[4]*2^204. // - // Write h = pq + r with 0 <= r < p. We want to compute r = h mod p. + // Write h = pq + r with 0 <= r < p. // - // Since h < 2^255, q = 0 or 1, with q = 0 when h < p and q = 1 when h >= p. + // We want to compute r = h mod p. + // + // If h < 2*p = 2^256 - 38, + // then q = 0 or 1, + // + // with q = 0 when h < p + // and q = 1 when h >= p. // // Notice that h >= p <==> h + 19 >= p + 19 <==> h + 19 >= 2^255. // Therefore q can be computed as the carry bit of h + 19. + // First, reduce the limbs to ensure h < 2*p. + let mut limbs = FieldElement64::reduce(self.0).0; + let mut q = (limbs[0] + 19) >> 51; q = (limbs[1] + q) >> 51; q = (limbs[2] + q) >> 51; From 0e9b8e0a6e3a9c09226aaad434d1bd5858daf020 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Fri, 19 Jan 2018 18:23:41 -0800 Subject: [PATCH 7/8] Refactor multiply implementation to eliminate reduce() call --- src/backend/u64/field.rs | 91 ++++++++++++++++++++++++++++++---------- 1 file changed, 69 insertions(+), 22 deletions(-) diff --git a/src/backend/u64/field.rs b/src/backend/u64/field.rs index 398e48e..87f9eeb 100644 --- a/src/backend/u64/field.rs +++ b/src/backend/u64/field.rs @@ -108,7 +108,19 @@ impl<'a, 'b> Mul<&'b FieldElement64> for &'a FieldElement64 { let a: &[u64; 5] = &self.0; let b: &[u64; 5] = &_rhs.0; - // 64-bit precomputations to avoid 128-bit multiplications + // Precondition: assume input limbs a[i], b[i] are bounded as + // + // a[i], b[i] < 2^(51 + b) + // + // where b is a real parameter measuring the "bit excess" of the limbs. + + // 64-bit precomputations to avoid 128-bit multiplications. + // + // This fits into a u64 whenever 51 + b + lg(19) < 64. + // + // Since 51 + b + lg(19) < 51 + 4.25 + b + // = 55.25 + b, + // this fits if b < 8.75. let b1_19 = b[1] * 19; let b2_19 = b[2] * 19; let b3_19 = b[3] * 19; @@ -121,34 +133,69 @@ impl<'a, 'b> Mul<&'b FieldElement64> for &'a FieldElement64 { let mut c3: u128 = m(a[3],b[0]) + m(a[2],b[1]) + m(a[1],b[2]) + m(a[0],b[3]) + m(a[4],b4_19); let mut c4: u128 = m(a[4],b[0]) + m(a[3],b[1]) + m(a[2],b[2]) + m(a[1],b[3]) + m(a[0],b[4]); - // Now c[i] < 2^2b * (1+i + (4-i)*19) < 2^(2b + lg(1+4*19)) < 2^(2b + 6.27) - // where b is the bitlength of the input limbs. - - // The carry (c[i] >> 51) fits into a u64 iff 2b+6.27 < 64+51 iff b <= 54. - // After the first carry pass, all c[i] fit into u64. + // How big are the c[i]? We have + // + // c[i] < 2^(102 + 2*b) * (1+i + (4-i)*19) + // < 2^(102 + lg(1 + 4*19) + 2*b) + // < 2^(108.27 + 2*b) + // + // The carry (c[i] >> 51) fits into a u64 when + // 108.27 + 2*b - 51 < 64 + // 2*b < 6.73 + // b < 3.365. + // + // So we require b < 3 to ensure this fits. debug_assert!(a[0] < (1 << 54)); debug_assert!(b[0] < (1 << 54)); debug_assert!(a[1] < (1 << 54)); debug_assert!(b[1] < (1 << 54)); debug_assert!(a[2] < (1 << 54)); debug_assert!(b[2] < (1 << 54)); debug_assert!(a[3] < (1 << 54)); debug_assert!(b[3] < (1 << 54)); debug_assert!(a[4] < (1 << 54)); debug_assert!(b[4] < (1 << 54)); - // The 128-bit output limbs are stored in two 64-bit registers - // (low/high part). By rebinding the names after carrying, we - // inform LLVM that the values have shrunk, so it can - // efficiently allocate registers. - let low_51_bit_mask = (1u64 << 51) - 1; - c1 += (c0 >> 51) as u128; - let mut c0: u64 = (c0 as u64) & low_51_bit_mask; - c2 += (c1 >> 51) as u128; - let c1: u64 = (c1 as u64) & low_51_bit_mask; - c3 += (c2 >> 51) as u128; - let c2: u64 = (c2 as u64) & low_51_bit_mask; - c4 += (c3 >> 51) as u128; - let c3: u64 = (c3 as u64) & low_51_bit_mask; - c0 += ((c4 >> 51) as u64) * 19; - let c4: u64 = (c4 as u64) & low_51_bit_mask; + // Casting to u64 and back tells the compiler that the carry is + // bounded by 2^64, so that the addition is a u128 + u64 rather + // than u128 + u128. - FieldElement64::reduce([c0,c1,c2,c3,c4]) + const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1; + let mut out = [0u64; 5]; + + c1 += ((c0 >> 51) as u64) as u128; + out[0] = (c0 as u64) & LOW_51_BIT_MASK; + + c2 += ((c1 >> 51) as u64) as u128; + out[1] = (c1 as u64) & LOW_51_BIT_MASK; + + c3 += ((c2 >> 51) as u64) as u128; + out[2] = (c2 as u64) & LOW_51_BIT_MASK; + + c4 += ((c3 >> 51) as u64) as u128; + out[3] = (c3 as u64) & LOW_51_BIT_MASK; + + let carry: u64 = (c4 >> 51) as u64; + out[4] = (c4 as u64) & LOW_51_BIT_MASK; + + // To see that this does not overflow, we need out[0] + carry * 19 < 2^64. + // + // c4 < a0*b4 + a1*b3 + a2*b2 + a3*b1 + a4*b0 + (carry from c3) + // < 5*(2^(51 + b) * 2^(51 + b)) + (carry from c3) + // < 2^(102 + 2*b + lg(5)) + 2^64. + // + // When b < 3 we get + // + // c4 < 2^110.33 so that carry < 2^59.33 + // + // so that + // + // out[0] + carry * 19 < 2^51 + 19 * 2^59.33 < 2^63.58 + // + // and there is no overflow. + out[0] = out[0] + carry * 19; + + // Now out[1] < 2^51 + 2^(64 -51) = 2^51 + 2^13 < 2^(51 + eps). + out[1] += out[0] >> 51; + out[0] &= LOW_51_BIT_MASK; + + // Now out[i] < 2^(51 + eps) for all i. + FieldElement64(out) } } From 86fd06db0050a177a3248ebd54553e637ce23fc6 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Wed, 24 Jan 2018 12:42:57 -0800 Subject: [PATCH 8/8] s/eps/epsilon/ for clarity --- src/backend/u64/field.rs | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/backend/u64/field.rs b/src/backend/u64/field.rs index 87f9eeb..20cb6da 100644 --- a/src/backend/u64/field.rs +++ b/src/backend/u64/field.rs @@ -190,11 +190,11 @@ impl<'a, 'b> Mul<&'b FieldElement64> for &'a FieldElement64 { // and there is no overflow. out[0] = out[0] + carry * 19; - // Now out[1] < 2^51 + 2^(64 -51) = 2^51 + 2^13 < 2^(51 + eps). + // Now out[1] < 2^51 + 2^(64 -51) = 2^51 + 2^13 < 2^(51 + epsilon). out[1] += out[0] >> 51; out[0] &= LOW_51_BIT_MASK; - // Now out[i] < 2^(51 + eps) for all i. + // Now out[i] < 2^(51 + epsilon) for all i. FieldElement64(out) } } @@ -246,7 +246,7 @@ impl FieldElement64 { FieldElement64([2251799813685228, 2251799813685247, 2251799813685247, 2251799813685247, 2251799813685247]) } - /// Given 64-bit input limbs, reduce to enforce the bound 2^(51 + eps). + /// Given 64-bit input limbs, reduce to enforce the bound 2^(51 + epsilon). #[inline(always)] fn reduce(mut limbs: [u64; 5]) -> FieldElement64 { const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1; @@ -499,11 +499,11 @@ impl FieldElement64 { // and there is no overflow. a[0] = a[0] + carry * 19; - // Now a[1] < 2^51 + 2^(64 -51) = 2^51 + 2^13 < 2^(51 + eps). + // Now a[1] < 2^51 + 2^(64 -51) = 2^51 + 2^13 < 2^(51 + epsilon). a[1] += a[0] >> 51; a[0] &= LOW_51_BIT_MASK; - // Now all a[i] < 2^(51 + eps) and a = self^(2^k). + // Now all a[i] < 2^(51 + epsilon) and a = self^(2^k). k = k - 1; if k == 0 {