testing/bigcalc.py

"""Big number calculator — arbitrary precision integer/decimal arithmetic.

All operations use Python's built-in big integers scaled by 10^n for decimal tracking.
No floating point at any intermediate step for +, -, *, /, pow (integer exponent).

Operations: add, sub, mul, div, pow, sqrt, square, fact, mod, gcd, lcm
"""

import math
import re


def parse_number(s: str):
    """Parse a number string into (sign: int, int_part: str, frac_part: str).

    Supports:
      - Integers: 42, -123
      - Decimals: 3.14, -0.001
      - Scientific notation: 2.5e3, 1.23E-4, -5e10

    Returns None on invalid input.
    """
    s = (s or '').strip()
    if not s:
        return None

    sign = 1
    if s[0] == '-':
        sign = -1
        s = s[1:]
    elif s[0] == '+':
        s = s[1:]
    if not s or s == '.':
        return None

    # Scientific notation
    exp = 0
    m = re.search(r'[eE]', s)
    if m:
        e_str = s[m.end():]
        if not e_str or (e_str[0] not in '+-' and not e_str[0].isdigit()):
            return None
        try:
            exp = int(e_str)
        except ValueError:
            return None
        s = s[:m.start()]
    if not s or s == '.':
        return None

    # Split integer and fractional parts
    if '.' in s:
        int_part, frac_part = s.split('.', 1)
        int_part = int_part.lstrip('0') or '0'
        frac_part = frac_part.rstrip('0')
    else:
        int_part = s.lstrip('0') or '0'
        frac_part = ''

    if not int_part.isdigit() or (frac_part and not frac_part.isdigit()):
        return None

    # Apply exponent (shift decimal point)
    if exp > 0:
        move = min(exp, len(frac_part))
        int_part += frac_part[:move]
        frac_part = frac_part[move:]
        exp -= move
        if exp > 0:
            int_part += '0' * exp
            exp = 0
    elif exp < 0:
        shift = -exp  # number of places to shift left
        # Move as many digits from int_part as possible
        take = min(shift, len(int_part))
        frac_part = int_part[-take:] + frac_part
        int_part = int_part[:-take] or '0'
        shift -= take
        # Remaining shift: prepend zeros to frac_part
        if shift > 0:
            frac_part = '0' * shift + frac_part

    int_part = int_part.lstrip('0') or '0'
    frac_part = frac_part.rstrip('0')
    return (sign, int_part, frac_part)


def to_bigint(n, target_dec):
    """Convert parsed number to Python int scaled by 10^target_dec."""
    sign, int_part, frac_part = n
    s = int_part + (frac_part or '').ljust(target_dec, '0')
    s = s.lstrip('0') or '0'
    val = int(s)
    return -val if sign == -1 else val


def from_bigint(val, dec_places, precision):
    """Format a bigint with decimal tracking back to string.

    Rounds to `precision` digits after the decimal point.
    """
    if val == 0:
        return '0'
    negative = val < 0
    s = str(-val) if negative else str(val)

    if dec_places <= 0:
        if dec_places < 0:
            s += '0' * (-dec_places)
        return ('-' if negative else '') + s

    # Pad with leading zeros
    while len(s) <= dec_places:
        s = '0' + s
    ins = len(s) - dec_places
    int_part = s[:ins] or '0'
    frac_part = s[ins:]

    # Round
    if len(frac_part) > precision:
        round_digit = int(frac_part[precision])
        frac_part = frac_part[:precision]
        if round_digit >= 5:
            new_s = str(int(int_part + frac_part) + 1)
            if len(new_s) > len(int_part) + len(frac_part):
                int_part = new_s[:len(new_s) - len(frac_part)]
                frac_part = new_s[len(new_s) - len(frac_part):]
            else:
                int_part = new_s[:len(new_s) - len(frac_part)] or '0'
                frac_part = new_s[len(new_s) - len(frac_part):]

    frac_part = frac_part.rstrip('0')
    result = int_part
    if frac_part:
        result += '.' + frac_part
    return ('-' if negative else '') + result


# ========== String-based integer helpers for very large numbers ==========

