+//! C Interfaces and Implementations: Techniques for Creating Reusable Software [Hanson 1996-08-30]

+//!

+

+use std::cmp::Ordering;

+#[cfg(test)]

+mod test;

+/// An extended-iprecision unsigned integer is represented in base by an

+/// array of n digits, least significant digit first. Most XP functions take n as

+/// an argument along with source and destination Ts; it is SAD for n<1 or

+/// for n not to be the length of the corresponding Ts. It is SAD to pass a

+/// null T or a T that is too small to any XP function.

+pub type XP = [u8];

+

+macro_rules! base {

+ () => {

+ 256

+ };

+}

+/// sets z[0..n-1] to u mod 2⁸ⁿ and returns u/2⁸ⁿ.

+pub const fn from_int(n: u32, z: &mut XP, mut u: u128) -> u128 {

+ let mut i = 0;

+ let n = n as usize;

+ while {

+ z[i] = (u % base!()) as u8;

+ i += 1;

+

+ u /= base!();

+ u > 0 && i < n

+ } {}

+ while i < n {

+ z[i] = 0;

+ i += 1;

+ }

+ u

+}

+/// returns x mod (ULONG_MAX+1).

+pub const fn to_int(n: u32, x: &XP) -> u128 {

+ let mut u = 0;

+ let mut i = size_of::<u128>().min(n as usize);

+

+ while i > 0 {

+ i -= 1;

+ u = base!() * u + x[i] as u128;

+ }

+ u

+}

+/// returns the length of x; that is, the index plus one of the most signifi-

+/// cant nonzero digit in x[0..n-1].

+pub const fn len(mut n: u32, x: &XP) -> u32 {

+ while (n > 1 && x[n as usize - 1] == 0) {

+ n -= 1;

+ }

+ n

+}

+macro_rules! u1 {

+ ($x:expr) => {

+ match $x {

+ 0 => false,

+ 1 => true,

+ _ => panic!("not u1"),

+ }

+ };

+}

+

+/// sets z[0..n-1] to x + y + carry and returns the carry out of z[n-1].

+pub const fn add(n: u32, z: &mut XP, x: &XP, y: &XP, carry: bool) -> bool {

+ let mut carry = carry as usize;

+ let mut i = 0;

+ while i < n as usize {

+ carry += x[i] as usize + y[i] as usize;

+ z[i] = (carry % base!()) as u8;

+ carry /= base!();

+

+ i += 1;

+ }

+ u1!(carry)

+}

+/// sets z[0..n-1] to x − y − borrow and returns the borrow into z[n-1].

+pub const fn sub(

+ n: u32,

+ z: &mut XP,

+ x: &XP,

+ y: &XP,

+ borrow: bool,

+) -> bool {

+ let mut i = 0;

+ let mut borrow = borrow as u32;

+ while i < n as usize {

+ let mut d = (x[i] as u32 + base!()) - borrow - y[i] as u32;

+ z[i] = (d % base!()) as u8;

+ borrow = 1 - (d / base!());

+ i += 1;

+ }

+ u1!(borrow)

+}

+#[doc(hidden)]

+pub const unsafe fn subp(

+ n: u32,

+ z: *mut u8,

+ x: *const u8,

+ y: *const u8,

+ borrow: bool,

+) -> bool {

+ let mut i = 0;

+ let mut borrow = borrow as u32;

+ while i < n as usize {

+ let mut d =

+ (*x.add(i) as u32 + base!()) - borrow - *y.add(i) as u32;

+ *z.add(i) = (d % base!()) as u8;

+ borrow = 1 - (d / base!());

+ i += 1;

+ }

+ u1!(borrow)

+}

+/// sets z[0..n-1] to x + y, where y is a single digit, and returns the carry-out of z[n-1].

+pub const fn sum(n: u32, z: &mut XP, x: &XP, mut y: u8) -> u8 {

+ unsafe { sum_p(n, z.as_mut_ptr(), x.as_ptr(), y) }

+}

+#[doc(hidden)]

