pidigits.R

# Constructors and string conversion functions

zero_mag = integer(0)
zero = c(0L, zero_mag)  # ( = 0L )
one_mag = c(1L)
one = c(1L, one_mag)
ten_mag = c(10L)
ten = c(1L, ten_mag)

elem_max = 10000
elem_digits = as.integer(log10(elem_max))

signum <- function(v) v[[1]]
mag <- function(v) v[2:length(v)]

str_to_mag <- function(s) {
    strip_leading_zeros <- function(s) {
        for (i in 1:nchar(s))
            if (substr(s, i, i) != '0')
                return(substr(s, i, nchar(s)))
        return("")
    }

    len = nchar(s <- strip_leading_zeros(s))
    if (len == 0)
        return(zero_mag)
    mod <- len %% elem_digits
    if (mod > 0) chunks <- substr(s, 1, mod)
    else chunks <- character(0)
    if (1 + mod < len) {
        chunk_inds_begin <- seq(1 + mod, len, elem_digits)
        chunk_inds_end <- chunk_inds_begin + elem_digits - 1
        chunks <- c(chunks, substring(s, chunk_inds_begin, chunk_inds_end))
    }
    return(sapply(chunks, as.integer))
}
hi =
bigint <- function(s, i=NA) {
    if (!is.na(i)) {
        s <- as.character(i)
        #x <- abs(i)
        #sgn <- if (i < 0L) -1L else 1L
        #while (i >= max_elem) {
        #}
     }
#    } else {
        if (substr(s, 1, 1) == '-') {
            mag <- str_to_mag(substr(s, 2L, nchar(s)))
            sgn <- -1L
        } else {
            mag <- str_to_mag(s)
            sgn <- 1L
        }
#    }
    return(c(if (length(mag)) sgn else 0L, mag))
}

to_int <- function(x) {
    ret <- 0L
    if ((x_len = length(x)) > 1L)
        for (i in length(x):2)
            ret <- ret * elem_max + x[[i]]
    return(x[[1]] * ret)
}

mag_to_str <- function(x) {
    len = length(x)
    if (len == 1) return(as.character(x[[1]]))
    xs <- sapply(x[2:length(x)], function(e) zeropad(e, 4))
    x1 <- x[[1]]
    return(paste(collapse="", sep="", c(
                x[[1]],
                sapply(x[2:length(x)], function(e) zeropad(e, elem_digits)))))
}

to_str <- function(x) {
    if (x[[1]] == 0L)
        return("0")
    return(paste(sep="", if (x[[1]] < 0L) '-' else '', mag_to_str(mag(x))))
}

# Comparison functions

cmp_elem <- function(x, y) (x > y) - (x < y)
cmp_mag <- function(x, y) {
    x_len = length(x)
    y_len = length(y)
    if (x_len < y_len) return(-1L)
    if (x_len > y_len) return(1L)
    for (i in 1:x_len)
        if (c <- cmp_elem(x[[i]], y[[i]]))
	    return(c)
    return(0L)
}
cmp <- function(x, y) {
    x_sign = x[[1]]
    y_sign = y[[1]]
    if (x_sign == y_sign) {
        if (x[[1]] == 1L) return(cmp_mag(mag(x), mag(y)))
        if (x[[1]] == -1L) return(cmp_mag(mag(y), mag(x)))
        return(0L)
    }
    return((x_sign > y_sign) - (x_sign < y_sign))
}

eq <- function(x, y) cmp(x, y) == 0L
ne <- function(x, y) cmp(x, y) != 0L
le <- function(x, y) cmp(x, y) <= 0L
lt <- function(x, y) cmp(x, y) < 0L
ge <- function(x, y) cmp(x, y) >= 0L
gt <- function(x, y) cmp(x, y) > 0L

# Arithmetic operations