def _str_cmp(a, b):
    """Compare two positive digit strings."""
    if len(a) != len(b):
        return 1 if len(a) > len(b) else -1
    if a == b:
        return 0
    return 1 if a > b else -1


def _str_add(a, b):
    """Add two positive integer strings."""
    carry = 0
    result = []
    i, j = len(a) - 1, len(b) - 1
    while i >= 0 or j >= 0 or carry:
        da = int(a[i]) if i >= 0 else 0
        db = int(b[j]) if j >= 0 else 0
        s = da + db + carry
        result.append(str(s % 10))
        carry = s // 10
        i -= 1; j -= 1
    return ''.join(reversed(result))


def _str_sub(a, b):
    """Subtract b from a (|a| >= |b|), both positive digit strings."""
    borrow = 0
    result = []
    i, j = len(a) - 1, len(b) - 1
    while i >= 0:
        da = int(a[i])
        db = int(b[j]) if j >= 0 else 0
        d = da - db - borrow
        if d < 0:
            d += 10
            borrow = 1
        else:
            borrow = 0
        result.append(str(d))
        i -= 1; j -= 1
    return ''.join(reversed(result)).lstrip('0') or '0'


def _str_mul(a, b):
    """Multiply two positive integer strings."""
    if a == '0' or b == '0':
        return '0'
    res = [0] * (len(a) + len(b))
    for i in range(len(a) - 1, -1, -1):
        for j in range(len(b) - 1, -1, -1):
            p = int(a[i]) * int(b[j]) + res[i + j + 1]
            res[i + j + 1] = p % 10
            res[i + j] += p // 10
    return ''.join(map(str, res)).lstrip('0') or '0'


def _str_divmod(dividend, divisor):
    """Divide two positive integer strings. Returns (quotient, remainder)."""
    if divisor == '0':
        return None
    quot = []
    rem = ''
    for ch in dividend:
        rem += ch
        rem = rem.lstrip('0') or '0'
        count = 0
        while _str_cmp(rem, divisor) >= 0:
            rem = _str_sub(rem, divisor)
            count += 1
        quot.append(str(count))
    q = ''.join(quot).lstrip('0') or '0'
    return (q, rem)


def _gcd_str(a, b):
    """GCD of two positive integer strings."""
    while b != '0':
        _, rem = _str_divmod(a, b)
        a, b = b, rem
        if a is None:
            return None
    return a


