Euler: Common Routines

Common Routines

Notes

In a number of solutions, I pass along a vector of prime numbers. I generally pass along the first 10,000 primes which are contained in this file which is a slightly modified version from The Primes Page. To read it in as a vector, issue

primes <- as.vector(t(as.matrix(read.csv('primes_10000.txt', header=F, sep='', skip=5))))

Routines


IPI = "Infinite Precision Integer" (format A)
IPIB = "Infinite Precision Integer" (format B)

integer to integer
numDivisors(n):  number of divisors of natural number n
powerIt(x,p,n):  calculates integer x raised to a natural power p and retaining the last n digits.
integer to IPI
expand0(a):  converts positive integer a to IPI

IPI to integer
cmpIt(x,y): compares x to y
ucmpIt(x,y): compares |x| to |y|
signIt(x): returns -1 for negative numbers, 1 otherwise

IPI to IPI
lastDigits(x, n):  return last n digits of x
addIt(x,y): x + y
multiplyIt(x,y):  x*y
powerIt(x, n, digits):  calculate x^n and keep last digits

integer to list
expand(a, b, string, base): calculate a/b in base base
findFactors(n, primes):  factor n given a list of primes

integer/IPIB to IPIB
addIt2(x, y, str, block=8): calculate x+y using blocks of size block digits

IPIB to integer
digitsIt2(x, block): returns number of digits

IPIB to string
stringIt(x, block): converts IPIB format number into a string

IPI format

  1. Integer is converted to a vector of single digit integers ordered from most significant to least significant.  For example, 256 --> c(2, 5, 6).
  2. Leading zeros are acceptable.
  3. Negative numbers are indicated by having the first nonzero digit be negative.

IPIB format

  1. Integer is converted to a vector of integers (of up to a specified number of digits) ordered from most significant to least significant. Defaults to integers of 8 digits.  For example, 5,123,456,789 --> c(51, 23456789).
  2. Leading zeros are acceptable.
  3. Negative numbers are indicated by having the first nonzero digit be negative.


cmpIt <- function(x, y) {
    # x: integer vector
    # y: integer vector
    # Returns -1 if xy

    if (identical(x,y)) {
        return (0)
    }
   
    sx <- sign(x[1])
    sy <- sign(y[1])
    if (sx == sy)  {
        lenx <- length(x)
        leny <- length(y)
        if (lenx == leny) {
            d <- x - y
            return (sx*sign(d[which(d != 0)[1]]))
        }
        if (lenx < leny) {
            return (-sx)
        }
        else {
            return (sx)
        }
    }
    return ((sx - sy)/2)
}

ucmpIt <- function(x, y) {
    # x: integer vector
    # y: integer vector
    # Returns -1 if |x|<|y|, 0 if equal, 1 if |x|>|y|

    if (identical(x,y)) {
        return (0)
    }

    lenx <- length(x)
    leny <- length(y)
    if (lenx == leny) {
        d <- abs(x) - abs(y)
        s <- d[which(d != 0)[1]]
        if (is.na(s)) {
            return (0)
        }
        else {
            return (s)
        }
    }
    if (lenx < leny) {
        return (-1)
    }
    else {
        return (1)
    }
}
   
signIt <- function(x) {
    # if x is zero, then it is positive
    # else returns the first non-zero sign

    if (is.na(match(-1, sign(x)))) { 1 } else { -1}
}