+pub const unsafe fn sum_p(

+ n: u32,

+ z: *mut u8,

+ x: *const u8,

+ mut y: u8,

+) -> u8 {

+ let mut i = 0;

+ let mut y = y as u32;

+ while i < n as usize {

+ y += *x.add(i) as u32;

+ *z.add(i) = (y % base!()) as u8;

+ y /= base!();

+ i += 1;

+ }

+ y as u8

+}

+/// sets z[0..n-1] to x − y, where y is a single digit, and returns the borrow into z[n-1].

+pub const fn diff(n: u32, z: &mut XP, x: &XP, y: u8) -> u8 {

+ let mut y = y as u32;

+ let mut i = 0;

+ while i < n as usize {

+ let d = (x[i] as u32 + base!()) - y;

+ z[i] = (d % base!()) as u8;

+ y = 1 - d / base!();

+ i += 1;

+ }

+ y as u8

+}

+/// sets z[0..n-1] to !x + carry, where carry is zero or one, and returns

+/// the carry-out of z[n-1].

+pub const fn neg(n: u32, z: &mut XP, x: &XP, carry: bool) -> bool {

+ let mut carry = carry as u32;

+ let mut i = 0;

+ while i < n as usize {

+ carry += (!x[i]) as u32;

+ z[i] = (carry % base!()) as u8;

+ carry /= base!();

+ i += 1;

+ }

+ u1!(carry)

+}

+/// adds x[0..n-1]•y[0..m-1] to z[0..n+m-1] and returns the carry-out of

+/// z[n+m-1]. If z=0, Xmul computes x•y. It is SAD for z to be the same T as x or y.

+pub const fn mul(z: &mut XP, n: u32, x: &XP, m: u32, y: &XP) -> u32 {

+ let (mut i, mut carryout) = (0, 0);

+ while i < n as usize {

+ let mut carry = 0;

+ let mut j = 0;

+ while j < m as usize {

+ carry += x[i as usize] as u32 * y[j as usize] as u32

+ + z[i + j] as u32;

+ z[i + j] = (carry % base!()) as u8;

+ carry /= base!();

+ j += 1;

+ }

+ while j < (n as usize + m as usize - i) {

+ carry += z[i + j] as u32;

+ z[i + j] = (carry % base!()) as u8;

+ carry /= base!();

+ j += 1;

+ }

+ assert!(carry == 0 || carry == 1);

+ carryout |= carry;

+ i += 1;

+ }

+ carryout

+}

+

+const unsafe fn quot_p(n: u32, z: *mut u8, x: *const u8, y: u8) -> u8 {

+ let mut i = n as isize - 1;

+ let mut carry = 0;

+

+ while i >= 0 {

+ carry = carry * base!() + (*x.add(i as _)) as u32;

+ *z.add(i as _) = (carry / y as u32) as _;

+ carry %= y as u32;

+ i -= 1;

+ }

+ carry as u8

+}

+/// sets z[0..n-1] to x/y, where y is a single digit, and returns x mod y. It

+/// is SAD for y=0

+pub const fn quot(n: u32, z: &mut XP, x: &XP, y: u8) -> u8 {

+ unsafe { quot_p(n, z.as_mut_ptr(), x.as_ptr(), y) }

+}

+/// XP_product sets z[0..n-1] to x•y and returns the carry-out of the

+/// most significant digit; the carry can be as large as 2^8- XP_quotient

+/// sets z[0..n-1] to x/y and returns the remainder, x mod y; the remainder

+/// can be as large as . For XP_product and XP_quotient, y must not

+/// exceed − 1.

+pub const fn product(n: u32, z: &mut XP, x: &XP, y: u8) -> u8 {

+ unsafe { product_p(n, z.as_mut_ptr(), x.as_ptr(), y) }

+}

+const unsafe fn product_p(n: u32, z: *mut u8, x: *const u8, y: u8) -> u8 {

+ let mut i = 0;

+ let mut carry = 0;

+ while i < n as usize {

+ carry += *x.add(i) as u32 * y as u32;

+ *z.add(i) = (carry % base!()) as u8;

+ carry /= base!();

+ i += 1;

+ }

+ carry.try_into().ok().unwrap()

+}

