diff --git a/src/ITensorNumericalAnalysis.jl b/src/ITensorNumericalAnalysis.jl
index a235e02..4df1a0a 100644
--- a/src/ITensorNumericalAnalysis.jl
+++ b/src/ITensorNumericalAnalysis.jl
@@ -28,7 +28,8 @@ export AbstractIndexMap,
   calculate_ind_values,
   dimension,
   dimensions,
-  grid_points
+  grid_points,
+  rand_p
 export IndsNetworkMap,
   continuous_siteinds,
   complex_continuous_siteinds,
diff --git a/src/IndexMaps/abstractindexmap.jl b/src/IndexMaps/abstractindexmap.jl
index 916731b..c6cb9cd 100644
--- a/src/IndexMaps/abstractindexmap.jl
+++ b/src/IndexMaps/abstractindexmap.jl
@@ -2,7 +2,7 @@ using Base: Base
 using Dictionaries: Dictionary, set!
 using ITensors: ITensors, Index, dim
 using ITensorNetworks: IndsNetwork, vertex_data
-
+using Random: AbstractRNG, default_rng
 abstract type AbstractIndexMap{VB,VD} end
 
 #These functions need to be defined on the concrete type for implementation
@@ -153,3 +153,9 @@ function grid_points(imap::AbstractIndexMap, d::Int)
   L = length(dimension_inds(imap, d))
   return grid_points(imap, base^L, d)
 end
+
+function rand_p(rng::AbstractRNG, imap::AbstractIndexMap)
+  return [rand_p(rng, imap, i) for i in 1:dimension(imap)]
+end
+
+rand_p(imap::AbstractIndexMap, args...) = rand_p(default_rng(), imap, args...)
diff --git a/src/IndexMaps/complexindexmap.jl b/src/IndexMaps/complexindexmap.jl
index eda4bc2..c631265 100644
--- a/src/IndexMaps/complexindexmap.jl
+++ b/src/IndexMaps/complexindexmap.jl
@@ -2,6 +2,7 @@ using Base: Base
 using Dictionaries: Dictionaries, Dictionary, set!
 using ITensors: ITensors, Index, dim, hastags
 using ITensorNetworks: IndsNetwork, vertex_data
+using Random: AbstractRNG
 
 struct ComplexIndexMap{VB,VD,VR} <: AbstractIndexMap{VB,VD}
   index_digit::VB
@@ -135,10 +136,30 @@ function grid_points(imap::ComplexIndexMap, N::Int, d::Int)
   dims = dim.(dimension_inds(imap, d))
   @assert all(y -> y == first(dims), dims)
   base = float(first(dims))
-  Lre, Lim = length(real_indices(imap, d)), length(imag_indices(imap, d))
+  Lre, Lim = length(real_indices(imap, d)), length(imaginary_indices(imap, d))
   are, aim = round(base^Lre / N), round(base^Lim / N)
   grid_points = [
     i * (are / base^Lre) + im * j * (aim / base^Lim) for i in 0:(N + 1) for j in 0:(N + 1)
   ]
   return filter(x -> real(x) < 1 && imag(x) < 1, grid_points)
 end
+
+""" 
+  Picks a random grid point from `imap` given a dimension
+"""
+function rand_p(rng::AbstractRNG, imap::ComplexIndexMap, d::Integer)
+  dims = dim.(dimension_inds(imap, d))
+  @assert all(y -> y == first(dims), dims)
+  base = float(first(dims))
+  Lre, Lim = length(real_indices(imap, d)), length(imaginary_indices(imap, d))
+
+  # generate a random bitstring of length Lre, convert to decimal
+  bitstring_real = rand(rng, [j for j in 0:(base - 1)], Lre)
+  real_part = sum(((i, b),) -> b * ((base * 1.0)^-i), enumerate(bitstring_real))
+
+  # generate a random bitstring of length Lim, convert to decimal
+  bitstring_imag = rand(rng, [j for j in 0:(base - 1)], Lim)
+  imag_part = sum(((i, b),) -> b * ((base * 1.0)^-i), enumerate(bitstring_imag))
+
+  return real_part + im * imag_part
+end
diff --git a/src/IndexMaps/realindexmap.jl b/src/IndexMaps/realindexmap.jl
index 51063e9..0cf1e19 100644
--- a/src/IndexMaps/realindexmap.jl
+++ b/src/IndexMaps/realindexmap.jl
@@ -2,6 +2,7 @@ using Base: Base
 using Dictionaries: Dictionaries, Dictionary, set!
 using ITensors: ITensors, Index, dim
 using ITensorNetworks: IndsNetwork, vertex_data
