diff --git a/CHANGELOG.md b/CHANGELOG.md
index 98cd436..4f825de 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,20 @@
 
 Entries are listed in reverse chronological order.
 
+## 1.2.3
+
+* Fix an issue identified by a Quarkslab audit (and Jack Grigg), where manually
+  constructing unreduced `Scalar` values, as needed for X/Ed25519, and then
+  performing scalar/scalar arithmetic could compute incorrect results.
+* Switch to upstream Rust intrinsics for the IFMA backend now that they exist in
+  Rust and don't need to be defined locally.
+* Ensure that the NAF computation works correctly, even for parameters never
+  used elsewhere in the codebase.
+* Minor refactoring to EdwardsPoint decompression.
+* Fix broken links in documentation.
+* Fix compilation on nightly broken due to changes to the `#[doc(include)]` path
+  root (not quite correctly done in 1.2.2).
+
 ## 1.2.2
 
 * Fix a typo in an internal doc-comment.
diff --git a/Cargo.toml b/Cargo.toml
index 6243498..864e7df 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "curve25519-dalek"
-version = "1.2.2"
+version = "1.2.3"
 authors = ["Isis Lovecruft <isis@patternsinthevoid.net>",
            "Henry de Valence <hdevalence@hdevalence.ca>"]
 readme = "README.md"
diff --git a/build.rs b/build.rs
index 2d020e1..0ffa5fc 100644
--- a/build.rs
+++ b/build.rs
@@ -1,13 +1,6 @@
-#![cfg_attr(
-    all(feature = "simd_backend", target_feature = "avx512ifma"),
-    feature(simd_ffi)
-)]
-#![cfg_attr(
-    all(feature = "simd_backend", target_feature = "avx512ifma"),
-    feature(link_llvm_intrinsics)
-)]
 #![cfg_attr(all(feature = "alloc", not(feature = "std")), feature(alloc))]
 #![cfg_attr(feature = "nightly", feature(doc_cfg))]
+#![cfg_attr(feature = "simd_backend", feature(stdsimd))]
 #![allow(unused_variables)]
 #![allow(non_snake_case)]
 #![allow(dead_code)]
diff --git a/src/backend/vector/avx2/mod.rs b/src/backend/vector/avx2/mod.rs
index a499e57..a1a21eb 100644
--- a/src/backend/vector/avx2/mod.rs
+++ b/src/backend/vector/avx2/mod.rs
@@ -19,7 +19,7 @@
 // missing).
 #![cfg_attr(
     all(feature = "nightly", feature = "stage2_build"),
-    doc(include = "../../docs/avx2-notes.md")
+    doc(include = "../../../../docs/avx2-notes.md")
 )]
 
 pub(crate) mod field;
diff --git a/src/backend/vector/ifma/field.rs b/src/backend/vector/ifma/field.rs
index d220a5a..80bd86f 100644
--- a/src/backend/vector/ifma/field.rs
+++ b/src/backend/vector/ifma/field.rs
@@ -14,12 +14,18 @@ use packed_simd::{u64x4, IntoBits};
 
 use backend::serial::u64::field::FieldElement51;
 