+/// sets q[0..n-1] to x[0..n-1]/y[0..m-1], r[0..m-1] to x[0..n-1] mod

+/// y[0..m-1]. If y=0, panics.

+/// tmp must hold at least n+m+2 digits. It is unexpected

+/// for q or r to be one of x or y, for q and r to be the same T, or for

+/// tmp to be too small

+pub const fn div(

+ mut n: u32,

+ q: &mut XP,

+ x: &XP,

+ mut m: u32,

+ y: &XP,

+ r: &mut XP,

+ tmp: &mut XP,

+) {

+ let nx = n;

+ let my = m;

+ n = len(n, x);

+ m = len(m, y);

+ if m == 1 {

+ if (y[0] == 0) {

+ panic!()

+ }

+ r[0] = quot(nx, q, x, y[0] as _);

+ let mut i = 0;

+ while i < my as usize - 1 {

+ r[1 + i] = 0;

+ i += 1;

+ }

+ } else if m > n {

+ let mut i = 0;

+ while i < nx as usize {

+ q[i] = 0;

+ i += 1;

+ }

+ r[..n as usize].copy_from_slice(&x[..n as usize]);

+ let mut i = 0;

+ while i < my - n {

+ r[(n + i) as usize] = 0;

+ i += 1;

+ }

+ } else {

+ let (rem, dq) = tmp.split_at_mut(n as usize + 1);

+ // let rem = &mut tmp[..n as usize];

+ // let dq = &mut tmp[n as usize + 1..];

+ assert!(2 <= m && m <= n);

+ rem[..n as usize].copy_from_slice(&x[..n as usize]);

+ rem[n as usize] = 0;

+ let mut k = (n - m) as i32;

+ while k >= 0 {

+ let k_ = k as u32;

+ let qk = {

+ let i = 0;

+ assert!(2 <= m && m <= k_ + m && k_ + m <= n);

+

+ let mut qk = {

+ let km = (k_ + m) as usize;

+ let y2 = y[m as usize - 1] as u64 * base!() as u64

+ + y[m as usize - 2] as u64;

+ let r3 = rem[km] as u64 * (base!() * base!())

+ + rem[km - 1] as u64 * base!()

+ + rem[km - 2] as u64;

+ let mut qk = r3 / y2;

+ if (qk >= base!()) {

+ qk = base!() - 1

+ };

+ qk as u8

+ };

+ // qk y[m-2..m-1]/rem[k+m-2..k+m] 314

+ dq[m as usize] = product(m, dq, y, qk) as _;

+ let mut i = m as i32;

+ // ch

+ while i > 0 && rem[(i + k) as usize] == dq[i as usize] {

+ i -= 1;

+ }

+

+ if (rem[(i + k) as usize] < dq[i as usize]) {

+ qk -= 1;

+ dq[m as usize] = product(m, dq, y, qk) as _;

+ }

+ qk

+ };

+ q[k as usize] = qk;

+ {

+ assert!(0 <= k && k <= k + m as i32);

+ let borrow = unsafe {

+ subp(

+ m + 1,

+ rem.as_mut_ptr().add(k as _),

+ rem.as_ptr().add(k as _),

+ dq.as_ptr(),

+ false,

+ )

+ };

+ assert!(!borrow);

+ }

+ k -= 1;

+ }

+ r[..m as usize].copy_from_slice(&rem[..m as usize]);

+ let mut i = n - m + 1;

+ while i < nx {

+ q[i as usize] = 0;

+ i += 1;

+ }

+ let mut i = m;

+ while i < my {

+ r[i as usize] = 0;

+ i += 1;

+ }

+ r[..m as usize].copy_from_slice(&rem[..m as usize]);

+ // memcpy(r, rem, m);

+ }

+}

+

+/// sets z[0..n-1] to x[0..m-1] shifted left by s bits, and fills the vacated

+/// bits with fill, which must be zero or one. It is SADk for s<0.