addIt <- function(x, y) {

    sx <- signIt(x)
    sy <- signIt(y)
   
    flip <- 1
    if (sx != sy) {
        if (ucmpIt(x, y) < 0) {
            tmp <- x
            x <- y
            y <- tmp
        }
        flip <- sign(x[1])
        x <- abs(x)
        y <- -abs(y) 
    }
    else if (sx < 0) {
        x <- abs(x)
        y <- abs(y)
       
        flip <- -1
    }

    lenx <- length(x)
    leny <- length(y)
  
    if (lenx != leny) {
        if (lenx > leny) {
            y <- c(rep(0, lenx - leny), y)
        }
        else {
            x <- c(rep(0, leny - lenx), x)
        }
    }
  
    lenx <- length(x)
    carry <- 0
    i <- lenx
    #cat(x, y, "\n")
    while (i > 0) {      
        p <- x[i] + y[i] + carry
        x[i] <- p%%10
        carry <- trunc(p/10)
        #cat(p, x[i], carry, "\n")
        if (p < 0) {
            carry <- -1
        }
        i <- i - 1
    }
    if (carry > 0) {
        x <- c(carry, x)
    }
    flip*x
}

 multiplyIt <- function(x, y) {
    # x: vector (eg, c(1,2,8) represents the number 128)
    # y: vector
    # Runtime is  O(length(x)*length(y))

    lenx <- length(x)
    leny <- length(y)
    z <- rep(0, lenx + leny - 1)
  
    sx <- signIt(x)
    sy <- signIt(y)

    for (i in 1:lenx) {
       for (j in 1:leny) {
          z[i+j-1] <- z[i+j-1] + abs(x[i]*y[j])
       }
    } 

    i <- length(z)
    carry <- 0;
    while (i > 0) {
        z[i] <- z[i] + carry
        if (z[i] > 9) {
            carry <- floor(z[i]/10)
            z[i] <- z[i] %% 10
        }
        else {
            carry <- 0
        }
        i <- i - 1
    }
    if (carry > 0) {
        # expand z. carry is always less than 10.
        z <- c(carry, z)
    }
    (sx*sy)*z         
}

lastDigits <- function(x, n)
{
    if (n == 0) {
        return (x)
    }

    lenx <- length(x)
    x[max(1, lenx - n +1):lenx]
}

powerIt <- function(x, n, y=0) {
    # x: number
    # n: positive integer
    # y: number of last digits to retain. If y=0, then it retains all
   
    p <- c(1)
   
    while (n > 0) {
        if (n %% 2 == 1) {
            p <- lastDigits(multiplyIt(p, x), y)
        }
        x <- lastDigits(multiplyIt(x, x), y)
       
        n <- floor(n/2)
    }
    lastDigits(p, y)
}

expand <- function(a, b, collapse=F, base=10) {
    # a: positive numerator
    # b: positive denominator
    # collapse: if T, will return a string. if F, returns a vector of characters
    # base: numeric base. Allows up to base 32 (does not use I or O)
    #
    # returns M.N;P where M.N represents decimal and ;P represents the terminating repeated
    #               sequence. For example, 1/3 (base 10) = 0.;3

    baseSymb <- c(0,1,2,3,4,5,6,7,8,9,
                  'A','B','C','D','E','F','G','H','J','K',
                  'L','M','N','P','Q','R','S','T','U','V',
                  'W','X');

    remset <- c()   # remainder
    v <- c()
   
    # whole number conversion
 
    q <- trunc(a/b)
    repeat {
       r <- q %% base
       q <- trunc(q/base)
      
       v <- c(baseSymb[r +1], v)
       if (q == 0) {
           break
       }
    }

    # decimal conversion

    r <- a %% b
    if (r != 0) {
        v <- c(v, ".")
        offset <- length(v) 
       
        repeat {
           if (r == 0) {
               break
           }
           idx <- which(remset == r)
           if (length(idx) > 0) {
               v <- c(v[1:(idx+offset-1)], ";", v[-(1:(idx+offset-1))])
               break
           }
           remset <- c(remset, r)
         
           q <- trunc(r*base/b)   # quotient
           r <- (r*base) %% b       # remainder
         
           v <- c(v, baseSymb[q +1])   
        }
    }
    if (collapse) {
        paste(v, collapse="")
    }
    else {
        v
    }
}

