use super::addition::{__add2, add2};

use super::subtraction::sub2;

use super::{biguint_from_vec, cmp_slice, BigUint, IntDigits};

use crate::big_digit::{self, BigDigit, DoubleBigDigit};

use crate::Sign::{self, Minus, NoSign, Plus};

use crate::{BigInt, UsizePromotion};

use core::cmp::Ordering;

use core::iter::Product;

use core::ops::{Mul, MulAssign};

use num_traits::{CheckedMul, FromPrimitive, One, Zero};

#[inline]

pub(super) fn mac_with_carry(

    a: BigDigit,

    b: BigDigit,

    c: BigDigit,

    acc: &mut DoubleBigDigit,

) -> BigDigit {

    *acc += DoubleBigDigit::from(a);

    *acc += DoubleBigDigit::from(b) * DoubleBigDigit::from(c);

    let lo = *acc as BigDigit;

    *acc >>= big_digit::BITS;

    lo

}

#[inline]

fn mul_with_carry(a: BigDigit, b: BigDigit, acc: &mut DoubleBigDigit) -> BigDigit {

    *acc += DoubleBigDigit::from(a) * DoubleBigDigit::from(b);

    let lo = *acc as BigDigit;

    *acc >>= big_digit::BITS;

    lo

}

/// Three argument multiply accumulate:

/// acc += b * c

fn mac_digit(acc: &mut [BigDigit], b: &[BigDigit], c: BigDigit) {

    if c == 0 {

        return;

    }

    let mut carry = 0;

    let (a_lo, a_hi) = acc.split_at_mut(b.len());

    for (a, &b) in a_lo.iter_mut().zip(b) {

        *a = mac_with_carry(*a, b, c, &mut carry);

    }

    let (carry_hi, carry_lo) = big_digit::from_doublebigdigit(carry);

    let final_carry = if carry_hi == 0 {

        __add2(a_hi, &[carry_lo])

    } else {

        __add2(a_hi, &[carry_hi, carry_lo])

    };

    assert_eq!(final_carry, 0, "carry overflow during multiplication!");

}

fn bigint_from_slice(slice: &[BigDigit]) -> BigInt {

    BigInt::from(biguint_from_vec(slice.to_vec()))

}

/// Three argument multiply accumulate:

/// acc += b * c

#[allow(clippy::many_single_char_names)]