add_mag <- function(x, y) {
    # if x is shorter, swap the two vectors
    if (length(x) < length(y)) {
        tmp <- x; x <- y; y <- tmp
    }
    x_len = length(x)
    y_len = length(y)

    x_index <- x_len
    y_index <- y_len
    result <- integer(x_len)
    sum <- 0L

    # add common parts of both numbers
    while (y_index > 0L) {
        sum <- (x[[x_index]] %% elem_max + y[[y_index]] %% elem_max +
                sum %/% elem_max)  # TODO shift
        result[[x_index]] <- sum %% elem_max
	x_index <- x_index - 1L
        y_index <- y_index - 1L
    }

    # copy remainder of the longer number while carry propagation is required
    carry <- (sum >= elem_max)
    while (x_index > 0L & carry) {
        carry <- (result[[x_index]] <- x[[x_index]] + 1L) == 0L
        x_index <- x_index - 1L
    }

    # copy remainder of the longer number
    while (x_index > 0L) {
        result[[x_index]] <- x[[x_index]]
        x_index <- x_index - 1L
    }

    # grow result if necessary
    if (carry)
        return(c(0x01L, result))
    return(result)
}

add <- function(x, y) {
    if (y[[1]] == 0L)
        return(x)
    if (x[[1]] == 0L)
        return(y)
    if (x[[1]] == y[[1]])
        return(c(x[[1]], add_mag(mag(x), mag(y))))

    c <- cmp_mag(mag(x), mag(y))
    if (c == 0L)
        return(zero)
    if (c > 0L)
        result_mag <- sub_mag(mag(x), mag(y))
    else
        result_mag <- sub_mag(mag(y), mag(x))
    return(c(sign_prod(c, x[[1]]),
             strip_leading_zero_elems(result_mag)))
}

sub_mag <- function(big, little) {
    big_len = length(big)
    little_len = length(little)
    result <- integer(big_len)
    big_index <- big_len
    little_index <- little_len
    difference <- 0L

    # subtract common parts of both numbers
    while (little_index > 0L) {
        difference <- (big[[big_index]] - little[[little_index]] +
                       if (difference < 0L) -1L else 0L)
        result[[big_index]] <- difference %% elem_max
        big_index <- big_index - 1L
        little_index <- little_index - 1L
    }

    # subtract remainder of longer number while borrow propagates
    borrow <- if (difference < 0L) -1 else 0L
    while (big_index > 0L && borrow) {
        borrow <- (result[[big_index]] <- big[[big_index]] - 1L) == -1L
        big_index <- big_index - 1L
    }

    # copy remainder of the longer number
    while (big_index > 0L) {
        result[[big_index]] <- big[[big_index]]
        big_index <- big_index - 1L
    }

    return(result)
}

sub <- function(x, y) {
    if (y[[1]] == 0L)
        return(x)
    if (x[[1]] == 0L)
        return(negate(y))
    if (x[[1]] != y[[1]])
        return(c(x[[1]], add_mag(mag(x), mag(y))))

    c <- cmp_mag(mag(x), mag(y))
    if (c == 0L)
        return(zero)
    if (c > 0L)
        result_mag <- sub_mag(mag(x), mag(y))
    else
        result_mag <- sub_mag(mag(y), mag(x))
    return(c(sign_prod(c, x[[1]]),
             strip_leading_zero_elems(result_mag)))
}

negate <- function(x) c(-x[[1]], x[2:length(x)])

multiply_mag <- function(x, y) {
    x_len = length(x)
    y_len = length(y)
    x_start = x_len  # remove
    y_start = y_len  # remove
    c <- integer(x_len + y_len)

    carry <- 0L
    k <- y_start + x_start
    for (j in seq(y_start, 1L, -1L)) {
        product = y[[j]] * x[[x_start]] + carry
        c[k] <- product %% elem_max
        carry <- product %/% elem_max
        k <- k - 1L
    }
    c[x_start] <- carry

    if (x_len <= 1)
        return(c)

    for (i in seq(x_start - 1L, 1L, -1L)) {
        carry <- 0
        k <- y_start + i
        for (j in seq(y_start, 1L, -1L)) {
            product = y[[j]] * x[[i]] + c[[k]] + carry
            c[[k]] <- product %% elem_max
            carry <- product %/% elem_max
            k <- k - 1L
        }
        c[[i]] <- carry %% elem_max
    }
    return(c)
}

mul <- function(x, y) {
    if (y[[1]] == 0L || x[[1]] == 0)
        return(zero)

    return(c(sign_prod(x[[1]], y[[1]]),
            strip_leading_zero_elems(multiply_mag(mag(x), mag(y)))))
}