findFactors <- function(n, primes) {
   
    pri <- c()  # prime
    deg <- c()  # degree

    m <- findInterval(n, primes)
    if (n == 1 || n == primes[m]) {
        pri <- c(n)
        deg <- c(1)
    }
    else {
        m <- findInterval(trunc(sqrt(n)), primes)
       
        x <- n
        i <- 1
        for (i in 1:m) {
            p <- primes[i]
            d <- 0
            while ( (x %% p == 0) && x > 0 ) {
                d <- d+1
                x <- x/p
            }
            if (d > 0) {
                pri <- c(pri, p)
                deg <- c(deg, d)

                if (x == 1) {
                    break
                }
            }
        }
        if (x != 1) {
            pri <- c(pri, x)
            if(log(x) >= 2*log(primes[length(primes)])) {
                 deg <- c(deg, 0)
            }
            else {
                 deg <- c(deg, 1)
            }
        }
    }
    list(pri=pri, deg=deg)
}

numDivisors <- function(n) {
   cnt <- 2
   m <- trunc(sqrt(n))
   if (m == 1) {
       return (cnt)
   }

   if (n %% m == 0) {
        cnt <- cnt + 1
   }
   d <- 2:(m-1)
   cnt <- cnt + 2*length(which (n %% d == 0))

   return (cnt)
}
expand0 <- function(n) {
   as.numeric(strsplit(as.character(n),"")[[1]])
}
lastDigits <- function(x, n)
{
    if (n == 0) {
        return (x)
    }

    lenx <- length(x)
    x[max(1, lenx - n +1):lenx]
}

powerIt <- function(x, n, y=0) {
    # x: number
    # n: number
    # y: number of last digits to retain. If y=0, then it retains all
   
    p <- c(1)
    x <- expand0(x)
   
    while (n > 0) {
        if (n %% 2 == 1) {
            p <- lastDigits(multiplyIt(p, x), y)
        }
        x <- lastDigits(multiplyIt(x, x), y)
       
        n <- floor(n/2)
    }
    lastDigits(p, y)
}

 stringIt <- function(x, block) {
    # converts *x* to a string
    # block: blocksize. If negative, will NOT reverse the string
 
    r <- sign(block)
    block <- abs(block)
    if (block == 1) {
         if (r < 0) {
             z <- paste(x, collapse="")
         }
         else {
             z <- paste(rev(x), collapse="")
         }
    }
    else {
        lenx <- length(x)
        x <- c(sapply(x[-lenx],
                      function(z) formatC(z, format="d", width=block, flag="0")),
               x[lenx])
        if (r < 0) {
            z <- paste(x, collapse="")
        }
        else {
            z <- paste(rev(x), collapse="")
        }
    }
    z
}

 digitsIt2 <- function(x, block) {
    # returns the number of digits represented by *x*

    if (!identical(x, 0)) {
        return ((length(x)-1)*block + 1 + floor(log10(x[length(x)])))
    }
    1
}
   
addIt2 <- function(x, y, str=T, block=8) {
    # if str=F, returns a vector of integers where each vector entry
    # has at most *block* digits except possibly for the last entry which
    # has at least 1 digit.   

    if (is.character(x)) {
        xp <- splitIt(x, block)
    }
    else {
        xp <- x
    }
    if (is.character(y)) {
        yp <- splitIt(y, block)
    }
    else {
        yp <- y
    }
   
    carry <- 0
    plen <- max(length(xp), length(yp))
    xp <- c(xp, rep(0, plen - length(xp)))
    yp <- c(yp, rep(0, plen - length(yp)))
    for (i in 1:plen) {
        p <- xp[i] + yp[i] + carry
        xp[i] <- p%%10^block
        carry <- trunc(p/10^block)
    }
    if (carry > 0) {
        xp[plen+1] <- carry
    }
 
    if (str) {
        xp <- stringIt(xp, block)
    }
    xp
}

No comments:

Post a Comment

Subscribe to: Comments (Atom)