+using Random: AbstractRNG
 
 struct RealIndexMap{VB,VD} <: AbstractIndexMap{VB,VD}
   index_digit::VB
@@ -75,12 +76,102 @@ function calculate_ind_values(imap::RealIndexMap, xs::Vector, dims::Vector{Int})
   return ind_to_ind_value_map
 end
 
-function grid_points(imap::RealIndexMap, N::Int, d::Int)
+""" 
+    grid_points(imap, N, d, span; exact_grid=true, enforced=[])
+
+  Gives `N` grid points from a given dimension of `imap` within a specified range.
+
+  # Arguments
+  - `imap::RealIndexMap`: An IndexMap specifying the structure of the TN being used
+  - `N::Int`: The number of grid points requested.
+  - `d::Int` The index of the dimension of `imap` requested
+  - `span::AbstractVector{<:Number}`: A two element number vector [a,b] with 0≤a<b≤1. The right endpoint of this span is not included as a gridpoint in the output.
+  - `exact_grid::Bool`: Flag specifying whether the function should give exact grid points 
+  (note: using `exact_grid=true` may cause less than `N` grid points to be returned)
+  - `enforced::AbstractVector{<:Number}`: A list of points that we want to enforce to show up in the grid.
+"""
+function grid_points(
+  imap::RealIndexMap,
+  N::Int,
+  d::Int;
+  span::AbstractVector{<:Number}=[0, 1],
+  exact_grid::Bool=true,
+  enforced=[],
+)
+  if length(span) != 2 || span[1] >= span[2] || span[1] < 0 || span[2] > 1
+    throw(
+      "expected a two-element vector [a,b] with 0≤a<b≤1 as input span, instead found $span"
+    )
+  end
+  dims = dim.(dimension_inds(imap, d))
+  @assert all(y -> y == first(dims), dims)
+  base = first(dims)
+  L = length(dimension_inds(imap, d))
+
+  # generate grid_points within the span (exclusive of right endpoint)
+  grid_points = collect(
+    range(
+      Float64(span[1]),
+      Float64(span[2]) - (Float64(span[2]) - Float64(span[1])) / N;
+      length=N,
+    ),
+  )
+
+  # if exact_grid is true, round to exact gridpoints
+  if exact_grid
+    grid_points = round_to_nearest_exact_point.(grid_points, (L,), (base,))
+    grid_points = unique(grid_points)
+  end
+
+  # add enforced points
+  if !isempty(enforced)
+    grid_points = sort(unique(vcat(grid_points, enforced_points)))
+  end
+
+  return grid_points
+end
+
+function round_to_nearest_exact_point(point::Float64, L::Int, base::Int=2)
+  return round(point * Float64(base)^min(63, L)) / Float64(base)^min(63, L)
+end
+#TODO: avoid using 2.0^min(63,L) to prevent overflow
+
+function grid_points(imap::RealIndexMap, d::Int; kwargs...)
+  dims = dim.(dimension_inds(imap, d))
+  @assert all(y -> y == first(dims), dims)
+  base = dims[d]
+  L = length(dimension_inds(imap, d))
+  return grid_points(imap, base^L, d; kwargs...)
+end
+
+grid_points(imap::RealIndexMap; kwargs...) = grid_points(imap, 1; kwargs...)
+
+#multi-dimensional grid_points
+function grid_points(imap::RealIndexMap, Ns::Vector{Int}, dims::Vector{Int}; kwargs...)
+  if length(Ns) != length(dims)
+    throw("length of Ns and dims do not match!")
+  end
+  coords = [grid_points(imap, pair[1], pair[2]; kwargs...) for pair in zip(Ns, dims)]
+  gp = Base.Iterators.product(coords...)
+  return collect.(gp)
+end
+
+function grid_points(imap::RealIndexMap, dims::Vector{Int}; kwargs...)
+  coords = [grid_points(imap, d; kwargs...) for d in dims]
+  gp = Base.Iterators.product(coords...)
+  return collect.(gp)
+end
+
+""" 
+  Picks a random grid point from `imap` given a dimension
+"""
+function rand_p(rng::AbstractRNG, imap::RealIndexMap, d::Integer)
   dims = dim.(dimension_inds(imap, d))
   @assert all(y -> y == first(dims), dims)
