diff --git a/Project.toml b/Project.toml
index c51ef7b..09d78a6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ repo = "https://github.com/JuliaAlgebra/MultivariateBases.jl.git"
 version = "0.2.1"
 
 [deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 
diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl
index f22122b..dc84aed 100644
--- a/src/MultivariateBases.jl
+++ b/src/MultivariateBases.jl
@@ -29,6 +29,8 @@ export generators,
     reccurence_third_coef,
     reccurence_deno_coef
 include("fixed.jl")
+
+import LinearAlgebra
 include("orthogonal.jl")
 include("hermite.jl")
 include("laguerre.jl")
diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index 86ccc87..83dac46 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -28,6 +28,17 @@ function degree_one_univariate_polynomial(
     MA.@rewrite(variable + 0)
 end
 
+function _scalar_product_function(::Type{<:ChebyshevBasisFirstKind}, i::Int)
+    if i == 0
+        return π
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return (π / 2^i) * prod(n+1:i) / factorial(n)
+    end
+end
+
 """
     struct ChebyshevBasisSecondKind{P} <: AbstractChebyshevBasis{P}
         polynomials::Vector{P}
@@ -45,3 +56,14 @@ function degree_one_univariate_polynomial(
 )
     MA.@rewrite(2variable + 0)
 end
+
+function _scalar_product_function(::Type{<:ChebyshevBasisSecondKind}, i::Int)
+    if i == 0
+        return π / 2
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return π / (2^(i + 1)) * prod(n+2:i) / factorial(n)
+    end
+end
diff --git a/src/hermite.jl b/src/hermite.jl
index e7322c7..cf64476 100644
--- a/src/hermite.jl
+++ b/src/hermite.jl
@@ -30,6 +30,17 @@ function degree_one_univariate_polynomial(
     MA.@rewrite(1variable)
 end
 
+function _scalar_product_function(::Type{<:ProbabilistsHermiteBasis}, i::Int)
+    if i == 0
+        return √(2 * π)
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return (√(2 * π) / (2^n)) * prod(n+1:2*n)
+    end
+end
+
 """
     struct PhysicistsHermiteBasis{P} <: AbstractHermiteBasis{P}
         polynomials::Vector{P}
@@ -48,3 +59,14 @@ function degree_one_univariate_polynomial(
 )
     MA.@rewrite(2variable)
 end
+
+function _scalar_product_function(::Type{<:PhysicistsHermiteBasis}, i::Int)
+    if i == 0
+        return √(π)
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return (√(π) / (2^i)) * prod(n+1:i)
+    end
+end
diff --git a/src/laguerre.jl b/src/laguerre.jl
index e731b70..fcec105 100644
--- a/src/laguerre.jl
+++ b/src/laguerre.jl
@@ -26,3 +26,7 @@ function degree_one_univariate_polynomial(
 )
     MA.@rewrite(1 - variable)
 end
+
+function _scalar_product_function(::Type{<:LaguerreBasis}, i::Int)
+    return factorial(i)
+end
diff --git a/src/legendre.jl b/src/legendre.jl
index 0582307..bb9b8e1 100644
--- a/src/legendre.jl
+++ b/src/legendre.jl
@@ -35,3 +35,11 @@ function degree_one_univariate_polynomial(
 )
     MA.@rewrite(variable + 0)
 end
+
+function _scalar_product_function(::Type{<:LegendreBasis}, i::Int)
+    if isodd(i)
+        return 0
+    else
+        return 2 / (i + 1)
+    end
+end
diff --git a/src/monomial.jl b/src/monomial.jl
index 9a38316..ee2dfc6 100644
--- a/src/monomial.jl
+++ b/src/monomial.jl
@@ -91,6 +91,10 @@ function MP.polynomial(Q::AbstractMatrix, mb::MonomialBasis, T::Type)
     return MP.polynomial(Q, mb.monomials, T)
 end
 
+function MP.coefficients(p, basis::MonomialBasis)
+    return MP.coefficients(p, basis.monomials)
+end
+
 function MP.coefficients(p, ::Type{<:MonomialBasis})
     return MP.coefficients(p)
 end
diff --git a/src/orthogonal.jl b/src/orthogonal.jl
index 4b3c7e4..c40eff6 100644
--- a/src/orthogonal.jl
+++ b/src/orthogonal.jl
@@ -158,3 +158,67 @@ function basis_covering_monomials(
         MP.monomial_vector(collect(m)),
     )
 end
+
+function _scalar_product_function(
+    ::Type{<:AbstractMultipleOrthogonalBasis},
+    i::Int,
+) end
+
+function LinearAlgebra.dot(
+    p,
+    q,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return _integral(p * q, basis_type)
+end
+
+function _integral(
+    p::Number,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return p * _scalar_product_function(basis_type, 0)
+end
+
+function _integral(
+    p::MP.AbstractVariable,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return _scalar_product_function(basis_type, 1)
+end
+
+function _integral(
+    p::MP.AbstractMonomial,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return prod([
+        _scalar_product_function(basis_type, i) for i in MP.exponents(p)
+    ])
+end
+
+function _integral(
+    p::MP.AbstractTerm,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return MP.coefficient(p) * _integral(MP.monomial(p), basis_type)
+end
+
+function _integral(
+    p::MP.AbstractPolynomial,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return sum([_integral(t, basis_type) for t in MP.terms(p)])
+end
+
+function MP.coefficients(
+    p,
+    basis::AbstractMultipleOrthogonalBasis;
+    check = true,
+)
+    B = typeof(basis)
+    coeffs = [
+        LinearAlgebra.dot(p, el, B) / LinearAlgebra.dot(el, el, B) for
+        el in basis
+    ]
+    idx = findall(c -> !isapprox(c, 0; atol = 1e-10), coeffs)
+    return coeffs[idx]
+end
diff --git a/src/scaled.jl b/src/scaled.jl
index c4e634d..b3483e9 100644
--- a/src/scaled.jl
+++ b/src/scaled.jl
@@ -95,3 +95,6 @@ unscale_coef(t::MP.AbstractTerm) = MP.coefficient(t) / scaling(MP.monomial(t))
 function MP.coefficients(p, ::Type{<:ScaledMonomialBasis})
     return unscale_coef.(MP.terms(p))
 end
+function MP.coefficients(p, basis::ScaledMonomialBasis)
+    return MP.coefficients(p, basis.monomials) ./ scaling.(MP.monomials(p))
+end
diff --git a/test/Project.toml b/test/Project.toml
index 97b28c0..04cebf8 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,5 @@
 [deps]
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/chebyshev.jl b/test/chebyshev.jl
index 512c714..1616923 100644
--- a/test/chebyshev.jl
+++ b/test/chebyshev.jl
@@ -7,14 +7,50 @@ using MultivariateBases
         x -> [1, x, 2x^2 - 1, 4x^3 - 3x, 8x^4 - 8x^2 + 1],
         true,
     )
+    univ_orthogonal_test(
+        ChebyshevBasis,
+        i -> π * (1 / 2 + (0^i) / 2);
+        atol = 1e-12,
+    )
     orthogonal_test(
         ChebyshevBasisSecondKind,
         x -> [1, 2x, 4x^2 - 1, 8x^3 - 4x, 16x^4 - 12x^2 + 1],
         true,
     )
+
+    univ_orthogonal_test(ChebyshevBasisSecondKind, i -> π / 2; atol = 1e-12)
 end
 
 @testset "API degree = $degree" for degree in 0:3
     api_test(ChebyshevBasis, degree)
     api_test(ChebyshevBasisSecondKind, degree)
 end
+
+@testset "Coefficients" begin
+    coefficient_test(
+        ChebyshevBasis,
+        reverse([
+            0.0625,
+            0.0625,
+            0.0625,
+            -0.25,
+            0.0625,
+            -0.3125,
+            -0.3125,
+            0.625,
+        ]),
+    )
+    coefficient_test(
+        ChebyshevBasisSecondKind,
+        reverse([
+            0.015625,
+            0.015625,
+            0.015625,
+            -0.09375,
+            0.015625,
+            -0.109375,
+            -0.109375,
+            0.875,
+        ]),
+    )
+end
diff --git a/test/fixed.jl b/test/fixed.jl
index d63c574..a0455b0 100644
--- a/test/fixed.jl
+++ b/test/fixed.jl
@@ -72,3 +72,11 @@ end
     p = @inferred polynomial(zeros(Int, 0, 0), basis, Int)
     @test iszero(p)
 end
+
+@testset "Enumerate" begin
+    monos = [1, x, y^2]
+    basis = FixedPolynomialBasis(monos)
+    for (i, e) in enumerate(basis)
+        @test e == monos[i]
+    end
+end
diff --git a/test/hermite.jl b/test/hermite.jl
index 38222ca..12f9fdb 100644
--- a/test/hermite.jl
+++ b/test/hermite.jl
@@ -7,14 +7,41 @@ using MultivariateBases
         x -> [1, x, x^2 - 1, x^3 - 3x, x^4 - 6x^2 + 3],
         true,
     )
+    univ_orthogonal_test(ProbabilistsHermiteBasis, i -> √(2 * π) * factorial(i))
     orthogonal_test(
         PhysicistsHermiteBasis,
         x -> [1, 2x, 4x^2 - 2, 8x^3 - 12x, 16x^4 - 48x^2 + 12],
         true,
     )
+    univ_orthogonal_test(
+        PhysicistsHermiteBasis,
+        i -> √(π) * factorial(i) * 2^i;
+        atol = 1e-12,
+    ) #precision issue
 end
 
 @testset "API degree = $degree" for degree in 0:3
     api_test(ProbabilistsHermiteBasis, degree)
     api_test(PhysicistsHermiteBasis, degree)
 end
+
+@testset "Coefficients" begin
+    coefficient_test(
+        ProbabilistsHermiteBasis,
+        [4, 6, 6, 1, 9, 1, 1, 1];
+        atol = 1e-12,
+    )
+    coefficient_test(
+        PhysicistsHermiteBasis,
+        reverse([
+            0.015625,
+            0.015625,
+            0.03125,
+            0.1875,
+            0.03125,
+            0.1875,
+            0.1875,
+            1.0,
+        ]),
+    )
+end
diff --git a/test/laguerre.jl b/test/laguerre.jl
index cdad6e0..76c3c75 100644
--- a/test/laguerre.jl
+++ b/test/laguerre.jl
@@ -13,8 +13,39 @@ using MultivariateBases
         ],
         false,
     )
+
+    univ_orthogonal_test(LaguerreBasis, i -> 1; atol = 1e-12) # there are precision issues
 end
 
 @testset "API degree = $degree" for degree in 0:3
     api_test(LaguerreBasis, degree)
 end
+
+@testset "Coefficients" begin
+    coefficient_test(
+        LaguerreBasis,
+        reverse([
+            48.0,
+            48.0,
+            -96.0,
+            -192.0,
+            -192.0,
+            -96.0,
+            48.0,
+            384.0,
+            564.0,
+            384.0,
+            48.0,
+            -192.0,
+            -744.0,
+            -744.0,
+            -192.0,
+            324.0,
+            720.0,
+            324.0,
+            -264.0,
+            -264.0,
+            85.0,
+        ]),
+    )
+end
diff --git a/test/legendre.jl b/test/legendre.jl
index f72fbdb..46e3c20 100644
--- a/test/legendre.jl
+++ b/test/legendre.jl
@@ -7,8 +7,26 @@ using MultivariateBases
         x -> [1, x, (3x^2 - 1) / 2, (5x^3 - 3x) / 2, (35x^4 - 30x^2 + 3) / 8],
         true,
     )
+
+    univ_orthogonal_test(LegendreBasis, i -> 2 / (2 * i + 1))
 end
 
 @testset "API degree = $degree" for degree in 0:3
     api_test(LegendreBasis, degree)
 end
+
+@testset "Coefficients" begin
+    coefficient_test(
+        LegendreBasis,
+        reverse([
+            0.1523809523809524,
+            0.1523809523809524,
+            0.0761904761904762,
+            -0.5714285714285714,
+            0.0761904761904762,
+            -0.3428571428571428,
+            -0.3428571428571428,
+            0.8,
+        ]),
+    )
+end
diff --git a/test/monomial.jl b/test/monomial.jl
index a0e6f2a..3743427 100644
--- a/test/monomial.jl
+++ b/test/monomial.jl
@@ -66,3 +66,15 @@ end
     p = @inferred polynomial(zeros(Int, 0, 0), basis, Int)
     @test iszero(p)
 end
+
+@testset "Enumerate" begin
+    monos = [1, y, x]
+    basis = MonomialBasis(monos)
+    for (i, e) in enumerate(basis)
+        @test e == monos[i]
+    end
+end
+
+@testset "Coefficients" begin
+    coefficient_test(MonomialBasis, [1, -3, 1, 1])
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index fc9afee..3bd703b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,7 @@
 using Test
 
 using MultivariateBases
-
+using LinearAlgebra
 using DynamicPolynomials
 
 function api_test(B::Type{<:AbstractPolynomialBasis}, degree)
@@ -27,6 +27,25 @@ function api_test(B::Type{<:AbstractPolynomialBasis}, degree)
     end
 end
 
+function univ_orthogonal_test(
+    B::Type{<:AbstractMultipleOrthogonalBasis},
+    univ::Function;
+    kwargs...,
+)
+    @polyvar x
+    basis = maxdegree_basis(B, [x], 4)
+    for i in eachindex(basis)
+        @test isapprox(
+            dot(basis[i], basis[i], B),
+            univ(maxdegree(basis[i]));
+            kwargs...,
+        )
+        for j in 1:i-1
+            @test isapprox(dot(basis[i], basis[j], B), 0.0; kwargs...)
+        end
+    end
+end
+
 function orthogonal_test(
     B::Type{<:AbstractMultipleOrthogonalBasis},
     univ::Function,
@@ -68,6 +87,15 @@ function orthogonal_test(
     end
 end
 
+function coefficient_test(B::Type{<:AbstractPolynomialBasis}, coefs; kwargs...)
+    @polyvar x y
+    p = x^4 * y^2 + x^2 * y^4 - 3 * x^2 * y^2 + 1
+    basis = basis_covering_monomials(B, monomials(p))
+    cc = coefficients(p, basis)
+    @test isapprox(coefs, cc; kwargs...)
+    @test isapprox(p, polynomial(cc, basis); kwargs...)
+end
+
 @testset "Monomial" begin
     include("monomial.jl")
 end
diff --git a/test/scaled.jl b/test/scaled.jl
index bec5e74..ddd8c95 100644
--- a/test/scaled.jl
+++ b/test/scaled.jl
@@ -29,3 +29,15 @@ end
 @testset "API degree = $degree" for degree in 0:3
     api_test(ScaledMonomialBasis, degree)
 end
+
+@testset "Enumerate" begin
+    monos = [1, y, x]
+    basis = ScaledMonomialBasis(monos)
+    for (i, e) in enumerate(basis)
+        @test e == monos[i]
+    end
+end
+
+@testset "Coefficients" begin
+    coefficient_test(ScaledMonomialBasis, [1, -√3 / √2, 1 / √15, 1 / √15])
+end