inverse elligator works, but at what cost?!?!!?

aux/ristretto/ristretto.sage

@@ -53,6 +53,8 @@ def isqrt(x,exn=InvalidEncodingException("Not on curve")):
     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)"""
     if x==0: return True,0
@@ -325,8 +327,8 @@ class Decaf_1_1_Point(QuotientEdwardsPoint):
             den = x*y
             isr = isqrt(num*(a-d)*den^2)
     
-            iden = isr * den * self.isoMagic
-            inum = isr * num
+            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
             
             if negative(iden*inum*self.i*t^2*(d-a)) != toggle_rotation:
                 iden,inum = inum,iden
@@ -386,18 +388,29 @@ class Decaf_1_1_Point(QuotientEdwardsPoint):
                     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
+                        #print
+                        #print toggle_rotation,toggle_altx,toggle_s
+                        #print m1
+                        #print m12
                     
                     
-                        if self == self.__class__() and self.cofactor == 4:
-                            # Hacks for identity!
-                            if toggle_altx: m12 = 1
-                            elif toggle_s: m1 = 1
-                            elif toggle_r: continue
-                            ## BOTH???
+                        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)
@@ -406,7 +419,7 @@ class Decaf_1_1_Point(QuotientEdwardsPoint):
                         ok,sr = isqrt_i(rnum*rden*self.qnr)
                         if not ok: continue
                         sr *= rnum
-                        print "Works! %d %x" % (swap,sr)
+                        #print "Works! %d %x" % (swap,sr)
                     
                         if negative(sr) != toggle_r: sr = -sr
                         ret = self.gfToBytes(sr)
@@ -681,12 +694,12 @@ def test(cls,n):
         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)
+#test(Ed25519Point,100)
+#test(NegEd25519Point,100)
+#test(IsoEd25519Point,100)
+#test(IsoEd448Point,100)
+#test(TwistedEd448GoldilocksPoint,100)
+#test(Ed448GoldilocksPoint,100)
         
    
 def testElligator(cls,n):
@@ -710,12 +723,12 @@ def testElligator(cls,n):
             pass # TODO
         
 
-testElligator(Ed25519Point,100)
-testElligator(NegEd25519Point,100)
-testElligator(IsoEd25519Point,100)
-testElligator(IsoEd448Point,100)
-testElligator(Ed448GoldilocksPoint,100)
-testElligator(TwistedEd448GoldilocksPoint,100)
+#testElligator(Ed25519Point,100)
+#testElligator(NegEd25519Point,100)
+#testElligator(IsoEd25519Point,100)
+#testElligator(IsoEd448Point,100)
+#testElligator(Ed448GoldilocksPoint,100)
+#testElligator(TwistedEd448GoldilocksPoint,100)
 
 def gangtest(classes,n):
     print "Gang test",[cls.__name__ for cls in classes]
@@ -743,5 +756,5 @@ 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)
+#gangtest([IsoEd448Point,TwistedEd448GoldilocksPoint,Ed448GoldilocksPoint],100)
+#gangtest([Ed25519Point,IsoEd25519Point],100)

src/GENERATED/c/curve25519/decaf.c

@@ -48,7 +48,8 @@ static const scalar_t point_scalarmul_adjustment = {{{
 
 const uint8_t decaf_x25519_base_point[DECAF_X25519_PUBLIC_BYTES] = { 0x09 };
 
-static const gf RISTRETTO_ISOMAGIC = {{{
+#define RISTRETTO_FACTOR DECAF_255_RISTRETTO_FACTOR
+const gf RISTRETTO_FACTOR = {{{
     0x0fdaa805d40ea, 0x2eb482e57d339, 0x007610274bc58, 0x6510b613dc8ff, 0x786c8905cfaff
 }}};
 
@@ -157,7 +158,7 @@ void API_NS(deisogenize) (
     gf_mulw(t2,t1,-1-TWISTED_D); /* -x^2 * (a-d) * num */
     gf_isr(t1,t2);    /* t1 = isr */
     gf_mul(t2,t1,t3); /* t2 = ratio */
-    gf_mul(t4,t2,RISTRETTO_ISOMAGIC);
+    gf_mul(t4,t2,RISTRETTO_FACTOR);
     mask_t negx = gf_lobit(t4) ^ toggle_altx;
     gf_cond_neg(t2, negx);
     gf_mul(t3,t2,p->z);
@@ -183,12 +184,12 @@ void API_NS(deisogenize) (
     gf_mulw(t1,t4,-1-TWISTED_D);
     gf_isr(t4,t1);         /* isqrt(num*(a-d)*den^2) */
     gf_mul(t1,t2,t4);
-    gf_mul(t2,t1,RISTRETTO_ISOMAGIC); /* t2 = "iden" in ristretto.sage */
+    gf_mul(t2,t1,RISTRETTO_FACTOR); /* t2 = "iden" in ristretto.sage */
     gf_mul(t1,t3,t4);                 /* t1 = "inum" in ristretto.sage */
 
     /* Calculate altxy = iden*inum*i*t^2*(d-a) */
     gf_mul(t3,t1,t2);
-    gf_mul_qnr(t4,t3);
+    gf_mul_i(t4,t3);
     gf_mul(t3,t4,p->t);
     gf_mul(t4,t3,p->t);
     gf_mulw(t3,t4,TWISTED_D+1);      /* iden*inum*i*t^2*(d-1) */
@@ -196,10 +197,10 @@ void API_NS(deisogenize) (
     
     /* Rotate if altxy is negative */
     gf_cond_swap(t1,t2,rotate);
-    gf_mul_qnr(t4,p->x);
+    gf_mul_i(t4,p->x);
     gf_cond_sel(t4,p->y,t4,rotate);  /* t4 = "fac" = ix if rotate, else y */
     
-    gf_mul_qnr(t5,RISTRETTO_ISOMAGIC); /* t5 = imi */
+    gf_mul_i(t5,RISTRETTO_FACTOR); /* t5 = imi */
     gf_mul(t3,t5,t2);                /* iden * imi */
     gf_mul(t2,t5,t1);
     gf_mul(t5,t2,p->t);              /* "altx" = iden*imi*t */
@@ -258,20 +259,20 @@ decaf_error_t API_NS(point_decode) (
     gf_add(tmp2,tmp2,tmp2);        /* 2*s*isr*den */
     gf_mul(tmp,tmp2,isr);          /* 2*s*isr^2*den */
     gf_mul(p->x,tmp,num);          /* 2*s*isr^2*den*num */
-    gf_mul(tmp,tmp2,RISTRETTO_ISOMAGIC); /* 2*s*isr*den*magic */
+    gf_mul(tmp,tmp2,RISTRETTO_FACTOR); /* 2*s*isr*den*magic */
     gf_cond_neg(p->x,gf_lobit(tmp)); /* flip x */
     
 #if COFACTOR==8
     /* Additionally check y != 0 and x*y*isomagic nonegative */
     succ &= ~gf_eq(p->y,ZERO);
     gf_mul(tmp,p->x,p->y);
-    gf_mul(tmp2,tmp,RISTRETTO_ISOMAGIC);
+    gf_mul(tmp2,tmp,RISTRETTO_FACTOR);
     succ &= ~gf_lobit(tmp2);
 #endif
 
 #if IMAGINE_TWIST
     gf_copy(tmp,p->x);
-    gf_mul_qnr(p->x,tmp);
+    gf_mul_i(p->x,tmp);
 #endif
 
     /* Fill in z and t */
@@ -1077,9 +1078,9 @@ void API_NS(point_mul_by_cofactor_and_encode_like_eddsa) (
         gf_mul ( y, u, t ); // (x^2+y^2)(2z^2-y^2+x^2)
         gf_mul ( u, z, t );
         gf_copy( z, u );
-        gf_mul ( u, x, RISTRETTO_ISOMAGIC );
+        gf_mul ( u, x, RISTRETTO_FACTOR );
 #if IMAGINE_TWIST
-        gf_mul_qnr( x, u );
+        gf_mul_i( x, u );
 #else
 #error "... probably wrong"
         gf_copy( x, u );
@@ -1090,7 +1091,7 @@ void API_NS(point_mul_by_cofactor_and_encode_like_eddsa) (
     {
         API_NS(point_double)(q,q);
         API_NS(point_double)(q,q);
-        gf_mul_qnr(x, q->x);
+        gf_mul_i(x, q->x);
         gf_copy(y, q->y);
         gf_copy(z, q->z);
     }
@@ -1188,8 +1189,8 @@ decaf_error_t API_NS(point_decode_like_eddsa_and_ignore_cofactor) (
         gf_sqr ( p->x, p->z );
         gf_add ( p->z, p->x, p->x );
         gf_sub ( c, p->z, p->t ); // 2z^2 - y^2 + x^2
-        gf_div_qnr ( a, c );
-        gf_mul ( c, a, RISTRETTO_ISOMAGIC );
+        gf_div_i ( a, c );
+        gf_mul ( c, a, RISTRETTO_FACTOR );
         gf_mul ( p->x, b, p->t); // (2xy)(y^2-x^2)
         gf_mul ( p->z, p->t, c ); // (y^2-x^2)sd(2z^2 - y^2 + x^2)
         gf_mul ( p->y, d, c ); // (y^2+x^2)sd(2z^2 - y^2 + x^2)

src/GENERATED/c/curve25519/elligator.c

@@ -21,6 +21,10 @@
 #define IMAGINE_TWIST 1
 #define COFACTOR 8
 static const int EDWARDS_D = -121665;
+
+#define RISTRETTO_FACTOR DECAF_255_RISTRETTO_FACTOR
+extern const gf RISTRETTO_FACTOR;
+
 /* End of template stuff */
 
 extern mask_t API_NS(deisogenize) (
@@ -145,8 +149,23 @@ API_NS(invert_elligator_nonuniform) (
     API_NS(deisogenize)(a,b,c,p,sgn_s,sgn_altx,sgn_ed_T);
     
     mask_t is_identity = gf_eq(p->t,ZERO);
+#if COFACTOR==4
     gf_cond_sel(b,b,ONE,is_identity & sgn_altx);
     gf_cond_sel(c,c,ONE,is_identity & sgn_s &~ sgn_altx);
+#elif IMAGINE_TWIST
+    /* Terrible, terrible special casing due to lots of 0/0 is deisogenize
+     * Basically we need to generate -D and +- i*RISTRETTO_FACTOR
+     */
+    gf_mul_i(a,RISTRETTO_FACTOR);
+    gf_cond_sel(b,b,ONE,is_identity);
+    gf_cond_neg(a,sgn_altx);
+    gf_cond_sel(c,c,a,is_identity & sgn_ed_T);
+    gf_cond_sel(c,c,ZERO,is_identity & ~sgn_ed_T);
+    gf_mulw(a,ONE,-EDWARDS_D);
+    gf_cond_sel(c,c,a,is_identity & ~sgn_ed_T &~ sgn_altx);
+#else
+#error "Different special-casing goes here!"
+#endif
     
 #if IMAGINE_TWIST
     gf_mulw(a,b,-EDWARDS_D);
@@ -169,7 +188,8 @@ API_NS(invert_elligator_nonuniform) (
 #endif
     
     gf_cond_neg(b, sgn_r0^gf_lobit(b));
-    succ &= ~(gf_eq(b,ZERO) & sgn_r0);
+    /* Eliminate duplicate values for identity ... */
+    succ &= ~(gf_eq(b,ZERO) & (sgn_r0 | sgn_s));
     // #if COFACTOR == 8
     //     succ &= ~(is_identity & sgn_ed_T); /* NB: there are no preimages of rotated identity. */
     // #endif

src/GENERATED/c/ed448goldilocks/decaf.c

@@ -48,7 +48,8 @@ static const scalar_t point_scalarmul_adjustment = {{{
 
 const uint8_t decaf_x448_base_point[DECAF_X448_PUBLIC_BYTES] = { 0x05 };
 
-static const gf RISTRETTO_ISOMAGIC = {{{
+#define RISTRETTO_FACTOR DECAF_448_RISTRETTO_FACTOR
+const gf RISTRETTO_FACTOR = {{{
     0x42ef0f45572736, 0x7bf6aa20ce5296, 0xf4fd6eded26033, 0x968c14ba839a66, 0xb8d54b64a2d780, 0x6aa0a1f1a7b8a5, 0x683bf68d722fa2, 0x22d962fbeb24f7
 }}};
 
@@ -157,7 +158,7 @@ void API_NS(deisogenize) (
     gf_mulw(t2,t1,-1-TWISTED_D); /* -x^2 * (a-d) * num */
     gf_isr(t1,t2);    /* t1 = isr */
     gf_mul(t2,t1,t3); /* t2 = ratio */
-    gf_mul(t4,t2,RISTRETTO_ISOMAGIC);
+    gf_mul(t4,t2,RISTRETTO_FACTOR);
     mask_t negx = gf_lobit(t4) ^ toggle_altx;
     gf_cond_neg(t2, negx);
     gf_mul(t3,t2,p->z);
@@ -183,12 +184,12 @@ void API_NS(deisogenize) (
     gf_mulw(t1,t4,-1-TWISTED_D);
     gf_isr(t4,t1);         /* isqrt(num*(a-d)*den^2) */
     gf_mul(t1,t2,t4);
-    gf_mul(t2,t1,RISTRETTO_ISOMAGIC); /* t2 = "iden" in ristretto.sage */
+    gf_mul(t2,t1,RISTRETTO_FACTOR); /* t2 = "iden" in ristretto.sage */
     gf_mul(t1,t3,t4);                 /* t1 = "inum" in ristretto.sage */
 
     /* Calculate altxy = iden*inum*i*t^2*(d-a) */
     gf_mul(t3,t1,t2);
-    gf_mul_qnr(t4,t3);
+    gf_mul_i(t4,t3);
     gf_mul(t3,t4,p->t);
     gf_mul(t4,t3,p->t);
     gf_mulw(t3,t4,TWISTED_D+1);      /* iden*inum*i*t^2*(d-1) */
@@ -196,10 +197,10 @@ void API_NS(deisogenize) (
     
     /* Rotate if altxy is negative */
     gf_cond_swap(t1,t2,rotate);
-    gf_mul_qnr(t4,p->x);
+    gf_mul_i(t4,p->x);
     gf_cond_sel(t4,p->y,t4,rotate);  /* t4 = "fac" = ix if rotate, else y */
     
-    gf_mul_qnr(t5,RISTRETTO_ISOMAGIC); /* t5 = imi */
+    gf_mul_i(t5,RISTRETTO_FACTOR); /* t5 = imi */
     gf_mul(t3,t5,t2);                /* iden * imi */
     gf_mul(t2,t5,t1);
     gf_mul(t5,t2,p->t);              /* "altx" = iden*imi*t */
@@ -258,20 +259,20 @@ decaf_error_t API_NS(point_decode) (
     gf_add(tmp2,tmp2,tmp2);        /* 2*s*isr*den */
     gf_mul(tmp,tmp2,isr);          /* 2*s*isr^2*den */
     gf_mul(p->x,tmp,num);          /* 2*s*isr^2*den*num */
-    gf_mul(tmp,tmp2,RISTRETTO_ISOMAGIC); /* 2*s*isr*den*magic */
+    gf_mul(tmp,tmp2,RISTRETTO_FACTOR); /* 2*s*isr*den*magic */
     gf_cond_neg(p->x,gf_lobit(tmp)); /* flip x */
     
 #if COFACTOR==8
     /* Additionally check y != 0 and x*y*isomagic nonegative */
     succ &= ~gf_eq(p->y,ZERO);
     gf_mul(tmp,p->x,p->y);
-    gf_mul(tmp2,tmp,RISTRETTO_ISOMAGIC);
+    gf_mul(tmp2,tmp,RISTRETTO_FACTOR);
     succ &= ~gf_lobit(tmp2);
 #endif
 
 #if IMAGINE_TWIST
     gf_copy(tmp,p->x);
-    gf_mul_qnr(p->x,tmp);
+    gf_mul_i(p->x,tmp);
 #endif
 
     /* Fill in z and t */
@@ -1077,9 +1078,9 @@ void API_NS(point_mul_by_cofactor_and_encode_like_eddsa) (
         gf_mul ( y, u, t ); // (x^2+y^2)(2z^2-y^2+x^2)
         gf_mul ( u, z, t );
         gf_copy( z, u );
-        gf_mul ( u, x, RISTRETTO_ISOMAGIC );
+        gf_mul ( u, x, RISTRETTO_FACTOR );
 #if IMAGINE_TWIST
-        gf_mul_qnr( x, u );
+        gf_mul_i( x, u );
 #else
 #error "... probably wrong"
         gf_copy( x, u );
@@ -1090,7 +1091,7 @@ void API_NS(point_mul_by_cofactor_and_encode_like_eddsa) (
     {
         API_NS(point_double)(q,q);
         API_NS(point_double)(q,q);
-        gf_mul_qnr(x, q->x);
+        gf_mul_i(x, q->x);
         gf_copy(y, q->y);
         gf_copy(z, q->z);
     }
@@ -1188,8 +1189,8 @@ decaf_error_t API_NS(point_decode_like_eddsa_and_ignore_cofactor) (
         gf_sqr ( p->x, p->z );
         gf_add ( p->z, p->x, p->x );
         gf_sub ( c, p->z, p->t ); // 2z^2 - y^2 + x^2
-        gf_div_qnr ( a, c );
-        gf_mul ( c, a, RISTRETTO_ISOMAGIC );
+        gf_div_i ( a, c );
+        gf_mul ( c, a, RISTRETTO_FACTOR );
         gf_mul ( p->x, b, p->t); // (2xy)(y^2-x^2)
         gf_mul ( p->z, p->t, c ); // (y^2-x^2)sd(2z^2 - y^2 + x^2)
         gf_mul ( p->y, d, c ); // (y^2+x^2)sd(2z^2 - y^2 + x^2)

src/GENERATED/c/ed448goldilocks/elligator.c

@@ -21,6 +21,10 @@
 #define IMAGINE_TWIST 0
 #define COFACTOR 4
 static const int EDWARDS_D = -39081;
+
+#define RISTRETTO_FACTOR DECAF_448_RISTRETTO_FACTOR
+extern const gf RISTRETTO_FACTOR;
+
 /* End of template stuff */
 
 extern mask_t API_NS(deisogenize) (
@@ -145,8 +149,23 @@ API_NS(invert_elligator_nonuniform) (
     API_NS(deisogenize)(a,b,c,p,sgn_s,sgn_altx,sgn_ed_T);
     
     mask_t is_identity = gf_eq(p->t,ZERO);
+#if COFACTOR==4
     gf_cond_sel(b,b,ONE,is_identity & sgn_altx);
     gf_cond_sel(c,c,ONE,is_identity & sgn_s &~ sgn_altx);
+#elif IMAGINE_TWIST
+    /* Terrible, terrible special casing due to lots of 0/0 is deisogenize
+     * Basically we need to generate -D and +- i*RISTRETTO_FACTOR
+     */
+    gf_mul_i(a,RISTRETTO_FACTOR);
+    gf_cond_sel(b,b,ONE,is_identity);
+    gf_cond_neg(a,sgn_altx);
+    gf_cond_sel(c,c,a,is_identity & sgn_ed_T);
+    gf_cond_sel(c,c,ZERO,is_identity & ~sgn_ed_T);
+    gf_mulw(a,ONE,-EDWARDS_D);
+    gf_cond_sel(c,c,a,is_identity & ~sgn_ed_T &~ sgn_altx);
+#else
+#error "Different special-casing goes here!"
+#endif
     
 #if IMAGINE_TWIST
     gf_mulw(a,b,-EDWARDS_D);
@@ -169,7 +188,8 @@ API_NS(invert_elligator_nonuniform) (
 #endif
     
     gf_cond_neg(b, sgn_r0^gf_lobit(b));
-    succ &= ~(gf_eq(b,ZERO) & sgn_r0);
+    /* Eliminate duplicate values for identity ... */
+    succ &= ~(gf_eq(b,ZERO) & (sgn_r0 | sgn_s));
     // #if COFACTOR == 8
     //     succ &= ~(is_identity & sgn_ed_T); /* NB: there are no preimages of rotated identity. */
     // #endif

src/include/field.h

@@ -103,5 +103,10 @@ static DECAF_INLINE void gf_div_qnr(gf_s *__restrict__ out, const gf x) {
 #endif
 }
 
+#if P_MOD_8 == 5
+#define gf_mul_i gf_mul_qnr
+#define gf_div_i gf_div_qnr
+#endif
+
 
 #endif // __GF_H__

src/per_curve/decaf.tmpl.c

@@ -37,7 +37,8 @@ static const scalar_t point_scalarmul_adjustment = {{{
 
 const uint8_t decaf_x$(gf_shortname)_base_point[DECAF_X$(gf_shortname)_PUBLIC_BYTES] = { $(ser(mont_base,8)) };
 
-static const gf RISTRETTO_ISOMAGIC = {{{
+#define RISTRETTO_FACTOR $(C_NS)_RISTRETTO_FACTOR
+const gf RISTRETTO_FACTOR = {{{
     $(ser(msqrt(d-1 if imagine_twist else -d,modulus,lo_bit_clear=True),gf_lit_limb_bits))
 }}};
 
@@ -146,7 +147,7 @@ void API_NS(deisogenize) (
     gf_mulw(t2,t1,-1-TWISTED_D); /* -x^2 * (a-d) * num */
     gf_isr(t1,t2);    /* t1 = isr */
     gf_mul(t2,t1,t3); /* t2 = ratio */
-    gf_mul(t4,t2,RISTRETTO_ISOMAGIC);
+    gf_mul(t4,t2,RISTRETTO_FACTOR);
     mask_t negx = gf_lobit(t4) ^ toggle_altx;
     gf_cond_neg(t2, negx);
     gf_mul(t3,t2,p->z);
@@ -172,12 +173,12 @@ void API_NS(deisogenize) (
     gf_mulw(t1,t4,-1-TWISTED_D);
     gf_isr(t4,t1);         /* isqrt(num*(a-d)*den^2) */
     gf_mul(t1,t2,t4);
-    gf_mul(t2,t1,RISTRETTO_ISOMAGIC); /* t2 = "iden" in ristretto.sage */
+    gf_mul(t2,t1,RISTRETTO_FACTOR); /* t2 = "iden" in ristretto.sage */
     gf_mul(t1,t3,t4);                 /* t1 = "inum" in ristretto.sage */
 
     /* Calculate altxy = iden*inum*i*t^2*(d-a) */
     gf_mul(t3,t1,t2);
-    gf_mul_qnr(t4,t3);
+    gf_mul_i(t4,t3);
     gf_mul(t3,t4,p->t);
     gf_mul(t4,t3,p->t);
     gf_mulw(t3,t4,TWISTED_D+1);      /* iden*inum*i*t^2*(d-1) */
@@ -185,10 +186,10 @@ void API_NS(deisogenize) (
     
     /* Rotate if altxy is negative */
     gf_cond_swap(t1,t2,rotate);
-    gf_mul_qnr(t4,p->x);
+    gf_mul_i(t4,p->x);
     gf_cond_sel(t4,p->y,t4,rotate);  /* t4 = "fac" = ix if rotate, else y */
     
-    gf_mul_qnr(t5,RISTRETTO_ISOMAGIC); /* t5 = imi */
+    gf_mul_i(t5,RISTRETTO_FACTOR); /* t5 = imi */
     gf_mul(t3,t5,t2);                /* iden * imi */
     gf_mul(t2,t5,t1);
     gf_mul(t5,t2,p->t);              /* "altx" = iden*imi*t */
@@ -247,20 +248,20 @@ decaf_error_t API_NS(point_decode) (
     gf_add(tmp2,tmp2,tmp2);        /* 2*s*isr*den */
     gf_mul(tmp,tmp2,isr);          /* 2*s*isr^2*den */
     gf_mul(p->x,tmp,num);          /* 2*s*isr^2*den*num */
-    gf_mul(tmp,tmp2,RISTRETTO_ISOMAGIC); /* 2*s*isr*den*magic */
+    gf_mul(tmp,tmp2,RISTRETTO_FACTOR); /* 2*s*isr*den*magic */
     gf_cond_neg(p->x,gf_lobit(tmp)); /* flip x */
     
 #if COFACTOR==8
     /* Additionally check y != 0 and x*y*isomagic nonegative */
     succ &= ~gf_eq(p->y,ZERO);
     gf_mul(tmp,p->x,p->y);
-    gf_mul(tmp2,tmp,RISTRETTO_ISOMAGIC);
+    gf_mul(tmp2,tmp,RISTRETTO_FACTOR);
     succ &= ~gf_lobit(tmp2);
 #endif
 
 #if IMAGINE_TWIST
     gf_copy(tmp,p->x);
-    gf_mul_qnr(p->x,tmp);
+    gf_mul_i(p->x,tmp);
 #endif
 
     /* Fill in z and t */
@@ -1066,9 +1067,9 @@ void API_NS(point_mul_by_cofactor_and_encode_like_eddsa) (
         gf_mul ( y, u, t ); // (x^2+y^2)(2z^2-y^2+x^2)
         gf_mul ( u, z, t );
         gf_copy( z, u );
-        gf_mul ( u, x, RISTRETTO_ISOMAGIC );
+        gf_mul ( u, x, RISTRETTO_FACTOR );
 #if IMAGINE_TWIST
-        gf_mul_qnr( x, u );
+        gf_mul_i( x, u );
 #else
 #error "... probably wrong"
         gf_copy( x, u );
@@ -1079,7 +1080,7 @@ void API_NS(point_mul_by_cofactor_and_encode_like_eddsa) (
     {
         API_NS(point_double)(q,q);
         API_NS(point_double)(q,q);
-        gf_mul_qnr(x, q->x);
+        gf_mul_i(x, q->x);
         gf_copy(y, q->y);
         gf_copy(z, q->z);
     }
@@ -1177,8 +1178,8 @@ decaf_error_t API_NS(point_decode_like_eddsa_and_ignore_cofactor) (
         gf_sqr ( p->x, p->z );
         gf_add ( p->z, p->x, p->x );
         gf_sub ( c, p->z, p->t ); // 2z^2 - y^2 + x^2
-        gf_div_qnr ( a, c );
-        gf_mul ( c, a, RISTRETTO_ISOMAGIC );
+        gf_div_i ( a, c );
+        gf_mul ( c, a, RISTRETTO_FACTOR );
         gf_mul ( p->x, b, p->t); // (2xy)(y^2-x^2)
         gf_mul ( p->z, p->t, c ); // (y^2-x^2)sd(2z^2 - y^2 + x^2)
         gf_mul ( p->y, d, c ); // (y^2+x^2)sd(2z^2 - y^2 + x^2)

src/per_curve/elligator.tmpl.c

@@ -10,6 +10,10 @@
 #define IMAGINE_TWIST $(imagine_twist)
 #define COFACTOR $(cofactor)
 static const int EDWARDS_D = $(d);
+
+#define RISTRETTO_FACTOR $(C_NS)_RISTRETTO_FACTOR
+extern const gf RISTRETTO_FACTOR;
+
 /* End of template stuff */
 
 extern mask_t API_NS(deisogenize) (
@@ -134,8 +138,23 @@ API_NS(invert_elligator_nonuniform) (
     API_NS(deisogenize)(a,b,c,p,sgn_s,sgn_altx,sgn_ed_T);
     
     mask_t is_identity = gf_eq(p->t,ZERO);
+#if COFACTOR==4
     gf_cond_sel(b,b,ONE,is_identity & sgn_altx);
     gf_cond_sel(c,c,ONE,is_identity & sgn_s &~ sgn_altx);
+#elif IMAGINE_TWIST
+    /* Terrible, terrible special casing due to lots of 0/0 is deisogenize
+     * Basically we need to generate -D and +- i*RISTRETTO_FACTOR
+     */
+    gf_mul_i(a,RISTRETTO_FACTOR);
+    gf_cond_sel(b,b,ONE,is_identity);
+    gf_cond_neg(a,sgn_altx);
+    gf_cond_sel(c,c,a,is_identity & sgn_ed_T);
+    gf_cond_sel(c,c,ZERO,is_identity & ~sgn_ed_T);
+    gf_mulw(a,ONE,-EDWARDS_D);
+    gf_cond_sel(c,c,a,is_identity & ~sgn_ed_T &~ sgn_altx);
+#else
+#error "Different special-casing goes here!"
+#endif
     
 #if IMAGINE_TWIST
     gf_mulw(a,b,-EDWARDS_D);
@@ -158,7 +177,8 @@ API_NS(invert_elligator_nonuniform) (
 #endif
     
     gf_cond_neg(b, sgn_r0^gf_lobit(b));
-    succ &= ~(gf_eq(b,ZERO) & sgn_r0);
+    /* Eliminate duplicate values for identity ... */
+    succ &= ~(gf_eq(b,ZERO) & (sgn_r0 | sgn_s));
     // #if COFACTOR == 8
     //     succ &= ~(is_identity & sgn_ed_T); /* NB: there are no preimages of rotated identity. */
     // #endif