diff --git a/.travis.yml b/.travis.yml
index a77e76ffe2354..e3c4e6620dd0d 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -5,12 +5,12 @@ compiler:
 notifications:
     email: false
 before_install:
-    - BUILDOPTS="LLVM_CONFIG=llvm-config-3.2 USE_QUIET=0 USE_LIB64=0"; for lib in LLVM ZLIB SUITESPARSE ARPACK BLAS FFTW LAPACK GMP PCRE LIBUNWIND READLINE GRISU OPENLIBM RMATH LIBUV; do export BUILDOPTS="$BUILDOPTS USE_SYSTEM_$lib=1"; done
+    - BUILDOPTS="LLVM_CONFIG=llvm-config-3.2 USE_QUIET=0 USE_LIB64=0"; for lib in LLVM ZLIB SUITESPARSE ARPACK BLAS FFTW LAPACK GMP MPFR PCRE LIBUNWIND READLINE GRISU OPENLIBM RMATH LIBUV; do export BUILDOPTS="$BUILDOPTS USE_SYSTEM_$lib=1"; done
     - sudo apt-get update -qq -y
     - sudo apt-get install zlib1g-dev
     - sudo add-apt-repository ppa:staticfloat/julia-deps -y
     - sudo apt-get update -qq -y
-    - sudo apt-get install patchelf gfortran llvm-3.2-dev libsuitesparse-dev libncurses5-dev libopenblas-dev liblapack-dev libarpack2-dev libfftw3-dev libgmp-dev libpcre3-dev libunwind7-dev libreadline-dev libdouble-conversion-dev libopenlibm-dev librmath-dev libuv-julia-dev -y
+    - sudo apt-get install patchelf gfortran llvm-3.2-dev libsuitesparse-dev libncurses5-dev libopenblas-dev liblapack-dev libarpack2-dev libfftw3-dev libgmp-dev libpcre3-dev libunwind7-dev libreadline-dev libdouble-conversion-dev libopenlibm-dev librmath-dev libuv-julia-dev libmpfr-dev -y
 script:
     - make $BUILDOPTS PREFIX=/tmp/julia install
     - cd .. && mv julia julia2
