diff --git a/src/itensornetworks_elementary_functions.jl b/src/itensornetworks_elementary_functions.jl
index ae6392a..0fe1011 100644
--- a/src/itensornetworks_elementary_functions.jl
+++ b/src/itensornetworks_elementary_functions.jl
@@ -17,14 +17,15 @@ default_k_value() = 1.0
 default_nterms() = 20
 default_dimension() = 1
 
-"""Construct the product state representation of the function f(x) = const."""
+"""Build a representation of the function f(x,y,z,...) = c, with flexible choice of linkdim"""
 function const_itensornetwork(
-  s::IndsNetwork, bit_map; c::Union{Float64,ComplexF64}=default_c_value()
+  s::IndsNetwork, bit_map; c::Union{Float64,ComplexF64}=default_c_value(), linkdim::Int64=1
 )
-  ψ = copy_tensor_network(s)
+  ψ = random_tensornetwork(s; link_space=linkdim)
   inv_L = Float64(1.0 / nv(s))
   for v in vertices(ψ)
-    ψ[v] *= c^inv_L
+    virt_inds = setdiff(inds(ψ[v]), Index[only(s[v])])
+    ψ[v] = c^inv_L * c_tensor(only(s[v]), virt_inds)
   end
 
   return ITensorNetworkFunction(ψ, bit_map)
@@ -216,7 +217,7 @@ function polynomial_itensornetwork(
   root_vertices_dim = filter(v -> v ∈ dimension_vertices, leaf_vertices(s_tree))
   @assert !isempty(root_vertices_dim)
   root_vertex = first(root_vertices_dim)
-  ψ = copy_tensor_network(s_tree; linkdim=n + 1)
+  ψ = const_itensornetwork(s_tree, bit_map; linkdim=n + 1)
   #Place the Q_n tensors, making sure we get the right index pointing towards the root
   for v in dimension_vertices
     siteindex = s_tree[v][]
@@ -245,7 +246,7 @@ function polynomial_itensornetwork(
     end
   end
 
-  return ITensorNetworkFunction(ψ, bit_map)
+  return ψ
 end
 
 function random_itensornetwork(s::IndsNetwork, bit_map; kwargs...)
diff --git a/src/itensornetworksutils.jl b/src/itensornetworksutils.jl
index 9d4fe15..05ac429 100644
--- a/src/itensornetworksutils.jl
+++ b/src/itensornetworksutils.jl
@@ -1,7 +1,7 @@
 using ITensors: Index, dim, inds
 using ITensorNetworks: random_tensornetwork, IndsNetwork
 
-"""Build the order L tensor corresponding to fx(x): x ∈ [0,1]."""
+"""Build the order L tensor corresponding to fx(x): x ∈ [0,1], default decomposition is binary"""
 function build_full_rank_tensor(L::Int64, fx::Function; base::Int64=2)
   inds = [Index(base, "$i") for i in 1:L]
   dims = Tuple([base for i in 1:L])
@@ -15,10 +15,10 @@ function build_full_rank_tensor(L::Int64, fx::Function; base::Int64=2)
   return ITensor(array, inds)
 end
 
+"""Build the tensor C such that C_{phys_ind, virt_inds...} = delta_{virt_inds...}"""
 function c_tensor(phys_ind::Index, virt_inds::Vector)
   inds = vcat(phys_ind, virt_inds)
   @assert allequal(dim.(virt_inds))
-  #Build tensor to be delta on inds and independent of phys_ind
   T = ITensor(0.0, inds...)
   for i in 1:dim(phys_ind)
     for j in 1:dim(first(virt_inds))
@@ -29,13 +29,3 @@ function c_tensor(phys_ind::Index, virt_inds::Vector)
 
   return T
 end
-
-function copy_tensor_network(s::IndsNetwork; linkdim::Int64=1)
-  tn = random_tensornetwork(s; link_space=linkdim)
-  for v in vertices(tn)
-    virt_inds = setdiff(inds(tn[v]), Index[only(s[v])])
-    tn[v] = c_tensor(only(s[v]), virt_inds)
-  end
-
-  return tn
-end