-#[allow(improper_ctypes)]
-extern "C" {
-    #[link_name = "llvm.x86.avx512.vpmadd52l.uq.256"]
-    fn madd52lo(z: u64x4, x: u64x4, y: u64x4) -> u64x4;
-    #[link_name = "llvm.x86.avx512.vpmadd52h.uq.256"]
-    fn madd52hi(z: u64x4, x: u64x4, y: u64x4) -> u64x4;
+/// A wrapper around `vpmadd52luq` that works on `u64x4`.
+#[inline(always)]
+unsafe fn madd52lo(z: u64x4, x: u64x4, y: u64x4) -> u64x4 {
+    use core::arch::x86_64::_mm256_madd52lo_epu64;
+    _mm256_madd52lo_epu64(z.into_bits(), x.into_bits(), y.into_bits()).into_bits()
+}
+
+/// A wrapper around `vpmadd52huq` that works on `u64x4`.
+#[inline(always)]
+unsafe fn madd52hi(z: u64x4, x: u64x4, y: u64x4) -> u64x4 {
+    use core::arch::x86_64::_mm256_madd52hi_epu64;
+    _mm256_madd52hi_epu64(z.into_bits(), x.into_bits(), y.into_bits()).into_bits()
 }
 
 /// A vector of four field elements in radix 2^51, with unreduced coefficients.
@@ -203,11 +209,7 @@ use subtle::ConditionallySelectable;
 
 impl ConditionallySelectable for F51x4Reduced {
     #[inline]
-    fn conditional_select(
-        a: &F51x4Reduced,
-        b: &F51x4Reduced,
-        choice: Choice,
-    ) -> F51x4Reduced {
+    fn conditional_select(a: &F51x4Reduced, b: &F51x4Reduced, choice: Choice) -> F51x4Reduced {
         let mask = (-(choice.unwrap_u8() as i64)) as u64;
         let mask_vec = u64x4::splat(mask);
         F51x4Reduced([
diff --git a/src/backend/vector/ifma/mod.rs b/src/backend/vector/ifma/mod.rs
index 6be5e72..f1bb22a 100644
--- a/src/backend/vector/ifma/mod.rs
+++ b/src/backend/vector/ifma/mod.rs
@@ -9,7 +9,7 @@
 
 #![cfg_attr(
     all(feature = "nightly", feature = "stage2_build"),
-    doc(include = "../../docs/ifma-notes.md")
+    doc(include = "../../../../docs/ifma-notes.md")
 )]
 
 pub mod field;
diff --git a/src/backend/vector/mod.rs b/src/backend/vector/mod.rs
index 1b42da3..95b5446 100644
--- a/src/backend/vector/mod.rs
+++ b/src/backend/vector/mod.rs
@@ -19,7 +19,7 @@
 // missing).
 #![cfg_attr(
     all(feature = "nightly", feature = "stage2_build"),
-    doc(include = "../../docs/parallel-formulas.md")
+    doc(include = "../../../docs/parallel-formulas.md")
 )]
 
 #[cfg(not(any(target_feature = "avx2", target_feature = "avx512ifma", rustdoc)))]
diff --git a/src/edwards.rs b/src/edwards.rs
index c60069c..e03230b 100644
--- a/src/edwards.rs
+++ b/src/edwards.rs
@@ -84,7 +84,7 @@
 //! successful decompression of a compressed point, or else by
 //! operations on other (valid) `EdwardsPoint`s.
 //!
-//! [curve_models]: https://doc-internal.dalek.rs/curve25519_dalek/curve_models/index.html
+//! [curve_models]: https://doc-internal.dalek.rs/curve25519_dalek/backend/serial/curve_models/index.html
 
 // We allow non snake_case names because coordinates in projective space are
 // traditionally denoted by the capitalisation of their respective
@@ -190,11 +190,10 @@ impl CompressedEdwardsY {
 
         if is_valid_y_coord.unwrap_u8() != 1u8 { return None; }
 
-        // Flip the sign of X if it's not correct
+         // FieldElement::sqrt_ratio_i always returns the nonnegative square root,
+         // so we negate according to the supplied sign bit.
         let compressed_sign_bit = Choice::from(self.as_bytes()[31] >> 7);
-        let    current_sign_bit = X.is_negative();
-
-        X.conditional_negate(current_sign_bit ^ compressed_sign_bit);
+        X.conditional_negate(compressed_sign_bit);
 
         Some(EdwardsPoint{ X: X, Y: Y, Z: Z, T: &X * &Y })
     }
diff --git a/src/field.rs b/src/field.rs
index e154f24..58dd60a 100644
--- a/src/field.rs
+++ b/src/field.rs
@@ -154,6 +154,9 @@ impl FieldElement {
             acc = &acc * input;
         }
 
+	// acc is nonzero iff all inputs are nonzero
+        assert_eq!(acc.is_zero().unwrap_u8(), 0);
+
         // Compute the inverse of all products
         acc = acc.invert();
 
diff --git a/src/lib.rs b/src/lib.rs
index 07c3bef..f1ae938 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -9,17 +9,11 @@
 // - Henry de Valence <hdevalence@hdevalence.ca>
 
 #![no_std]
-#![cfg_attr(
-    any(
-        all(feature = "simd_backend", target_feature = "avx512ifma"),
-        all(feature = "nightly", rustdoc)
-    ),
-    feature(simd_ffi, link_llvm_intrinsics)
-)]
 #![cfg_attr(feature = "nightly", feature(test))]
 #![cfg_attr(all(feature = "alloc", not(feature = "std")), feature(alloc))]
 #![cfg_attr(feature = "nightly", feature(external_doc))]
 #![cfg_attr(feature = "nightly", feature(doc_cfg))]
+#![cfg_attr(feature = "simd_backend", feature(stdsimd))]
 // Refuse to compile if documentation is missing, but only on nightly.
 //
 // This means that missing docs will still fail CI, but means we can use
diff --git a/src/scalar.rs b/src/scalar.rs
index 575fefe..b5253ac 100644
--- a/src/scalar.rs
+++ b/src/scalar.rs
@@ -297,7 +297,7 @@ 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();
+        *self = *self + _rhs;
     }
 }
 