-  base = float(first(dims))
+  base = dims[d]
   L = length(dimension_inds(imap, d))
-  a = round(base^L / N)
-  grid_points = [i * (a / base^L) for i in 0:(N + 1)]
-  return filter(x -> x < 1, grid_points)
+  # generate a random bitstring of length L, converts to decimal
+  bitstring = rand(rng, [j for j in 0:(base - 1)], L)
+  decimal = sum(((i, b),) -> b * ((base * 1.0)^-i), enumerate(bitstring))
+  return decimal
 end
diff --git a/src/indsnetworkmap.jl b/src/indsnetworkmap.jl
index ac61c0f..f222aba 100644
--- a/src/indsnetworkmap.jl
+++ b/src/indsnetworkmap.jl
@@ -5,6 +5,7 @@ using NamedGraphs.GraphsExtensions: rem_vertex
 using ITensors: ITensors
 using ITensorNetworks:
   ITensorNetworks, AbstractIndsNetwork, IndsNetwork, data_graph, underlying_graph
+using Random: AbstractRNG
 
 struct IndsNetworkMap{V,I,IN<:IndsNetwork{V,I},IM<:AbstractIndexMap} <:
        AbstractIndsNetwork{V,I}
@@ -95,6 +96,7 @@ for f in [
   :calculate_ind_values,
   :calculate_p,
   :grid_points,
+  :rand_p,
   :index_value_to_scalar,
   :index_values_to_scalars,
 ]