fn mac3(mut acc: &mut [BigDigit], mut b: &[BigDigit], mut c: &[BigDigit]) {

    // Least-significant zeros have no effect on the output.

    if let Some(&0) = b.first() {

        if let Some(nz) = b.iter().position(|&d| d != 0) {

            b = &b[nz..];

            acc = &mut acc[nz..];

        } else {

            return;

        }

    }

    if let Some(&0) = c.first() {

        if let Some(nz) = c.iter().position(|&d| d != 0) {

            c = &c[nz..];

            acc = &mut acc[nz..];

        } else {

            return;

        }

    }

    let acc = acc;

    let (x, y) = if b.len() < c.len() { (b, c) } else { (c, b) };

    // We use four algorithms for different input sizes.

    //

    // - For small inputs, long multiplication is fastest.

    // - If y is at least least twice as long as x, split using Half-Karatsuba.

    // - Next we use Karatsuba multiplication (Toom-2), which we have optimized

    //   to avoid unnecessary allocations for intermediate values.

    // - For the largest inputs we use Toom-3, which better optimizes the

    //   number of operations, but uses more temporary allocations.

    //

    // The thresholds are somewhat arbitrary, chosen by evaluating the results

    // of `cargo bench --bench bigint multiply`.

    if x.len() <= 32 {

        // Long multiplication:

        for (i, xi) in x.iter().enumerate() {

            mac_digit(&mut acc[i..], y, *xi);

        }

    } else if x.len() * 2 <= y.len() {

        // Karatsuba Multiplication for factors with significant length disparity.

        //

        // The Half-Karatsuba Multiplication Algorithm is a specialized case of

        // the normal Karatsuba multiplication algorithm, designed for the scenario

        // where y has at least twice as many base digits as x.

        //

        // In this case y (the longer input) is split into high2 and low2,

        // at m2 (half the length of y) and x (the shorter input),

        // is used directly without splitting.

        //

        // The algorithm then proceeds as follows:

        //

        // 1. Compute the product z0 = x * low2.

        // 2. Compute the product temp = x * high2.

        // 3. Adjust the weight of temp by adding m2 (* NBASE ^ m2)

        // 4. Add temp and z0 to obtain the final result.

        //

        // Proof:

        //

        // The algorithm can be derived from the original Karatsuba algorithm by

        // simplifying the formula when the shorter factor x is not split into

        // high and low parts, as shown below.

        //

        // Original Karatsuba formula:

        //

        //     result = (z2 * NBASE ^ (m2 × 2)) + ((z1 - z2 - z0) * NBASE ^ m2) + z0

        //

        // Substitutions:

        //

        //     low1 = x

        //     high1 = 0

        //

        // Applying substitutions:

        //

        //     z0 = (low1 * low2)

        //        = (x * low2)

        //

        //     z1 = ((low1 + high1) * (low2 + high2))

        //        = ((x + 0) * (low2 + high2))

        //        = (x * low2) + (x * high2)

        //

        //     z2 = (high1 * high2)

        //        = (0 * high2)

        //        = 0

        //

        // Simplified using the above substitutions:

        //

        //     result = (z2 * NBASE ^ (m2 × 2)) + ((z1 - z2 - z0) * NBASE ^ m2) + z0

        //            = (0 * NBASE ^ (m2 × 2)) + ((z1 - 0 - z0) * NBASE ^ m2) + z0

        //            = ((z1 - z0) * NBASE ^ m2) + z0

        //            = ((z1 - z0) * NBASE ^ m2) + z0

        //            = (x * high2) * NBASE ^ m2 + z0

        let m2 = y.len() / 2;

        let (low2, high2) = y.split_at(m2);

        // (x * high2) * NBASE ^ m2 + z0

        mac3(acc, x, low2);

        mac3(&mut acc[m2..], x, high2);

    } else if x.len() <= 256 {

        // Karatsuba multiplication:

        //

        // The idea is that we break x and y up into two smaller numbers that each have about half

        // as many digits, like so (note that multiplying by b is just a shift):

        //

        // x = x0 + x1 * b

        // y = y0 + y1 * b

        //

        // With some algebra, we can compute x * y with three smaller products, where the inputs to

        // each of the smaller products have only about half as many digits as x and y:

        //

        // x * y = (x0 + x1 * b) * (y0 + y1 * b)

        //

        // x * y = x0 * y0

        //       + x0 * y1 * b

        //       + x1 * y0 * b

        //       + x1 * y1 * b^2

        //

        // Let p0 = x0 * y0 and p2 = x1 * y1:

        //

        // x * y = p0

        //       + (x0 * y1 + x1 * y0) * b

        //       + p2 * b^2

        //

        // The real trick is that middle term:

        //

        //         x0 * y1 + x1 * y0

        //

        //       = x0 * y1 + x1 * y0 - p0 + p0 - p2 + p2

        //

        //       = x0 * y1 + x1 * y0 - x0 * y0 - x1 * y1 + p0 + p2

        //

        // Now we complete the square:

        //

        //       = -(x0 * y0 - x0 * y1 - x1 * y0 + x1 * y1) + p0 + p2

        //

        //       = -((x1 - x0) * (y1 - y0)) + p0 + p2

        //

        // Let p1 = (x1 - x0) * (y1 - y0), and substitute back into our original formula:

        //

        // x * y = p0

        //       + (p0 + p2 - p1) * b

        //       + p2 * b^2

        //

        // Where the three intermediate products are:

        //

        // p0 = x0 * y0

        // p1 = (x1 - x0) * (y1 - y0)

        // p2 = x1 * y1

        //

        // In doing the computation, we take great care to avoid unnecessary temporary variables

        // (since creating a BigUint requires a heap allocation): thus, we rearrange the formula a

        // bit so we can use the same temporary variable for all the intermediate products:

        //

        // x * y = p2 * b^2 + p2 * b

        //       + p0 * b + p0

        //       - p1 * b

        //

        // The other trick we use is instead of doing explicit shifts, we slice acc at the

        // appropriate offset when doing the add.

        // When x is smaller than y, it's significantly faster to pick b such that x is split in

        // half, not y:

        let b = x.len() / 2;

        let (x0, x1) = x.split_at(b);

        let (y0, y1) = y.split_at(b);

        // We reuse the same BigUint for all the intermediate multiplies and have to size p

        // appropriately here: x1.len() >= x0.len and y1.len() >= y0.len():

        let len = x1.len() + y1.len() + 1;

        let mut p = BigUint { data: vec![0; len] };

        // p2 = x1 * y1

        mac3(&mut p.data, x1, y1);

        // Not required, but the adds go faster if we drop any unneeded 0s from the end:

        p.normalize();

        add2(&mut acc[b..], &p.data);

        add2(&mut acc[b * 2..], &p.data);

        // Zero out p before the next multiply:

        p.data.truncate(0);

        p.data.resize(len, 0);

        // p0 = x0 * y0

        mac3(&mut p.data, x0, y0);

        p.normalize();

        add2(acc, &p.data);

        add2(&mut acc[b..], &p.data);

        // p1 = (x1 - x0) * (y1 - y0)

        // We do this one last, since it may be negative and acc can't ever be negative:

        let (j0_sign, j0) = sub_sign(x1, x0);

        let (j1_sign, j1) = sub_sign(y1, y0);

        match j0_sign * j1_sign {

            Plus => {

                p.data.truncate(0);

                p.data.resize(len, 0);

                mac3(&mut p.data, &j0.data, &j1.data);

                p.normalize();

                sub2(&mut acc[b..], &p.data);

            }

            Minus => {

                mac3(&mut acc[b..], &j0.data, &j1.data);

            }

            NoSign => (),

        }

    } else {

        // Toom-3 multiplication:

        //

        // Toom-3 is like Karatsuba above, but dividing the inputs into three parts.

        // Both are instances of Toom-Cook, using `k=3` and `k=2` respectively.

        //

        // The general idea is to treat the large integers digits as

        // polynomials of a certain degree and determine the coefficients/digits

        // of the product of the two via interpolation of the polynomial product.

        let i = y.len() / 3 + 1;

        let x0_len = Ord::min(x.len(), i);

        let x1_len = Ord::min(x.len() - x0_len, i);

        let y0_len = i;

        let y1_len = Ord::min(y.len() - y0_len, i);

        // Break x and y into three parts, representating an order two polynomial.

        // t is chosen to be the size of a digit so we can use faster shifts

        // in place of multiplications.

        //

        // x(t) = x2*t^2 + x1*t + x0

        let x0 = bigint_from_slice(&x[..x0_len]);

        let x1 = bigint_from_slice(&x[x0_len..x0_len + x1_len]);

        let x2 = bigint_from_slice(&x[x0_len + x1_len..]);

        // y(t) = y2*t^2 + y1*t + y0

        let y0 = bigint_from_slice(&y[..y0_len]);

        let y1 = bigint_from_slice(&y[y0_len..y0_len + y1_len]);

        let y2 = bigint_from_slice(&y[y0_len + y1_len..]);

        // Let w(t) = x(t) * y(t)

        //

        // This gives us the following order-4 polynomial.

        //

        // w(t) = w4*t^4 + w3*t^3 + w2*t^2 + w1*t + w0

        //

        // We need to find the coefficients w4, w3, w2, w1 and w0. Instead

        // of simply multiplying the x and y in total, we can evaluate w

        // at 5 points. An n-degree polynomial is uniquely identified by (n + 1)

        // points.

        //

        // It is arbitrary as to what points we evaluate w at but we use the

        // following.

        //

        // w(t) at t = 0, 1, -1, -2 and inf

        //

        // The values for w(t) in terms of x(t)*y(t) at these points are:

        //

        // let a = w(0)   = x0 * y0

        // let b = w(1)   = (x2 + x1 + x0) * (y2 + y1 + y0)

        // let c = w(-1)  = (x2 - x1 + x0) * (y2 - y1 + y0)

        // let d = w(-2)  = (4*x2 - 2*x1 + x0) * (4*y2 - 2*y1 + y0)

        // let e = w(inf) = x2 * y2 as t -> inf

        // x0 + x2, avoiding temporaries

        let p = &x0 + &x2;

        // y0 + y2, avoiding temporaries

        let q = &y0 + &y2;

        // x2 - x1 + x0, avoiding temporaries

        let p2 = &p - &x1;

        // y2 - y1 + y0, avoiding temporaries

        let q2 = &q - &y1;

        // w(0)

        let r0 = &x0 * &y0;

        // w(inf)

        let r4 = &x2 * &y2;

        // w(1)

        let r1 = (p + x1) * (q + y1);

        // w(-1)

        let r2 = &p2 * &q2;

        // w(-2)

        let r3 = ((p2 + x2) * 2 - x0) * ((q2 + y2) * 2 - y0);

        // Evaluating these points gives us the following system of linear equations.

        //

        //  0  0  0  0  1 | a

        //  1  1  1  1  1 | b

        //  1 -1  1 -1  1 | c

        // 16 -8  4 -2  1 | d

        //  1  0  0  0  0 | e

        //

        // The solved equation (after gaussian elimination or similar)

        // in terms of its coefficients:

        //

        // w0 = w(0)

        // w1 = w(0)/2 + w(1)/3 - w(-1) + w(-2)/6 - 2*w(inf)

        // w2 = -w(0) + w(1)/2 + w(-1)/2 - w(inf)

        // w3 = -w(0)/2 + w(1)/6 + w(-1)/2 - w(-2)/6 + 2*w(inf)

        // w4 = w(inf)

        //

        // This particular sequence is given by Bodrato and is an interpolation

        // of the above equations.

        let mut comp3: BigInt = (r3 - &r1) / 3u32;

        let mut comp1: BigInt = (r1 - &r2) >> 1;

        let mut comp2: BigInt = r2 - &r0;

        comp3 = ((&comp2 - comp3) >> 1) + (&r4 << 1);

        comp2 += &comp1 - &r4;

        comp1 -= &comp3;

        // Recomposition. The coefficients of the polynomial are now known.

        //

        // Evaluate at w(t) where t is our given base to get the result.

        //

        //     let bits = u64::from(big_digit::BITS) * i as u64;

        //     let result = r0

        //         + (comp1 << bits)

        //         + (comp2 << (2 * bits))

        //         + (comp3 << (3 * bits))

        //         + (r4 << (4 * bits));

        //     let result_pos = result.to_biguint().unwrap();

        //     add2(&mut acc[..], &result_pos.data);

        //

        // But with less intermediate copying:

        for (j, result) in [&r0, &comp1, &comp2, &comp3, &r4].iter().enumerate().rev() {

            match result.sign() {

                Plus => add2(&mut acc[i * j..], result.digits()),

                Minus => sub2(&mut acc[i * j..], result.digits()),

                NoSign => {}

            }

        }

    }

}

