//! Low-order (quadratic and cubic) Bézier curves.
// NOTE: Most info from https://pomax.github.io/bezierinfo

use num_traits::{Zero, real::Real};
use crate::ops::*;
use std::ops::*;
use std::ops::Add;
use crate::vec::repr_c::{
    Vec3 as CVec3,
    Vec4 as CVec4,
};

// WISH: OOBBs from beziers
// WISH: "Tracing a curve at fixed distance intervals"
// WISH: Line-curve intersection (Especially straight horizontal and straight vertical)

macro_rules! bezier_impl_any {
    (3 $Bezier:ident $Point:ident) => {

        bezier_impl_any!{$Bezier $Point}

        impl<T: Real> $Bezier<T> {
            /// Gets the Axis-Aligned Bounding Box for this curve.
            pub fn aabb(self) -> Aabb<T> {
                let (min_x, max_x) = self.x_bounds();
                let (min_y, max_y) = self.y_bounds();
                let (min_z, max_z) = self.z_bounds();
                Aabb {
                    min: Vec3::new(min_x, min_y, min_z),
                    max: Vec3::new(max_x, max_y, max_z),
                }
            }
            /// Returns this curve, flipping the `y` coordinate of each of its points.
            pub fn flipped_z(self) -> Self {
                self.into_vector().map(|mut p| {p.z = -p.z; p}).into()
            }
            /// Flips the `x` coordinate of all points of this curve.
            pub fn flip_z(&mut self) {
                *self = self.flipped_z();
            }
        }
        impl<T> Mul<$Bezier<T>> for Rows3<T> where T: Real + MulAdd<T,T,Output=T> {
            type Output = $Bezier<T>;
            fn mul(self, rhs: $Bezier<T>) -> $Bezier<T> {
                rhs.into_vector().map(|p| self * p).into()
            }
        }
        impl<T> Mul<$Bezier<T>> for Cols3<T> where T: Real + MulAdd<T,T,Output=T> {
            type Output = $Bezier<T>;
            fn mul(self, rhs: $Bezier<T>) -> $Bezier<T> {
                rhs.into_vector().map(|p| self * p).into()
            }
        }
        impl<T> Mul<$Bezier<T>> for Rows4<T> where T: Real + MulAdd<T,T,Output=T> {
            type Output = $Bezier<T>;
            fn mul(self, rhs: $Bezier<T>) -> $Bezier<T> {
                rhs.into_vector().map(|p| self.mul_point(p).into()).into()
            }
        }
        impl<T> Mul<$Bezier<T>> for Cols4<T> where T: Real + MulAdd<T,T,Output=T> {
            type Output = $Bezier<T>;
            fn mul(self, rhs: $Bezier<T>) -> $Bezier<T> {
                rhs.into_vector().map(|p| self.mul_point(p).into()).into()
            }
        }
    };
    (2 $Bezier:ident $Point:ident) => {

        bezier_impl_any!{$Bezier $Point}

        impl<T> Mul<$Bezier<T>> for Rows2<T> where T: Real + MulAdd<T,T,Output=T> {
            type Output = $Bezier<T>;
            fn mul(self, rhs: $Bezier<T>) -> $Bezier<T> {
                rhs.into_vector().map(|p| self * p).into()
            }
        }
        impl<T> Mul<$Bezier<T>> for Cols2<T> where T: Real + MulAdd<T,T,Output=T> {
            type Output = $Bezier<T>;
            fn mul(self, rhs: $Bezier<T>) -> $Bezier<T> {
                rhs.into_vector().map(|p| self * p).into()
            }
        }
        impl<T> Mul<$Bezier<T>> for Rows3<T> where T: Real + MulAdd<T,T,Output=T> {
            type Output = $Bezier<T>;
            fn mul(self, rhs: $Bezier<T>) -> $Bezier<T> {
                rhs.into_vector().map(|p| self.mul_point_2d(p).into()).into()
            }
        }
        impl<T> Mul<$Bezier<T>> for Cols3<T> where T: Real + MulAdd<T,T,Output=T> {
            type Output = $Bezier<T>;
            fn mul(self, rhs: $Bezier<T>) -> $Bezier<T> {
                rhs.into_vector().map(|p| self.mul_point_2d(p).into()).into()
            }
        }
    };
    ($Bezier:ident $Point:ident) => {
        impl<T: Real> $Bezier<T> {
            /// Evaluates the normalized tangent at interpolation factor `t`.
            pub fn normalized_tangent(self, t: T) -> Point<T> where T: Add<T, Output=T> {
                self.evaluate_derivative(t).normalized()
            }
            // WISH: better length approximation estimations (e.g see https://math.stackexchange.com/a/61796)
            /// Approximates the curve's length by subdividing it into step_count+1 segments.
            pub fn length_by_discretization(self, step_count: u16) -> T
                where T: Add<T, Output=T> + From<u16>
            {
                let mut length = T::zero();
                let mut prev_point = self.evaluate(T::zero());
                for i in 1..(step_count+2) {
                    let t = <T as From<u16>>::from(i)/(<T as From<u16>>::from(step_count)+T::one());
                    let next_point = self.evaluate(t);
                    length = length + (next_point - prev_point).magnitude();
                    prev_point = next_point;
                }
                length
            }

            /// Gets the Axis-Aligned Bounding Rectangle for this curve.
            ///
            /// On 3D curves, this discards the `z` values.
            pub fn aabr(self) -> Aabr<T> {
                let (min_x, max_x) = self.x_bounds();
                let (min_y, max_y) = self.y_bounds();
                Aabr {
                    min: Vec2::new(min_x, min_y),
                    max: Vec2::new(max_x, max_y),
                }
            }
            /// Returns this curve, flipping the `x` coordinate of each of its points.
            pub fn flipped_x(self) -> Self {
                self.into_vector().map(|mut p| {p.x = -p.x; p}).into()
            }
            /// Returns this curve, flipping the `y` coordinate of each of its points.
            pub fn flipped_y(self) -> Self {
                self.into_vector().map(|mut p| {p.y = -p.y; p}).into()
            }
            /// Flips the `x` coordinate of all points of this curve.
            pub fn flip_x(&mut self) {
                *self = self.flipped_x();
            }
            /// Flips the `y` coordinate of all points of this curve.
            pub fn flip_y(&mut self) {
                *self = self.flipped_y();
            }

            // TODO: Test this! binary_search_point_by_steps
            /// Searches for the point lying on this curve that is closest to `p`.
            ///
            /// `steps` is the number of points to sample in the curve for the "broad phase"
            /// that takes place before the binary search.
            ///
            /// `epsilon` denotes the desired precision for the result. The higher it is, the
            /// sooner the algorithm will finish, but the result would be less satisfactory.
            ///
            /// # Panics
            /// Panics if `epsilon` is less than or equal to `T::epsilon()`.  
            /// `epsilon` must be positive and not approximately equal to zero.
            pub fn binary_search_point_by_steps(self, p: Point<T>, steps: u16, epsilon: T) -> (T, Point<T>) 
                where T: Add<T, Output=T> + From<u16>
            {
                let steps_f = <T as From<u16>>::from(steps);
                let it = (0..steps).map(|i| {
                    let i = <T as From<u16>>::from(i);
                    let t = i / steps_f;
                    (t, self.evaluate(t))
                });
                // half_interval = 1/(2*steps)
                let h = (steps_f + steps_f).recip();
                self.binary_search_point(p, it, h, epsilon)
            }
            // TODO: Test this! binary_search_point
            /// Searches for the point lying on this curve that is closest to `p`.
            ///
            /// For an example usage, see the source code of `binary_search_point_by_steps()`.
            ///
            /// `coarse` is an iterator over pairs of `(interpolation_value, point)` that are
            /// assumed to, together, represent a discretization of the curve.  
            /// This parameter is used for a "broad phase" - the point yielded by `coarse` that is
            /// closest to `p` is the starting point for the binary search.  
            /// `coarse` may very well yield a single pair; Also, it was designed so that,
            /// if you already have the values handy, there is no need to recompute them.  
            /// This function doesn't panic if `coarse` yields no element, but then it would be
            /// very unlikely for the result to be satisfactory.
            ///
            /// `half_interval` is the starting value for the half of the binary search interval.
            ///
            /// `epsilon` denotes the desired precision for the result. The higher it is, the
            /// sooner the algorithm will finish, but the result would be less satisfactory.
            ///
            /// # Panics
            /// Panics if `epsilon` is less than or equal to `T::epsilon()`.  
            /// `epsilon` must be positive and not approximately equal to zero.
            pub fn binary_search_point<I>(self, p: Point<T>, coarse: I, half_interval: T, epsilon: T) -> (T, Point<T>)
                where T: Add<T, Output=T>, I: IntoIterator<Item=(T, Point<T>)>
            {
                debug_assert!(epsilon > T::epsilon());
                let mut t = T::one();
                let mut pt = self.end;
                let mut d = pt.distance_squared(p);
                for (t_, pt_) in coarse {
                    let d_ = pt_.distance_squared(p);
                    if d_ < d {
                        d = d_; pt = pt_; t = t_;
                    }
                }
                let mut h = half_interval;
                while h >= epsilon {
                    let (p1, p2) = (self.evaluate(t-h), self.evaluate(t+h));
                    let (d1, d2) = (p.distance_squared(p1), p.distance_squared(p2));
                    if d1 < d || d2 < d {
                        if d1 < d2 {
                            d = d1; pt = p1; t = t - h;
                        } else {
                            d = d2; pt = p2; t = t + h;
                        }
                        continue;
                    }
                    h = h / (T::one() + T::one());
                }
                (t, pt)
            }
        }
    };
}

