def modular_geometric_sum(x, n, mod):
    """ Compute a_n = (1 + a^1 + ... + a^{n-1}) % mod 
    using that 
        a_{2n} = ((x_n + 1) * a_n) % mod
        a_{2n+1} = (x^{2n} + a_{2n}) % mod
    """
    if n == 1:
        return 1 % mod
    elif n % 2 == 0:
        return ((pow(x, n // 2, mod) + 1) * modular_geometric_sum(x, n // 2, mod)) % mod
    else:
        return (pow(x, n - 1, mod) + modular_geometric_sum(x, n - 1, mod)) % mod