diff --git a/LICENSE.md b/LICENSE.md
index 7f58931c0a54b..fb06837f22f72 100644
--- a/LICENSE.md
+++ b/LICENSE.md
@@ -58,6 +58,7 @@ External libraries, if used, include their own licenses:
 - [LAPACK](http://netlib.org/lapack/LICENSE.txt)
 - [LIBUNWIND](http://git.savannah.gnu.org/gitweb/?p=libunwind.git;a=blob_plain;f=LICENSE;hb=master)
 - [LLVM](http://llvm.org/releases/3.0/LICENSE.TXT)
+- [MPFR](http://www.mpfr.org/mpfr-current/mpfr.html#Copying)
 - [OPENBLAS](https://raw.github.com/xianyi/OpenBLAS/master/LICENSE)
 - [PCRE](http://www.pcre.org/licence.txt)
 - [READLINE](http://cnswww.cns.cwru.edu/php/chet/readline/rltop.html)
diff --git a/Make.inc b/Make.inc
index 8b781055b59f3..99c380f50ab19 100644
--- a/Make.inc
+++ b/Make.inc
@@ -163,6 +163,7 @@ USE_SYSTEM_BLAS=0
 USE_SYSTEM_LAPACK=0
 USE_SYSTEM_FFTW=0
 USE_SYSTEM_GMP=0
+USE_SYSTEM_MPFR=0
 USE_SYSTEM_ARPACK=0
 USE_SYSTEM_SUITESPARSE=0
 USE_SYSTEM_ZLIB=0
diff --git a/Makefile b/Makefile
index 37c1777a3bda7..b5b4bd7c39c3c 100644
--- a/Makefile
+++ b/Makefile
@@ -66,7 +66,7 @@ JL_PRIVATE_LIBS = amd arpack camd ccolamd cholmod colamd \
                   fftw3 fftw3f fftw3_threads fftw3f_threads \
                   gmp grisu openlibm openlibm-extras pcre \
                   random Rmath spqr suitesparse_wrapper \
-                  umfpack z openblas
+                  umfpack z openblas mpfr
 
 PREFIX ?= julia-$(JULIA_COMMIT)
 install:
diff --git a/README.md b/README.md
index 0d756c88326ef..19a6dba48f3a0 100644
--- a/README.md
+++ b/README.md
@@ -152,6 +152,7 @@ Julia uses the following external libraries, which are automatically downloaded
 - **[FFTW]**                — library for computing fast Fourier transforms very quickly and efficiently.
 - **[PCRE]**                — Perl-compatible regular expressions library.
 - **[GMP]**                 — the GNU multiple precision arithmetic library, needed for bigint support.
+- **[MPFR]**                — the GNU multiple precision floating point library, needed for arbitrary precision floating point support.
 - **[double-conversion]**   — efficient number-to-text conversion.
 - **[Rmath]**               — basic RNGs and distributions.
 
@@ -179,6 +180,7 @@ Julia uses the following external libraries, which are automatically downloaded
 [FemtoLisp]:    https://github.com/JeffBezanson/femtolisp
 [readline]:     http://cnswww.cns.cwru.edu/php/chet/readline/rltop.html
 [GMP]:          http://gmplib.org/
+[MPFR]:         http://www.mpfr.org/
 [double-conversion]: http://double-conversion.googlecode.com/
 [Rmath]:        http://cran.r-project.org/doc/manuals/R-admin.html#The-standalone-Rmath-library
 [libuv]:        https://github.com/JuliaLang/libuv
diff --git a/base/bigfloat.jl b/base/bigfloat.jl
deleted file mode 100644
index d32c693626dd7..0000000000000
--- a/base/bigfloat.jl
+++ /dev/null
@@ -1,133 +0,0 @@
-BigFloat_clear(mpf::Vector{Int32}) = ccall((:__gmpf_clear, :libgmp), Void, (Ptr{Void},), mpf)
-
-immutable BigFloat <: FloatingPoint
-    mpf::Vector{Int32}
-    function BigFloat() 
-        z = Array(Int32, 6)
-        ccall((:__gmpf_init,:libgmp), Void, (Ptr{Void},), z)
-        b = new(z)
-        finalizer(b.mpf, BigFloat_clear)
-        return b
-    end
-end
-
-function BigFloat(x::BigFloat)
-    z = BigFloat()
-    ccall((:__gmpf_set, :libgmp), Void, (Ptr{Void}, Ptr{Void}), z.mpf, x.mpf)
-    return z
-end
-
-function BigFloat(x::String)
-    z = BigFloat()
-    err = ccall((:__gmpf_set_str, :libgmp), Int32, (Ptr{Void}, Ptr{Uint8}, Int32), z.mpf, bytestring(x), 0)
-    if err != 0; error("Invalid input"); end
-    return z
-end
-
-function BigFloat(x::Float64)
-    z = BigFloat()
-    ccall((:__gmpf_set_d, :libgmp), Void, (Ptr{Void}, Float64), z.mpf, x)
-    return z
-end
-
-function BigFloat(x::Uint)
-    z = BigFloat()
-    ccall((:__gmpf_set_ui, :libgmp), Void, (Ptr{Void}, Uint), z.mpf, x)
-    return z
-end
-
-function BigFloat(x::Int)
-    z = BigFloat()
-    ccall((:__gmpf_set_si, :libgmp), Void, (Ptr{Void}, Int), z.mpf, x)
-    return z
-end
-
-function BigFloat(x::BigInt)
-    z = BigFloat()
-    ccall((:__gmpf_set_z, :libgmp), Void, (Ptr{Void}, Ptr{Void}), z.mpf, x.mpz)
-    return z
-end
-
-BigFloat(x::Bool) = BigFloat(uint(x))
-BigFloat(x::Signed) = BigFloat(int(x))
-BigFloat(x::Unsigned) = BigFloat(uint(x))
-#BigFloat(x::Int128) = BigFloat(BigInt(x))
-#BigFloat(x::Uint128) = BigFloat(BigInt(x))
-if WORD_SIZE == 32
-    BigFloat(x::Int64) = BigFloat(string(x))
-    BigFloat(x::Uint64) = BigFloat(BigInt(x))
-end
-BigFloat(x::Float32) = BigFloat(float64(x))
-BigFloat(x::Rational) = BigFloat(num(x)) / BigFloat(den(x))
-
-convert(::Type{BigFloat}, x::Rational) = BigFloat(x) # to resolve ambiguity
-convert(::Type{BigFloat}, x::Real) = BigFloat(x)
-
-convert(::Type{Float64}, x::BigFloat) = ccall((:__gmpf_get_d,:libgmp), Float64, (Ptr{Void},), x.mpf)
-convert(::Type{FloatingPoint}, x::BigInt) = BigFloat(x)
-
-promote_rule{T<:Union(Integer,FloatingPoint)}(::Type{BigFloat}, ::Type{T}) = BigFloat
-promote_rule{T<:FloatingPoint}(::Type{BigInt},::Type{T}) = BigFloat
-
-
-# mpf doesn't have inf or nan
-isnan(x::BigFloat) = false
-isinf(x::BigFloat) = false
-
-# Binary ops
-for (fJ, fC) in ((:+,:add), (:-,:sub), (:*,:mul), (:/,:div))
-    @eval begin 
-        function ($fJ)(x::BigFloat, y::BigFloat)
-            z = BigFloat()
-            ccall(($(string(:__gmpf_,fC)),:libgmp), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}), z.mpf, x.mpf, y.mpf)
-            return z
-        end
-    end
-end
-
-function -(x::BigFloat)
-    z = BigFloat()
-    ccall((:__gmpf_neg, :libgmp), Void, (Ptr{Void}, Ptr{Void}), z.mpf, x.mpf)
-    return z
-end
-
-function cmp(x::BigFloat, y::BigFloat)
-    ccall((:__gmpf_cmp, :libgmp), Int32, (Ptr{Void}, Ptr{Void}), x.mpf, y.mpf)
-end
-
-for f in (:sqrt, :ceil, :floor, :trunc)
-    @eval begin
-        function ($f)(x::BigFloat)
-            z = BigFloat()
-            ccall(($(string(:__gmpf_,f)), :libgmp), Void, (Ptr{Void}, Ptr{Void}), z.mpf, x.mpf)
-            return z
-        end
-    end
-end
-
-function ^(x::BigFloat, y::Uint)
-    z = BigFloat()
-    ccall((:__gmpf_pow_ui, :libgmp), Void, (Ptr{Void}, Ptr{Void}, Uint), z.mpf, x.mpf, y)
-    return z
-end
-
-^(x::Float32, y::BigInt) = BigFloat(x)^y
-^(x::Float64, y::BigInt) = BigFloat(x)^y
-
-==(x::BigFloat, y::BigFloat) = cmp(x,y) == 0
-<=(x::BigFloat, y::BigFloat) = cmp(x,y) <= 0
->=(x::BigFloat, y::BigFloat) = cmp(x,y) >= 0
-<(x::BigFloat, y::BigFloat) = cmp(x,y) < 0
->(x::BigFloat, y::BigFloat) = cmp(x,y) > 0
-
-function string(x::BigFloat)
-    lng = 128
-    for i = 1:2
-        z = Array(Uint8, lng)
-        lng = ccall((:__gmp_snprintf,:libgmp), Int32, (Ptr{Uint8}, Uint, Ptr{Uint8}, Ptr{Void}...), z, lng, "%.Fe", x.mpf)
-        if lng < 128 || i == 2; return bytestring(convert(Ptr{Uint8}, z[1:lng])); end
-    end
-end
-
-show(io::IO, b::BigFloat) = print(io, string(b))
-showcompact(io::IO, b::BigFloat) = print(io, string(b))
diff --git a/base/bigint.jl b/base/bigint.jl
index 206cb2b809bbe..5f09497512ccd 100644
--- a/base/bigint.jl
+++ b/base/bigint.jl
@@ -1,39 +1,33 @@
-BigInt_clear(mpz::Vector{Int32}) = ccall((:__gmpz_clear, :libgmp), Void, (Ptr{Void},), mpz)
-
-immutable BigInt <: Integer
-    mpz::Vector{Int32}
+type BigInt <: Integer
+    alloc::Cint
+    size::Cint
+    d::Ptr{Void}
     function BigInt() 
-        z = Array(Int32, 4)
-        ccall((:__gmpz_init,:libgmp), Void, (Ptr{Void},), z)
-        b = new(z)
-        finalizer(b.mpz, BigInt_clear)
+        b = new(zero(Cint), zero(Cint), C_NULL)
+        ccall((:__gmpz_init,:libgmp), Void, (Ptr{BigInt},), &b)
+        finalizer(b, BigInt_clear)
         return b
     end
 end
+BigInt_clear(mpz::BigInt) = ccall((:__gmpz_clear, :libgmp), Void, (Ptr{BigInt},), &mpz)
 
-function BigInt(x::BigInt)
-    z = BigInt()
-    ccall((:__gmpz_set, :libgmp), Void, (Ptr{Void}, Ptr{Void}), z.mpz, x.mpz)
-    return z
-end
-
+BigInt(x::BigInt) = x
 function BigInt(x::String)
     z = BigInt()
-    err = ccall((:__gmpz_set_str, :libgmp), Int32, (Ptr{Void}, Ptr{Uint8}, Int32), z.mpz, bytestring(x), 0)
+    err = ccall((:__gmpz_set_str, :libgmp), Int32, (Ptr{BigInt}, Ptr{Uint8}, Int32), &z, bytestring(x), 0)
     if err != 0; error("Invalid input"); end
     return z
 end
 
 function BigInt(x::Int)
     z = BigInt()
-    ccall((:__gmpz_set_si, :libgmp), Void, (Ptr{Void}, Int), z.mpz, x)
+    ccall((:__gmpz_set_si, :libgmp), Void, (Ptr{BigInt}, Clong), &z, x)
     return z
 end
 
 function BigInt(x::Uint)
     z = BigInt()
-    ccall((:__gmpz_set_ui, :libgmp), Void,
-        (Ptr{Void}, Uint), z.mpz, x)
+    ccall((:__gmpz_set_ui, :libgmp), Void,(Ptr{BigInt}, Culong), &z, x)
     return z
 end
 
@@ -50,10 +44,10 @@ end
 convert{T<:Integer}(::Type{BigInt}, x::T) = BigInt(x)
 
 convert(::Type{Int}, n::BigInt) =
-    ccall((:__gmpz_get_si, :libgmp), Int, (Ptr{Void},), n.mpz)
+    convert(Int, ccall((:__gmpz_get_si, :libgmp), Clong, (Ptr{BigInt},), &n))
 
 convert(::Type{Uint}, n::BigInt) =
-    ccall((:__gmpz_get_ui, :libgmp), Uint, (Ptr{Void},), n.mpz)
+    convert(Uint, ccall((:__gmpz_get_ui, :libgmp), Culong, (Ptr{BigInt},), &n))
 
 promote_rule{T<:Integer}(::Type{BigInt}, ::Type{T}) = BigInt
 
@@ -65,18 +59,62 @@ for (fJ, fC) in ((:+, :add), (:-,:sub), (:*, :mul),
     @eval begin
         function ($fJ)(x::BigInt, y::BigInt)
             z = BigInt()
-            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}), z.mpz, x.mpz, y.mpz)
+            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &x, &y)
             return z
         end
     end
 end
 
+# Basic arithmetic without promotion
+function +(x::BigInt, c::Culong)
+    z = BigInt()
+    ccall((:__gmpz_add_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
+    return z
+end
++(c::Culong, x::BigInt) = x + c
++(c::Unsigned, x::BigInt) = x + convert(Culong, c)
++(x::BigInt, c::Unsigned) = x + convert(Culong, c)
++(x::BigInt, c::Signed) = c < 0 ? -(x, convert(Culong, -c)) : x + convert(Culong, c)
++(c::Signed, x::BigInt) = c < 0 ? -(x, convert(Culong, -c)) : x + convert(Culong, c)
+
+function -(x::BigInt, c::Culong)
+    z = BigInt()
+    ccall((:__gmpz_sub_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
+    return z
+end
+function -(c::Culong, x::BigInt)
+    z = BigInt()
+    ccall((:__gmpz_ui_sub, :libgmp), Void, (Ptr{BigInt}, Culong, Ptr{BigInt}), &z, c, &x)
+    return z
+end
+-(x::BigInt, c::Unsigned) = -(x, convert(Culong, c))
+-(c::Unsigned, x::BigInt) = -(convert(Culong, c), x)
+-(x::BigInt, c::Signed) = c < 0 ? +(x, convert(Culong, -c)) : -(x, convert(Culong, c))
+-(c::Signed, x::BigInt) = c < 0 ? -(x + convert(Culong, -c)) : -(convert(Culong, c), x)
+
+function *(x::BigInt, c::Culong)
+    z = BigInt()
+    ccall((:__gmpz_mul_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
+    return z
+end
+*(c::Culong, x::BigInt) = x * c
+*(c::Unsigned, x::BigInt) = x * convert(Culong, c)
+*(x::BigInt, c::Unsigned) = x * convert(Culong, c)
+function *(x::BigInt, c::Clong)
+    z = BigInt()
+    ccall((:__gmpz_mul_si, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Clong), &z, &x, c)
+    return z
+end
+*(c::Clong, x::BigInt) = x * c
+*(x::BigInt, c::Signed) = x * convert(Clong, c)
+*(c::Signed, x::BigInt) = x * convert(Clong, c)
+
 # unary ops
 for (fJ, fC) in ((:-, :neg), (:~, :com))
     @eval begin
         function ($fJ)(x::BigInt)
             z = BigInt()
-            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{Void}, Ptr{Void}), z.mpz, x.mpz)
+            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}), &z, &x)
             return z
         end
     end
@@ -84,7 +122,7 @@ end
 
 function <<(x::BigInt, c::Uint)
     z = BigInt()
-    ccall((:__gmpz_mul_2exp, :libgmp), Void, (Ptr{Void}, Ptr{Void}, Uint), z.mpz, x.mpz, c)
+    ccall((:__gmpz_mul_2exp, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
     return z
 end
 <<(x::BigInt, c::Int32)   = c<0 ? throw(DomainError()) : x<<uint(c)
@@ -92,7 +130,7 @@ end
 
 function >>(x::BigInt, c::Uint)
     z = BigInt()
-    ccall((:__gmpz_fdiv_q_2exp, :libgmp), Void, (Ptr{Void}, Ptr{Void}, Uint), z.mpz, x.mpz, c)
+    ccall((:__gmpz_fdiv_q_2exp, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
     return z
 end
 >>(x::BigInt, c::Int32)   = c<0 ? throw(DomainError()) : x>>uint(c)
@@ -101,23 +139,23 @@ end
 function divrem(x::BigInt, y::BigInt)
     z1 = BigInt()
     z2 = BigInt()
-    ccall((:__gmpz_tdiv_qr, :libgmp), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}), z1.mpz, z2.mpz, x.mpz, y.mpz)
+    ccall((:__gmpz_tdiv_qr, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &(z1.mpz), &(z2.mpz), &x, &y)
     z1, z2
 end
 
 function cmp(x::BigInt, y::BigInt)
-    ccall((:__gmpz_cmp, :libgmp), Int32, (Ptr{Void}, Ptr{Void}), x.mpz, y.mpz)
+    ccall((:__gmpz_cmp, :libgmp), Int32, (Ptr{BigInt}, Ptr{BigInt}), &x, &y)
 end
 
 function sqrt(x::BigInt)
     z = BigInt()
-    ccall((:__gmpz_sqrt, :libgmp), Void, (Ptr{Void}, Ptr{Void}), z.mpz, x.mpz)
+    ccall((:__gmpz_sqrt, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}), &z, &x)
     return z
 end
 
 function ^(x::BigInt, y::Uint)
     z = BigInt()
-    ccall((:__gmpz_pow_ui, :libgmp), Void, (Ptr{Void}, Ptr{Void}, Uint), z.mpz, x.mpz, y)
+    ccall((:__gmpz_pow_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, y)
     return z
 end
 
@@ -139,8 +177,8 @@ function gcdx(a::BigInt, b::BigInt)
     s = BigInt()
     t = BigInt()
     ccall((:__gmpz_gcdext, :libgmp), Void,
-        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}),
-        g, s, t, a.mpz, b.mpz)
+        (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
+        &g, &s, &t, &a, &b)
     BigInt(g), BigInt(s), BigInt(t)
 end
 
@@ -152,14 +190,14 @@ function factorial(bn::BigInt)
     end
     z = BigInt()
     ccall((:__gmpz_fac_ui, :libgmp), Void,
-        (Ptr{Void}, Uint), z.mpz, n)
+        (Ptr{BigInt}, Culong), &z, n)
     return z
 end
 
 function binomial(n::BigInt, k::Uint)
     z = BigInt()
     ccall((:__gmpz_bin_ui, :libgmp), Void,
-        (Ptr{Void}, Ptr{Void}, Uint), z.mpz, n.mpz, k)
+        (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &n, k)
     return z
 end
 binomial(n::BigInt, k::Integer) = k<0 ? throw(DomainError()) : binomial(n, uint(k))
@@ -173,7 +211,7 @@ binomial(n::BigInt, k::Integer) = k<0 ? throw(DomainError()) : binomial(n, uint(
 function string(x::BigInt)
     lng = ndigits(x) + 2
     z = Array(Uint8, lng)
-    lng = ccall((:__gmp_snprintf,:libgmp), Int32, (Ptr{Uint8}, Uint, Ptr{Uint8}, Ptr{Void}...), z, lng, "%Zd", x.mpz)
+    lng = ccall((:__gmp_snprintf,:libgmp), Int32, (Ptr{Uint8}, Culong, Ptr{Uint8}, Ptr{BigInt}...), z, lng, "%Zd", &x)
     return bytestring(convert(Ptr{Uint8}, z[1:lng]))
 end
 
@@ -181,4 +219,4 @@ function show(io::IO, x::BigInt)
     print(io, string(x))
 end
 
-ndigits(x::BigInt) = ccall((:__gmpz_sizeinbase,:libgmp), Uint, (Ptr{Void}, Int32), x.mpz, 10)
+ndigits(x::BigInt) = ccall((:__gmpz_sizeinbase,:libgmp), Culong, (Ptr{BigInt}, Int32), &x, 10)
diff --git a/base/exports.jl b/base/exports.jl
index 5f3dfc083e154..5194645f31f43 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -9,6 +9,7 @@ export
     LibRandom,
     Random,
     Math,
+    MPFR,
     GMP,
     Sort,
     Test,
@@ -820,6 +821,16 @@ export
     randn,
     srand,
 
+# precision
+    exp10,
+    get_precision,
+    get_bigfloat_precision,
+    set_bigfloat_precision,
+    with_bigfloat_precision,
+    get_bigfloat_rounding,
+    set_bigfloat_rounding,
+    with_bigfloat_rounding,
+
 # statistics
     cor,
     cov,
diff --git a/base/floatfuncs.jl b/base/floatfuncs.jl
index 6eb85d352fe7f..4bcc6a6aa9591 100644
--- a/base/floatfuncs.jl
+++ b/base/floatfuncs.jl
@@ -31,6 +31,10 @@ maxintfloat() = maxintfloat(Float64)
 
 integer_valued(x::FloatingPoint) = (trunc(x)==x)&isfinite(x)
 
+## precision, as defined by the effective number of bits in the mantissa ##
+get_precision(::Float32) = 24
+get_precision(::Float64) = 53
+
 num2hex(x::Float32) = hex(box(Uint32,unbox(Float32,x)),8)
 num2hex(x::Float64) = hex(box(Uint64,unbox(Float64,x)),16)
 
diff --git a/base/gmp.jl b/base/gmp.jl
index f0dae0bc6429b..be33590a5c93b 100644
--- a/base/gmp.jl
+++ b/base/gmp.jl
@@ -1,6 +1,6 @@
 module GMP
 
-export BigInt, BigFloat
+export BigInt
 
 import
     Base.(*),
@@ -43,7 +43,12 @@ import
     Base.string,
     Base.trunc
 
+type mpz_struct
+    alloc::Cint
+    size::Cint
+    d::Ptr{Void}
+end
+
 include("bigint.jl")
-include("bigfloat.jl")
 
 end # module
diff --git a/base/mpfr.jl b/base/mpfr.jl
new file mode 100644
index 0000000000000..fdc4e8ecd1e69
--- /dev/null
+++ b/base/mpfr.jl
@@ -0,0 +1,664 @@
+module MPFR
+
+export
+    # Types
+    BigFloat,
+    # Functions
+    exp10,
+    get_bigfloat_precision,
+    set_bigfloat_precision,
+    with_bigfloat_precision,
+    set_bigfloat_rounding,
+    get_bigfloat_rounding,
+    with_bigfloat_rounding
+
+import
+    Base: (*), +, -, /, <, <=, ==, >, >=, ^, besselj, besselj0, besselj1,
+        bessely, bessely0, bessely1, ceil, cmp, convert, copysign, exp, exp2,
+        exponent, factorial, floor, hypot, integer_valued, iround, isfinite,
+        isinf, isnan, ldexp, log, log2, log10, max, min, mod, modf, nextfloat,
+        prevfloat, promote_rule, rem, round, show, showcompact, sum, sqrt,
+        string, trunc, get_precision,
+    # import trigonometric functions
+        sin, cos, tan, sec, csc, cot, acos, asin, atan, cosh, sinh, tanh,
+        sech, csch, coth, acosh, asinh, atanh, atan2
+
+const ROUNDING_MODE = [0]
+const DEFAULT_PRECISION = [256]
+
+# Rounding modes
+const RoundToNearest = 0
+const RoundToZero = 1
+const RoundUp = 2
+const RoundDown = 3
+const RoundAwayZero = 4
+
+# Basic type and initialization definitions
+
+type BigFloat <: FloatingPoint
+    prec::Clong
+    sign::Cint
+    exp::Clong
+    d::Ptr{Void}
+    function BigFloat()
+        N = get_bigfloat_precision()
+        z = new(zero(Clong), zero(Cint), zero(Clong), C_NULL)
+        ccall((:mpfr_init2,:libmpfr), Void, (Ptr{BigFloat}, Clong), &z, N)
+        finalizer(z, MPFR_clear)
+        return z
+    end
+    # Not recommended for general use
+    function BigFloat(prec::Clong, sign::Cint, exp::Clong, d::Ptr{Void})
+        new(prec, sign, exp, d)
+    end
+end
+MPFR_clear(mpfr::BigFloat) = ccall((:mpfr_clear, :libmpfr), Void, (Ptr{BigFloat},), &mpfr)
+
+BigFloat(x::BigFloat) = x
+
+for (fJ, fC) in ((:si,:Int), (:ui,:Uint), (:d,:Float64))
+    @eval begin
+        function BigFloat(x::($fC))
+            z = BigFloat()
+            ccall(($(string(:mpfr_set_,fJ)), :libmpfr), Int32, (Ptr{BigFloat}, ($fC), Int32), &z, x, ROUNDING_MODE[end])
+            return z
+        end
+    end
+end
+
+function BigFloat(x::BigInt)
+    z = BigFloat()
+    ccall((:mpfr_set_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, ROUNDING_MODE[end])   
+    return z
+end
+
+function BigFloat(x::String, base::Int)
+    z = BigFloat()
+    err = ccall((:mpfr_set_str, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{Uint8}, Int32, Int32), &z, x, base, ROUNDING_MODE[end])
+    if err != 0; error("Invalid input"); end
+    return z
+end
+BigFloat(x::String) = BigFloat(x, 10)
+
+
+BigFloat(x::Bool) = BigFloat(uint(x))
+BigFloat(x::Signed) = BigFloat(int(x))
+BigFloat(x::Unsigned) = BigFloat(uint(x))
+if WORD_SIZE == 32
+    BigFloat(x::Int64) = BigFloat(string(x))
+    BigFloat(x::Uint64) = BigFloat(BigInt(x))
+end
+BigFloat(x::Float32) = BigFloat(float64(x))
+BigFloat(x::Rational) = BigFloat(num(x)) / BigFloat(den(x))
+
+convert(::Type{BigFloat}, x::Rational) = BigFloat(x) # to resolve ambiguity
+convert(::Type{BigFloat}, x::Real) = BigFloat(x)
+
+
+convert(::Type{Int64}, x::BigFloat) = int64(convert(Clong, x))
+convert(::Type{Int32}, x::BigFloat) = int32(convert(Clong, x))
+convert(::Type{Clong}, x::BigFloat) = integer_valued(x) ?
+    ccall((:mpfr_get_si,:libmpfr), Clong, (Ptr{BigFloat}, Int32), &x, ROUNDING_MODE[end]) :
+    throw(InexactError())
+convert(::Type{Uint32}, x::BigFloat) = uint64(convert(Culong, x))
+convert(::Type{Uint32}, x::BigFloat) = uint32(convert(Culong, x))
+convert(::Type{Culong}, x::BigFloat) = integer_valued(x) ?
+    ccall((:mpfr_get_ui,:libmpfr), Culong, (Ptr{BigFloat}, Int32), &x, ROUNDING_MODE[end]) :
+    throw(InexactError())
+function convert(::Type{BigInt}, x::BigFloat) 
+    if integer_valued(x)
+        z = BigInt()
+        ccall((:mpfr_get_z,:libmpfr), Int32, (Ptr{BigInt}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+        return z
+    else
+        throw(InexactError())
+    end
+end
+convert(::Type{Float64}, x::BigFloat) = ccall((:mpfr_get_d,:libmpfr), Float64, (Ptr{BigFloat},), &x)
+convert(::Type{Float32}, x::BigFloat) = ccall((:mpfr_get_flt,:libmpfr), Float32, (Ptr{BigFloat},), &x)
+
+promote_rule{T<:Real}(::Type{BigFloat}, ::Type{T}) = BigFloat
+
+promote_rule{T<:FloatingPoint}(::Type{BigInt},::Type{T}) = BigFloat
+promote_rule{T<:FloatingPoint}(::Type{BigFloat},::Type{T}) = BigFloat
+
+# Basic arithmetic without promotion
+# Unsigned addition
+function +(x::BigFloat, c::Culong)
+    z = BigFloat()
+    ccall((:mpfr_add_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
++(c::Culong, x::BigFloat) = x + c
++(c::Unsigned, x::BigFloat) = x + convert(Culong, c)
++(x::BigFloat, c::Unsigned) = x + convert(Culong, c)
+
+# Signed addition
+function +(x::BigFloat, c::Clong)
+    z = BigFloat()
+    ccall((:mpfr_add_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
++(c::Clong, x::BigFloat) = x + c
++(x::BigFloat, c::Signed) = x + convert(Clong, c)
++(c::Signed, x::BigFloat) = x + convert(Clong, c)
+
+# Float64 addition
+function +(x::BigFloat, c::Float64)
+    z = BigFloat()
+    ccall((:mpfr_add_d, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Float64, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
++(c::Float64, x::BigFloat) = x + c
++(c::Float32, x::BigFloat) = x + convert(Float64, c)
++(x::BigFloat, c::Float32) = x + convert(Float64, c)
+
+# BigInt addition
+function +(x::BigFloat, c::BigInt)
+    z = BigFloat()
+    ccall((:mpfr_add_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, ROUNDING_MODE[end])
+    return z
+end
++(c::BigInt, x::BigFloat) = x + c
+
+# Unsigned subtraction
+function -(x::BigFloat, c::Culong)
+    z = BigFloat()
+    ccall((:mpfr_sub_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
+function -(c::Culong, x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_ui_sub, :libmpfr), Int32, (Ptr{BigFloat}, Culong, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
+    return z
+end
+-(x::BigFloat, c::Unsigned) = -(x, convert(Culong, c))
+-(c::Unsigned, x::BigFloat) = -(convert(Culong, c), x)
+
+# Signed subtraction
+function -(x::BigFloat, c::Clong)
+    z = BigFloat()
+    ccall((:mpfr_sub_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
+function -(c::Clong, x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_si_sub, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
+    return z
+end
+-(x::BigFloat, c::Signed) = -(x, convert(Clong, c))
+-(c::Signed, x::BigFloat) = -(convert(Clong, c), x)
+
+# Float64 subtraction
+function -(x::BigFloat, c::Float64)
+    z = BigFloat()
+    ccall((:mpfr_sub_d, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Float64, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
+function -(c::Float64, x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_d_sub, :libmpfr), Int32, (Ptr{BigFloat}, Float64, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
+    return z
+end
+-(x::BigFloat, c::Float32) = -(x, convert(Float64, c))
+-(c::Float32, x::BigFloat) = -(convert(Float64, c), x)
+
+# BigInt subtraction
+function -(x::BigFloat, c::BigInt)
+    z = BigFloat()
+    ccall((:mpfr_sub_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, ROUNDING_MODE[end])
+    return z
+end
+function -(c::BigInt, x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_z_sub, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigInt}, Ptr{BigFloat}, Int32), &z, &c, &x, ROUNDING_MODE[end])
+    return z
+end
+
+# Unsigned multiplication
+function *(x::BigFloat, c::Culong)
+    z = BigFloat()
+    ccall((:mpfr_mul_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
+*(c::Culong, x::BigFloat) = x * c
+*(c::Unsigned, x::BigFloat) = x * convert(Culong, c)
+*(x::BigFloat, c::Unsigned) = x * convert(Culong, c)
+
+# Signed multiplication
+function *(x::BigFloat, c::Clong)
+    z = BigFloat()
+    ccall((:mpfr_mul_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
+*(c::Clong, x::BigFloat) = x * c
+*(x::BigFloat, c::Signed) = x * convert(Clong, c)
+
+# Float64 multiplication
+function *(x::BigFloat, c::Float64)
+    z = BigFloat()
+    ccall((:mpfr_mul_d, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Float64, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
+*(c::Float64, x::BigFloat) = x * c
+*(c::Float32, x::BigFloat) = x * convert(Float64, c)
+*(x::BigFloat, c::Float32) = x * convert(Float64, c)
+
+# BigInt multiplication
+*(c::Signed, x::BigFloat) = x * convert(Clong, c)
+function *(x::BigFloat, c::BigInt)
+    z = BigFloat()
+    ccall((:mpfr_mul_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, ROUNDING_MODE[end])
+    return z
+end
+*(c::BigInt, x::BigFloat) = x * c
+
+# Unsigned division
+function /(x::BigFloat, c::Culong)
+    z = BigFloat()
+    ccall((:mpfr_div_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
+function /(c::Culong, x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_ui_div, :libmpfr), Int32, (Ptr{BigFloat}, Culong, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
+    return z
+end
+/(x::BigFloat, c::Unsigned) = /(x, convert(Culong, c))
+/(c::Unsigned, x::BigFloat) = /(convert(Culong, c), x)
+
+# Signed division
+function /(x::BigFloat, c::Clong)
+    z = BigFloat()
+    ccall((:mpfr_div_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
+function /(c::Clong, x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_si_div, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
+    return z
+end
+/(x::BigFloat, c::Signed) = /(x, convert(Clong, c))
+/(c::Signed, x::BigFloat) = /(convert(Clong, c), x)
+
+# Float64 division
+function /(x::BigFloat, c::Float64)
+    z = BigFloat()
+    ccall((:mpfr_div_d, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Float64, Int32), &z, &x, c, ROUNDING_MODE[end])
+    return z
+end
+function /(c::Float64, x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_d_div, :libmpfr), Int32, (Ptr{BigFloat}, Float64, Ptr{BigFloat}, Int32), &z, c, &x, ROUNDING_MODE[end])
+    return z
+end
+/(x::BigFloat, c::Float32) = /(x, convert(Float64, c))
+/(c::Float32, x::BigFloat) = /(convert(Float64, c), x)
+
+# BigInt division
+function /(x::BigFloat, c::BigInt)
+    z = BigFloat()
+    ccall((:mpfr_div_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &c, ROUNDING_MODE[end])
+    return z
+end
+
+# Basic operations
+for (fJ, fC) in ((:+,:add), (:-,:sub), (:*,:mul), (:/,:div), (:^, :pow))
+    @eval begin 
+        function ($fJ)(x::BigFloat, y::BigFloat)
+            z = BigFloat()
+            ccall(($(string(:mpfr_,fC)),:libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
+            return z
+        end
+    end
+end
+
+function -(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_neg, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function cmp(x::BigFloat, y::BigFloat)
+    ccall((:mpfr_cmp, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y)
+end
+
+function sqrt(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_sqrt, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    if isnan(z)
+        throw(DomainError())
+    end
+    return z
+end
+
+for f in (:ceil, :floor, :trunc)
+    @eval begin
+        function ($f)(x::BigFloat)
+            z = BigFloat()
+            ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &z, &x)
+            return z
+        end
+    end
+end
+
+function ^(x::BigFloat, y::Unsigned)
+    z = BigFloat()
+    ccall((:mpfr_pow_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, y, ROUNDING_MODE[end])
+    return z
+end
+
+function ^(x::BigFloat, y::Signed)
+    z = BigFloat()
+    ccall((:mpfr_pow_si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, y, ROUNDING_MODE[end])
+    return z
+end
+
+function ^(x::BigFloat, y::BigInt)
+    z = BigFloat()
+    ccall((:mpfr_pow_z, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigInt}, Int32), &z, &x, &y, ROUNDING_MODE[end])
+    return z
+end
+
+function exp(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_exp, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function exp2(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_exp2, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function exp10(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_exp10, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function ldexp(x::BigFloat, n::Clong)
+    z = BigFloat()
+    ccall((:mpfr_mul_2si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, n, ROUNDING_MODE[end])
+    return z
+end
+function ldexp(x::BigFloat, n::Culong)
+    z = BigFloat()
+    ccall((:mpfr_mul_2ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32), &z, &x, n, ROUNDING_MODE[end])
+    return z
+end
+ldexp(x::BigFloat, n::Signed) = ldexp(x, convert(Clong, n))
+ldexp(x::BigFloat, n::Unsigned) = ldexp(x, convert(Culong, n))
+
+function besselj0(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_j0, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function besselj1(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_j1, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function besselj(n::Integer, x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_jn, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, n, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function bessely0(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_y0, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function bessely1(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_y1, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function bessely(n::Integer, x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_yn, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, n, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function factorial(x::BigFloat)
+    if x < 0 || !integer_valued(x)
+        throw(DomainError())
+    end
+    ui = convert(Culong, x)
+    z = BigFloat()
+    ccall((:mpfr_fac_ui, :libmpfr), Int32, (Ptr{BigFloat}, Culong, Int32), &z, ui, ROUNDING_MODE[end])
+    return z
+end
+
+function hypot(x::BigFloat, y::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_hypot, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
+    return z
+end
+
+function log(x::BigFloat)
+    if x < 0
+        throw(DomainError())
+    end
+    z = BigFloat()
+    ccall((:mpfr_log, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function log2(x::BigFloat)
+    if x < 0
+        throw(DomainError())
+    end
+    z = BigFloat()
+    ccall((:mpfr_log2, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function log10(x::BigFloat)
+    if x < 0
+        throw(DomainError())
+    end
+    z = BigFloat()
+    ccall((:mpfr_log10, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function max(x::BigFloat, y::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_max, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
+    return z
+end
+
+function min(x::BigFloat, y::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_min, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
+    return z
+end
+
+function modf(x::BigFloat)
+    if isinf(x)
+        return (BigFloat(NaN), x)
+    end
+    zint = BigFloat()
+    zfloat = BigFloat()
+    ccall((:mpfr_modf, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &zint, &zfloat, &x, ROUNDING_MODE[end])
+    return (zfloat, zint)
+end
+
+function rem(x::BigFloat, y::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_remainder, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
+    return z
+end
+
+# function sum{T<:BigFloat}(arr::AbstractArray{T})
+#     z = BigFloat()
+#     n = length(arr)
+#     ptrarr = [pointer(&x) for x in arr]
+#     ccall((:mpfr_sum, :libmpfr), Int32,
+#         (Ptr{BigFloat}, Ptr{Void}, Culong, Int32), 
+#         &z, ptrarr, n, ROUNDING_MODE[1])
+#     return z
+# end
+
+# Trigonometric functions
+# Every NaN is thrown as an error, and it follows somewhat closely
+# the Base functions behavior
+for f in (:sin,:cos,:tan,:sec,:csc,:cot,:acos,:asin,:atan,
+        :cosh,:sinh,:tanh,:sech,:csch,:coth,:acosh,:asinh,:atanh)
+    @eval begin
+        function ($f)(x::BigFloat)
+            z = BigFloat()
+            ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+            if isnan(z)
+                throw(DomainError())
+            end
+            return z
+        end
+    end
+end
+
+function atan2(y::BigFloat, x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_atan2, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &y, &x, ROUNDING_MODE[end])
+    return z
+end
+
+# Utility functions
+==(x::BigFloat, y::BigFloat) = ccall((:mpfr_equal_p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y) != 0
+<=(x::BigFloat, y::BigFloat) = ccall((:mpfr_lessequal_p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y) != 0
+>=(x::BigFloat, y::BigFloat) = ccall((:mpfr_greaterequal_p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y) != 0
+<(x::BigFloat, y::BigFloat) = ccall((:mpfr_less_p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y) != 0
+>(x::BigFloat, y::BigFloat) = ccall((:mpfr_greater_p, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}), &x, &y) != 0
+
+function get_precision(x::BigFloat)
+    return ccall((:mpfr_get_prec, :libmpfr), Clong, (Ptr{BigFloat},), &x)
+end
+
+get_bigfloat_precision() = DEFAULT_PRECISION[end]
+function set_bigfloat_precision(x::Int)
+    if x < 2
+        throw(DomainError())
+    end
+    DEFAULT_PRECISION[end] = x
+end
+
+get_bigfloat_rounding() = ROUNDING_MODE[end]
+function set_bigfloat_rounding(x::Int)
+    if x < 0 || x > 4
+        throw(DomainError())
+    end
+    ROUNDING_MODE[end] = x
+end
+
+function copysign(x::BigFloat, y::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_copysign, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
+    return z
+end
+
+function exponent(x::BigFloat)
+    if x == 0 || !isfinite(x)
+        throw(DomainError())
+    end
+    # The '- 1' is to make it work as Base.exponent
+    return ccall((:mpfr_get_exp, :libmpfr), Clong, (Ptr{BigFloat},), &x) - 1
+end
+
+function integer_valued(x::BigFloat)
+    return ccall((:mpfr_integer_p, :libmpfr), Int32, (Ptr{BigFloat},), &x) != 0
+end
+
+function iround(x::BigFloat)
+    fits = ccall((:mpfr_fits_slong_p, :libmpfr), Int32, (Ptr{BigFloat}, Int32), &x, ROUNDING_MODE[end])
+    if fits != 0
+        return ccall((:mpfr_get_si, :libmpfr), Clong, (Ptr{BigFloat}, Int32), &x, ROUNDING_MODE[end])
+    end
+    z = BigInt()
+    ccall((:mpfr_get_z, :libmpfr), Int32, (Ptr{BigInt}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+for (f,t) in ((:uint8,:Uint8), (:uint16,:Uint16), (:uint32,:Uint32),
+              (:int64,:Int64), (:uint64,:Uint64),
+              # requires int128/uint128(::Type{BigInt}) support
+              # (:int128,:Int128), (:uint128,:Uint128),
+              (:unsigned,:Uint), (:uint,:Uint))
+    @eval iround(::Type{$t}, x::BigFloat) = ($f)(iround(x))
+end
+
+isfinite(x::BigFloat) = !isinf(x)
+function isinf(x::BigFloat)
+    return ccall((:mpfr_inf_p, :libmpfr), Int32, (Ptr{BigFloat},), &x) != 0
+end
+
+function isnan(x::BigFloat)
+    return ccall((:mpfr_nan_p, :libmpfr), Int32, (Ptr{BigFloat},), &x) != 0
+end
+
+function nextfloat(x::BigFloat)
+   z = copy(x)
+   ccall((:mpfr_nextabove, :libmpfr), Int32, (Ptr{BigFloat},), &z) != 0
+   return z
+end
+
+function prevfloat(x::BigFloat)
+    z = copy(x)
+   ccall((:mpfr_nextbelow, :libmpfr), Int32, (Ptr{BigFloat},), &z) != 0
+   return z
+end
+
+function copy(x::BigFloat)
+    z = BigFloat()
+    ccall((:mpfr_set, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[end])
+    return z
+end
+
+function with_bigfloat_precision(f::Function, precision::Integer)
+    old_precision = get_bigfloat_precision()
+    set_bigfloat_precision(precision)
+    ret = f()
+    set_bigfloat_precision(old_precision)
+    return ret
+end
+
+function with_bigfloat_rounding(f::Function, rounding::Integer)
+    old_rounding = get_bigfloat_rounding()
+    set_bigfloat_rounding(rounding)
+    ret = f()
+    set_bigfloat_rounding(old_rounding)
+    return ret
+end
+
+function round(x::BigFloat, prec::Integer)
+    if prec < 1
+        throw(DomainError())
+    end
+    prec = int(ceil(log2(10^prec)))
+    z = BigFloat(x)
+    ccall((:mpfr_prec_round, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Int32), &z, prec, ROUNDING_MODE[end])
+    return z
+end
+
+function string(x::BigFloat)
+    lng = 128
+    for i = 1:2
+        z = Array(Uint8, lng)
+        lng = ccall((:mpfr_snprintf,:libmpfr), Int32, (Ptr{Uint8}, Culong, Ptr{Uint8}, Ptr{BigFloat}...), z, lng, "%.Re", &x)
+        if lng < 128 || i == 2
+            return bytestring(convert(Ptr{Uint8}, z[1:lng]))
+        end
+    end
+end
+
+show(io::IO, b::BigFloat) = print(io, string(b) * " with $(get_precision(b)) bits of precision")
+showcompact(io::IO, b::BigFloat) = print(io, string(b))
+
+end #module
diff --git a/base/sysimg.jl b/base/sysimg.jl
index 5c904b2f79488..25e4a21bc3577 100644
--- a/base/sysimg.jl
+++ b/base/sysimg.jl
@@ -172,6 +172,8 @@ importall .DSP
 # BigInts and BigFloats
 include("gmp.jl")
 importall .GMP
+include("mpfr.jl")
+importall .MPFR
 
 # deprecated functions
 include("deprecated.jl")
diff --git a/deps/.gitignore b/deps/.gitignore
index ac2b0f8b7af17..86ee7cfaca893 100644
--- a/deps/.gitignore
+++ b/deps/.gitignore
@@ -15,6 +15,7 @@
 /libunwind-*
 /lighttpd-*
 /llvm-*
+/mpfr-*
 /openblas-*
 /patchelf-*
 /pcre-*
diff --git a/deps/Makefile b/deps/Makefile
index e05965117d0d8..493dcdfc78d7c 100644
--- a/deps/Makefile
+++ b/deps/Makefile
@@ -21,7 +21,7 @@ ifeq ($(OS),WINNT)
 CONFIGURE_COMMON += LDFLAGS=-Wl,--stack,8388608
 endif
 
-#autoconf configure-driven scripts: llvm readline pcre arpack fftw unwind gmp patchelf
+#autoconf configure-driven scripts: llvm readline pcre arpack fftw unwind gmp mpfr patchelf
 #custom configure-driven script: zlib
 #custom Makefile rules: openlibm Rmath double-conversion random suitesparse-wrapper suitesparse lapack openblas uv
 
@@ -102,6 +102,10 @@ endif
 #STAGE1_DEPS += zlib
 #endif
 
+ifeq ($(USE_SYSTEM_MPFR), 0)
+STAGE1_DEPS += mpfr
+endif
+
 ifeq ($(USE_SYSTEM_ARPACK), 0)
 STAGE2_DEPS += arpack
 endif
@@ -134,7 +138,7 @@ install: $(addprefix install-, $(LIBS))
 cleanall: $(addprefix clean-, $(LIBS))
 distclean: $(addprefix distclean-, $(LIBS))
 	rm -rf $(BUILD)
-getall: get-llvm get-readline get-uv get-pcre get-double-conversion get-openlibm get-random get-openblas get-fftw get-suitesparse get-unwind get-gmp get-zlib get-patchelf
+getall: get-llvm get-readline get-uv get-pcre get-double-conversion get-openlibm get-random get-openblas get-fftw get-suitesparse get-unwind get-gmp get-mpfr get-zlib get-patchelf
 
 ## PATHS ##
 DIRS = $(addprefix $(BUILD)/,lib include bin share etc)
@@ -1182,6 +1186,44 @@ compile-gmp: $(GMP_SRC_TARGET)
 check-gmp: gmp-$(GMP_VER)/checked
 install-gmp: $(GMP_OBJ_TARGET)
 
+## MPFR ##
+
+MPFR_SRC_TARGET = mpfr-$(MPFR_VER)/.libs/libmpfr.$(SHLIB_EXT)
+MPFR_OBJ_TARGET = $(BUILD)/lib/libmpfr.$(SHLIB_EXT)
+
+mpfr-$(MPFR_VER).tar.bz2:
+	$(WGET_DASH_O) $@ http://www.mpfr.org/mpfr-current/$@
+mpfr-$(MPFR_VER)/configure: mpfr-$(MPFR_VER).tar.bz2
+	tar jxf $<
+	touch -c $@
+mpfr-$(MPFR_VER)/config.status: mpfr-$(MPFR_VER)/configure
+	cd mpfr-$(MPFR_VER) && \
+	./configure $(CONFIGURE_COMMON) F77= --enable-shared --disable-static
+	touch -c $@
+$(MPFR_SRC_TARGET): mpfr-$(MPFR_VER)/config.status
+	$(MAKE) -C mpfr-$(MPFR_VER) $(LIBTOOL_CCLD)
+mpfr-$(MPFR_VER)/checked: $(MPFR_SRC_TARGET)
+ifeq ($(OS),$(BUILD_OS))
+	$(MAKE) -C mpfr-$(MPFR_VER) $(LIBTOOL_CCLD) check
+endif
+	echo 1 > $@
+$(MPFR_OBJ_TARGET): $(MPFR_SRC_TARGET) mpfr-$(MPFR_VER)/checked
+	$(MAKE) -C mpfr-$(MPFR_VER) $(LIBTOOL_CCLD) install
+	$(INSTALL_NAME_CMD)libmpfr.dylib $@
+	touch -c $@
+
+clean-mpfr:
+	-$(MAKE) -C mpfr-$(MPFR_VER) clean
+	-rm -f $(MPFR_OBJ_TARGET)
+distclean-mpfr:
+	-rm -rf mpfr-$(MPFR_VER).tar.bz2 mpfr-$(MPFR_VER)
+
+get-mpfr: mpfr-$(MPFR_VER).tar.bz2
+configure-mpfr: mpfr-$(MPFR_VER)/config.status
+compile-mpfr: $(MPFR_SRC_TARGET)
+check-mpfr: mpfr-$(MPFR_VER)/checked
+install-mpfr: $(MPFR_OBJ_TARGET)
+
 ## ZLIB ##
 
 ifeq ($(OS),WINNT)
diff --git a/deps/Versions.make b/deps/Versions.make
index a1854bb3c9c7d..0480610a58c94 100644
--- a/deps/Versions.make
+++ b/deps/Versions.make
@@ -11,6 +11,7 @@ FFTW_VER = 3.3.3
 SUITESPARSE_VER = 4.2.0
 UNWIND_VER = 1.1
 GMP_VER=5.1.1
+MPFR_VER=3.1.2
 ZLIB_VER = 1.2.7
 PATCHELF_VER = 0.6
 GIT_VER = 1.8.2.1
diff --git a/test/bigfloat.jl b/test/bigfloat.jl
deleted file mode 100644
index d5efd69705a57..0000000000000
--- a/test/bigfloat.jl
+++ /dev/null
@@ -1,81 +0,0 @@
-tol = 1e-12
-
-a = BigFloat("12.34567890121")
-b = BigFloat("12.34567890122")
-
-@test typeof(a+1e-11) == BigFloat
-@test_approx_eq_eps a+1e-11 b tol
-@test !(b == a)
-@test b > a
-@test b >= a
-@test !(b < a)
-@test !(b <= a)
-
-c = BigFloat("24.69135780242")
-@test typeof(a * 2) == BigFloat
-@test_approx_eq_eps a*2 c tol
-@test_approx_eq_eps (c-a) a tol
-
-
-d = BigFloat("-24.69135780242")
-@test typeof(d) == BigFloat
-@test_approx_eq_eps d+c 0 tol
-
-#Multiple calls for sanity check, since we're doing direct memory manipulation
-@test string(a) == "1.234567890121e+01"
-@test string(b) == "1.234567890122e+01"
-@test string(c) == "2.469135780242e+01"
-@test string(d) == "-2.469135780242e+01"
-
-@test_approx_eq_eps (BigFloat(3)/BigFloat(2)) BigFloat(1.5) tol
-
-@test typeof(BigFloat(typemax(Int8))) == BigFloat
-@test typeof(BigFloat(typemax(Int16))) == BigFloat
-@test typeof(BigFloat(typemax(Int32))) == BigFloat
-@test typeof(BigFloat(typemax(Int64))) == BigFloat
-#@test typeof(BigFloat(typemax(Int128))) == BigFloat
-
-@test typeof(BigFloat(true)) == BigFloat
-@test typeof(BigFloat(typemax(Uint8))) == BigFloat
-@test typeof(BigFloat(typemax(Uint16))) == BigFloat
-@test typeof(BigFloat(typemax(Uint32))) == BigFloat
-@test typeof(BigFloat(typemax(Uint64))) == BigFloat
-#@test typeof(BigFloat(typemax(Uint128))) == BigFloat
-
-@test typeof(BigFloat(realmax(Float32))) == BigFloat
-@test typeof(BigFloat(realmax(Float64))) == BigFloat
-
-@test typeof(BigFloat(BigInt(1))) == BigFloat
-@test typeof(BigFloat(BigFloat(1))) == BigFloat
-
-@test typeof(BigFloat(1//1)) == BigFloat
-@test typeof(BigFloat(one(Rational{BigInt}))) == BigFloat
-
-f = BigFloat("1234567890.123")
-g = BigFloat("1234567891.123")
-
-tol = 1e-3
-
-@test_approx_eq_eps f+int8(1) g tol
-@test_approx_eq_eps f+int16(1) g tol
-@test_approx_eq_eps f+int32(1) g tol
-@test_approx_eq_eps f+int64(1) g tol
-#@test_approx_eq_eps f+int128(1) g tol
-
-@test_approx_eq_eps f+true g tol
-@test_approx_eq_eps f+uint8(1) g tol
-@test_approx_eq_eps f+uint16(1) g tol
-@test_approx_eq_eps f+uint32(1) g tol
-@test_approx_eq_eps f+uint64(1) g tol
-#@test_approx_eq_eps f+uint128(1) g tol
-
-@test_approx_eq_eps f+BigInt(1) g tol
-
-@test_approx_eq_eps f+1f0 g tol
-@test_approx_eq_eps f+1e0 g tol
-
-@test_approx_eq_eps f+BigFloat(1) g tol
-
-@test_approx_eq_eps f+(1//1) g tol
-
-@test_approx_eq_eps f+one(Rational{BigInt}) g tol
diff --git a/test/bigint.jl b/test/bigint.jl
index f8966d3884cc1..6fe5daf607ac8 100644
--- a/test/bigint.jl
+++ b/test/bigint.jl
@@ -54,18 +54,96 @@ end
 
 @test typeof(BigInt(BigInt(1))) == BigInt
 
+
+# Signed addition
 @test a+int8(1) == b
 @test a+int16(1) == b
 @test a+int32(1) == b
 @test a+int64(1) == b
-#@test a+int128(1) == b
-
+@test int8(1)+ a == b
+@test int16(1)+a == b
+@test int32(1)+a == b
+@test int64(1)+a == b
+@test b+int8(-1) == a
+@test b+int16(-1) == a
+@test b+int32(-1) == a
+@test b+int64(-1) == a
+@test int8(-1)+ b == a
+@test int16(-1)+b == a
+@test int32(-1)+b == a
+@test int64(-1)+b == a
+
+# Unsigned addition
 @test a+true == b
 @test a+uint8(1) == b
 @test a+uint16(1) == b
 @test a+uint32(1) == b
 @test a+uint64(1) == b
-#@test a+uint128(1) == b
+@test true+a == b
+@test uint8(1)+ a == b
+@test uint16(1)+a == b
+@test uint32(1)+a == b
+@test uint64(1)+a == b
+
+# Signed subtraction
+@test b-int8(1) == a
+@test b-int16(1) == a
+@test b-int32(1) == a
+@test b-int64(1) == a
+@test int8(1)- b == -a
+@test int16(1)-b == -a
+@test int32(1)-b == -a
+@test int64(1)-b == -a
+@test a-int8(-1) == b
+@test a-int16(-1) == b
+@test a-int32(-1) == b
+@test a-int64(-1) == b
+@test int8(-1)- a == -b
+@test int16(-1)-a == -b
+@test int32(-1)-a == -b
+@test int64(-1)-a == -b
+
+# Unsigned subtraction
+@test b-true == a
+@test b-uint8(1) == a
+@test b-uint16(1) == a
+@test b-uint32(1) == a
+@test b-uint64(1) == a
+@test true-b == -a
+@test uint8(1)- b == -a
+@test uint16(1)-b == -a
+@test uint32(1)-b == -a
+@test uint64(1)-b == -a
+
+# Signed multiplication
+@test a*int8(1) == a
+@test a*int16(1) == a
+@test a*int32(1) == a
+@test a*int64(1) == a
+@test int8(1)* a == a
+@test int16(1)*a == a
+@test int32(1)*a == a
+@test int64(1)*a == a
+@test a*int8(-1) == -a
+@test a*int16(-1) == -a
+@test a*int32(-1) == -a
+@test a*int64(-1) == -a
+@test int8(-1)* a == -a
+@test int16(-1)*a == -a
+@test int32(-1)*a == -a
+@test int64(-1)*a == -a
+
+# Unsigned multiplication
+@test a*true == a
+@test a*uint8(1) == a
+@test a*uint16(1) == a
+@test a*uint32(1) == a
+@test a*uint64(1) == a
+@test true*a == a
+@test uint8(1)* a == a
+@test uint16(1)*a == a
+@test uint32(1)*a == a
+@test uint64(1)*a == a
 
 @test a+BigInt(1) == b
 
diff --git a/test/mpfr.jl b/test/mpfr.jl
new file mode 100644
index 0000000000000..5a4eb20da27cc
--- /dev/null
+++ b/test/mpfr.jl
@@ -0,0 +1,606 @@
+# constructors
+with_bigfloat_precision(53) do
+    x = BigFloat()
+    x = BigFloat(12)
+end
+x = BigFloat(12)
+y = BigFloat(x)
+@test_approx_eq x y
+y = BigFloat(0xc)
+@test_approx_eq x y
+y = BigFloat(12.)
+@test_approx_eq x y
+y = BigFloat(BigInt(12))
+@test_approx_eq x y
+y = BigFloat(BigFloat(12))
+@test_approx_eq x y
+y = BigFloat("12")
+@test_approx_eq x y
+y = BigFloat(float32(12.))
+@test_approx_eq x y
+y = BigFloat(12//1)
+@test_approx_eq x y
+
+# +
+x = BigFloat(12)
+y = BigFloat(30)
+@test x + y == BigFloat(42)
+
+# -
+x = BigFloat(12)
+y = BigFloat(-30)
+@test x - y == BigFloat(42)
+
+# *
+x = BigFloat(6)
+y = BigFloat(9)
+@test x * y != BigFloat(42)
+@test x * y == BigFloat(54)
+
+# /
+x = BigFloat(9)
+y = BigFloat(6)
+@test x / y == BigFloat(9/6)
+
+# < / > / <= / >=
+x = BigFloat(12)
+y = BigFloat(42)
+z = BigFloat(30)
+@test y > x
+@test y >= x
+@test y > z
+@test y >= z
+@test x < y
+@test x <= y
+@test z < y
+@test z <= y
+@test y - x >= z
+@test y - x <= z
+@test !(x >= z)
+@test !(y <= z)
+
+# rounding modes
+with_bigfloat_precision(4) do
+    # default mode is round to nearest
+    down, up =  with_bigfloat_rounding(MPFR.RoundToNearest) do
+        BigFloat("0.0938"), BigFloat("0.102")
+    end
+    with_bigfloat_rounding(MPFR.RoundDown) do
+        @test BigFloat(0.1) == down
+        @test BigFloat(0.1) != up
+    end
+    with_bigfloat_rounding(MPFR.RoundUp) do
+        @test BigFloat(0.1) != down
+        @test BigFloat(0.1) == up
+    end
+end
+
+# ^
+x = BigFloat(12)
+y = BigFloat(4)
+@test x^y == BigFloat(20736)
+
+# ceil
+x = BigFloat(12.042)
+@test BigFloat(13) == ceil(x)
+
+# copysign
+x = BigFloat(1)
+y = BigFloat(-1)
+@test copysign(x, y) == y
+@test copysign(y, x) == x
+
+# isfinite / isinf
+x = BigFloat(Inf)
+y = BigFloat(1)
+@test isinf(x) == true
+@test isinf(y) == false
+@test isfinite(x) == false
+@test isinf(x) == true
+
+# isnan
+x = BigFloat(NaN)
+y = BigFloat(1)
+@test isnan(x) == true
+@test isnan(y) == false
+
+# convert to
+@test convert(BigFloat, 1//2) == BigFloat("0.5")
+@test convert(BigFloat, 0.5) == BigFloat("0.5")
+@test convert(BigFloat, 40) == BigFloat("40")
+@test convert(BigFloat, float32(0.5)) == BigFloat("0.5")
+
+# convert from
+@test convert(Float64, BigFloat(0.5)) == 0.5
+@test convert(Float32, BigFloat(0.5)) == float32(0.5)
+
+# exponent
+x = BigFloat(0)
+@test_fails exponent(x)
+x = BigFloat(Inf)
+@test_fails exponent(x)
+x = BigFloat(15.674)
+@test exponent(x) == exponent(15.674)
+
+# nextfloat/prevfloat should be immutable
+x = 12.
+y = BigFloat(x)
+@test x == y
+nextfloat(y)
+@test x == y
+prevfloat(y)
+@test x == y
+
+# sqrt DomainError
+@test_fails sqrt(BigFloat(-1))
+
+# precision
+old_precision = get_bigfloat_precision()
+x = BigFloat(0)
+@test get_precision(x) == old_precision
+set_bigfloat_precision(256)
+x = BigFloat(0)
+@test get_precision(x) == 256
+set_bigfloat_precision(old_precision)
+z = with_bigfloat_precision(240) do
+    z = x + 20
+    return z
+end
+@test float(z) == 20.
+@test get_precision(z) == 240
+x = BigFloat(12)
+@test get_precision(x) == old_precision
+@test_fails set_bigfloat_precision(1)
+
+# integer_valued
+@test integer_valued(BigFloat(12))
+@test !integer_valued(BigFloat(12.12))
+
+# nextfloat / prevfloat
+with_bigfloat_precision(53) do
+    x = BigFloat(12.12)
+    @test BigFloat(nextfloat(12.12)) == nextfloat(x)
+    @test BigFloat(prevfloat(12.12)) == prevfloat(x)
+end
+@test isnan(nextfloat(BigFloat(NaN)))
+@test isnan(prevfloat(BigFloat(NaN)))
+
+# comparisons
+x = BigFloat(1)
+y = BigFloat(-1)
+z = BigFloat(NaN)
+ipl = BigFloat(Inf)
+imi = BigFloat(-Inf)
+@test x > y
+@test x >= y
+@test x >= x
+@test y < x
+@test y <= x
+@test y <= y
+@test x < ipl
+@test x <= ipl
+@test x > imi
+@test x >= imi
+@test imi == imi
+@test ipl == ipl
+@test imi < ipl
+@test z != z
+@test !(z == z)
+@test !(z <= z)
+@test !(z < z)
+@test !(z >= z)
+@test !(z > z)
+
+# modf
+x = BigFloat(12)
+y = BigFloat(0.5)
+@test modf(x+y) == (y, x)
+x = BigFloat(NaN)
+@test map(isnan, modf(x)) == (true, true)
+x = BigFloat(Inf)
+y = modf(x)
+@test (isnan(y[1]), isinf(y[2])) == (true, true)
+
+# rem
+with_bigfloat_precision(53) do
+    x = BigFloat(2)
+    y = BigFloat(1.67)
+    @test rem(x,y) == rem(2, 1.67)
+    y = BigFloat(NaN)
+    @test isnan(rem(x,y))
+    @test isnan(rem(y,x))
+    y = BigFloat(Inf)
+    @test rem(x,y) == x
+    @test isnan(rem(y,x))
+end
+
+# min/max
+x = BigFloat(4)
+y = BigFloat(2)
+@test max(x,y) == x
+@test min(x,y) == y
+y = BigFloat(NaN)
+@test max(x,y) == x
+@test min(x,y) == x
+@test isnan(max(y,y))
+@test isnan(min(y,y))
+
+# sum
+x = BigFloat(1)
+y = BigFloat(2)
+z = BigFloat(3)
+w = BigFloat(4)
+@test sum([x,y,z,w]) == BigFloat(10)
+big_array = ones(BigFloat, 100)
+@test sum(big_array) == BigFloat(100)
+
+# promotion
+# the array converts everyone to the DEFAULT_PRECISION!
+x = BigFloat(12)
+y = with_bigfloat_precision(60) do
+    BigFloat(42)
+end
+@test [x,y] == [BigFloat(12), BigFloat(42)]
+
+# log / log2 / log10
+with_bigfloat_precision(53) do
+x = BigFloat(42)
+    @test log(x) == log(42)
+    @test isinf(log(BigFloat(0)))
+    @test_fails log(BigFloat(-1))
+    @test log2(x) == log2(42)
+    @test isinf(log2(BigFloat(0)))
+    @test_fails log2(BigFloat(-1))
+    @test log10(x) == log10(42)
+    @test isinf(log10(BigFloat(0)))
+    @test_fails log10(BigFloat(-1))
+end
+
+# exp / exp2 / exp10
+with_bigfloat_precision(53) do
+    x = BigFloat(10)
+    @test exp(x) == exp(10)
+    @test exp2(x) == 1024
+    @test exp10(x) == 10000000000
+end
+
+# convert to integer types
+x = BigFloat(12.1)
+y = BigFloat(42)
+@test_fails convert(Int32, x)
+@test_fails convert(Int64, x)
+@test_fails convert(BigInt, x)
+@test_fails convert(Uint32, x)
+@test_fails convert(Uint32, x)
+@test convert(Int32, y) == 42
+@test convert(Int64, y) == 42
+@test convert(BigInt, y) == 42
+@test convert(Uint32, y) == 42
+@test convert(Uint32, y) == 42
+
+# iround
+x = BigFloat(42.42)
+y = with_bigfloat_precision(256) do
+    BigFloat("9223372036854775809.2324")
+end
+z = BigInt("9223372036854775809")
+@test iround(x) == 42
+@test iround(y) == z
+@test typeof(iround(Uint8, x)) == Uint8 && iround(Uint8, x) == 0x2a
+@test typeof(iround(Uint16, x)) == Uint16 && iround(Uint16, x) == 0x2a
+@test typeof(iround(Uint32, x)) == Uint32 && iround(Uint32, x) == 0x2a
+@test typeof(iround(Uint64, x)) == Uint64 && iround(Uint64, x) == 0x2a
+@test typeof(iround(Int64, x)) == Int64 && iround(Int64, x) == 42
+@test typeof(iround(Int, x)) == Int && iround(Int, x) == 42
+@test typeof(iround(Uint, x)) == Uint && iround(Uint, x) == 0x2a
+
+# factorial
+with_bigfloat_precision(256) do
+    x = BigFloat(42)
+    @test factorial(x) == factorial(BigInt(42))
+    x = BigFloat(10)
+    @test factorial(x) == factorial(10)
+    @test_fails factorial(BigFloat(-1))
+    @test_fails factorial(BigFloat(331.3))
+end
+
+# bessel functions
+with_bigfloat_precision(53) do
+    @test_approx_eq besselj(4, BigFloat(2)) besselj(4, 2.)
+    @test_approx_eq besselj0(BigFloat(2))  besselj0(2.)
+    @test_approx_eq besselj1(BigFloat(2))  besselj1(2.)
+    @test_approx_eq bessely(4, BigFloat(2))  bessely(4, 2.)
+    @test_approx_eq bessely0(BigFloat(2))  bessely0(2.)
+    @test_approx_eq bessely1(BigFloat(2))  bessely1(2.)
+end
+
+# trigonometric functions
+with_bigfloat_precision(53) do
+    for f in (:sin,:cos,:tan,:sec,:csc,:cot,:acos,:asin,:atan,
+            :cosh,:sinh,:tanh,:sech,:csch,:coth,:asinh),
+        j in (-1., -0.5, -0.25, .25, .5, 1.)
+        @eval begin
+            @test_approx_eq ($f)(BigFloat($j)) ($f)($j)
+        end
+    end
+    for f in (:acos,:asin,:acosh,:atanh),
+        j in (-2, -1.5)
+        @eval begin
+            @test_fails ($f)(BigFloat($j))
+        end
+    end
+    for f in (:sin,:cos,:tan,:sec,:csc,:cot,:cosh,:sinh,:tanh,
+            :sech,:csch,:coth,:acosh,:asinh),
+        j in (1., 1.5, 1.9)
+        @eval begin
+            @test_approx_eq ($f)(BigFloat($j)) ($f)($j)
+        end
+    end
+    for j in (.25, .5)
+        @test_approx_eq atanh(BigFloat(j)) atanh(j)
+    end
+end
+
+# hypot
+@test hypot(BigFloat(3), BigFloat(4)) == 5
+
+# atan2
+with_bigfloat_precision(53) do
+    @test isequal(atan2(12,2), atan2(BigFloat(12), BigFloat(2)))
+end
+
+# ldexp
+with_bigfloat_precision(53) do
+    @test ldexp(BigFloat(24.5), 72) == ldexp(24.5, 72)
+    @test ldexp(BigFloat(24.5), int16(72)) == ldexp(24.5, 72)
+    @test ldexp(BigFloat(24.5), -72) == ldexp(24.5, -72)
+    @test ldexp(BigFloat(24.5), int16(-72)) == ldexp(24.5, -72)
+    @test ldexp(BigFloat(24.5), uint(72)) == ldexp(24.5, 72)
+    @test ldexp(BigFloat(24.5), 0x48) == ldexp(24.5, 72)
+end
+
+# basic arithmetic
+# Signed addition
+a = BigFloat("123456789012345678901234567890")
+b = BigFloat("123456789012345678901234567891")
+@test a+int8(1) == b
+@test a+int16(1) == b
+@test a+int32(1) == b
+@test a+int64(1) == b
+@test int8(1)+ a == b
+@test int16(1)+a == b
+@test int32(1)+a == b
+@test int64(1)+a == b
+@test b+int8(-1) == a
+@test b+int16(-1) == a
+@test b+int32(-1) == a
+@test b+int64(-1) == a
+@test int8(-1)+ b == a
+@test int16(-1)+b == a
+@test int32(-1)+b == a
+@test int64(-1)+b == a
+
+# Unsigned addition
+@test a+true == b
+@test a+uint8(1) == b
+@test a+uint16(1) == b
+@test a+uint32(1) == b
+@test a+uint64(1) == b
+@test true+a == b
+@test uint8(1)+ a == b
+@test uint16(1)+a == b
+@test uint32(1)+a == b
+@test uint64(1)+a == b
+
+# Float64 addition
+@test a + 1.0f0 == b
+@test 1.0f0 + a == b
+@test a + 1.0 == b
+@test 1.0 + a == b
+
+# BigInt addition
+@test a + BigInt(1) == b
+@test BigInt(1) + a == b
+
+# Signed subtraction
+@test b-int8(1) == a
+@test b-int16(1) == a
+@test b-int32(1) == a
+@test b-int64(1) == a
+@test int8(1)- b == -a
+@test int16(1)-b == -a
+@test int32(1)-b == -a
+@test int64(1)-b == -a
+@test a-int8(-1) == b
+@test a-int16(-1) == b
+@test a-int32(-1) == b
+@test a-int64(-1) == b
+@test int8(-1)- a == -b
+@test int16(-1)-a == -b
+@test int32(-1)-a == -b
+@test int64(-1)-a == -b
+
+# Unsigned subtraction
+@test b-true == a
+@test b-uint8(1) == a
+@test b-uint16(1) == a
+@test b-uint32(1) == a
+@test b-uint64(1) == a
+@test true-b == -a
+@test uint8(1)- b == -a
+@test uint16(1)-b == -a
+@test uint32(1)-b == -a
+@test uint64(1)-b == -a
+
+# Float64 subtraction
+@test b - 1.0f0 == a
+@test 1.0f0 - b == -a
+@test b - 1.0 == a
+@test 1.0 - b == -a
+
+# BigInt subtraction
+@test b - BigInt(1) == a
+@test BigInt(1) - b == -a
+
+# Signed multiplication
+@test a*int8(1) == a
+@test a*int16(1) == a
+@test a*int32(1) == a
+@test a*int64(1) == a
+@test int8(1)* a == a
+@test int16(1)*a == a
+@test int32(1)*a == a
+@test int64(1)*a == a
+@test a*int8(-1) == -a
+@test a*int16(-1) == -a
+@test a*int32(-1) == -a
+@test a*int64(-1) == -a
+@test int8(-1)* a == -a
+@test int16(-1)*a == -a
+@test int32(-1)*a == -a
+@test int64(-1)*a == -a
+
+# Unsigned multiplication
+@test a*true == a
+@test a*uint8(1) == a
+@test a*uint16(1) == a
+@test a*uint32(1) == a
+@test a*uint64(1) == a
+@test true*a == a
+@test uint8(1)* a == a
+@test uint16(1)*a == a
+@test uint32(1)*a == a
+@test uint64(1)*a == a
+
+# Float64 multiplication
+@test a * 1.0f0 == a
+@test 1.0f0 * a == a
+@test a * 1.0 == a
+@test 1.0 * a == a
+
+# BigInt multiplication
+@test a * BigInt(1) == a
+@test BigInt(1) * a == a
+
+# Signed division
+c = BigInt("61728394506172839450617283945")
+# d = 2^200
+d = BigFloat("1606938044258990275541962092341162602522202993782792835301376")
+f = BigFloat("6.223015277861141707144064053780124240590252168721167133101116614789698834035383e-61")
+
+@test a/int8(2) == c
+@test a/int16(2) == c
+@test a/int32(2) == c
+@test a/int64(2) == c
+@test int8(1)/ d == f
+@test int16(1)/d == f
+@test int32(1)/d == f
+@test int64(1)/d == f
+@test a/int8(-2) == -c
+@test a/int16(-2) == -c
+@test a/int32(-2) == -c
+@test a/int64(-2) == -c
+@test int8(-1)/ d == -f
+@test int16(-1)/d == -f
+@test int32(-1)/d == -f
+@test int64(-1)/d == -f
+
+# Unsigned division
+@test a/true == a
+@test a/uint8(2) == c
+@test a/uint16(2) == c
+@test a/uint32(2) == c
+@test a/uint64(2) == c
+@test true/d == f
+@test uint8(1)/ d == f
+@test uint16(1)/d == f
+@test uint32(1)/d == f
+@test uint64(1)/d == f
+
+# Float64 division
+@test a / 2.0f0 == c
+@test 1.0f0 / d == f
+@test a / 2.0 == c
+@test 1.0 / d == f
+
+# BigInt division
+@test a / BigInt(2) == c
+
+# old tests
+tol = 1e-12
+
+a = BigFloat("12.34567890121")
+b = BigFloat("12.34567890122")
+
+@test_approx_eq_eps a+1e-11 b tol
+@test !(b == a)
+@test b > a
+@test b >= a
+@test !(b < a)
+@test !(b <= a)
+
+c = BigFloat("24.69135780242")
+@test typeof(a * 2) == BigFloat
+@test_approx_eq_eps a*2 c tol
+@test_approx_eq_eps (c-a) a tol
+
+
+d = BigFloat("-24.69135780242")
+@test typeof(d) == BigFloat
+@test_approx_eq_eps d+c 0 tol
+
+@test_approx_eq_eps (BigFloat(3)/BigFloat(2)) BigFloat(1.5) tol
+
+@test typeof(BigFloat(typemax(Int8))) == BigFloat
+@test typeof(BigFloat(typemax(Int16))) == BigFloat
+@test typeof(BigFloat(typemax(Int32))) == BigFloat
+@test typeof(BigFloat(typemax(Int64))) == BigFloat
+#@test typeof(BigFloat(typemax(Int128))) == BigFloat
+
+@test typeof(BigFloat(true)) == BigFloat
+@test typeof(BigFloat(typemax(Uint8))) == BigFloat
+@test typeof(BigFloat(typemax(Uint16))) == BigFloat
+@test typeof(BigFloat(typemax(Uint32))) == BigFloat
+@test typeof(BigFloat(typemax(Uint64))) == BigFloat
+#@test typeof(BigFloat(typemax(Uint128))) == BigFloat
+
+@test typeof(BigFloat(realmax(Float32))) == BigFloat
+@test typeof(BigFloat(realmax(Float64))) == BigFloat
+
+@test typeof(BigFloat(BigInt(1))) == BigFloat
+@test typeof(BigFloat(BigFloat(1))) == BigFloat
+
+@test typeof(BigFloat(1//1)) == BigFloat
+@test typeof(BigFloat(one(Rational{BigInt}))) == BigFloat
+
+f = BigFloat("1234567890.123")
+g = BigFloat("1234567891.123")
+
+tol = 1e-3
+
+@test_approx_eq_eps f+int8(1) g tol
+@test_approx_eq_eps f+int16(1) g tol
+@test_approx_eq_eps f+int32(1) g tol
+@test_approx_eq_eps f+int64(1) g tol
+#@test_approx_eq_eps f+int128(1) g tol
+
+@test_approx_eq_eps f+true g tol
+@test_approx_eq_eps f+uint8(1) g tol
+@test_approx_eq_eps f+uint16(1) g tol
+@test_approx_eq_eps f+uint32(1) g tol
+@test_approx_eq_eps f+uint64(1) g tol
+#@test_approx_eq_eps f+uint128(1) g tol
+
+@test_approx_eq_eps f+BigInt(1) g tol
+
+@test_approx_eq_eps f+1f0 g tol
+@test_approx_eq_eps f+1e0 g tol
+
+@test_approx_eq_eps f+BigFloat(1) g tol
+
+@test_approx_eq_eps f+(1//1) g tol
+
+@test_approx_eq_eps f+one(Rational{BigInt}) g tol
+
+# new tests
+
diff --git a/test/runtests.jl b/test/runtests.jl
index 3e9f27f688eba..f180696d7afda 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,8 +3,8 @@ testnames = ["core", "keywordargs", "numbers", "strings", "unicode",
              "linalg", "blas", "fft", "dsp", "sparse", "bitarray",
              "random", "math", "functional", "bigint", "sorting",
              "statistics", "spawn", "parallel", "priorityqueue",
-             "arpack", "bigfloat", "file", "perf", "suitesparse", "version",
-             "pollfd"]
+             "arpack", "file", "perf", "suitesparse", "version",
+             "pollfd", "mpfr"]
 
 # Disabled: "complex"