macro_rules! bezier_impl_quadratic_axis {
    ($QuadraticBezier:ident $Point:ident ($x_s:expr) $x:ident $x_inflection:ident $x_min:ident $x_max:ident $x_bounds:ident) => {
        impl<T: Real> $QuadraticBezier<T> {
            /// Returns the evaluation factor that gives an inflection point along the
            #[doc=$x_s]
            /// axis, if any.
            // Code in part taken from `lyon` crate, geom.
            // Also explained at https://pomax.github.io/bezierinfo/#extremities
            pub fn $x_inflection(self) -> Option<T> {
                let div = self.start.$x - (self.ctrl.$x + self.ctrl.$x) + self.end.$x;
                if div.abs() <= T::epsilon() {
                    return None;
                }
                let t = (self.start.$x - self.ctrl.$x) / div;
                if T::zero() <= t && t <= T::one() {
                    return Some(t);
                }
                return None;
            }
            /// Returns the evaluation factor that gives the point on the curve which
            #[doc=$x_s]
            /// coordinate is the minimum.
            pub fn $x_min(self) -> T {
                if let Some(t) = self.$x_inflection() {
                    let p = self.evaluate(t);
                    if p.$x < self.start.$x && p.$x < self.end.$x {
                        return t;
                    }
                }
                if self.start.$x < self.end.$x { T::zero() } else { T::one() }
            }
            /// Returns the evaluation factor that gives the point on the curve which
            #[doc=$x_s]
            /// coordinate is the maximum.
            pub fn $x_max(self) -> T {
                if let Some(t) = self.$x_inflection() {
                    let p = self.evaluate(t);
                    if p.$x > self.start.$x && p.$x > self.end.$x {
                        return t;
                    }
                }
                if self.start.$x > self.end.$x { T::zero() } else { T::one() }
            }
            /// Returns the evaluation factors that give the points on the curve which
            #[doc=$x_s]
            /// coordinates are the respective minimum and maximum.
            pub fn $x_bounds(self) -> (T, T) {
                // PERF: We don't need to compute $x_inflections twice!
                (self.$x_min(), self.$x_max())
            }
        }
    };
}

