Arbitrary Base Conversion

Python

# divides longnum src (in base src_base) by divisor
# returns a pair of (longnum dividend, remainder)
def divmod_long(src, src_base, divisor):
  dividend=[]
  remainder=0
  for d in src:
    (e, remainder) = divmod(d + remainder * src_base, divisor)
    if dividend or e: dividend += [e]
  return (dividend, remainder)

def convert(src, src_base, dst_base):
  result = []
  while src:
    (src, remainder) = divmod_long(src, src_base, dst_base)
    result = [remainder] + result
  return result

Here's a Haskell solution

import Data.List
import Control.Monad

type Numeral = (Int, [Int])

swap              ::  (a,b) -> (b,a)
swap (x,y)        =   (y,x)

unfoldl           ::  (b -> Maybe (b,a)) -> b -> [a]
unfoldl f         =   reverse . unfoldr (fmap swap . f)

normalize         ::  Numeral -> Numeral
normalize (r,ds)  =   (r, dropWhile (==0) ds)

divModLongInt            ::  Numeral -> Int -> (Numeral,Int)
divModLongInt (r,dd) dv  =   let  divDigit c d = swap ((c*r+d) `divMod` dv)
                                  (remainder, quotient) = mapAccumR divDigit 0 (reverse dd)
                             in   (normalize (r,reverse quotient), remainder)

changeRadixLongInt       ::  Numeral -> Int -> Numeral
changeRadixLongInt n r'  =   (r', unfoldl produceDigit n)
  where  produceDigit  (_,[])   =  Nothing
         produceDigit  x        =  Just (divModLongInt x r')

changeRadix :: [Int] -> Int -> Int -> [Int]
changeRadix digits origBase newBase = snd $ changeRadixLongInt (origBase,digits) newBase

doLine line = let [(digits,rest0)] = reads line
                  [(origBase,rest1)] = reads rest0
                  [(newBase,rest2)] = reads rest1
              in show $ changeRadix digits origBase newBase

main = interact (unlines . map doLine . lines)

And running the tests from the question:

$ ./a.out 
[1,2,3,4] 16 16
[1,2,3,4]
[1,0] 10 100
[10]
[41, 15, 156, 123, 254, 156, 141, 2, 24] 256 16
[2,9,0,15,9,12,7,11,15,14,9,12,8,13,0,2,1,8]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] 2 10
[1,2,3,7,9,4,0,0,3,9,2,8,5,3,8,0,2,7,4,8,9,9,1,2,4,2,2,3]
[41, 42, 43] 256 36
[1,21,29,22,3]

R

Handles many thousands of elements^* in less than a minute.

addb <- function(v1,v2,b) {
    ml <- max(length(v1),length(v2))
    v1 <- c(rep(0, ml-length(v1)),v1)
    v2 <- c(rep(0, ml-length(v2)),v2)
    v1 = v1 + v2
    resm = v1%%b
    resd = c(floor(v1/b),0)
    while (any(resd != 0)) {
        v1 = c(0,resm) + resd
        resm = v1%%b
        resd = c(floor(v1/b),0)
    }
    while (v1[1] == 0) v1 = v1[-1]
    return(v1)
}

redb <- function(v,b) {
    return (addb(v,0,b))
}

mm = rbind(1)

mulmat <- function(fromb, tob, n) {
    if (dim(mm)[2] >= n) return(mm)
    if (n == 1) return(1)
    newr = addb(mulmat(fromb,tob,n-1) %*% rep(fromb-1,n-1), 1, tob)
    newm = mulmat(fromb,tob,n-1)
    while (is.null(dim(newm)) || dim(newm)[1] < length(newr)) newm = rbind(0,newm)
    mm <<-  cbind(newr, newm)
    return(mm)
}

dothelocomotion <- function(fromBase, toBase, v) {
    mm  <<- rbind(1)
    return(redb(mulmat(fromBase, toBase, length(v)) %*% v, toBase))
}

^{* for >500 elements you have to raise the default recursion level or do not reset the mm matrix on dothelocomotion()}

Examples:

v1 = c(41, 15, 156, 123, 254, 156, 141, 2, 24)
dothelocomotion(256,16,v1)
2  9  0 15  9 12  7 11 15 14  9 12  8 13  0  2  1  8

dothelocomotion(256,36,c(41,42,43))
1 21 29 22  3

dothelocomotion(2,10, rep(1,90))
1 2 3 7 9 4 0 0 3 9 2 8 5 3 8 0 2 7 4 8 9 9 1 2 4 2 2 3