fn mul3(x: &[BigDigit], y: &[BigDigit]) -> BigUint {

    let len = x.len() + y.len() + 1;

    let mut prod = BigUint { data: vec![0; len] };

    mac3(&mut prod.data, x, y);

    prod.normalized()

}

fn scalar_mul(a: &mut BigUint, b: BigDigit) {

    match b {

        0 => a.set_zero(),

        1 => {}

        _ => {

            if b.is_power_of_two() {

                *a <<= b.trailing_zeros();

            } else {

                let mut carry = 0;

                for a in a.data.iter_mut() {

                    *a = mul_with_carry(*a, b, &mut carry);

                }

                if carry != 0 {

                    a.data.push(carry as BigDigit);

                }

            }

        }

    }

}

fn sub_sign(mut a: &[BigDigit], mut b: &[BigDigit]) -> (Sign, BigUint) {

    // Normalize:

    if let Some(&0) = a.last() {

        a = &a[..a.iter().rposition(|&x| x != 0).map_or(0, |i| i + 1)];

    }

    if let Some(&0) = b.last() {

        b = &b[..b.iter().rposition(|&x| x != 0).map_or(0, |i| i + 1)];

    }

    match cmp_slice(a, b) {

        Ordering::Greater => {

            let mut a = a.to_vec();

            sub2(&mut a, b);

            (Plus, biguint_from_vec(a))

        }

        Ordering::Less => {

            let mut b = b.to_vec();

            sub2(&mut b, a);

            (Minus, biguint_from_vec(b))

        }

        Ordering::Equal => (NoSign, BigUint::ZERO),

    }

}