macro_rules! bezier_impl_cubic_axis {
    ($CubicBezier:ident $Point:ident ($x_s:expr) $x:ident $x_inflections:ident $x_min:ident $x_max:ident $x_bounds:ident) => {
        impl<T: Real> $CubicBezier<T> {
            /// Returns the evaluation factor that gives an inflection point along the
            #[doc=$x_s]
            /// axis, if any.
            // Code in part taken from `lyon` crate, geom.
            // Also explained at https://pomax.github.io/bezierinfo/#extremities
            pub fn $x_inflections(self) -> Option<(T, Option<T>)> {
                // See www.faculty.idc.ac.il/arik/quality/appendixa.html for an explanation
                // The derivative of a cubic bezier curve is a curve representing a second degree polynomial function
                // f(x) = a * x² + b * x + c such as :
                let two = T::one() + T::one();
                let three = two + T::one();
                let four = three + T::one();
                let six = three + three;
                let a = three * (self.end.$x - three * self.ctrl1.$x + three * self.ctrl0.$x - self.start.$x);
                let b = six * (self.ctrl1.$x - two * self.ctrl0.$x + self.start.$x);
                let c = three * (self.ctrl0.$x - self.start.$x);

                // If the derivative is a linear function
                if a.abs() <= T::epsilon() {
                    return if b.abs() <= T::epsilon() {
                        // If the derivative is a constant function
                        if c.abs() <= T::epsilon() {
                            Some((T::zero(), None))
                        } else {
                            None
                        }
                    } else {
                        Some((-c / b, None))
                    };
                }

                // Wants to use IsBetween01, but that would annoyingly propagate the trait bound.
                let is_between01 = |t| { T::zero() < t && t < T::one() };

                let discriminant = b * b - four * a * c;

                // There is no Real solution for the equation
                if discriminant < T::zero() {
                    return None;
                }

                // There is one Real solution for the equation
                if discriminant.abs() <= T::epsilon() {
                    let t = -b / (a + a);
                    return if is_between01(t) {
                        Some((t, None))
                    } else {
                        None
                    };
                }

                // There are two Real solutions for the equation
                let discriminant_sqrt = discriminant.sqrt();

                let first_extremum = (-b - discriminant_sqrt) / (a + a);
                let second_extremum = (-b + discriminant_sqrt) / (a + a);

                if is_between01(first_extremum) {
                    if is_between01(second_extremum) {
                        Some((first_extremum, Some(second_extremum)))
                    } else {
                        Some((first_extremum, None))
                    }
                } else {
                    if is_between01(second_extremum) {
                        Some((second_extremum, None))
                    } else {
                        None
                    }
                }
            }
            /// Returns the evaluation factor that gives the point on the curve which
            #[doc=$x_s]
            /// coordinate is the minimum.
            pub fn $x_min(self) -> T {
                if let Some((t1, t2)) = self.$x_inflections() {
                    let p1 = self.evaluate(t1);
                    if let Some(t2) = t2 {
                        let p2 = self.evaluate(t2);
                        if (p1.$x < self.start.$x && p1.$x < self.end.$x)
                        || (p2.$x < self.start.$x && p2.$x < self.end.$x)
                        {
                            return if p1.$x < p2.$x { t1 } else { t2 };
                        }
                    } else {
                        if p1.$x < self.start.$x && p1.$x < self.end.$x {
                            return t1;
                        }
                    }
                }
                if self.start.$x < self.end.$x { T::zero() } else { T::one() }
            }
            /// Returns the evaluation factor that gives the point on the curve which
            #[doc=$x_s]
            /// coordinate is the maximum.
            pub fn $x_max(self) -> T {
                if let Some((t1, t2)) = self.$x_inflections() {
                    let p1 = self.evaluate(t1);
                    if let Some(t2) = t2 {
                        let p2 = self.evaluate(t2);
                        if (p1.$x > self.start.$x && p1.$x > self.end.$x)
                        || (p2.$x > self.start.$x && p2.$x > self.end.$x)
                        {
                            return if p1.$x > p2.$x { t1 } else { t2 };
                        }
                    } else {
                        if p1.$x > self.start.$x && p1.$x > self.end.$x {
                            return t1;
                        }
                    }
                }
                if self.start.$x > self.end.$x { T::zero() } else { T::one() }
            }

            /// Returns the evaluation factors that give the points on the curve which
            #[doc=$x_s]
            /// coordinates are the respective minimum and maximum.
            pub fn $x_bounds(self) -> (T, T) {
                // PERF: We don't need to compute $x_inflections twice!
                (self.$x_min(), self.$x_max())
            }
        }
    };
}
macro_rules! bezier_impl_quadratic {
    ($(#[$attrs:meta])* 3 $QuadraticBezier:ident $CubicBezier:ident $Point:ident $LineSegment:ident) => {
        bezier_impl_quadratic!{$(#[$attrs])* $QuadraticBezier $CubicBezier $Point $LineSegment}
        bezier_impl_quadratic_axis!{$QuadraticBezier $Point ("Z") z z_inflection min_z max_z z_bounds}
        bezier_impl_any!(3 $QuadraticBezier $Point);
    };
    ($(#[$attrs:meta])* 2 $QuadraticBezier:ident $CubicBezier:ident $Point:ident $LineSegment:ident) => {
        bezier_impl_quadratic!{$(#[$attrs])* $QuadraticBezier $CubicBezier $Point $LineSegment}
        bezier_impl_any!(2 $QuadraticBezier $Point);
    };
    ($(#[$attrs:meta])* $QuadraticBezier:ident $CubicBezier:ident $Point:ident $LineSegment:ident) => {

        type Point<T> = $Point<T>;
        
        $(#[$attrs])*
        #[derive(Debug, Default, Copy, Clone, Hash, PartialEq, Eq, /*PartialOrd, Ord*/)]
        #[cfg_attr(feature="serde", derive(Serialize, Deserialize))]
        ///
        /// Also note that quadratic Bézier curves are quite bad at approximating circles.
        /// See [the relevant section of "A Primer on Bézier Curves"](https://pomax.github.io/bezierinfo/#circles)
        /// for an explanation.
        pub struct $QuadraticBezier<T> {
            /// Starting point of the curve.
            pub start: Point<T>,
            /// Control point of the curve.
            pub ctrl: Point<T>, 
            /// End point of the curve.
            pub end: Point<T>,
        }
        
        impl<T: Real> $QuadraticBezier<T> {
            /// Evaluates the position of the point lying on the curve at interpolation factor `t`.
            ///
            /// This is one of the most important Bézier curve operations,
            /// because, in one way or another, it is used to render a curve
            /// to the screen.
            /// The common use case is to successively evaluate a curve at a set of values
            /// that range from 0 to 1, to approximate the curve as an array of
            /// line segments which are then rendered.
            pub fn evaluate(self, t: T) -> Point<T> {
                let l = T::one();
                let two = l+l;
                self.start*(l-t)*(l-t) + self.ctrl*two*(l-t)*t + self.end*t*t
            }
            /// Evaluates the derivative tangent at interpolation factor `t`, which happens to give
            /// a non-normalized tangent vector.
            ///
            /// See also `normalized_tangent()`.
            pub fn evaluate_derivative(self, t: T) -> Point<T> {
                let l = T::one();
                let n = l+l;
                (self.ctrl-self.start)*(l-t)*n + (self.end-self.ctrl)*t*n
            }
            /// Returns the constant matrix M such that,
            /// given `T = [1, t*t, t*t*t]` and `P` the vector of control points,
            /// `dot(T * M, P)` evalutes the Bezier curve at 't'.
            ///
            /// This function name is arguably dubious.
            pub fn matrix() -> Mat3<T> {
                let zero = T::zero();
                let one = T::one();
                let two = one+one;
                Mat3 {
                    rows: CVec3::new(
                        Vec3::new( one,  zero, zero),
                        Vec3::new(-two,  two, zero),
                        Vec3::new( one, -two, one),
                    )
                }
            }
            /// Splits this quadratic Bézier curve into two curves, at interpolation factor `t`.
            // NOTE that some computations may be reused, but the compiler can
            // reason about these. Clarity wins here IMO.
            pub fn split(self, t: T) -> [Self; 2] {
                let l = T::one();
                let two = l+l;
                let first = $QuadraticBezier {
                    start: self.start,
                    ctrl:  self.ctrl*t - self.start*(t-l),
                    end:   self.end*t*t - self.ctrl*two*t*(t-l) + self.start*(t-l)*(t-l),
                };
                let second = $QuadraticBezier {
                    start: self.end*t*t - self.ctrl*two*t*(t-l) + self.start*(t-l)*(t-l),
                    ctrl:  self.end*t - self.ctrl*(t-l),
                    end:   self.end,
                };
                [first, second]
            }
            /// Elevates this curve into a cubic Bézier curve.
            pub fn into_cubic(self) -> $CubicBezier<T> {
                self.into()
            }
        }
        impl<T> $QuadraticBezier<T> {
            /// Gets this curve reversed, i.e swaps `start` with `end`.
            pub fn reversed(mut self) -> Self {
                self.reverse();
                self
            }
            /// Reverses this curve, i.e swaps `start` with `end`.
            pub fn reverse(&mut self) {
                std::mem::swap(&mut self.start, &mut self.end);
            }
            // Convenience for this module
            pub(crate) fn into_vector(self) -> Vec3<Point<T>> {
                self.into_vec3()
            }
            /// Converts this curve into a `Vec3` of points.
            pub fn into_vec3(self) -> Vec3<Point<T>> {
                self.into()
            }
            /// Converts this curve into a tuple of points.
            pub fn into_tuple(self) -> (Point<T>, Point<T>, Point<T>) {
                self.into_vec3().into_tuple()
            }
            /// Converts this curve into an array of points.
            pub fn into_array(self) -> [Point<T>; 3] {
                self.into_vec3().into_array()
            }
        }
        
        impl<T> From<Vec3<Point<T>>> for $QuadraticBezier<T> {
            fn from(v: Vec3<Point<T>>) -> Self {
                $QuadraticBezier {
                    start: v.x, 
                    ctrl: v.y, 
                    end: v.z
                }
            }
        }
        impl<T> From<$QuadraticBezier<T>> for Vec3<Point<T>> {
            fn from(v: $QuadraticBezier<T>) -> Self {
                Vec3::new(v.start, v.ctrl, v.end)
            }
        }

        impl<T: Real> From<$LineSegment<T>> for $QuadraticBezier<T> {
            fn from(line_segment: $LineSegment<T>) -> Self {
                let ctrl = (line_segment.start + line_segment.end) / (T::one() + T::one());
                Self {
                    start: line_segment.start, 
                    ctrl, 
                    end:   line_segment.end
                }
            }
        }
        impl<T: Real> From<Range<Point<T>>> for $QuadraticBezier<T> {
            fn from(range: Range<Point<T>>) -> Self {
                Self::from($LineSegment::from(range))
            }
        }
        
        bezier_impl_quadratic_axis!{$QuadraticBezier $Point ("X") x x_inflection min_x max_x x_bounds}
        bezier_impl_quadratic_axis!{$QuadraticBezier $Point ("Y") y y_inflection min_y max_y y_bounds}
    }
}

// NOTE: The reason to split 2D-3D conversion into two macros is
// to make sure that each is displayed in the correct documentation page.
macro_rules! bezier_impl_2d_into_3d {
    ($Bezier2:ident $Bezier3:ident) => {
        impl<T: Zero> $Bezier2<T> {
            /// Converts this 2D curve to a 3D one, setting the `z` elements to zero.
            pub fn into_3d(self) -> $Bezier3<T> {
                self.into()
            }
        }
        impl<T: Zero> From<$Bezier2<T>> for $Bezier3<T> {
            fn from(c: $Bezier2<T>) -> Self {
                c.into_vector().map(Into::into).into()
            }
        }
    };
}
macro_rules! bezier_impl_3d_into_2d {
    ($Bezier3:ident $Bezier2:ident) => {
        impl<T> $Bezier3<T> {
            /// Converts this 3D curve to a 2D one, dropping the `z` elements.
            pub fn into_2d(self) -> $Bezier2<T> { 
                self.into() 
            }
        }
        impl<T> From<$Bezier3<T>> for $Bezier2<T> {
            fn from(c: $Bezier3<T>) -> Self {
                c.into_vector().map(Into::into).into()
            }
        }
    };
}

macro_rules! bezier_impl_cubic {
    ($(#[$attrs:meta])* 3 $QuadraticBezier:ident $CubicBezier:ident $Point:ident $LineSegment:ident) => {
        bezier_impl_cubic!{$(#[$attrs])* $QuadraticBezier $CubicBezier $Point $LineSegment}
        bezier_impl_cubic_axis!{$CubicBezier $Point ("Z") z z_inflections min_z max_z z_bounds}
        bezier_impl_any!(3 $CubicBezier $Point);
    };
    ($(#[$attrs:meta])* 2 $QuadraticBezier:ident $CubicBezier:ident $Point:ident $LineSegment:ident) => {
        bezier_impl_cubic!{$(#[$attrs])* $QuadraticBezier $CubicBezier $Point $LineSegment}
        bezier_impl_any!(2 $CubicBezier $Point);
    };
    ($(#[$attrs:meta])* $QuadraticBezier:ident $CubicBezier:ident $Point:ident $LineSegment:ident) => {

        type Point<T> = $Point<T>;
        
        $(#[$attrs])*
        #[derive(Debug, Default, Copy, Clone, Hash, PartialEq, Eq, /*PartialOrd, Ord*/)]
        #[cfg_attr(feature="serde", derive(Serialize, Deserialize))]
        pub struct $CubicBezier<T> {
            /// Starting point of the curve.
            pub start: Point<T>, 
            /// First control point of the curve, associated with `start`.
            pub ctrl0: Point<T>, 
            /// Second control point of the curve, associated with `end`.
            pub ctrl1: Point<T>,
            /// End point of the curve.
            pub end: Point<T>,
        }

        impl<T: Real> $CubicBezier<T> {
            /// Evaluates the position of the point lying on the curve at interpolation factor `t`.
            ///
            /// This is one of the most important Bézier curve operations,
            /// because, in one way or another, it is used to render a curve
            /// to the screen.
            /// The common use case is to successively evaluate a curve at a set of values
            /// that range from 0 to 1, to approximate the curve as an array of
            /// line segments which are then rendered.
            pub fn evaluate(self, t: T) -> Point<T> {
                let l = T::one();
                let three = l+l+l;
                self.start*(l-t)*(l-t)*(l-t) + self.ctrl0*three*(l-t)*(l-t)*t + self.ctrl1*three*(l-t)*t*t + self.end*t*t*t
            }
            /// Evaluates the derivative tangent at interpolation factor `t`, which happens to give
            /// a non-normalized tangent vector.
            ///
            /// See also `normalized_tangent()`.
            pub fn evaluate_derivative(self, t: T) -> Point<T> {
                let l = T::one();
                let n = l+l+l;
                let two = l+l;
                (self.ctrl0-self.start)*(l-t)*(l-t)*n + (self.ctrl1-self.ctrl0)*two*(l-t)*t*n + (self.end-self.ctrl1)*t*t*n
            }
            /// Returns the constant matrix M such that,
            /// given `T = [1, t*t, t*t*t, t*t*t*t]` and `P` the vector of control points,
            /// `dot(T * M, P)` evalutes the Bezier curve at 't'.
            ///
            /// This function name is arguably dubious.
            pub fn matrix() -> Mat4<T> {
                let zero = T::zero();
                let one = T::one();
                let three = one+one+one;
                let six = three + three;
                Mat4 {
                    rows: CVec4::new(
                        Vec4::new( one,  zero,  zero, zero),
                        Vec4::new(-three,  three,  zero, zero),
                        Vec4::new( three, -six,  three, zero),
                        Vec4::new(-one,  three, -three, one),
                    )
                }
            }
            /// Splits this cubic Bézier curve into two curves, at interpolation factor `t`.
            // NOTE that some computations may be reused, but the compiler can
            // reason about these. Clarity wins here IMO.
            pub fn split(self, t: T) -> [Self; 2] {
                let l = T::one();
                let two = l+l;
                let three = l+l+l;
                let first = $CubicBezier {
                    start: self.start,
                    ctrl0: self.ctrl0*t - self.start*(t-l),
                    ctrl1: self.ctrl1*t*t - self.ctrl0*two*t*(t-l) + self.start*(t-l)*(t-l),
                    end:   self.end*t*t*t - self.ctrl1*three*t*t*(t-l) + self.ctrl0*three*t*(t-l)*(t-l) - self.start*(t-l)*(t-l)*(t-l),
                };
                let second = $CubicBezier {
                    start: self.end*t*t*t - self.ctrl1*three*t*t*(t-l) + self.ctrl0*three*t*(t-l)*(t-l) - self.start*(t-l)*(t-l)*(t-l),
                    ctrl0: self.end*t*t - self.ctrl1*two*t*(t-l) + self.ctrl0*(t-l)*(t-l),
                    ctrl1: self.end*t - self.ctrl1*(t-l),
                    end:   self.end,
                };
                [first, second]
            }
            /// Gets the cubic Bézier curve that approximates a unit quarter circle.
            ///
            /// You can build a good-looking circle out of 4 curves by applying
            /// symmetries to this curve.
            pub fn unit_quarter_circle() -> Self {
                let (two, three) = (T::one()+T::one(), T::one()+T::one()+T::one());
                let coeff = (two+two)*(two.sqrt()-T::one())/three;
                Self {
                    start: Vec2::unit_x().into(),
                    ctrl0: (Vec2::unit_x() + Vec2::unit_y()*coeff).into(),
                    ctrl1: (Vec2::unit_x()*coeff + Vec2::unit_y()).into(),
                    end: Vec2::unit_y().into(),
                }
            }
            /// Gets the 4 cubic Bézier curves that, used together, approximate a unit quarter circle.
            ///
            /// The returned tuple is `(north-east, north-west, south-west, south-east)`.
            pub fn unit_circle() -> [Self; 4] {
                let a = Self::unit_quarter_circle();
                let b = a.flipped_x();
                let c = b.flipped_y();
                let d = a.flipped_y();
                [a, b, c, d]
            }
        }
        impl<T> $CubicBezier<T> {
            /// Gets this curve reversed, i.e swaps `start` with `end` and `ctrl0` with `ctrl1`.
            pub fn reversed(mut self) -> Self {
                self.reverse();
                self
            }
            /// Reverses this curve, i.e swaps `start` with `end` and `ctrl0` with `ctrl1`.
            pub fn reverse(&mut self) {
                std::mem::swap(&mut self.start, &mut self.end);
                std::mem::swap(&mut self.ctrl0, &mut self.ctrl1);
            }
            // Convenience for this module
            pub(crate) fn into_vector(self) -> Vec4<Point<T>> {
                self.into_vec4()
            }
            /// Converts this curve into a `Vec4` of points.
            pub fn into_vec4(self) -> Vec4<Point<T>> {
                self.into()
            }
            /// Converts this curve into a tuple of points.
            pub fn into_tuple(self) -> (Point<T>, Point<T>, Point<T>, Point<T>) {
                self.into_vec4().into_tuple()
            }
            /// Converts this curve into an array of points.
            pub fn into_array(self) -> [Point<T>; 4] {
                self.into_vec4().into_array()
            }
        }
        
        impl<T> From<Vec4<Point<T>>> for $CubicBezier<T> {
            fn from(v: Vec4<Point<T>>) -> Self {
                $CubicBezier {
                    start: v.x, 
                    ctrl0: v.y,
                    ctrl1: v.z, 
                    end: v.w
                }
            }
        }
        impl<T> From<$CubicBezier<T>> for Vec4<Point<T>> {
            fn from(v: $CubicBezier<T>) -> Self {
                Vec4::new(v.start, v.ctrl0, v.ctrl1, v.end)
            }
        }
        impl<T: Real + Lerp<T,Output=T>> From<$LineSegment<T>> for $CubicBezier<T> 
        {
            fn from(line_segment: $LineSegment<T>) -> Self {
                let three = T::one() + T::one() + T::one();
                let t = three.recip();
                let ctrl0 = Lerp::lerp_unclamped(line_segment.start, line_segment.end, t);
                let ctrl1 = Lerp::lerp_unclamped(line_segment.start, line_segment.end, t+t);
                Self {
                    start: line_segment.start, 
                    ctrl0, 
                    ctrl1, 
                    end:   line_segment.end
                }
            }
        }
        impl<T: Real + Lerp<T, Output=T>> From<Range<Point<T>>> for $CubicBezier<T> {
            fn from(range: Range<Point<T>>) -> Self {
                Self::from($LineSegment::from(range))
            }
        }
        impl<T: Real> From<$QuadraticBezier<T>> for $CubicBezier<T> {
            fn from(b: $QuadraticBezier<T>) -> Self {
                let three = T::one() + T::one() + T::one();
                $CubicBezier {
                    start: b.start,
                    ctrl0: (b.start + b.ctrl + b.ctrl) / three,
                    ctrl1: (b.end + b.ctrl + b.ctrl) / three,
                    end: b.end,
                }
            }
        }
        
        bezier_impl_cubic_axis!{$CubicBezier $Point ("X") x x_inflections min_x max_x x_bounds}
        bezier_impl_cubic_axis!{$CubicBezier $Point ("Y") y y_inflections min_y max_y y_bounds}
    }
}

macro_rules! impl_all_beziers {
    ($mod:ident) => {
        use  crate::vec::$mod::{Vec3, Vec4, Vec2};
        use  crate::mat::$mod::row_major::{Mat2 as Rows2, Mat3 as Rows3, Mat4 as Rows4};
        use  crate::mat::$mod::column_major::{Mat2 as Cols2, Mat3 as Cols3, Mat4 as Cols4};
        use crate::geom::$mod::{LineSegment2, LineSegment3, Aabr, Aabb};
        use self::Rows4 as Mat4;
        use self::Rows3 as Mat3;

        mod quadratic_bezier2 {
            use super::*;
            bezier_impl_quadratic!{
                /// A 2D Bézier curve with one control point.
                ///
                /// 2x2 and 3x3 matrices can be multiplied by a Bézier curve to transform all of its points.
                2 QuadraticBezier2 CubicBezier2 Vec2 LineSegment2
            }
            bezier_impl_2d_into_3d!{QuadraticBezier2 QuadraticBezier3}
        }

        mod quadratic_bezier3 {
            use super::*;
            bezier_impl_quadratic!{
                /// A 3D Bézier curve with one control point.
                ///
                /// 3x3 and 4x4 matrices can be multiplied by a Bézier curve to transform all of its points.
                3 QuadraticBezier3 CubicBezier3 Vec3 LineSegment3
            }
            bezier_impl_3d_into_2d!{QuadraticBezier3 QuadraticBezier2}
        }


        mod cubic_bezier2 {
            use super::*;
            bezier_impl_cubic!{
                /// A 2D Bézier curve with two control points.
                ///
                /// 2x2 and 3x3 matrices can be multiplied by a Bézier curve to transform all of its points.
                2 QuadraticBezier2 CubicBezier2 Vec2 LineSegment2
            }
            bezier_impl_2d_into_3d!{CubicBezier2 CubicBezier3}
        }


        mod cubic_bezier3 {
            use super::*;
            bezier_impl_cubic!{
                /// A 3D Bézier curve with two control points.
                ///
                /// 3x3 and 4x4 matrices can be multiplied by a Bézier curve to transform all of its points.
                3 QuadraticBezier3 CubicBezier3 Vec3 LineSegment3
            }
            bezier_impl_3d_into_2d!{CubicBezier3 CubicBezier2}
        }

        pub use quadratic_bezier2::*;
        pub use quadratic_bezier3::*;
        pub use cubic_bezier2::*;
        pub use cubic_bezier3::*;
    };
}

#[cfg(all(nightly, feature="repr_simd"))]
pub mod repr_simd {
    //! Bézier curve structs that use `#[repr(simd)]` vectors.
    use super::*;
    impl_all_beziers!{repr_simd}
}
pub mod repr_c {
    //! Bézier curve structs that use `#[repr(C)]` vectors.
    use super::*;
    impl_all_beziers!{repr_c}
}

pub use self::repr_c::*;

#[cfg(test)]
mod tests {
    use super::*;
    use crate::vec::{Vec2, Vec3};

    macro_rules! test {
        ($Bezier:ident $bezier:ident $Vec:ident) => {
            mod $bezier {
                use super::*;
                #[test] fn lerp_from_line_segment() {
                    let count = 32;
                    let t_iter = (0..(count+1)).into_iter().map(|i| i as f32 / (count as f32));
                    let l = || $Vec::<f32>::unit_x() .. $Vec::<f32>::unit_y();
                    let c = $Bezier::from(l());
                    for t in t_iter {
                        assert_relative_eq!(c.evaluate(t), Lerp::lerp_unclamped_precise(l().start, l().end, t))
                    }
                }
            }
        };
    }
    test!{QuadraticBezier2 quadratic2 Vec2}
    test!{QuadraticBezier3 quadratic3 Vec3}
    test!{CubicBezier2 cubic2 Vec2}
    test!{CubicBezier3 cubic3 Vec3}
}