def _bigint_sqrt(n):
    """Integer square root using Newton's method (floor)."""
    if n < 0:
        return None
    if n == 0:
        return 0
    x = n
    y = (x + 1) // 2
    while y < x:
        x = y
        y = (y + n // y) // 2
    return x


# ========== Core operations ==========

def add(x_str, y_str, precision=20):
    """X + Y"""
    x = parse_number(x_str)
    y = parse_number(y_str)
    if x is None or y is None:
        return None
    dec = max(len(x[2]), len(y[2]))
    a = to_bigint(x, dec)
    b = to_bigint(y, dec)
    return from_bigint(a + b, dec, precision)


def sub(x_str, y_str, precision=20):
    """X - Y"""
    x = parse_number(x_str)
    y = parse_number(y_str)
    if x is None or y is None:
        return None
    dec = max(len(x[2]), len(y[2]))
    a = to_bigint(x, dec)
    b = to_bigint(y, dec)
    return from_bigint(a - b, dec, precision)


def mul(x_str, y_str, precision=20):
    """X × Y"""
    x = parse_number(x_str)
    y = parse_number(y_str)
    if x is None or y is None:
        return None
    dec_x = len(x[2])
    dec_y = len(y[2])
    a = to_bigint(x, 0)
    b = to_bigint(y, 0)
    return from_bigint(a * b, dec_x + dec_y, precision)


def div(x_str, y_str, precision=20):
    """X / Y"""
    x = parse_number(x_str)
    y = parse_number(y_str)
    if x is None or y is None:
        return None
    if y[1] == '0' and not y[2]:
        return None  # division by zero
    dec_x = len(x[2])
    dec_y = len(y[2])
    scale = dec_y + precision
    a = to_bigint(x, scale)
    b = to_bigint(y, 0)
    if b == 0:
        return None
    result = a // b
    final_dec = scale - dec_x
    sign_mul = 1 if x[0] == y[0] else -1
    return from_bigint(result * sign_mul, final_dec, precision)


def pow_int(x_str, y_str, precision=20):
    """X^Y (Y as integer). Returns None if result too large."""
    x = parse_number(x_str)
    y = parse_number(y_str)
    if x is None or y is None:
        return None
    if y[0] == -1:
        return '0'  # negative exponent → 0 for integer precision
    if y[2]:
        # Decimal exponent — fall back to float approximation
        return pow_float(x_str, y_str, precision)

    # Integer exponent via binary exponentiation on strings
    exp_str = y[1]
    if exp_str == '0':
        return '1'

    if len(exp_str) > 10:
        return None  # exponent too large

    x_sign = x[0]
    x_dec = len(x[2])
    x_str_clean = x[1] + (x[2] or '')

    result = '1'
    base = x_str_clean
    e = exp_str

    while e != '0':
        if int(e[-1]) % 2 == 1:
            # result = result * base (tracking decimals)
            res_clean = result.lstrip('0') or '0'
            result = _str_mul(res_clean, base)
            # Determine the resulting scale and apply it
            res_scale = 0
            base_scale = len(e)  # rough tracking
        # base = base * base
        base = _str_mul(base, base)
        e = _str_divmod(e, '2')[0]

    return result


def pow_float(x_str, y_str, precision=20):
    """X^Y via float approximation (used when Y is decimal)."""
    x = parse_number(x_str)
    y = parse_number(y_str)
    if x is None or y is None:
        return None
    x_val = float(x[0] * int(x[1] + (x[2] or '0')[:15])) / (10 ** len(x[2]))
    y_val = float(y[0] * int(y[1] + (y[2] or '0')[:15])) / (10 ** len(y[2]))
    if x_val <= 0:
        return None
    result = x_val ** y_val
    if not math.isfinite(result):
        return None
    return format(float_result_to_str(result, precision))


def sqrt(x_str, precision=20):
    """√X"""
    x = parse_number(x_str)
    if x is None:
        return None
    if x[0] == -1:
        return None
    x_dec = len(x[2])
    scale = precision * 2
    a = to_bigint(x, scale)
    r = _bigint_sqrt(a)
    if r is None:
        return None
    return from_bigint(r, precision, precision)


def square(x_str, precision=20):
    """X²"""
    return mul(x_str, x_str, precision)


def fact(x_str, precision=0):
    """X! (integer factorial)"""
    x = parse_number(x_str)
    if x is None:
        return None
    if x[0] == -1:
        return None
    if x[2]:
        return None  # integer required
    n = x[1]
    if n in ('0', '1'):
        return '1'

    n_int = int(n)
    if n_int > 100000:
        return None  # too large

    result = '1'
    for i in range(2, n_int + 1):
        result = _str_mul(result, str(i))
    return result


def mod(x_str, y_str, precision=0):
    """X MOD Y"""
    x = parse_number(x_str)
    y = parse_number(y_str)
    if x is None or y is None:
        return None
    if y[1] == '0' and not y[2]:
        return None
    a = int(x[0] * int(x[1]))
    b = int(y[0] * int(y[1]))
    if b == 0:
        return None
    return str(a % b)


def gcd(x_str, y_str, precision=0):
    """GCD of X and Y"""
    x = parse_number(x_str)
    y = parse_number(y_str)
    if x is None or y is None:
        return None
    a = x[1].lstrip('0') or '0'
    b = y[1].lstrip('0') or '0'
    if a == '0' and b == '0':
        return '0'
    if a == '0':
        return b
    if b == '0':
        return a
    return _gcd_str(a, b)


def lcm(x_str, y_str, precision=0):
    """LCM of X and Y"""
    x = parse_number(x_str)
    y = parse_number(y_str)
    if x is None or y is None:
        return None
    a = x[1].lstrip('0') or '0'
    b = y[1].lstrip('0') or '0'
    if a == '0' or b == '0':
        return '0'
    g = _gcd_str(a, b)
    if g is None:
        return None
    prod = _str_mul(a, b)
    q, _ = _str_divmod(prod, g)
    return q