@@ -104,6 +106,14 @@ for f in [
     end
   end
 end
+# Forward RNG functionality
+for f in [:rand_p]
+  @eval begin
+    function $f(rng::AbstractRNG, inm::IndsNetworkMap, args...; kwargs...)
+      return $f(rng, indexmap(inm), args...; kwargs...)
+    end
+  end
+end
 
 #Functions on indsnetwork that don't seem to autoforward
 function ITensors.inds(inm::IndsNetworkMap, args...; kwargs...)
diff --git a/test/test_indexmaps.jl b/test/test_indexmaps.jl
index 4344833..3c17745 100644
--- a/test/test_indexmaps.jl
+++ b/test/test_indexmaps.jl
@@ -2,7 +2,7 @@ using Test
 using ITensorNumericalAnalysis
 
 using NamedGraphs: vertices
-using NamedGraphs.NamedGraphGenerators: named_grid
+using NamedGraphs.NamedGraphGenerators: named_grid, named_comb_tree
 using NamedGraphs.GraphsExtensions: is_tree
 using ITensors: siteinds, inds
 using Dictionaries: Dictionary
@@ -70,3 +70,136 @@ Random.seed!(1234)
     xyzvals ≈ xyzvals_approx
   end
 end
+
+@testset "grid_points tests" begin
+  @testset "test grid_points irregular span" begin
+    L = 16
+    base = 2
+    g = named_comb_tree((2, L ÷ 2))
+    s = continuous_siteinds(g; map_dimension=2)
+
+    #first set -- irregular span
+    N = 125
+    a = 0.12
+    b = 0.95
+
+    # test that giving an irregular span still works and it obeys the constraints
+    test_gridpoints = grid_points(s, N, 1; span=[a, b])
+    @test length(test_gridpoints) == N
+    @test test_gridpoints[1] >= a
+    @test test_gridpoints[end] < b
+  end
+
+  @testset "test grid_points spacing" begin
+    L = 16
+    g = named_comb_tree((2, L ÷ 2))
+    s = continuous_siteinds(g; map_dimension=2)
+
+    #second set
+    N = 32
+    a = 0.25
+    b = 0.5
+    test2 = grid_points(s, N, 1; span=[a, b])
+    @test length(test2) == N
+    @test test2[1] >= a
+    @test test2[end] < b
+
+    # test that the spacing between all points is the same
+    spacing_equal = true
+    gap = (test2[2] - test2[1])
+    for i in 2:(N - 1)
+      if !(test2[i + 1] - test2[i] ≈ gap)
+        spacing_equal = false
+      end
+    end
+    @test spacing_equal
+  end
+
+  @testset "test grid_points exact_grid" begin
+    L = 16
+    g = named_comb_tree((2, L ÷ 2))
+    s = continuous_siteinds(g; map_dimension=2)
+    a = 0.25
+    b = 0.6
+
+    test3 = grid_points(s, 1; span=[a, b], exact_grid=true)
+
+    # tests that the number of points generated matches the number of exact gridpoints that fall in the range [a,b).
+    @test length(test3) == round(b * 2^(L ÷ 2) - 0.5) - round(a * 2^(L ÷ 2)) + 1
+    @test test3[1] >= a
+    @test test3[end] < b
+  end
+
+  @testset "test grid_points high precision" begin
+    #fourth set -- very large L
+    L = 140
+    base = 2
+    g = named_comb_tree((2, L ÷ 2))
+    s = continuous_siteinds(g; map_dimension=2)
+    n_grid = 16
+    @test length(grid_points(s, n_grid, 1)) == n_grid
+  end
+end
+
+@testset "test rand_p" begin
+  @testset "test rand_p for realindexmaps" begin
+    #test the rand_p() function to see if it succeeds on large L
+    L = 140
+    g = named_comb_tree((2, L ÷ 2))
+    s = continuous_siteinds(g; map_dimension=2)
+    ψ = cos_itn(s; dim=1) * cos_itn(s; dim=2)
+
+    #test the use of seedeed randomness
+    rng = Random.Xoshiro(42)
+    rand_gridpoint = rand_p(rng, s)
+    x1 = real(ITensorNumericalAnalysis.evaluate(ψ, rand_gridpoint))
+    @test x1 >= -1 && x1 <= 1 # check to make sure ψ can be evaluated at these points
+
+    rand_gridpoint1 = rand_p(rng, s, 1)
+    rand_gridpoint2 = rand_p(rng, s, 2)
+    x2 = real(ITensorNumericalAnalysis.evaluate(ψ, [rand_gridpoint1, rand_gridpoint2]))
+    @test x2 >= -1 && x2 <= 1 # check to make sure ψ can be evaluated at these points
+
+    #test the use of default rng
+    default_rng_gridpoint1 = rand_p(s)
+    y = real(ITensorNumericalAnalysis.evaluate(ψ, default_rng_gridpoint1))
+    @test y >= -1 && y <= 1 # check to make sure ψ can be evaluated at these points
+
+    # same things but with smaller L as well
+    L = 12
+    g = named_comb_tree((2, L ÷ 2))
+    s = continuous_siteinds(g; map_dimension=2)
+    ψ = cos_itn(s; dim=1) * cos_itn(s; dim=2)
+
+    #test the use of seedeed randomness
+    rng = Random.Xoshiro(42)
+    rand_gridpoint = rand_p(rng, s)
+    x1 = real(ITensorNumericalAnalysis.evaluate(ψ, rand_gridpoint))
+    @test x1 >= -1 && x1 <= 1 # check to make sure ψ can be evaluated at these points
+
+    rand_gridpoint1 = rand_p(rng, s, 1)
+    rand_gridpoint2 = rand_p(rng, s, 2)
+    x2 = real(ITensorNumericalAnalysis.evaluate(ψ, [rand_gridpoint1, rand_gridpoint2]))
+    @test x2 >= -1 && x2 <= 1 # check to make sure ψ can be evaluated at these points
+
+    #test the use of default rng
+    default_rng_gridpoint1 = rand_p(s)
+    y = real(ITensorNumericalAnalysis.evaluate(ψ, default_rng_gridpoint1))
+    @test y >= -1 && y <= 1 # check to make sure ψ can be evaluated at these points
+  end
+
+  @testset "test rand_p for complexindexmaps" begin
+    L = 16
+    g = named_comb_tree((2, L ÷ 2))
+    vs = collect(vertices(g))
+    split = div(length(vs), 2)
+    real_vs = [vs[1:split]]
+    imag_vs = [vs[(split + 1):end]]
+    s = complex_continuous_siteinds(g, real_vs, imag_vs)
+    ψ = poly_itn(s, [0, 0, 1])
+
+    rand_gridpoint = rand_p(s)
+    x = ITensorNumericalAnalysis.evaluate(ψ, rand_gridpoint)
+    @test real(x) >= -1 && real(x) <= 1 && imag(x) >= -1 && imag(x) <= 1
+  end
+end