div_mag <- function(x_mag, y_mag) {
    if (length(y_mag) == 1L && y_mag == one_mag)
        return(x_mag)

    x <- c(1L, x_mag)
    y <- c(1L, y_mag)
    x_mag_log10 = log10_mag(x_mag); y_mag_log10 = log10_mag(y_mag)
    lo_log10 = x_mag_log10 - y_mag_log10 - (x_mag_log10 != y_mag_log10)
    hi_log10 <- x_mag_log10 - y_mag_log10 + 1L
    lo <- bigint_pow10(lo_log10)

    # try pruning hi > 10 or lo <= 10
    if (lo_log10 == 0L && hi_log10 > 1L) {
        lo10 = mul(lo, ten)
        if (gt(mul(lo10, y), x))
            hi <- lo10
        else {
            lo <- lo10
            hi <- bigint_pow10(hi_log10)
        }
    } else
      hi <- bigint_pow10(hi_log10)

    while (lt(lo, hi)) {
        mid <- div2(add(add(lo, hi), one))
        if (le(mul(mid, y), x))
            lo <- mid
        else
            hi <- sub(mid, one)
    }
    return(mag(lo))
}

div <- function(x, y) {
    # check if division by zero
    if (x[[1]] == 0L)
        return(zero)
    c <- cmp_mag(mag(x), mag(y))
    if (c == 0L)
        return(one)
    if (c < 0L)
        return(zero)
    return(c(sign_prod(x[[1]], y[[1]]), div_mag(mag(x), mag(y))))
}

# Helper arithmetic functions

div2_mag <- function(x) {
    x_len <- length(x)
    if (x_len == 1L)
        return(x %/% 2L)

    borrow <- (x[[1]] == 1)
    x_start <- borrow + 1L
    x_end <- x_len
    result_index <- 1L
    result <- integer(x_end - x_start + 1L)
    for (x_index in x_start:x_end) {
        d = x[[x_index]] + elem_max * borrow
        result[[result_index]] <- d %/% 2
        borrow <- d %% 2
        result_index <- result_index + 1L
    }
    return(result)
}

div2 <- function(x) {
    if (x[[1]] == 0L || (length(x) == 2L && x[[2]] == 1L))
        return(zero)
    return(c(x[[1]], div2_mag(mag(x))))
}

log10_mag <- function(m) elem_digits * (length(m) - 1L) + as.integer(log(m[[1]], 10))

bigint_pow10 <- function(n) c(1L, as.integer(10^(n %% elem_digits)), rep.int(0L, n %/% elem_digits))

# Misc functions
sign_prod <- function(x, y) (x == y) - (x != y)
strip_leading_zero_elems <- function(x) {
    for (i in 1:length(x))
        if (x[[i]] != 0L)
            return(x[i:length(x)])
    return(zero_mag)
}

zeropad <- function(s, n)
    paste(sep="", paste(collapse="", rep('0', max(0L, n - nchar(s)))), s)


# PIDIGITS program

pidigits <- function(args) {
    N = if (length(args)) as.integer(args[[1]]) else 100L
    ONE = one
    TWO = bigint("2")
    TEN = bigint("10")
    THREE = add(ONE, TWO)
    i <- k <- ns <- 0L
    k1 <- 1L
    a <- t <- u <- bigint("0")
    n <- d <- bigint("1")
    while (TRUE) {
        k <- k + 1L
        t <- mul(n, TWO)
        n <- mul(n, bigint(k))
        a <- add(a, t)
        k1 <- k1 + 2L
        k1_big <- bigint(i=k1)
        a <- mul(a, k1_big)
        d <- mul(d, k1_big)
        if (ge(a, n)) {
            n3a <- add(mul(n, THREE), a)
            t <- div(n3a, d)
	    td = mul(t, d)
            u <- add(sub(n3a, td), n)
            if (gt(d, u)) {
                ns <- ns * 10L + to_int(t)
                i <- i + 1L
                if (i %% 5L == 0L) {
		    cat(zeropad(as.character(ns), 5))
		    if (i %% 2L == 0L)
                       cat(sep="", "\t:", i, "\n")
                    ns = 0L
                }
                if (i >= N)
                    break
                a <- sub(a, td)  # TODO use td
                a <- mul(a, TEN)
                n <- mul(n, TEN)
            }
        }
    }
}

if (!exists("i_am_wrapper"))
    pidigits(commandArgs(trailingOnly=TRUE))