+pub const fn shl(

+ n: u32,

+ z: &mut XP,

+ m: u32,

+ x: &XP,

+ mut s: u32,

+ fill: u8,

+) {

+ assert!(matches!(fill, 0 | 0xff));

+ let mut j = n as i32 - 1;

+ let mut i = (if (n > m) { m - 1 } else { n - s / 8 - 1 }) as i32;

+ while j >= (m + s / 8) as _ {

+ z[j as usize] = 0;

+ j -= 1;

+ }

+ //const fill when

+ while i >= 0 {

+ z[j as usize] = x[i as usize];

+ j -= 1;

+ i -= 1;

+ }

+ while j >= 0 {

+ z[j as usize] = fill;

+ j -= 1;

+ }

+ s %= 8;

+ if s > 0 {

+ unsafe { product_p(n, z.as_mut_ptr(), z.as_ptr(), 1 << s) };

+ z[0] |= fill >> (8 - s);

+ }

+}

+//// shifts right; see [`shl`]. If n>m, the excess bits are treated as if

+/// they were equal to fill

+pub const fn shr(

+ n: u32,

+ z: &mut XP,

+ m: u32,

+ x: &XP,

+ s: u32,

+ fill: u8,

+) -> u8 {

+ assert!(matches!(fill, 0 | 0xff));

+ let mut i = s / 8;

+ let mut j = 0;

+ while i < m && j < n {

+ z[j as usize] = x[i as usize];

+ j += 1;

+ i += 1;

+ }

+ while j < n {

+ z[j as usize] = fill;

+ j += 1;

+ }

+ let s = s % 8;

+ if s > 0 {

+ let r = unsafe { quot_p(n, z.as_mut_ptr(), z.as_ptr(), 1 << s) };

+ z[n as usize - 1] |= fill << (8 - s);

+ r

+ } else {

+ 0

+ }

+}

+/// interprets str as an unsigned integer in base using z[0..n-1] as the

+/// initial value in the conversion, and returns the first nonzero carry-out

+/// of the conversion step. If end≠null, *end points to the character in

+/// str that terminated the scan or produced a nonzero carry.

+fn from_str(n: u32, z: &mut XP, base: u8, str: &[u8]) -> u8 {

+ assert!(base < 36);

+ let mut i = 0;

+ while i < str.len() {

+ let b = str[i];

+ let b = match b {

+ b'0'..=b'9' => b - b'0',

+ b'a'..=b'z' => b - b'a',

+ _ => return 255,

+ };

+ if b > base {

+ return 255;

+ }

+ let c = unsafe { product_p(n, z.as_mut_ptr(), z.as_ptr(), base) };

+ unsafe { sum_p(n, z.as_mut_ptr(), z.as_ptr(), b as _) };

+

+ if c != 0 {

+ return c;

+ }

+

+ i += 1;

+ }

+ 0

+}

+/// fills str[0..size-1] with the character representation of x in base,

+/// sets x to zero, and returns str. It is a panic for size to be

+/// too small, or for base<2 or base>36.

+pub const fn to_str<'a>(

+ put_in: &'a mut Vec<u8>,

+ base: u8,

+ mut n: u32,

+ x: &mut XP,

+) -> &'a [u8] {

+ let mut i = 0;

+ assert!(base >= 2 && base <= 36);

+ while {

+ let r = unsafe { quot_p(n, x.as_mut_ptr(), x.as_ptr(), base) };

+ put_in.push(b"0123456789abcdefghijklmnopqrstuvxyz"[r as usize]);

+ i += 1;

+ // str[i] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"[r];

+ while n > 1 && x[n as usize - 1] == 0 {

+ n -= 1;

+ }

+ n > 1 || x[0] != 0

+ } {}

+

+ put_in.as_mut_slice().reverse();

+ put_in.as_slice()

+}

+

+/// x[..n] <=> y[..n]

+pub const fn cmp(n: u32, x: &XP, y: &XP) -> Ordering {

+ let mut i = n as i32 - 1;

+ while (i > 0 && x[i as usize] == y[i as usize]) {

+ i -= 1

+ }

+ x[i as usize].cmp(&y[i as usize])

+}