macro_rules! impl_mul {

    ($(impl Mul<$Other:ty> for $Self:ty;)*) => {$(

        impl Mul<$Other> for $Self {

            type Output = BigUint;

            #[inline]

            fn mul(self, other: $Other) -> BigUint {

                match (&*self.data, &*other.data) {

                    // multiply by zero

                    (&[], _) | (_, &[]) => BigUint::ZERO,

                    // multiply by a scalar

                    (_, &[digit]) => self * digit,

                    (&[digit], _) => other * digit,

                    // full multiplication

                    (x, y) => mul3(x, y),

                }

            }

Unexecuted instantiation: <num_bigint::biguint::BigUint as core::ops::arith::Mul>::mul

Unexecuted instantiation: <num_bigint::biguint::BigUint as core::ops::arith::Mul<&num_bigint::biguint::BigUint>>::mul

Unexecuted instantiation: <&num_bigint::biguint::BigUint as core::ops::arith::Mul>::mul

Unexecuted instantiation: <&num_bigint::biguint::BigUint as core::ops::arith::Mul<num_bigint::biguint::BigUint>>::mul

        }

    )*}

}

impl_mul! {

    impl Mul<BigUint> for BigUint;

    impl Mul<BigUint> for &BigUint;

    impl Mul<&BigUint> for BigUint;

    impl Mul<&BigUint> for &BigUint;

}

