diff --git a/README.md b/README.md
index 897dceb..7863f5f 100644
--- a/README.md
+++ b/README.md
@@ -42,6 +42,7 @@ operations.
 
 The following types implement the MutableArithmetics API:
 * The API is implemented for `Base.BigInt` in `src/bigint.jl`.
+* The API is implemented for `Base.BigFloat` in `src/bigfloat.jl`.
 * The API is implemented for `Base.Array` in `src/linear_algebra.jl`.
 * The interface for multivariate polynomials [MultivariatePolynomials](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl)
   as well as its two implementations [DynamicPolynomials](https://github.com/JuliaAlgebra/DynamicPolynomials.jl)
diff --git a/src/MutableArithmetics.jl b/src/MutableArithmetics.jl
index 2ea43ce..5abbab5 100644
--- a/src/MutableArithmetics.jl
+++ b/src/MutableArithmetics.jl
@@ -51,6 +51,7 @@ end
 scaling_convert(T::Type, x::LinearAlgebra.UniformScaling) = convert(T, x.λ)
 scaling_convert(T::Type, x) = convert(T, x)
 include("bigint.jl")
+include("bigfloat.jl")
 include("linear_algebra.jl")
 include("sparse_arrays.jl")
 
diff --git a/src/bigfloat.jl b/src/bigfloat.jl
new file mode 100644
index 0000000..44ba20a
--- /dev/null
+++ b/src/bigfloat.jl
@@ -0,0 +1,78 @@
+mutability(::Type{BigFloat}) = IsMutable()
+mutable_copy(x::BigFloat) = deepcopy(x)
+
+@static if VERSION >= v"1.1.0-DEV.683"
+    const MPFRRoundingMode = Base.MPFR.MPFRRoundingMode
+else
+    const MPFRRoundingMode = Int32
+end
+
+# zero
+promote_operation(::typeof(zero), ::Type{BigFloat}) = BigFloat
+function _set_si!(x::BigFloat, value)
+    ccall((:mpfr_set_si, :libmpfr), Int32, (Ref{BigFloat}, Clong, MPFRRoundingMode), x, value, Base.MPFR.ROUNDING_MODE[])
+    return x
+end
+mutable_operate!(::typeof(zero), x::BigFloat) = _set_si!(x, 0)
+
+# one
+promote_operation(::typeof(one), ::Type{BigFloat}) = BigFloat
+mutable_operate!(::typeof(one), x::BigFloat) = _set_si!(x, 1)
+
+# +
+promote_operation(::typeof(+), ::Vararg{Type{BigFloat}, N}) where {N} = BigFloat
+function mutable_operate_to!(output::BigFloat, ::typeof(+), a::BigFloat, b::BigFloat)
+    ccall((:mpfr_add, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode), output, a, b, Base.MPFR.ROUNDING_MODE[])
+    return output
+end
+#function mutable_operate_to!(output::BigFloat, op::typeof(+), a::BigFloat, b::LinearAlgebra.UniformScaling)
+#    return mutable_operate_to!(output, op, a, b.λ)
+#end
+
+# *
+promote_operation(::typeof(*), ::Vararg{Type{BigFloat}, N}) where {N} = BigFloat
+function mutable_operate_to!(output::BigFloat, ::typeof(*), a::BigFloat, b::BigFloat)
+    ccall((:mpfr_mul, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode), output, a, b, Base.MPFR.ROUNDING_MODE[])
+    return output
+end
+
+function mutable_operate_to!(output::BigFloat, op::Union{typeof(*), typeof(+)},
+                             a::BigFloat, b::BigFloat, c::Vararg{BigFloat, N}) where N
+    mutable_operate_to!(output, op, a, b)
+    return mutable_operate!(op, output, c...)
+end
+function mutable_operate!(op::Function, x::BigFloat, args::Vararg{Any, N}) where N
+    mutable_operate_to!(x, op, x, args...)
+end
+
+# add_mul
+# Buffer to hold the product
+buffer_for(::typeof(add_mul), args::Vararg{Type{BigFloat}, N}) where {N} = BigFloat()
+function mutable_operate_to!(output::BigFloat, ::typeof(add_mul), x::BigFloat, y::BigFloat, z::BigFloat, args::Vararg{BigFloat, N}) where N
+    return mutable_buffered_operate_to!(BigFloat(), output, add_mul, x, y, z, args...)
+end
+
+function mutable_buffered_operate_to!(buffer::BigFloat, output::BigFloat, ::typeof(add_mul),
+                                      a::BigFloat, x::BigFloat, y::BigFloat, args::Vararg{BigFloat, N}) where N
+    mutable_operate_to!(buffer, *, x, y, args...)
+    return mutable_operate_to!(output, +, a, buffer)
+end
+function mutable_buffered_operate!(buffer::BigFloat, op::typeof(add_mul), x::BigFloat, args::Vararg{Any, N}) where N
+    return mutable_buffered_operate_to!(buffer, x, op, x, args...)
+end
+
+scaling_to_bigfloat(x::BigFloat) = x
+scaling_to_bigfloat(x::Number) = convert(BigFloat, x)
+scaling_to_bigfloat(J::LinearAlgebra.UniformScaling) = scaling_to_bigfloat(J.λ)
+function mutable_operate_to!(output::BigFloat, op::Union{typeof(+), typeof(*)}, args::Vararg{Scaling, N}) where N
+    return mutable_operate_to!(output, op, scaling_to_bigfloat.(args)...)
+end
+function mutable_operate_to!(output::BigFloat, op::typeof(add_mul), x::Scaling, y::Scaling, z::Scaling, args::Vararg{Scaling, N}) where N
+    return mutable_operate_to!(
+        output, op, scaling_to_bigfloat(x), scaling_to_bigfloat(y),
+        scaling_to_bigfloat(z), scaling_to_bigfloat.(args)...)
+end
+# Called for instance if `args` is `(v', v)` for a vector `v`.
+function mutable_operate_to!(output::BigFloat, op::typeof(add_mul), x, y, z, args::Vararg{Any, N}) where N
+    return mutable_operate_to!(output, +, x, *(y, z, args...))
+end
diff --git a/test/bigfloat.jl b/test/bigfloat.jl
new file mode 100644
index 0000000..b478702
--- /dev/null
+++ b/test/bigfloat.jl
@@ -0,0 +1,11 @@
+MA.Test.int_test(BigFloat)
+
+@testset "Allocation" begin
+    a = BigFloat(2)
+    b = BigFloat(3)
+    c = BigFloat(4)
+    alloc_test(() -> MA.add!(a, b), 0)
+    alloc_test(() -> MA.add_to!(c, a, b), 0)
+    alloc_test(() -> MA.mul!(a, b), 0)
+    alloc_test(() -> MA.mul_to!(c, a, b), 0)
+end
diff --git a/test/bigint.jl b/test/bigint.jl
index 1f2f202..ea544b9 100644
--- a/test/bigint.jl
+++ b/test/bigint.jl
@@ -6,4 +6,6 @@ MA.Test.int_test(BigInt)
     c = BigInt(4)
     alloc_test(() -> MA.add!(a, b), 0)
     alloc_test(() -> MA.add_to!(c, a, b), 0)
+    alloc_test(() -> MA.mul!(a, b), 0)
+    alloc_test(() -> MA.mul_to!(c, a, b), 0)
 end
diff --git a/test/rewrite.jl b/test/rewrite.jl
index 276e564..0f3d37c 100644
--- a/test/rewrite.jl
+++ b/test/rewrite.jl
@@ -40,6 +40,7 @@ using OffsetArrays
         (Int, true),
         (Float64, true),
         (BigInt, true),
+        (BigFloat, true),
         (DummyBigInt, false)
     ]
     @testset "Error" begin
diff --git a/test/runtests.jl b/test/runtests.jl
index 8415199..fdca76f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,6 +12,9 @@ end
 @testset "BigInt" begin
     include("bigint.jl")
 end
+@testset "BigFloat" begin
+    include("bigfloat.jl")
+end
 @testset "Broadcast" begin
     include("broadcast.jl")
 end