@@ -305,8 +305,17 @@ define_add_assign_variants!(LHS = Scalar, RHS = Scalar);
 
 impl<'a, 'b> Add<&'b Scalar> for &'a Scalar {
     type Output = Scalar;
+    #[allow(non_snake_case)]
     fn add(self, _rhs: &'b Scalar) -> Scalar {
-        UnpackedScalar::add(&self.unpack(), &_rhs.unpack()).pack()
+        // The UnpackedScalar::add function produces reduced outputs
+        // if the inputs are reduced.  However, these inputs may not
+        // be reduced -- they might come from Scalar::from_bits.  So
+        // after computing the sum, we explicitly reduce it mod l
+        // before repacking.
+        let sum = UnpackedScalar::add(&self.unpack(), &_rhs.unpack());
+        let sum_R = UnpackedScalar::mul_internal(&sum, &constants::R);
+        let sum_mod_l = UnpackedScalar::montgomery_reduce(&sum_R);
+        sum_mod_l.pack()
     }
 }
 
@@ -314,7 +323,7 @@ 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();
+        *self = *self - _rhs;
     }
 }
 
@@ -322,8 +331,18 @@ 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 {
-        UnpackedScalar::sub(&self.unpack(), &_rhs.unpack()).pack()
+    #[allow(non_snake_case)]
+    fn sub(self, rhs: &'b Scalar) -> Scalar {
+        // The UnpackedScalar::sub function requires reduced inputs
+        // and produces reduced output. However, these inputs may not
+        // be reduced -- they might come from Scalar::from_bits.  So
+        // we explicitly reduce the inputs.
+        let self_R = UnpackedScalar::mul_internal(&self.unpack(), &constants::R);
+        let self_mod_l = UnpackedScalar::montgomery_reduce(&self_R);
+        let rhs_R = UnpackedScalar::mul_internal(&rhs.unpack(), &constants::R);
+        let rhs_mod_l = UnpackedScalar::montgomery_reduce(&rhs_R);
+
+        UnpackedScalar::sub(&self_mod_l, &rhs_mod_l).pack()
     }
 }
 
@@ -331,8 +350,11 @@ define_sub_variants!(LHS = Scalar, RHS = Scalar, Output = Scalar);
 
 impl<'a> Neg for &'a Scalar {
     type Output = Scalar;
+    #[allow(non_snake_case)]
     fn neg(self) -> Scalar {
-        &Scalar::zero() - self
+        let self_R = UnpackedScalar::mul_internal(&self.unpack(), &constants::R);
+        let self_mod_l = UnpackedScalar::montgomery_reduce(&self_R);
+        UnpackedScalar::sub(&UnpackedScalar::zero(), &self_mod_l).pack() 
     }
 }
 
@@ -917,7 +939,7 @@ impl Scalar {
                 naf[pos] = window as i8;
             } else {
                 carry = 1;
-                naf[pos] = (window as i8) - (width as i8);
+                naf[pos] = (window as i8).wrapping_sub(width as i8);
             }
 
             pos += w;
@@ -1233,6 +1255,33 @@ mod test {
          0,0,0,11,0,0,0,0,0,15,0,0,0,0,0,-9,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,7,
          0,0,0,0,0,-15,0,0,0,0,0,15,0,0,0,0,15,0,0,0,0,15,0,0,0,0,0,1,0,0,0,0];
 