macro_rules! impl_mul_assign {

    ($(impl MulAssign<$Other:ty> for BigUint;)*) => {$(

        impl MulAssign<$Other> for BigUint {

            #[inline]

            fn mul_assign(&mut self, other: $Other) {

                match (&*self.data, &*other.data) {

                    // multiply by zero

                    (&[], _) => {},

                    (_, &[]) => self.set_zero(),

                    // multiply by a scalar

                    (_, &[digit]) => *self *= digit,

                    (&[digit], _) => *self = other * digit,

                    // full multiplication

                    (x, y) => *self = mul3(x, y),

                }

            }

Unexecuted instantiation: <num_bigint::biguint::BigUint as core::ops::arith::MulAssign<&num_bigint::biguint::BigUint>>::mul_assign

Unexecuted instantiation: <num_bigint::biguint::BigUint as core::ops::arith::MulAssign>::mul_assign

        }

    )*}

}

impl_mul_assign! {

    impl MulAssign<BigUint> for BigUint;

    impl MulAssign<&BigUint> for BigUint;

}

promote_unsigned_scalars!(impl Mul for BigUint, mul);

promote_unsigned_scalars_assign!(impl MulAssign for BigUint, mul_assign);

forward_all_scalar_binop_to_val_val_commutative!(impl Mul<u32> for BigUint, mul);

forward_all_scalar_binop_to_val_val_commutative!(impl Mul<u64> for BigUint, mul);

forward_all_scalar_binop_to_val_val_commutative!(impl Mul<u128> for BigUint, mul);

impl Mul<u32> for BigUint {

    type Output = BigUint;

    #[inline]

    fn mul(mut self, other: u32) -> BigUint {

        self *= other;

        self

    }

}

impl MulAssign<u32> for BigUint {

    #[inline]

    fn mul_assign(&mut self, other: u32) {

        scalar_mul(self, other as BigDigit);

    }

}

impl Mul<u64> for BigUint {

    type Output = BigUint;

    #[inline]

    fn mul(mut self, other: u64) -> BigUint {

        self *= other;

        self

    }

}

impl MulAssign<u64> for BigUint {

    cfg_digit!(

        #[inline]

        fn mul_assign(&mut self, other: u64) {

            if let Some(other) = BigDigit::from_u64(other) {

                scalar_mul(self, other);

            } else {

                let (hi, lo) = big_digit::from_doublebigdigit(other);

                *self = mul3(&self.data, &[lo, hi]);

            }

        }

        #[inline]

        fn mul_assign(&mut self, other: u64) {

            scalar_mul(self, other);

        }

    );

}

impl Mul<u128> for BigUint {

    type Output = BigUint;

    #[inline]

    fn mul(mut self, other: u128) -> BigUint {

        self *= other;

        self

    }

}

impl MulAssign<u128> for BigUint {

    cfg_digit!(

        #[inline]

        fn mul_assign(&mut self, other: u128) {

            if let Some(other) = BigDigit::from_u128(other) {

                scalar_mul(self, other);

            } else {

                *self = match super::u32_from_u128(other) {

                    (0, 0, c, d) => mul3(&self.data, &[d, c]),

                    (0, b, c, d) => mul3(&self.data, &[d, c, b]),

                    (a, b, c, d) => mul3(&self.data, &[d, c, b, a]),

                };

            }

        }

        #[inline]

        fn mul_assign(&mut self, other: u128) {

            if let Some(other) = BigDigit::from_u128(other) {

                scalar_mul(self, other);

            } else {

                let (hi, lo) = big_digit::from_doublebigdigit(other);

                *self = mul3(&self.data, &[lo, hi]);

            }

        }

    );

}

impl CheckedMul for BigUint {

    #[inline]

    fn checked_mul(&self, v: &BigUint) -> Option<BigUint> {

        Some(self.mul(v))

    }

}

impl_product_iter_type!(BigUint);

#[test]

fn test_sub_sign() {

    use crate::BigInt;

    use num_traits::Num;

    fn sub_sign_i(a: &[BigDigit], b: &[BigDigit]) -> BigInt {

        let (sign, val) = sub_sign(a, b);

        BigInt::from_biguint(sign, val)

    }

    let a = BigUint::from_str_radix("265252859812191058636308480000000", 10).unwrap();

    let b = BigUint::from_str_radix("26525285981219105863630848000000", 10).unwrap();

    let a_i = BigInt::from(a.clone());

    let b_i = BigInt::from(b.clone());

    assert_eq!(sub_sign_i(&a.data, &b.data), &a_i - &b_i);

    assert_eq!(sub_sign_i(&b.data, &a.data), &b_i - &a_i);

}