From de377290ee94d7e5089c4d245b90a979cb402b81 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Thu, 4 Jan 2018 13:43:43 -0800 Subject: [PATCH 01/13] 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 02/13] 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 03/13] 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 04/13] 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 05/13] 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 06/13] 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 07/13] 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 460355d5e06d590ab289ed063045793734cc4bd4 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Sat, 20 Jan 2018 02:30:33 +0000 Subject: [PATCH 08/13] Bump subtle dependency to 0.5.0. --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index f4a60ee..68e3127 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -40,7 +40,7 @@ version = "0.4" version = "0.7" [dependencies.subtle] -version = "^0.3" +version = "^0.5" default-features = false [dependencies.clear_on_drop] From 52d600db9500aaa372d7b4dc681148c18428b9c1 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Wed, 24 Jan 2018 10:58:22 -0800 Subject: [PATCH 09/13] Implement public-facing ops traits on all combos of &T/T The public-facing types with arithmetic operations are: - `Scalar`s - `ExtendedPoint`s - `RistrettoPoint`s For these types we define operators with all combinations of borrowed and non-borrowed inputs, to avoid forcing API consumers to write extra ampersands. Since all of the operations involved with these types are expensive relative to the cost of an unnecessary copy, this isn't a big deal. The `MontgomeryPoint` struct isn't included in the above because it's only useful for scalar multiplication. This commit is based on work by @UnlawfulMonad. --- build.rs | 5 ++ src/edwards.rs | 23 ++++++++- src/lib.rs | 4 ++ src/macros.rs | 123 +++++++++++++++++++++++++++++++++++++++++++++++ src/ristretto.rs | 22 +++++++++ src/scalar.rs | 19 ++++++++ 6 files changed, 195 insertions(+), 1 deletion(-) create mode 100644 src/macros.rs diff --git a/build.rs b/build.rs index 2c30435..6b6e652 100644 --- a/build.rs +++ b/build.rs @@ -26,6 +26,11 @@ extern crate serde; #[cfg(feature = "yolocrypto")] extern crate stdsimd; +// Macros come first! +#[path="src/macros.rs"] +#[macro_use] +mod macros; + // Public modules #[path="src/scalar.rs"] diff --git a/src/edwards.rs b/src/edwards.rs index 5fca4e9..4af3ec5 100644 --- a/src/edwards.rs +++ b/src/edwards.rs @@ -378,12 +378,16 @@ impl<'a, 'b> Add<&'b ExtendedPoint> for &'a ExtendedPoint { } } +define_add_variants!(LHS = ExtendedPoint, RHS = ExtendedPoint, Output = ExtendedPoint); + impl<'b> AddAssign<&'b ExtendedPoint> for ExtendedPoint { fn add_assign(&mut self, _rhs: &'b ExtendedPoint) { *self = (self as &ExtendedPoint) + _rhs; } } +define_add_assign_variants!(LHS = ExtendedPoint, RHS = ExtendedPoint); + impl<'a, 'b> Sub<&'b ExtendedPoint> for &'a ExtendedPoint { type Output = ExtendedPoint; fn sub(self, other: &'b ExtendedPoint) -> ExtendedPoint { @@ -391,12 +395,16 @@ impl<'a, 'b> Sub<&'b ExtendedPoint> for &'a ExtendedPoint { } } +define_sub_variants!(LHS = ExtendedPoint, RHS = ExtendedPoint, Output = ExtendedPoint); + impl<'b> SubAssign<&'b ExtendedPoint> for ExtendedPoint { fn sub_assign(&mut self, _rhs: &'b ExtendedPoint) { *self = (self as &ExtendedPoint) - _rhs; } } +define_sub_assign_variants!(LHS = ExtendedPoint, RHS = ExtendedPoint); + // ------------------------------------------------------------------------ // Negation // ------------------------------------------------------------------------ @@ -414,6 +422,14 @@ impl<'a> Neg for &'a ExtendedPoint { } } +impl Neg for ExtendedPoint { + type Output = ExtendedPoint; + + fn neg(self) -> ExtendedPoint { + -&self + } +} + // ------------------------------------------------------------------------ // Scalar multiplication // ------------------------------------------------------------------------ @@ -425,6 +441,11 @@ impl<'b> MulAssign<&'b Scalar> for ExtendedPoint { } } +define_mul_assign_variants!(LHS = ExtendedPoint, RHS = Scalar); + +define_mul_variants!(LHS = ExtendedPoint, RHS = Scalar, Output = ExtendedPoint); +define_mul_variants!(LHS = Scalar, RHS = ExtendedPoint, Output = ExtendedPoint); + impl<'a, 'b> Mul<&'b Scalar> for &'a ExtendedPoint { type Output = ExtendedPoint; /// Scalar multiplication: compute `scalar * self`. @@ -478,7 +499,7 @@ impl<'a, 'b> Mul<&'b ExtendedPoint> for &'a Scalar { /// For scalar multiplication of a basepoint, /// `EdwardsBasepointTable` is approximately 4x faster. fn mul(self, point: &'b ExtendedPoint) -> ExtendedPoint { - point * &self + point * self } } diff --git a/src/lib.rs b/src/lib.rs index 6c4ed32..0baf368 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -68,6 +68,10 @@ extern crate serde; #[cfg(all(test, feature = "serde"))] extern crate serde_cbor; +// Internal macros. Must come first! +#[macro_use] +pub(crate) mod macros; + //------------------------------------------------------------------------ // curve25519-dalek public modules //------------------------------------------------------------------------ diff --git a/src/macros.rs b/src/macros.rs new file mode 100644 index 0000000..448d32c --- /dev/null +++ b/src/macros.rs @@ -0,0 +1,123 @@ +// -*- mode: rust; -*- +// +// This file is part of curve25519-dalek. +// Copyright (c) 2016-2017 Isis Lovecruft, Henry de Valence +// See LICENSE for licensing information. +// +// Authors: +// - Isis Agora Lovecruft +// - Henry de Valence + +//! Internal macros. + +/// Define borrow and non-borrow variants of `Add`. +macro_rules! define_add_variants { + (LHS = $lhs:ty, RHS = $rhs:ty, Output = $out:ty) => { + impl<'b> Add<&'b $rhs> for $lhs { + type Output = $out; + fn add(self, rhs: &'b $rhs) -> $out { + &self + rhs + } + } + + impl<'a> Add<$rhs> for &'a $lhs { + type Output = $out; + fn add(self, rhs: $rhs) -> $out { + self + &rhs + } + } + + impl Add<$rhs> for $lhs { + type Output = $out; + fn add(self, rhs: $rhs) -> $out { + &self + &rhs + } + } + } +} + +/// Define non-borrow variants of `AddAssign`. +macro_rules! define_add_assign_variants { + (LHS = $lhs:ty, RHS = $rhs:ty) => { + impl AddAssign<$rhs> for $lhs { + fn add_assign(&mut self, rhs: $rhs) { + *self += &rhs; + } + } + } +} + +/// Define borrow and non-borrow variants of `Sub`. +macro_rules! define_sub_variants { + (LHS = $lhs:ty, RHS = $rhs:ty, Output = $out:ty) => { + impl<'b> Sub<&'b $rhs> for $lhs { + type Output = $out; + fn sub(self, rhs: &'b $rhs) -> $out { + &self - rhs + } + } + + impl<'a> Sub<$rhs> for &'a $lhs { + type Output = $out; + fn sub(self, rhs: $rhs) -> $out { + self - &rhs + } + } + + impl Sub<$rhs> for $lhs { + type Output = $out; + fn sub(self, rhs: $rhs) -> $out { + &self - &rhs + } + } + } +} + +/// Define non-borrow variants of `SubAssign`. +macro_rules! define_sub_assign_variants { + (LHS = $lhs:ty, RHS = $rhs:ty) => { + impl SubAssign<$rhs> for $lhs { + fn sub_assign(&mut self, rhs: $rhs) { + *self -= &rhs; + } + } + } +} + +/// Define borrow and non-borrow variants of `Mul`. +macro_rules! define_mul_variants { + (LHS = $lhs:ty, RHS = $rhs:ty, Output = $out:ty) => { + impl<'b> Mul<&'b $rhs> for $lhs { + type Output = $out; + fn mul(self, rhs: &'b $rhs) -> $out { + &self * rhs + } + } + + impl<'a> Mul<$rhs> for &'a $lhs { + type Output = $out; + fn mul(self, rhs: $rhs) -> $out { + self * &rhs + } + } + + impl Mul<$rhs> for $lhs { + type Output = $out; + fn mul(self, rhs: $rhs) -> $out { + &self * &rhs + } + } + } +} + +/// Define non-borrow variants of `MulAssign`. +macro_rules! define_mul_assign_variants { + (LHS = $lhs:ty, RHS = $rhs:ty) => { + impl MulAssign<$rhs> for $lhs { + fn mul_assign(&mut self, rhs: $rhs) { + *self *= &rhs; + } + } + } +} + diff --git a/src/ristretto.rs b/src/ristretto.rs index ee6260e..40d667b 100644 --- a/src/ristretto.rs +++ b/src/ristretto.rs @@ -925,12 +925,16 @@ impl<'a, 'b> Add<&'b RistrettoPoint> for &'a RistrettoPoint { } } +define_add_variants!(LHS = RistrettoPoint, RHS = RistrettoPoint, Output = RistrettoPoint); + impl<'b> AddAssign<&'b RistrettoPoint> for RistrettoPoint { fn add_assign(&mut self, _rhs: &RistrettoPoint) { *self = (self as &RistrettoPoint) + _rhs; } } +define_add_assign_variants!(LHS = RistrettoPoint, RHS = RistrettoPoint); + impl<'a, 'b> Sub<&'b RistrettoPoint> for &'a RistrettoPoint { type Output = RistrettoPoint; @@ -939,12 +943,16 @@ impl<'a, 'b> Sub<&'b RistrettoPoint> for &'a RistrettoPoint { } } +define_sub_variants!(LHS = RistrettoPoint, RHS = RistrettoPoint, Output = RistrettoPoint); + impl<'b> SubAssign<&'b RistrettoPoint> for RistrettoPoint { fn sub_assign(&mut self, _rhs: &RistrettoPoint) { *self = (self as &RistrettoPoint) - _rhs; } } +define_sub_assign_variants!(LHS = RistrettoPoint, RHS = RistrettoPoint); + impl<'a> Neg for &'a RistrettoPoint { type Output = RistrettoPoint; @@ -953,6 +961,14 @@ impl<'a> Neg for &'a RistrettoPoint { } } +impl Neg for RistrettoPoint { + type Output = RistrettoPoint; + + fn neg(self) -> RistrettoPoint { + -&self + } +} + impl<'b> MulAssign<&'b Scalar> for RistrettoPoint { fn mul_assign(&mut self, scalar: &'b Scalar) { let result = (self as &RistrettoPoint) * scalar; @@ -977,6 +993,12 @@ impl<'a, 'b> Mul<&'b RistrettoPoint> for &'a Scalar { } } +define_mul_assign_variants!(LHS = RistrettoPoint, RHS = Scalar); + +define_mul_variants!(LHS = RistrettoPoint, RHS = Scalar, Output = RistrettoPoint); +define_mul_variants!(LHS = Scalar, RHS = RistrettoPoint, Output = RistrettoPoint); + + /// Given a vector of (possibly secret) scalars and a vector of /// (possibly secret) points, compute `c_1 P_1 + ... + c_n P_n`. /// diff --git a/src/scalar.rs b/src/scalar.rs index 8d12b4e..64a59bc 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -186,6 +186,8 @@ impl<'b> MulAssign<&'b Scalar> for Scalar { } } +define_mul_assign_variants!(LHS = Scalar, RHS = Scalar); + impl<'a, 'b> Mul<&'b Scalar> for &'a Scalar { type Output = Scalar; fn mul(self, _rhs: &'b Scalar) -> Scalar { @@ -193,12 +195,16 @@ impl<'a, 'b> Mul<&'b Scalar> for &'a Scalar { } } +define_mul_variants!(LHS = Scalar, RHS = Scalar, Output = Scalar); + impl<'b> AddAssign<&'b Scalar> for Scalar { fn add_assign(&mut self, _rhs: &'b Scalar) { *self = UnpackedScalar::add(&self.unpack(), &_rhs.unpack()).pack(); } } +define_add_assign_variants!(LHS = Scalar, RHS = Scalar); + impl<'a, 'b> Add<&'b Scalar> for &'a Scalar { type Output = Scalar; fn add(self, _rhs: &'b Scalar) -> Scalar { @@ -206,12 +212,16 @@ impl<'a, 'b> Add<&'b Scalar> for &'a Scalar { } } +define_add_variants!(LHS = Scalar, RHS = Scalar, Output = Scalar); + impl<'b> SubAssign<&'b Scalar> for Scalar { fn sub_assign(&mut self, _rhs: &'b Scalar) { *self = UnpackedScalar::sub(&self.unpack(), &_rhs.unpack()).pack(); } } +define_sub_assign_variants!(LHS = Scalar, RHS = Scalar); + impl<'a, 'b> Sub<&'b Scalar> for &'a Scalar { type Output = Scalar; fn sub(self, _rhs: &'b Scalar) -> Scalar { @@ -219,6 +229,8 @@ impl<'a, 'b> Sub<&'b Scalar> for &'a Scalar { } } +define_sub_variants!(LHS = Scalar, RHS = Scalar, Output = Scalar); + impl<'a> Neg for &'a Scalar { type Output = Scalar; fn neg(self) -> Scalar { @@ -226,6 +238,13 @@ impl<'a> Neg for &'a Scalar { } } +impl<'a> Neg for Scalar { + type Output = Scalar; + fn neg(self) -> Scalar { + -&self + } +} + impl ConditionallyAssignable for Scalar { /// Conditionally assign another Scalar to this one. /// From 86fd06db0050a177a3248ebd54553e637ce23fc6 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Wed, 24 Jan 2018 12:42:57 -0800 Subject: [PATCH 10/13] 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 { From 7f3e82bab5dc9a53580fad719ab749c4b8b7855a Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Wed, 24 Jan 2018 21:25:03 +0000 Subject: [PATCH 11/13] Revert "Bump subtle dependency to 0.5.0." This reverts commit 460355d5e06d590ab289ed063045793734cc4bd4. --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index 68e3127..f4a60ee 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -40,7 +40,7 @@ version = "0.4" version = "0.7" [dependencies.subtle] -version = "^0.5" +version = "^0.3" default-features = false [dependencies.clear_on_drop] From 53908798a7b4fb1fd87ee22a5a791faf7c108dd1 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Wed, 24 Jan 2018 22:20:53 +0000 Subject: [PATCH 12/13] Remove Travis notifications from #curve25519-dalek channel. --- .travis.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 4a28c7e..a598e32 100644 --- a/.travis.yml +++ b/.travis.yml @@ -52,4 +52,3 @@ notifications: slack: rooms: - dalek-cryptography:Xxv9WotKYWdSoKlgKNqXiHoD#dalek-bots - - dalek-cryptography:Xxv9WotKYWdSoKlgKNqXiHoD#curve25519-dalek From 3b7399a06f6884e7a50beef514837a3c9213f476 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 25 Jan 2018 02:58:00 +0000 Subject: [PATCH 13/13] Bump curve25519-dalek version to 0.14.4. --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index f4a60ee..05ba1c4 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "curve25519-dalek" -version = "0.14.3" +version = "0.14.4" authors = ["Isis Lovecruft ", "Henry de Valence "] readme = "README.md"