+    static LARGEST_ED25519_S: Scalar = Scalar {
+        bytes: [
+            0xf8, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x7f,
+        ],
+    };
+
+    static CANONICAL_LARGEST_ED25519_S_PLUS_ONE: Scalar = Scalar {
+        bytes: [
+            0x7e, 0x34, 0x47, 0x75, 0x47, 0x4a, 0x7f, 0x97,
+            0x23, 0xb6, 0x3a, 0x8b, 0xe9, 0x2a, 0xe7, 0x6d,
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x0f,
+        ],
+    };
+
+    static CANONICAL_LARGEST_ED25519_S_MINUS_ONE: Scalar = Scalar {
+        bytes: [
+            0x7c, 0x34, 0x47, 0x75, 0x47, 0x4a, 0x7f, 0x97,
+            0x23, 0xb6, 0x3a, 0x8b, 0xe9, 0x2a, 0xe7, 0x6d,
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x0f,
+        ],
+    };
+
     #[test]
     fn fuzzer_testcase_reduction() {
         // LE bytes of 24519928653854221733733552434404946937899825954937634815
@@ -1264,13 +1313,42 @@ mod test {
     }
 
     #[test]
-    fn non_adjacent_form() {
+    fn non_adjacent_form_test_vector() {
         let naf = A_SCALAR.non_adjacent_form(5);
         for i in 0..256 {
             assert_eq!(naf[i], A_NAF[i]);
         }
     }
 
+    fn non_adjacent_form_iter(w: usize, x: &Scalar) {
+        let naf = x.non_adjacent_form(w);
+
+        // Reconstruct the scalar from the computed NAF
+        let mut y = Scalar::zero();
+        for i in (0..256).rev() {
+            y += y;
+            let digit = if naf[i] < 0 {
+                -Scalar::from((-naf[i]) as u64)
+            } else {
+                Scalar::from(naf[i] as u64)
+            };
+            y += digit;
+        }
+
+        assert_eq!(*x, y);
+    }
+
+    #[test]
+    fn non_adjacent_form_random() {
+        let mut rng = rand::thread_rng();
+        for _ in 0..1_000 {
+            let x = Scalar::random(&mut rng);
+            for w in &[5, 6, 7, 8] {
+                non_adjacent_form_iter(*w, &x);
+            }
+        }
+    }
+
     #[test]
     fn from_u64() {
         let val: u64 = 0xdeadbeefdeadbeef;
@@ -1293,6 +1371,82 @@ mod test {
         }
     }
 
+    #[test]
+    fn add_reduces() {
+        // Check that the addition works
+        assert_eq!(
+            (LARGEST_ED25519_S + Scalar::one()).reduce(),
+            CANONICAL_LARGEST_ED25519_S_PLUS_ONE
+        );
+        // Check that the addition reduces
+        assert_eq!(
+            LARGEST_ED25519_S + Scalar::one(),
+            CANONICAL_LARGEST_ED25519_S_PLUS_ONE
+        );
+    }
+
+    #[test]
+    fn sub_reduces() {
+        // Check that the subtraction works
+        assert_eq!(
+            (LARGEST_ED25519_S - Scalar::one()).reduce(),
+            CANONICAL_LARGEST_ED25519_S_MINUS_ONE
+        );
+        // Check that the subtraction reduces
+        assert_eq!(
+            LARGEST_ED25519_S - Scalar::one(),
+            CANONICAL_LARGEST_ED25519_S_MINUS_ONE
+        );
+    }
+
+    #[test]
+    fn quarkslab_scalar_overflow_does_not_occur() {
+        // Check that manually-constructing large Scalars with
+        // from_bits cannot produce incorrect results.
+        //
+        // The from_bits function is required to implement X/Ed25519,
+        // while all other methods of constructing a Scalar produce
+        // reduced Scalars.  However, this "invariant loophole" allows
+        // constructing large scalars which are not reduced mod l.
+        //
+        // This issue was discovered independently by both Jack
+        // "str4d" Grigg (issue #238), who noted that reduction was
+        // not performed on addition, and Laurent Grémy & Nicolas
+        // Surbayrole of Quarkslab, who noted that it was possible to
+        // cause an overflow and compute incorrect results.
+        //
+        // This test is adapted from the one suggested by Quarkslab.
+
+        let large_bytes = [
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x7f,
+        ];
+
+        let a = Scalar::from_bytes_mod_order(large_bytes);
+        let b = Scalar::from_bits(large_bytes);
+
+        assert_eq!(a, b.reduce());
+
+        let a_3 = a + a + a;
+        let b_3 = b + b + b;
+
+        assert_eq!(a_3, b_3);
+
+        let neg_a = -a;
+        let neg_b = -b;
+
+        assert_eq!(neg_a, neg_b);
+
+        let minus_a_3 = Scalar::zero() - a - a - a;
+        let minus_b_3 = Scalar::zero() - b - b - b;
+
+        assert_eq!(minus_a_3, minus_b_3);
+        assert_eq!(minus_a_3, -a_3);
+        assert_eq!(minus_b_3, -b_3);
+    }
+
     #[test]
     fn impl_add() {
         let two = Scalar::from(2u64);
diff --git a/src/window.rs b/src/window.rs
index 08f2d78..603760e 100644
--- a/src/window.rs
+++ b/src/window.rs
@@ -171,7 +171,7 @@ pub(crate) struct NafLookupTable8<T>(pub(crate) [T; 64]);
 impl<T: Copy> NafLookupTable8<T> {
     pub fn select(&self, x: usize) -> T {
         debug_assert_eq!(x & 1, 1);
-        debug_assert!(x < 256);
+        debug_assert!(x < 128);
 
         self.0[x / 2]
     }