diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 669750353..e8550798d 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -12,7 +12,7 @@ concurrency:
   cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: Julia ${{ matrix.version }} - ${{ matrix.arch }} - ${{ github.event_name }} - ExactPredicates ${{ matrix.preferences }} - ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
     env:
       JULIA_NUM_THREADS: ${{ matrix.julia-threads }}
@@ -21,11 +21,15 @@ jobs:
       matrix:
         version:
           - '1'
-          - '1.6'
+          - 'lts'
         os:
           - ubuntu-latest
         arch:
           - x64
+        preferences:
+          - 'default'
+          - 'true'
+          - 'false'
     steps:
       - name: Checkout repository for access
         uses: actions/checkout@v4
@@ -35,6 +39,7 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
+      - run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)'
 
       - name: Cache artifacts, etc.
         uses: julia-actions/cache@v2
@@ -44,6 +49,8 @@ jobs:
 
       - name: Run the tests
         uses: julia-actions/julia-runtest@v1
+        env:
+          USE_EXACTPREDICATES: ${{ matrix.preferences }}
 
       - name: Process the coverage
         if: always()
diff --git a/.gitignore b/.gitignore
index c34786994..f6b55de47 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,4 +6,5 @@ sw
 test/Manifest.toml
 docs/Manifest.toml 
 /docs/build/
-/.vscode
\ No newline at end of file
+/.vscode
+LocalPreferences.toml
\ No newline at end of file
diff --git a/NEWS.md b/NEWS.md
index 71957459b..8e26ee2ff 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,16 @@
 # Changelog
 
+## v.1.1.0
+
+- Added the option to disable ExactPredicates.jl using Preferences.jl. See [#131](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/131).
+- Added `DelauanyTriangulation.validate_triangulation` for validating triangulations. See [#131](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/131).
+- Fixed a bug with the currently unused `orient(p, q, r, s)` predicate. See [#131](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/131).
+- Added private functions `getz`, `_getz`, `getxyz`, and `_getxyz`. See [#131](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/131).
+
+## v1.0.5 
+
+- Disabled `deepcopy` on `PolygonTree`s and made it a no-op. See [#129](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/issues/129).
+
 ## v1.0.4
 
 Nothing breaking. Main changes:
diff --git a/Project.toml b/Project.toml
index 4cb4537bb..1ab26c5fc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,19 +1,21 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.0.5"
+version = "1.1.0-DEV"
 
 [deps]
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
+Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 EnumX = "1.0"
 ExactPredicates = "2.2"
-julia = "1"
 Random = "1"
 Test = "1"
+julia = "1"
+Preferences = "1.4"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/README.md b/README.md
index a51553662..553338b4d 100644
--- a/README.md
+++ b/README.md
@@ -113,8 +113,9 @@ When contributing in the form of a pull request, there are a few important featu
 2. **Document any new functions**. If your contribution involves any new functions, make sure they are well documented even if they are only for internal use.
 3. **You need to include appropriate tests for your contribution**. For small changes, simple tests are fine, but for larger changes that implement or change an algorithm, understanding what tests you need to write is crucial when you are implementing a new feature. For example, if you wanted to implement a new algorithm for inserting a curve into a triangulation, you need to make sure that you test things like (1) point sets in general position, (2) point sets with many cocircular points, (3) collinear edges, and so on. In particular, degeneracies are some of the most important things to test for. This is not simple to do if you are inexperienced, so do feel free to ask for guidance.
 4. **Your contribution should only do one thing**. To make sure that it is easy to track your changes, and to make it easier to track regressions across versions, please try to only make your contribution do one thing. For example, do not move around a bunch of files while at the same time implementing a function somewhere else - make two pull requests.
+5. **Add to NEWS.md**. Describe your change in the NEWS.md file.
 
-If you are just contributing something minimal, for example a typo, even a blank PR is fine so long as it is obvious what you have done.
+If you are just contributing something minimal, for example a typo, even a blank PR description is fine so long as it is obvious what you have done.
 
 ## Similar Packages
 
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index ac40fa4cd..973af2626 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "1618ed588410cd4dc664fc80ecbe5eaa4e31c80c"
+project_hash = "a5f97bd51af051946996607d64905eb7a7b0c1f9"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
@@ -58,9 +58,9 @@ version = "1.1.1"
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "ed2ec3c9b483842ae59cd273834e5b46206d6dda"
+git-tree-sha1 = "5c9b74c973181571deb6442d41e5c902e6b9f38e"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.11.0"
+version = "7.12.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
@@ -89,9 +89,9 @@ uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
 [[deps.Automa]]
 deps = ["PrecompileTools", "TranscodingStreams"]
-git-tree-sha1 = "588e0d680ad1d7201d4c6a804dcb1cd9cba79fbb"
+git-tree-sha1 = "014bc22d6c400a7703c0f5dc1fdc302440cf88be"
 uuid = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b"
-version = "1.0.3"
+version = "1.0.4"
 
 [[deps.AxisAlgorithms]]
 deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
@@ -146,10 +146,10 @@ uuid = "159f3aea-2a34-519c-b102-8c37f9878175"
 version = "1.0.5"
 
 [[deps.CairoMakie]]
-deps = ["CRC32c", "Cairo", "Colors", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "PrecompileTools"]
-git-tree-sha1 = "d69c7593fe9d7d617973adcbe4762028c6899b2c"
+deps = ["CRC32c", "Cairo", "Cairo_jll", "Colors", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "PrecompileTools"]
+git-tree-sha1 = "f84837ccd1411ba059bb0b752dab9c7f1b0b0826"
 uuid = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
-version = "0.11.11"
+version = "0.12.4"
 
 [[deps.Cairo_jll]]
 deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
@@ -175,9 +175,9 @@ weakdeps = ["SparseArrays"]
 
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
-git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73"
+git-tree-sha1 = "b8fe8546d52ca154ac556809e10c75e6e7430ac8"
 uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
-version = "0.7.4"
+version = "0.7.5"
 
 [[deps.ColorBrewer]]
 deps = ["Colors", "JSON", "Test"]
@@ -241,9 +241,9 @@ version = "1.1.1+0"
 
 [[deps.ConcurrentUtilities]]
 deps = ["Serialization", "Sockets"]
-git-tree-sha1 = "6cbbd4d241d7e6579ab354737f4dd95ca43946e1"
+git-tree-sha1 = "ea32b83ca4fefa1768dc84e504cc0a94fb1ab8d1"
 uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb"
-version = "2.4.1"
+version = "2.4.2"
 
 [[deps.ConstructionBase]]
 deps = ["LinearAlgebra"]
@@ -292,10 +292,10 @@ uuid = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
 version = "1.2.0"
 
 [[deps.DelaunayTriangulation]]
-deps = ["EnumX", "ExactPredicates", "Random"]
+deps = ["EnumX", "ExactPredicates", "Preferences", "Random"]
 path = ".."
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
-version = "1.0.4"
+version = "1.0.6"
 
 [[deps.DiffResults]]
 deps = ["StaticArraysCore"]
@@ -348,9 +348,9 @@ version = "0.9.3"
 
 [[deps.Documenter]]
 deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"]
-git-tree-sha1 = "5461b2a67beb9089980e2f8f25145186b6d34f91"
+git-tree-sha1 = "76deb8c15f37a3853f13ea2226b8f2577652de05"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "1.4.1"
+version = "1.5.0"
 
 [[deps.DocumenterTools]]
 deps = ["AbstractTrees", "Base64", "DocStringExtensions", "Documenter", "FileWatching", "Gumbo", "LibGit2", "OpenSSH_jll", "Sass"]
@@ -581,9 +581,9 @@ version = "1.3.14+0"
 
 [[deps.GridLayoutBase]]
 deps = ["GeometryBasics", "InteractiveUtils", "Observables"]
-git-tree-sha1 = "6f93a83ca11346771a93bbde2bdad2f65b61498f"
+git-tree-sha1 = "fc713f007cff99ff9e50accba6373624ddd33588"
 uuid = "3955a311-db13-416c-9275-1d80ed98e5e9"
-version = "0.10.2"
+version = "0.11.0"
 
 [[deps.Grisu]]
 git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
@@ -692,13 +692,11 @@ deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArr
 git-tree-sha1 = "88a101217d7cb38a7b481ccd50d21876e1d1b0e0"
 uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 version = "0.15.1"
+weakdeps = ["Unitful"]
 
     [deps.Interpolations.extensions]
     InterpolationsUnitfulExt = "Unitful"
 
-    [deps.Interpolations.weakdeps]
-    Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
-
 [[deps.IntervalArithmetic]]
 deps = ["CRlibm_jll", "MacroTools", "RoundingEmulator"]
 git-tree-sha1 = "433b0bb201cd76cb087b017e49244f10394ebe9c"
@@ -950,16 +948,16 @@ uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 version = "0.5.13"
 
 [[deps.Makie]]
-deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun"]
-git-tree-sha1 = "4d49c9ee830eec99d3e8de2425ff433ece7cc1bc"
+deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Dates", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun", "Unitful"]
+git-tree-sha1 = "57a1a2b3d12e04f9e9fb77d61cd12571d5541c5f"
 uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
-version = "0.20.10"
+version = "0.21.4"
 
 [[deps.MakieCore]]
-deps = ["Observables", "REPL"]
-git-tree-sha1 = "248b7a4be0f92b497f7a331aed02c1e9a878f46b"
+deps = ["ColorTypes", "GeometryBasics", "IntervalSets", "Observables"]
+git-tree-sha1 = "638bc817096742e8302f7b0b972ee5701fe00e97"
 uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
-version = "0.7.3"
+version = "0.8.3"
 
 [[deps.MappedArrays]]
 git-tree-sha1 = "2dab0221fe2b0f2cb6754eaa743cc266339f527e"
@@ -978,9 +976,9 @@ version = "0.1.2"
 
 [[deps.MathTeXEngine]]
 deps = ["AbstractTrees", "Automa", "DataStructures", "FreeTypeAbstraction", "GeometryBasics", "LaTeXStrings", "REPL", "RelocatableFolders", "UnicodeFun"]
-git-tree-sha1 = "96ca8a313eb6437db5ffe946c457a401bbb8ce1d"
+git-tree-sha1 = "1865d0b8a2d91477c8b16b49152a32764c7b1f5f"
 uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
-version = "0.5.7"
+version = "0.6.0"
 
 [[deps.MbedTLS]]
 deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"]
@@ -1050,9 +1048,9 @@ uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
 version = "0.5.5"
 
 [[deps.OffsetArrays]]
-git-tree-sha1 = "e64b4f5ea6b7389f6f046d13d4896a8f9c1ba71e"
+git-tree-sha1 = "1a27764e945a152f7ca7efa04de513d473e9542e"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.14.0"
+version = "1.14.1"
 weakdeps = ["Adapt"]
 
     [deps.OffsetArrays.extensions]
@@ -1500,9 +1498,9 @@ version = "0.1.1"
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "20833c5b7f7edf0e5026f23db7f268e4f23ec577"
+git-tree-sha1 = "eeafab08ae20c62c44c8399ccb9354a04b80db50"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.6"
+version = "1.9.7"
 weakdeps = ["ChainRulesCore", "Statistics"]
 
     [deps.StaticArrays.extensions]
@@ -1611,9 +1609,9 @@ uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
 version = "0.10.0"
 
 [[deps.TranscodingStreams]]
-git-tree-sha1 = "d73336d81cafdc277ff45558bb7eaa2b04a8e472"
+git-tree-sha1 = "60df3f8126263c0d6b357b9a1017bb94f53e3582"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
-version = "0.10.10"
+version = "0.11.0"
 weakdeps = ["Random", "Test"]
 
     [deps.TranscodingStreams.extensions]
@@ -1647,6 +1645,20 @@ git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
 uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
 version = "0.4.1"
 
+[[deps.Unitful]]
+deps = ["Dates", "LinearAlgebra", "Random"]
+git-tree-sha1 = "dd260903fdabea27d9b6021689b3cd5401a57748"
+uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
+version = "1.20.0"
+
+    [deps.Unitful.extensions]
+    ConstructionBaseUnitfulExt = "ConstructionBase"
+    InverseFunctionsUnitfulExt = "InverseFunctions"
+
+    [deps.Unitful.weakdeps]
+    ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+    InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
+
 [[deps.WoodburyMatrices]]
 deps = ["LinearAlgebra", "SparseArrays"]
 git-tree-sha1 = "c1a7aa6219628fcd757dede0ca95e245c5cd9511"
@@ -1760,10 +1772,10 @@ uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
 version = "1.6.43+1"
 
 [[deps.libsass_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "941afb93587dcec07f89e511057f5efc0bec6f0d"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "6044ffe7e7bf0602e2039dc747c3332a097ac74b"
 uuid = "47bcb7c8-5119-555a-9eeb-0afcc36cd728"
-version = "3.6.4+0"
+version = "3.6.6+0"
 
 [[deps.libsixel_jll]]
 deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "libpng_jll"]
diff --git a/docs/Project.toml b/docs/Project.toml
index bd7a5bf9b..9809d0605 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -10,6 +10,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
+Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 SimpleGraphs = "55797a34-41de-5266-9ec1-32ac4eb504d3"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
diff --git a/docs/make.jl b/docs/make.jl
index e18dd1787..8018f5bc5 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -146,6 +146,7 @@ const _PAGES = [
         "Geometrical Predicates" => "manual/predicates.md",
         "Triangulation Output" => "manual/triangulation_output.md",
         "Voronoi Tessellation Output" => "manual/voronoi_output.md",
+        "Disabling Exact Predicates" => "manual/disabling_ea.md",
     ],
     "API Reference" => [
         "Section Overview" => "api/overview.md",
diff --git a/docs/src/api/triangulation.md b/docs/src/api/triangulation.md
index a49e35873..2896d9f61 100644
--- a/docs/src/api/triangulation.md
+++ b/docs/src/api/triangulation.md
@@ -59,4 +59,5 @@ get_insertion_order
 get_vertices(::Triangulation)
 get_edges(::Triangulation)
 num_neighbours
+validate_triangulation
 ```
\ No newline at end of file
diff --git a/docs/src/extended/utils.md b/docs/src/extended/utils.md
index 3e00efb44..45ee8e486 100644
--- a/docs/src/extended/utils.md
+++ b/docs/src/extended/utils.md
@@ -28,4 +28,10 @@ Pages = ["src/predicates/boundaries_and_ghosts.jl"]
 
 ```@docs
 convert_certificate
+DefaultAdjacentValue
+𝒢
+GhostVertex
+ε
+∅
+USE_EXACTPREDICATES
 ```
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
index 63db4fd18..4f621609d 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -18,7 +18,7 @@ This is a package for computing Delaunay triangulations and Voronoi tessellation
 - Computation of the pole of inaccessibility.
 - The interface for defining geometric primitives is fully customisable.
 
-To ensure that the algorithms are robust, we use [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl) to define all geometric predicates in this package. Much of the work in this package is derived from the book [*Delaunay Mesh Generation* by Cheng, Dey, and Shewchuk (2013)](https://people.eecs.berkeley.edu/~jrs/meshbook.html).
+To ensure that the algorithms are robust, we use [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl) to define all geometric predicates in this package. (ExactPredicates.jl can be disabled, as described [here](manual/disabling_ea.md).) Much of the work in this package is derived from the book [*Delaunay Mesh Generation* by Cheng, Dey, and Shewchuk (2013)](https://people.eecs.berkeley.edu/~jrs/meshbook.html).
 
 # Documentation Structure 
 
@@ -71,6 +71,7 @@ When contributing in the form of a pull request, there are a few important featu
 2. **Document any new functions**. If your contribution involves any new functions, make sure they are well documented even if they are only for internal use.
 3. **You need to include appropriate tests for your contribution**. For small changes, simple tests are fine, but for larger changes that implement or change an algorithm, understanding what tests you need to write is crucial when you are implementing a new feature. For example, if you wanted to implement a new algorithm for inserting a curve into a triangulation, you need to make sure that you test things like (1) point sets in general position, (2) point sets with many cocircular points, (3) collinear edges, and so on. In particular, degeneracies are some of the most important things to test for. This is not simple to do if you are inexperienced, so do feel free to ask for guidance.
 4. **Your contribution should only do one thing**. To make sure that it is easy to track your changes, and to make it easier to track regressions across versions, please try to only make your contribution do one thing. For example, do not move around a bunch of files while at the same time implementing a function somewhere else - make two pull requests.
+5. **Add to NEWS.md**. Describe your change in the NEWS.md file.
 
 If you are just contributing something minimal, for example a typo, even a blank PR is fine so long as it is obvious what you have done.
 
diff --git a/docs/src/literate_tutorials/convex.jl b/docs/src/literate_tutorials/convex.jl
index 9f59b36a9..edffffb84 100644
--- a/docs/src/literate_tutorials/convex.jl
+++ b/docs/src/literate_tutorials/convex.jl
@@ -8,6 +8,7 @@ using DelaunayTriangulation
 using CairoMakie
 using ReferenceTests #src
 using Test #src
+using Preferences #src
 fig_path = joinpath(@__DIR__, "../figures") #src
 
 points = [
@@ -21,7 +22,7 @@ tri = triangulate_convex(points, 1:7)
 #-
 fig, ax, sc = triplot(tri)
 fig 
-@test_reference joinpath(fig_path, "convex_ex_1.png") fig #src
+load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "convex_ex_1.png") fig #src
 
 # This `tri` is our triangulation of the convex polygon. 
 # The first input is the set of points, and `S` defines 
@@ -44,7 +45,7 @@ tri = triangulate_convex(points, S)
 #-
 fig, ax, sc = triplot(tri)
 fig 
-@test_reference joinpath(fig_path, "convex_ex_2.png") fig #src
+load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "convex_ex_2.png") fig #src
 
 # Here is a comparison of the time it takes to triangulate this 
 # using `triangulate_convex` or `triangulate`.
diff --git a/docs/src/literate_tutorials/point_location.jl b/docs/src/literate_tutorials/point_location.jl
index 2bb0ecc51..89b0caca1 100644
--- a/docs/src/literate_tutorials/point_location.jl
+++ b/docs/src/literate_tutorials/point_location.jl
@@ -21,6 +21,7 @@ using StableRNGs
 using Test #src
 using ReferenceTests #src
 fig_path = joinpath(@__DIR__, "../figures") #src
+using Preferences #src
 
 points = [
     (-3.0, 6.0), (5.0, 1.0), (-5.0, 3.0), (2.0, -3.0),
@@ -78,7 +79,7 @@ V = jump_and_march(tri, q)
 fig, ax, sc = triplot(tri, show_ghost_edges=true)
 scatter!(ax, q)
 fig
-@test_reference joinpath(fig_path, "point_location_ex_1.png") fig #src
+load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "point_location_ex_1.png") fig #src
 
 # ## Region with concave boundaries and holes 
 # Now we give an example of point location for a reason with holes. Since the 
@@ -104,7 +105,7 @@ refine!(tri; max_area=0.01get_area(tri), rng);
 # To demonstrate this, see the following plot:
 fig, ax, sc = triplot(tri, show_ghost_edges=true)
 fig
-@test_reference joinpath(fig_path, "point_location_ex_2.png") fig #src
+load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "point_location_ex_2.png") fig #src
 
 # The ghost edges now intersect the boundary, which doesn't make sense, and creates difficulties.
 # Let us now demonstrate how the function still works here. We try finding the blue points shown below.
@@ -117,7 +118,7 @@ qs = [
 fig, ax, sc = triplot(tri, show_ghost_edges=false)
 scatter!(ax, qs, color=:blue, markersize=16)
 fig
-@test_reference joinpath(fig_path, "point_location_ex_3.png") fig #src
+load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "point_location_ex_3.png") fig #src
 
 # Now let's find the triangles.
 Vs = [jump_and_march(tri, q; rng) for q in qs]
@@ -181,7 +182,7 @@ qs = [
 fig, ax, sc = triplot(tri)
 scatter!(ax, qs, color=:blue, markersize=16)
 fig
-@test_reference joinpath(fig_path, "point_location_ex_4.png") fig by=psnr_equality(10) #src
+load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "point_location_ex_4.png") fig by=psnr_equality(10) #src
 
 # Here are the `jump_and_march` results. 
 Vs = [jump_and_march(tri, q; rng, concavity_protection=true) for q in qs]
diff --git a/docs/src/manual/disabling_ea.md b/docs/src/manual/disabling_ea.md
new file mode 100644
index 000000000..3098b4d7c
--- /dev/null
+++ b/docs/src/manual/disabling_ea.md
@@ -0,0 +1,100 @@
+```@meta 
+CurrentModule = DelaunayTriangulation
+```
+
+# Disabling Exact Predicates
+
+For performance reasons, you may find it useful to want to disable exact predicates using [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl). This can be easily done using a setup with [Preferences.jl](https://github.com/JuliaPackaging/Preferences.jl), but before you consider disabling exact predicates.jl, there are a few things to be aware of. If you just want to disable them without reading a lot of information warning you about the consequences, please skip to the end.
+
+## Why use exact predicates?
+
+Three great resources for understanding why we need exact predicates are
+
+1. Jonathan Shewchuk's paper on adaptive precision floating-point arithmetic [here](https://doi.org/10.1007/PL00009321).
+2. Jonathan Shewchuk's lecture notes on geometric robustness [here](https://people.eecs.berkeley.edu/~jrs/meshpapers/robnotes.pdf).
+3. [This paper](https://doi.org/10.1016/j.comgeo.2007.06.003) by Kettner et al. (2008) on some examples of issues with inexact arithmetic.
+
+We give a simple summary here. A big component of the algorithms used in this package are what are known as _geometric predicates_, some of these being:
+
+- `orient(p, q, r)`: Is `r` left, right, or on the line through `pq`?
+- `incircle(p, q, r, s)`: Is `s` inside, outside, or on the circle through `p`, `q`, and `r`?
+
+These predicates can be computed using determinants:
+
+```math 
+\begin{align*}
+O_{pqr} &:= \textrm{orient}(p, q, r) = \begin{vmatrix} p_x - r_x & p_y - r_y \\ q_x - r_x & q_y - r_y \end{vmatrix}, \\
+C_{pqrs} &:= \textrm{incircle}(p, q, r, s) = \begin{vmatrix} p_x - s_x & p_y - s_y & (p_x - s_x)^2 + (p_y - s_y)^2 \\ q_x - s_x & q_y - s_y & (q_x - s_x)^2 + (q_y - s_y)^2 \\ r_x - s_x & r_y - s_y & (r_x - s_x)^2 + (r_y - s_y)^2 \end{vmatrix}.
+\end{align*}
+```
+
+The signs of these determinants $O_{pqr}$ and $C_{pqrs}$ are used to determinant the answers to the above questions. In inexact arithmetic, it is common that the sign picked is wrong when the determinants are close to zero. The consequences of this can be catastrophic:
+
+1. The algorithms may hang or crash.
+2. The final triangulation may be completely invalid. For example, if a point is being added into a triangulation right onto an existing edge, then in exact arithmetic we would know to split the edge to the left and to the right. In inexact arithmetic, the point may be to the left of the edge but detected as being to the right of it, thus adding a triangle that crosses an edge.
+3. You may encounter `BoundsError`s from bad `Adjacent` queries where a triangle is expected to exist but doesn't.
+
+Another issue is due to the fact that floating point arithmetic is not associative. In exact arithmetic, we would expect for example that 
+
+```math 
+O_{pqr} = O_{qrp} = O_{rpq},
+```
+
+but this is not true in floating point arithmetic. This causes issues with consistency - a point may be found to be both left and right of an edge depending on the order of the points given to the `orient` predicate, inevitably leading to an invalid triangulation. With the use of exact predicates, this property is guaranteed to hold, ensuring that all the predicate results are consistent with each other. This has the following consequence: **Even if you think exact predicates are not necessary for you because none of your inputs are exact (for example), you still want them to guarantee consistency with predicates regardless of the input order**.
+
+## Will disabling exact predicates give me better performance?
+
+It is also not even the case that using inexact predicates will give you better performance than if you were to us exact predicates. ExactPredicates.jl uses clever filters to ensure that the fast path, meaning an inexact computation, is performed first. This means that, in most cases, computing the `orient` predicate would have the same speed as if you were to compute it using the determinant definition. 
+
+One case where you could see improvements would be if there were many collinear points, where the slow computations from exact predicates would actually be needed - but this is exactly the case where you expect inexact predicates to cause you problems! 
+
+You should always benchmark your problems to see if disabling exact predicates, if you choose to do, will actually give you better performance.
+
+## How do I disable exact predicates?
+
+If you still want to disable exact predicates, here is how you can do so. Before doing `using DelaunayTriangulation`, you can use 
+
+```julia
+using Preferences: set_preferences! 
+set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => false)
+using DelaunayTriangulation # only load after setting the preference 
+```
+
+The `set_preferences!` call will make a `LocalPreferences.toml` file in your directory that sets this preference. Once this file exists you are free to delete `Preferences.jl`. If you want to skip using Preferences.jl entirely, you can also just create `LocalPreferences.toml` in your working directory manually and put 
+
+```
+[DelaunayTriangulation]
+USE_EXACTPREDICATES = false
+```
+
+into it. If you later want to re-enable exact predicates, either delete the file or write `USE_EXACTPREDICATES = true` instead. This setup ensures that there is no slowdown in the package from checking if ExactPredicates.jl is being used at runtime, as it is all done at compile time instead.
+
+## Can I check if my computed triangulation is valid?
+
+When you are not using exact predicates, you may want to check if your computed triangulation is actually a valid Delaunay triangulation. We provide the function `DelaunayTriangulation.validate_triangulation` for this purpose. This functionality is quite slow to use and is not currently optimised or well-documented (contributions towards addressing these issues are welcome), but it will work. One important note is that this check does actually use predicates in certain areas, so this check is still not guaranteed to be 100% accurate. Here is an example of its use.
+
+```julia
+using DelaunayTriangulation
+tri = triangulate(rand(2, 50))
+DelaunayTriangulation.validate_triangulation(tri)
+```
+```julia 
+true
+```
+```julia
+T = first(each_solid_triangle(tri)) 
+DelaunayTriangulation.delete_triangle!(tri, T) # break the triangulation for this example
+DelaunayTriangulation.validate_triangulation(tri)
+```
+```julia
+The edge (12, 40) does not have two incident triangles.
+The edge (12, 40) appears as an edge in the graph but it and its reverse are not both a key of the adjacent map.
+
+false
+```
+```julia
+DelaunayTriangulation.validate_triangulation(tri; print_result = false)
+```
+```julia 
+false
+```
diff --git a/docs/src/manual/predicates.md b/docs/src/manual/predicates.md
index 4d77271e8..5d4e862f6 100644
--- a/docs/src/manual/predicates.md
+++ b/docs/src/manual/predicates.md
@@ -4,7 +4,7 @@ CurrentModule = DelaunayTriangulation
 
 # Geometrical Predicates
 
-This section discusses how geometrical predicates are defined in this package. The predicates in this package are primarily derived from those implemented in [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl). The choice of exact predicates is important for the robustness of the algorithms in this package. Without using exact predicates, you may quickly find issues such as infinite loops or errors in the algorithms, as discussed for example at the start of p.3 of [these notes](https://perso.uclouvain.be/jean-francois.remacle/LMECA2170/robnotes.pdf) by Shewchuk.
+This section discusses how geometrical predicates are defined in this package. The predicates in this package are primarily derived from those implemented in [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl). The choice of exact predicates is important for the robustness of the algorithms in this package. Without using exact predicates, you may quickly find issues such as infinite loops or errors in the algorithms, as discussed for example at the start of p.3 of [these notes](https://perso.uclouvain.be/jean-francois.remacle/LMECA2170/robnotes.pdf) by Shewchuk. If you do want to disable exact predicates, see [here](disabling_ea.md).
 
 ## Certificates 
 
@@ -24,4 +24,4 @@ DelaunayTriangulation.is_positively_oriented(flag)
 
 ## Predicates
 
-The predicates we define, as mentioned, are primarily derived from those in [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl), where we simply extend the [`orient`](@ref orient_predicate), [`incircle`](@ref incircle_predicate), [`parallelorder`](@ref parallelorder_predicate), [`sameside`](@ref sameside_predicate), and [`meet`](@ref meet_predicate) predicates, allowing us to predicates for, for example, [`triangle_orientation`](@ref) and [`point_position_relative_to_line`](@ref). Predicates for working with the boundary and ghost vertices are also implemented, for example, [`is_ghost_edge`](@ref is_ghost_edge) and [`is_boundary_node`](@ref).
\ No newline at end of file
+The predicates we define, as mentioned, are primarily derived from those in [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl), where we simply extend the [`orient`](@ref orient_predicate), [`incircle`](@ref incircle_predicate), [`parallelorder`](@ref parallelorder_predicate), [`sameside`](@ref sameside_predicate), and [`meet`](@ref meet_predicate) predicates, allowing us to define predicates for, for example, [`triangle_orientation`](@ref) and [`point_position_relative_to_line`](@ref). Predicates for working with the boundary and ghost vertices are also implemented, for example, [`is_ghost_edge`](@ref is_ghost_edge) and [`is_boundary_node`](@ref).
\ No newline at end of file
diff --git a/docs/src/math/clipped.md b/docs/src/math/clipped.md
index 3108fe181..3af108345 100644
--- a/docs/src/math/clipped.md
+++ b/docs/src/math/clipped.md
@@ -318,7 +318,7 @@ new_points[u] = p
 new_points[v] = q
 new_vertices[u] = u
 new_vertices[v] = v
-poly!(ax, [new_points; new_points[1]], color = :darkgreen, linewidth = 3)
+poly!(ax, [new_points; new_points[1]], color = :darkgreen)
 lines!(ax, [0.0, 5.0, 5.0, 0.0, 0.0], [-15.0, -15.0, 15.0, 15.0, -15.0], color = :red)
 resize_to_layout!(fig)
 ```
diff --git a/src/DelaunayTriangulation.jl b/src/DelaunayTriangulation.jl
index 4e2f6d752..491ed8f06 100644
--- a/src/DelaunayTriangulation.jl
+++ b/src/DelaunayTriangulation.jl
@@ -1,9 +1,43 @@
 module DelaunayTriangulation
 
+"""
+    DefaultAdjacentValue = 0 
+
+Default value used for representing an empty result 
+from an adjacency query.
+"""
 const DefaultAdjacentValue = 0
+
+"""
+    ∅ = DefaultAdjacentValue
+
+Alias for [`DefaultAdjacentValue`](@ref).
+"""
 const ∅ = DefaultAdjacentValue
+
+"""
+    GhostVertex = -1
+
+Number used for representing initial ghost vertices. 
+All other ghost vertices are derived from subtracting from 
+this number. See https://juliageometry.github.io/DelaunayTriangulation.jl/stable/manual/ghost_triangles/.
+"""
 const GhostVertex = -1
+
+"""
+    𝒢 = GhostVertex
+
+Alias for [`𝒢`](@ref).
+"""
 const 𝒢 = GhostVertex
+
+"""
+    ε = sqrt(eps(Float64))
+
+Number used as a tolerance in certain functions, e.g. 
+for mesh refinement when using [`check_precision`](@ref) to 
+avoid degenerate circumcenters.
+"""
 const ε = sqrt(eps(Float64))
 
 const INF_WARN = Ref(true)
@@ -15,6 +49,42 @@ By default, this warning is enabled.
 """
 toggle_inf_warn!() = (INF_WARN[] = !INF_WARN[])
 
+@static if VERSION ≥ v"1.6"
+    using Preferences
+end
+
+@static if VERSION ≥ v"1.6"
+    const USE_EXACTPREDICATES = @load_preference("USE_EXACTPREDICATES", true)::Bool
+else 
+    const USE_EXACTPREDICATES = true 
+end
+@doc """
+    USE_EXACTPREDICATES
+
+Whether to use ExactPredicates.jl for computing predicates. By default, 
+this is true, but a user can change this by defining a preference with Preferences.jl, i.e. 
+you could do the following 
+
+```julia-repl
+julia> using Preferences: set_preferences!
+
+julia> set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => false)
+
+julia> using DelaunayTriangulation # load only after setting the preference
+
+julia> DelaunayTriangulation.USE_EXACTPREDICATES
+false
+```
+
+You have set `USE_EXACTPREDICATES = $USE_EXACTPREDICATES`. 
+
+!!! note "Precision"
+
+    Even if you have disabled ExactPredicates.jl, the predicates 
+    are still computed in Float64 precision.
+"""
+USE_EXACTPREDICATES
+
 using ExactPredicates
 using EnumX
 using Random
@@ -115,6 +185,8 @@ include("algorithms/voronoi/clipped.jl")
 include("algorithms/voronoi/main.jl")
 include("algorithms/voronoi/unbounded.jl")
 
+include("validation.jl")
+
 export
     each_triangle,
     each_solid_triangle,
diff --git a/src/README.md b/src/README.md
deleted file mode 100644
index 8b1378917..000000000
--- a/src/README.md
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/src/geometric_primitives/points.jl b/src/geometric_primitives/points.jl
index dbc541035..d9fd72237 100644
--- a/src/geometric_primitives/points.jl
+++ b/src/geometric_primitives/points.jl
@@ -32,6 +32,13 @@ julia> gety(p)
 """
 gety(p) = p[2]
 
+"""
+    getz(p) -> Number 
+
+Get the z-coordinate of `p`.
+"""
+getz(p) = p[3]
+
 """
     getxy(p) -> NTuple{2, Number}
 
@@ -49,6 +56,13 @@ julia> getxy(p)
 """
 getxy(p) = (getx(p), gety(p))
 
+"""
+    getxyz(p) -> NTuple{3, Number}
+
+Given a three-dimensional `p`, returns its coordinates as a `Tuple`.
+"""
+getxyz(p) = (getx(p), gety(p), getz(p))
+
 """
     _getx(p) -> Float64
 
@@ -93,10 +107,17 @@ julia> DelaunayTriangulation._gety(p)
 """
 _gety(p) = Float64(gety(p))
 
+"""
+    _getz(p) -> Float64 
+
+Get the z-coordinate of `p` as a `Float64.`
+"""
+_getz(p) = Float64(getz(p))
+
 """
     _getxy(p) -> NTuple{2, Float64}
 
-Get the coordinates of `p` as a `Tuple`.
+Get the coordinates of `p` as a `Tuple` of `Float64`s.
 
 # Examples
 ```jldoctest
@@ -115,6 +136,13 @@ julia> DelaunayTriangulation._getxy(p)
 """
 _getxy(p) = (_getx(p), _gety(p))
 
+"""
+    _getxyz(p) -> NTuple{3, Float64}
+
+Given a three-dimemsional `p`, returns its coordinates as a `Tuple` of `Float64`s.
+"""
+_getxyz(p) = (_getx(p), _gety(p), _getz(p))
+
 @doc """
     getpoint(points, vertex) -> NTuple{2, Number}
 
diff --git a/src/predicates/exactpredicates_definitions.jl b/src/predicates/exactpredicates_definitions.jl
index f699c77fa..ea2824dc7 100644
--- a/src/predicates/exactpredicates_definitions.jl
+++ b/src/predicates/exactpredicates_definitions.jl
@@ -1,3 +1,78 @@
+#=
+Much of the code in the file is derived or taken from ExactPredicates.jl (https://github.com/lairez/ExactPredicates.jl). 
+The license for ExactPredicates.jl is printed below, taken from https://github.com/lairez/ExactPredicates.jl/blob/master/LICENSE.
+
+Copyright (c) 2019 Pierre Lairez
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+=#
+
+"""
+    ext(ux, uy, vx, vy)
+
+Computes `ExactPredicates.ext((ux, uy), (vx, vy))`, i.e. 
+returns `ux * vy - uy * vx`.
+"""
+ext(ux, uy, vx, vy) = ux * vy - uy * vx
+
+"""
+    inp(ux, uy, vx, vy)
+
+Computes `ExactPredicates.inp((ux, uy), (vx, vy))`, i.e.
+returns `ux * vx + uy * vy`.
+"""
+inp(ux, uy, vx, vy) = ux * vx + uy * vy
+
+"""
+    det(a, b, c, d)
+
+Computes `ExactPredicates.det(a, b, c, d)`,
+i.e. returns `a*d - b*c`.
+"""
+det(a, b, c, d) = a * d - b * c
+
+"""
+    det(a, b, c, d, e, f, g, h, i)
+
+Computes `ExactPredicates.det(a, b, c, d, e, f, g, h, i)`,
+i.e. returns the determinant of 
+```math
+\\det \\begin{bmatrix} a & b & c \\\\ d & e & f \\\\ g & h & i \\end{bmatrix}
+```
+"""
+det(a, b, c, d, e, f, g, h, i) = a * det(e, f, h, i) - d * det(b, c, h, i) + g * det(b, c, e, f)
+
+"""
+    opposite_signs(x, y) -> Bool
+
+From ExactPredicates.jl, returns `true` if `x` and `y` have opposite signs, and `false` otherwise.
+Assumes that `x` and `y` are in `[-1, 0, 1]`.
+"""
+opposite_signs(x, y) = xor(x, y) == -2
+
+"""
+    sgn(x) -> Int 
+
+Returns `Int(sign(x))`.
+""" 
+sgn(x) = Int(sign(x))
+
 """
     orient_predicate(p, q, r) -> Integer 
 
@@ -7,6 +82,10 @@ Returns `ExactPredicates.orient(p, q, r)`, in particular we return:
 - `0`: `(p, q, r)` is collinear / degenerate.
 - `-1`: `(p, q, r)` is negatively oriented.
 
+If ExactPredicates.jl has been disabled using Preferences.jl, the 
+determinant defining `orient` is evaluated numerically 
+without exact arithmetic.
+
 # Extended help 
 The orient predicate is defined by the determinant 
 
@@ -14,7 +93,21 @@ The orient predicate is defined by the determinant
 \\text{orient}(p, q, r) = \\text{sgn} \\det \\begin{bmatrix} p_x & p_y & 1 \\\\ q_x & q_y & 1 \\\\ r_x & r_y & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} p_x-r_x & p_y-r_y \\\\ q_x-r_x & q_y-r_y \\end{bmatrix}.
 ```
 """
-orient_predicate(p, q, r) = orient(_getxy(p), _getxy(q), _getxy(r))
+orient_predicate(::Any, ::Any, ::Any)
+
+@static if USE_EXACTPREDICATES
+    orient_predicate(p, q, r) = orient(_getxy(p), _getxy(q), _getxy(r))
+else
+    orient_predicate(p, q, r) = _orient_predicate(_getxy(p), _getxy(q), _getxy(r))
+end
+function _orient_predicate(p, q, r)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    rx, ry = _getxy(r)
+    ux, uy = px - rx, py - ry
+    vx, vy = qx - rx, qy - ry
+    return sgn(ext(ux, uy, vx, vy))
+end
 
 """
     orient_predicate(p, q, r, s) -> Integer
@@ -44,14 +137,34 @@ Here, a positively oriented tetrahedron `(p, q, r, s)` takes the form
                                'r       
                                  '\\.    
 
+If ExactPredicates.jl has been disabled using Preferences.jl, the 
+determinant defining `orient` is evaluated numerically 
+without exact arithmetic.
+
 # Extended help 
 The orient predicate is defined by the determinant 
 
 ```math 
-\\text{orient}(p, q, r, s) = \\text{sgn} \\det \\begin{bmatrix} p_x & p_y & p_y & 1 \\\\ q_x & q_y & q_r & 1 \\\\ r_x & r_y & r_z & 1 \\\\ s_x & s_y & s_z & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} p_x - s_x & p_y - s_y & p_z - s_y \\\\ q_x - s_x & q_y - s_y & q_z - s_z \\\\ r_x - s_x & r_y - s_y & r_z - s_z \\end{bmatrix}.
+\\text{orient}(p, q, r, s) = \\text{sgn} \\det \\begin{bmatrix} p_x & p_y & p_z & 1 \\\\ q_x & q_y & q_z & 1 \\\\ r_x & r_y & r_z & 1 \\\\ s_x & s_y & s_z & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} p_x - s_x & p_y - s_y & p_z - s_y \\\\ q_x - s_x & q_y - s_y & q_z - s_z \\\\ r_x - s_x & r_y - s_y & r_z - s_z \\end{bmatrix}.
 ```
 """
-orient_predicate(p, q, r, s) = orient(p, q, r, s)
+orient_predicate(::Any, ::Any, ::Any, ::Any)
+
+@static if USE_EXACTPREDICATES
+    orient_predicate(p, q, r, s) = orient(_getxyz(p), _getxyz(q), _getxyz(r), _getxyz(s))
+else
+    orient_predicate(p, q, r, s) = _orient_predicate(_getxyz(p), _getxyz(q), _getxyz(r), _getxyz(s))
+end
+function _orient_predicate(p, q, r, s)
+    px, py, pz = _getxyz(p)
+    qx, qy, qz = _getxyz(q)
+    rx, ry, rz = _getxyz(r)
+    sx, sy, sz = _getxyz(s)
+    a, b, c = px - sx, py - sy, pz - sz
+    d, e, f = qx - sx, qy - sy, qz - sz
+    g, h, i = rx - sx, ry - sy, rz - sz
+    return sgn(det(a, b, c, d, e, f, g, h, i))
+end
 
 """
     incircle_predicate(a, b, c, p) -> Integer
@@ -62,14 +175,38 @@ Returns `ExactPredicates.incircle(a, b, c, p)`, in particular we return:
 - `0`: If `p` is on the circle defined by `(a, b, c)`.
 - `-1`: If `p` is outside the circle defined by `(a, b, c)`.
 
+If ExactPredicates.jl has been disabled using Preferences.jl, the 
+determinant defining `incircle` is evaluated numerically 
+without exact arithmetic.
+
 # Extended help 
 The incircle predicate is defined by the determinant
 
 ```math
-\\text{incircle}(a, b, c, d) = \\text{sgn} \\det \\begin{bmatrix} a_x & a_y & a_x^2 + a_y^2 & 1 \\\\ b_x & b_y & b_x62 + b_y^2 & 1 \\\\ c_x & c_y & c_x^2 + c_y^2 & 1 \\\\ d_x & d_y & d_x^2 + d_y^2 & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} a_x - d_x & a_y - d_y & (a_x - d_x)^2 + (a_y - d_y)^2 \\\\ b_x - d_x & b_y - d_y & (b_x - d_x)^2 + (b_y - d_y)^2 \\\\ c_x - d_x & c_y - d_y & (c_x - d_x)^2 + (c_y - d_y)^2 \\end{bmatrix}.
+\\text{incircle}(a, b, c, d) = \\text{sgn} \\det \\begin{bmatrix} a_x & a_y & a_x^2 + a_y^2 & 1 \\\\ b_x & b_y & b_x^2 + b_y^2 & 1 \\\\ c_x & c_y & c_x^2 + c_y^2 & 1 \\\\ d_x & d_y & d_x^2 + d_y^2 & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} a_x - d_x & a_y - d_y & (a_x - d_x)^2 + (a_y - d_y)^2 \\\\ b_x - d_x & b_y - d_y & (b_x - d_x)^2 + (b_y - d_y)^2 \\\\ c_x - d_x & c_y - d_y & (c_x - d_x)^2 + (c_y - d_y)^2 \\end{bmatrix}.
 ```
 """
-incircle_predicate(a, b, c, p) = incircle(_getxy(a), _getxy(b), _getxy(c), _getxy(p))
+incircle_predicate
+
+@static if USE_EXACTPREDICATES
+    incircle_predicate(a, b, c, p) = incircle(_getxy(a), _getxy(b), _getxy(c), _getxy(p))
+else
+    incircle_predicate(a, b, c, p) = _incircle_predicate(_getxy(a), _getxy(b), _getxy(c), _getxy(p))
+end
+function _incircle_predicate(p, q, r, a)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    rx, ry = _getxy(r)
+    ax, ay = _getxy(a)
+    qpx, qpy = qx - px, qy - py
+    rpx, rpy = rx - px, ry - py
+    apx, apy = ax - px, ay - py
+    aqx, aqy = ax - qx, ay - qy
+    rqx, rqy = rx - qx, ry - qy
+    val1 = ext(qpx, qpy, apx, apy) * inp(rpx, rpy, rqx, rqy)
+    val2 = ext(qpx, qpy, rpx, rpy) * inp(apx, apy, aqx, aqy)
+    return sgn(val1 - val2)
+end
 
 """
     parallelorder_predicate(a, b, p, q) -> Integer
@@ -79,8 +216,27 @@ Returns `ExactPredicates.parallelorder(a, b, p, q)`, in particular we return:
 - `1`: `q` is closer to the line `(a, b)` than `p`.
 - `0`: `p` and `q` are equidistant from the line `(a, b)`.
 - `-1`: `p` is closer to the line `(a, b)` than `q`.
+
+If ExactPredicates.jl has been disabled using Preferences.jl, the 
+the predicates defining `parallelorder` are all evaluated numerically 
+without exact arithmetic.
 """
-parallelorder_predicate(a, b, p, q) = parallelorder(_getxy(a), _getxy(b), _getxy(p), _getxy(q))
+parallelorder_predicate
+
+@static if USE_EXACTPREDICATES
+    parallelorder_predicate(a, b, p, q) = parallelorder(_getxy(a), _getxy(b), _getxy(p), _getxy(q))
+else
+    parallelorder_predicate(a, b, p, q) = _parallelorder_predicate(_getxy(a), _getxy(b), _getxy(p), _getxy(q))
+end
+function _parallelorder_predicate(a, b, p, q)
+    ax, ay = _getxy(a)
+    bx, by = _getxy(b)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    δx, δy = bx - ax, by - ay
+    qpx, qpy = qx - px, qy - py
+    return sgn(ext(δx, δy, qpx, qpy))
+end
 
 """
     sameside_predicate(a, b, p) -> Integer
@@ -109,13 +265,6 @@ function sameside_predicate(a, b, p)
     end
 end
 
-"""
-    opposite_signs(x, y) -> Bool
-
-From ExactPredicates.jl, returns `true` if `x` and `y` have opposite signs, and `false` otherwise.
-"""
-opposite_signs(x, y) = xor(x, y) == -2
-
 """
     meet_predicate(p, q, a, b) 
 
@@ -124,6 +273,10 @@ Returns `ExactPredicates.meet(p, q, a, b)`, in particular we return:
 - `1`: The open line segments `(p, q)` and `(a, b)` meet in a single point.
 - `0`: The closed line segments `[p, q]` and `[a, b]` meet in one or several points. 
 - `-1`: Otherwise.
+
+If ExactPredicates.jl has been disabled using Preferences.jl, then all 
+the `orient` predicates computed for this result are evaluated numerically 
+without exact arithmetic.
 """
 function meet_predicate(p, q, a, b)
     pqa = orient_predicate(p, q, a)
@@ -146,4 +299,35 @@ function meet_predicate(p, q, a, b)
     else
         return 0
     end
-end
\ No newline at end of file
+end
+
+@static if USE_EXACTPREDICATES
+    ExactPredicates.Codegen.@genpredicate function angle_is_acute(p::2, q::2, r::2)
+        pr = p - r
+        qr = q - r
+        ExactPredicates.Codegen.group!(pr...)
+        ExactPredicates.Codegen.group!(qr...)
+        return pr[1] * qr[1] + pr[2] * qr[2]
+    end
+else
+    function angle_is_acute(p, q, r)
+        px, py = _getxy(p)
+        qx, qy = _getxy(q)
+        rx, ry = _getxy(r)
+        ux, uy = px - rx, py - ry
+        vx, vy = qx - rx, qy - ry
+        return sgn(inp(ux, uy, vx, vy))
+    end
+end
+
+@doc """
+    angle_is_acute(p, q, r)
+
+Tests if the angle opposite `(p, q)` in the triangle `(p, q, r)`, 
+meaning `∠prq`, is acute, returning:
+
+- `1`: `∠prq` is acute.
+- `0`: `∠prq` is a right angle.
+- `-1`: `∠prq` is obtuse.
+"""
+angle_is_acute
\ No newline at end of file
diff --git a/src/predicates/predicates.jl b/src/predicates/predicates.jl
index d004af529..d3351db2d 100644
--- a/src/predicates/predicates.jl
+++ b/src/predicates/predicates.jl
@@ -540,14 +540,6 @@ function point_position_relative_to_witness_plane(tri::Triangulation, i, j, k, 
     return convert_certificate(cert, Cert.Above, Cert.On, Cert.Below)
 end
 
-ExactPredicates.Codegen.@genpredicate function angle_is_acute(p::2, q::2, r::2)
-    pr = p - r
-    qr = q - r
-    ExactPredicates.Codegen.group!(pr...)
-    ExactPredicates.Codegen.group!(qr...)
-    return pr[1] * qr[1] + pr[2] * qr[2]
-end
-
 """
     opposite_angle(p, q, r) -> Certificate
 
diff --git a/test/triangulation_validation.jl b/src/validation.jl
similarity index 78%
rename from test/triangulation_validation.jl
rename to src/validation.jl
index ba1cea877..7fa78b6f0 100644
--- a/test/triangulation_validation.jl
+++ b/src/validation.jl
@@ -1,5 +1,3 @@
-# maybe this should belong inside src for users at some point?
-
 #=
 """
     abstract type AbstractTriangulationState 
@@ -19,6 +17,22 @@ abstract type AbstractTriangulationState end
 test_state(state::AbstractTriangulationState) = state.flag
 Base.show(io::IO, state::AbstractTriangulationState) = print(io, summary(state))
 
+function compare_edge_vectors(E1, E2)
+    E1s = sort_edge_vector(collect(E1))
+    E2s = sort_edge_vector(collect(E2))
+    return E1s == E2s
+end
+
+function sort_edge_vector(E)
+    sorted_E = similar(E)
+    for i in eachindex(E)
+        u, v = E[i]
+        e = (min(u, v), max(u, v))
+        sorted_E[i] = e
+    end
+    return sort(sorted_E)
+end
+
 struct TriangleOrientationState <: AbstractTriangulationState
     flag::Bool
     bad_triangle::NTuple{3,Int}
@@ -34,12 +48,12 @@ function Base.summary(state::TriangleOrientationState)
 end
 function test_triangle_orientation(tri)
     for T in each_solid_triangle(tri)
-        cert = DT.triangle_orientation(tri, T)
-        flag = DT.is_positively_oriented(cert)
-        orientation = DT.is_positively_oriented(cert) ? :positively : DT.is_negatively_oriented(cert) ? :negatively : :degenerately
+        cert = triangle_orientation(tri, T)
+        flag = is_positively_oriented(cert)
+        orientation = is_positively_oriented(cert) ? :positively : is_negatively_oriented(cert) ? :negatively : :degenerately
         !flag && return TriangleOrientationState(flag, Int.(triangle_vertices(T)), orientation)
     end
-    return TriangleOrientationState(true, (DT.∅, DT.∅, DT.∅), :positive)
+    return TriangleOrientationState(true, (∅, ∅, ∅), :positive)
 end
 
 struct DelaunayCriterionState <: AbstractTriangulationState
@@ -52,7 +66,7 @@ function Base.summary(state::DelaunayCriterionState)
     if test_state(state)
         return "All the triangles are Delaunay."
     else
-        if state.bad_vertex !== DT.∅
+        if state.bad_vertex !== ∅
             return "The Delaunay criterion does not hold for the triangle-vertex pair ($(state.bad_triangle), $(state.bad_vertex))."
         else
             return "The test of the Delaunay criterion failed as there was a BoundsError when testing the visibility."
@@ -62,23 +76,23 @@ end
 function test_delaunay_criterion(tri)
     try
         points = get_points(tri)
-        triangle_tree = DT.BoundaryRTree(points)
-        segment_tree = DT.BoundaryRTree(points)
-        failures = Tuple{DT.triangle_type(tri),DT.integer_type(tri)}[]
+        triangle_tree = BoundaryRTree(points)
+        segment_tree = BoundaryRTree(points)
+        failures = Tuple{triangle_type(tri),integer_type(tri)}[]
         for T in each_solid_triangle(tri)
-            i, j, k = DT.triangle_vertices(T)
+            i, j, k = triangle_vertices(T)
             p, q, r = get_point(tri, i, j, k)
-            c = DT.triangle_circumcenter(p, q, r)
-            cr = DT.triangle_circumradius(p, q, r)
+            c = triangle_circumcenter(p, q, r)
+            cr = triangle_circumradius(p, q, r)
             xmin, xmax = getx(c) - cr, getx(c) + cr
             ymin, ymax = gety(c) - cr, gety(c) + cr
             any(!isfinite, (xmin, xmax, ymin, ymax)) && continue
-            bbox = DT.BoundingBox(xmin, xmax, ymin, ymax)
-            dbx = DT.DiametralBoundingBox(bbox, (i, j)) # just store (i, j) and get k via get_adjacent 
+            bbox = BoundingBox(xmin, xmax, ymin, ymax)
+            dbx = DiametralBoundingBox(bbox, (i, j)) # just store (i, j) and get k via get_adjacent 
             insert!(triangle_tree.tree, dbx)
         end
         for e in each_segment(tri)
-            i, j = DT.edge_vertices(e)
+            i, j = edge_vertices(e)
             p, q = get_point(tri, i, j)
             px, py = getxy(p)
             qx, qy = getxy(q)
@@ -86,32 +100,32 @@ function test_delaunay_criterion(tri)
             xmax = max(px, qx)
             ymin = min(py, qy)
             ymax = max(py, qy)
-            bbox = DT.BoundingBox(xmin, xmax, ymin, ymax)
-            dbx = DT.DiametralBoundingBox(bbox, (i, j))
+            bbox = BoundingBox(xmin, xmax, ymin, ymax)
+            dbx = DiametralBoundingBox(bbox, (i, j))
             insert!(segment_tree.tree, dbx)
         end
         for r in shuffle(collect(each_solid_vertex(tri)))
             !isempty(failures) && break
-            intersects = DT.get_intersections(triangle_tree, r, cache_id=1) # can't use multithreading here
+            intersects = get_intersections(triangle_tree, r, cache_id=1) # can't use multithreading here
             for box in intersects
-                i, j = DT.get_edge(box)
+                i, j = get_edge(box)
                 k = get_adjacent(tri, i, j)
                 any(==(r), (i, j, k)) && continue
-                T = DT.construct_triangle(DT.triangle_type(tri), i, j, k)
-                cert = DT.point_position_relative_to_circumcircle(tri, T, r)
-                c = DT.triangle_centroid(get_point(tri, i, j, k)...)
-                if DT.is_inside(cert)
-                    DT.is_boundary_edge(tri, i, j) && DT.is_right(DT.point_position_relative_to_line(tri, i, j, r)) && continue # if it's outside of the domain relative to this edge, just continue
-                    DT.is_boundary_edge(tri, j, k) && DT.is_right(DT.point_position_relative_to_line(tri, j, k, r)) && continue
-                    DT.is_boundary_edge(tri, k, i) && DT.is_right(DT.point_position_relative_to_line(tri, k, i, r)) && continue
-                    for (i, j) in DT.triangle_edges(i, j, k)
-                        DT.contains_segment(tri, i, j) && continue # visibility is defined according to the relative interior of the simplex, which means that it's fine if a segment can see the vertex
+                T = construct_triangle(triangle_type(tri), i, j, k)
+                cert = point_position_relative_to_circumcircle(tri, T, r)
+                c = triangle_centroid(get_point(tri, i, j, k)...)
+                if is_inside(cert)
+                    is_boundary_edge(tri, i, j) && is_right(point_position_relative_to_line(tri, i, j, r)) && continue # if it's outside of the domain relative to this edge, just continue
+                    is_boundary_edge(tri, j, k) && is_right(point_position_relative_to_line(tri, j, k, r)) && continue
+                    is_boundary_edge(tri, k, i) && is_right(point_position_relative_to_line(tri, k, i, r)) && continue
+                    for (i, j) in triangle_edges(i, j, k)
+                        contains_segment(tri, i, j) && continue # visibility is defined according to the relative interior of the simplex, which means that it's fine if a segment can see the vertex
                         if rand() < 1 / 2 # just testing both
-                            cert = DT.test_visibility(tri, i, j, r, shift=0.01, attractor=c)
+                            cert = test_visibility(tri, i, j, r, shift=0.01, attractor=c)
                         else
                             cert = test_visibility(tri, segment_tree, i, j, r, c)
                         end
-                        flag = DT.is_visible(cert)
+                        flag = is_visible(cert)
                         flag && push!(failures, (T, r))
                         flag && break
                     end
@@ -119,35 +133,35 @@ function test_delaunay_criterion(tri)
             end
         end
         if isempty(failures)
-            return DelaunayCriterionState(true, (DT.∅, DT.∅, DT.∅), DT.∅)
+            return DelaunayCriterionState(true, (∅, ∅, ∅), ∅)
         else
             return DelaunayCriterionState(false, Int.(failures[1][1]), Int(failures[1][2]))
         end
     catch e
-        e isa BoundsError && return DelaunayCriterionState(false, (DT.∅, DT.∅, DT.∅), DT.∅)
+        e isa BoundsError && return DelaunayCriterionState(false, (∅, ∅, ∅), ∅)
         rethrow(e)
     end
 end
 function test_visibility(tri::Triangulation, segment_tree, i, j, k, centroid=nothing)
-    if DT.is_boundary_edge(tri, j, i)
+    if is_boundary_edge(tri, j, i)
         i, j = j, i
     end
     p, q, a = get_point(tri, i, j, k)
-    if DT.is_boundary_edge(tri, i, j)
-        side_e = DT.point_position_relative_to_line(p, q, a)
-        DT.is_right(side_e) && return DT.Cert.Invisible
+    if is_boundary_edge(tri, i, j)
+        side_e = point_position_relative_to_line(p, q, a)
+        is_right(side_e) && return Cert.Invisible
     end
     # Need to see if i or j is a boundary node without the other being a boundary node.
     # If this is the case, it's possible that one of them mistakenly sees k.
     for u in (i, j)
-        flag, g = DT.is_boundary_node(tri, u)
+        flag, g = is_boundary_node(tri, u)
         if flag
-            ℓ = DT.get_left_boundary_node(tri, u, g)
-            cert = DT.point_position_relative_to_line(tri, ℓ, u, k)
-            DT.is_right(cert) && return DT.Cert.Invisible
-            ℓ = DT.get_right_boundary_node(tri, u, g)
-            cert = DT.point_position_relative_to_line(tri, u, ℓ, k)
-            DT.is_right(cert) && return DT.Cert.Invisible
+            ℓ = get_left_boundary_node(tri, u, g)
+            cert = point_position_relative_to_line(tri, ℓ, u, k)
+            is_right(cert) && return Cert.Invisible
+            ℓ = get_right_boundary_node(tri, u, g)
+            cert = point_position_relative_to_line(tri, u, ℓ, k)
+            is_right(cert) && return Cert.Invisible
         end
     end
     flags = falses(10)
@@ -163,18 +177,17 @@ function test_visibility(tri::Triangulation, segment_tree, i, j, k, centroid=not
         mx, my = getxy(m)
         ax, ay = getxy(a)
         xmin, ymin, xmax, ymax = min(mx, ax), min(my, ay), max(mx, ax), max(my, ay)
-        bbox = DT.BoundingBox(xmin, xmax, ymin, ymax)
-        intersections = DT.get_intersections(segment_tree, bbox; cache_id=2)
+        bbox = BoundingBox(xmin, xmax, ymin, ymax)
+        intersections = get_intersections(segment_tree, bbox; cache_id=2)
         for box in intersections
-            u, v = DT.get_edge(box)
-            # !DT.edges_are_disjoint((i, j), (u, v)) && continue
+            u, v = get_edge(box)
             p′, q′ = get_point(tri, u, v)
-            cert = DT.line_segment_intersection_type(m, a, p′, q′)
-            flags[idx] = !DT.has_no_intersections(cert) && !DT.is_touching(cert)
+            cert = line_segment_intersection_type(m, a, p′, q′)
+            flags[idx] = !has_no_intersections(cert) && !is_touching(cert)
             flags[idx] && break
         end
     end
-    return all(flags) ? DT.Cert.Invisible : DT.Cert.Visible
+    return all(flags) ? Cert.Invisible : Cert.Visible
 end
 
 struct EdgesHaveTwoIncidentTrianglesState <: AbstractTriangulationState
@@ -191,22 +204,22 @@ function Base.summary(state::EdgesHaveTwoIncidentTrianglesState)
 end
 function test_each_edge_has_two_incident_triangles(tri)
     for e in each_edge(tri)
-        i, j = DT.edge_vertices(e)
+        i, j = edge_vertices(e)
         vᵢⱼ = get_adjacent(tri, i, j)
         vⱼᵢ = get_adjacent(tri, j, i)
-        if DT.is_boundary_edge(tri, j, i)
-            flag = DT.is_ghost_vertex(vᵢⱼ) && DT.edge_exists(vⱼᵢ)
+        if is_boundary_edge(tri, j, i)
+            flag = is_ghost_vertex(vᵢⱼ) && edge_exists(vⱼᵢ)
             !flag && return EdgesHaveTwoIncidentTrianglesState(flag, Int.(edge_vertices(e)))
-        elseif DT.is_boundary_edge(tri, i, j)
-            flag = DT.is_ghost_vertex(vⱼᵢ) && DT.edge_exists(vᵢⱼ)
+        elseif is_boundary_edge(tri, i, j)
+            flag = is_ghost_vertex(vⱼᵢ) && edge_exists(vᵢⱼ)
             !flag && return EdgesHaveTwoIncidentTrianglesState(flag, Int.(edge_vertices(e)))
         else
-            flag = DT.edge_exists(vᵢⱼ) && DT.edge_exists(vⱼᵢ)
+            flag = edge_exists(vᵢⱼ) && edge_exists(vⱼᵢ)
             !flag && return EdgesHaveTwoIncidentTrianglesState(flag, Int.(edge_vertices(e)))
         end
     end
-    DT.clear_empty_features!(tri)
-    return EdgesHaveTwoIncidentTrianglesState(true, (DT.∅, DT.∅))
+    clear_empty_features!(tri)
+    return EdgesHaveTwoIncidentTrianglesState(true, (∅, ∅))
 end
 
 struct AdjacentMapState <: AbstractTriangulationState
@@ -225,7 +238,7 @@ function Base.summary(state::AdjacentMapState)
 end
 function test_adjacent_map_matches_triangles(tri)
     for T in each_triangle(tri)
-        u, v, w = DT.triangle_vertices(T)
+        u, v, w = triangle_vertices(T)
         flag1 = get_adjacent(tri, u, v) == w
         !flag1 && return AdjacentMapState(flag1, Int.(triangle_vertices(T)), Int.((u, v)), get_adjacent(tri, u, v))
         flag2 = get_adjacent(tri, v, w) == u
@@ -233,8 +246,8 @@ function test_adjacent_map_matches_triangles(tri)
         flag3 = get_adjacent(tri, w, u) == v
         !flag3 && return AdjacentMapState(flag3, Int.(triangle_vertices(T)), Int.((w, u)), get_adjacent(tri, w, u))
     end
-    DT.clear_empty_features!(tri)
-    return AdjacentMapState(true, (DT.∅, DT.∅, DT.∅), (DT.∅, DT.∅), DT.∅)
+    clear_empty_features!(tri)
+    return AdjacentMapState(true, (∅, ∅, ∅), (∅, ∅), ∅)
 end
 
 struct Adjacent2VertexMapState <: AbstractTriangulationState
@@ -253,20 +266,20 @@ function Base.summary(state::Adjacent2VertexMapState)
     end
 end
 function test_adjacent2vertex_map_matches_triangles(tri)
-    E = DT.edge_type(tri)
+    E = edge_type(tri)
     for T in each_triangle(tri)
         u, v, w = triangle_vertices(T)
         for (u, v, w) in ((u, v, w), (v, w, u), (w, u, v))
-            vw = DT.construct_edge(E, v, w)
+            vw = construct_edge(E, v, w)
             Su = get_adjacent2vertex(tri, u)
-            flag = DT.contains_edge(vw, Su)
+            flag = contains_edge(vw, Su)
             if !flag
-                _Su = Set{NTuple{2,Int}}(Int.(DT.edge_vertices(e)) for e in each_edge(Su))
-                return Adjacent2VertexMapState(flag, Int.(DT.triangle_vertices(T)), Int.((v, w)), Int(u), _Su)
+                _Su = Set{NTuple{2,Int}}(Int.(edge_vertices(e)) for e in each_edge(Su))
+                return Adjacent2VertexMapState(flag, Int.(triangle_vertices(T)), Int.((v, w)), Int(u), _Su)
             end
         end
     end
-    return Adjacent2VertexMapState(true, (DT.∅, DT.∅, DT.∅), (DT.∅, DT.∅), DT.∅, Set{NTuple{2,Int}}())
+    return Adjacent2VertexMapState(true, (∅, ∅, ∅), (∅, ∅), ∅, Set{NTuple{2,Int}}())
 end
 
 struct AdjacentMapAdjacent2VertexMapState <: AbstractTriangulationState
@@ -286,14 +299,14 @@ function Base.summary(state::AdjacentMapAdjacent2VertexMapState)
 end
 function test_adjacent_map_matches_adjacent2vertex_map(tri)
     for (k, S) in get_adjacent2vertex(get_adjacent2vertex(tri))
-        _S = Set{NTuple{2,Int}}(Int.(DT.edge_vertices(e)) for e in each_edge(S))
+        _S = Set{NTuple{2,Int}}(Int.(edge_vertices(e)) for e in each_edge(S))
         for e in each_edge(S)
             flag = get_adjacent(tri, e) == k
             !flag && return AdjacentMapAdjacent2VertexMapState(flag, Int(k), get_adjacent(tri, e), _S, Int.(edge_vertices(e)))
         end
     end
-    DT.clear_empty_features!(tri)
-    return AdjacentMapAdjacent2VertexMapState(true, DT.∅, DT.∅, Set{NTuple{2,Int}}(), (DT.∅, DT.∅))
+    clear_empty_features!(tri)
+    return AdjacentMapAdjacent2VertexMapState(true, ∅, ∅, Set{NTuple{2,Int}}(), (∅, ∅))
 end
 
 struct Adjacent2VertexMapAdjacentMapState <: AbstractTriangulationState
@@ -320,14 +333,14 @@ function test_adjacent2vertex_map_matches_adjacent_map(tri)
         k = get_adjacent(tri, e)
         if haskey(adj2v, k)
             S = get_adjacent2vertex(tri, k)
-            flag = DT.contains_edge(e, S)
+            flag = contains_edge(e, S)
         else
             flag = false
         end
         !flag && return Adjacent2VertexMapAdjacentMapState(flag, Int(k), Int.(edge_vertices(e)), haskey(adj2v, k))
     end
-    DT.clear_empty_features!(tri)
-    return Adjacent2VertexMapAdjacentMapState(true, DT.∅, (DT.∅, DT.∅), true)
+    clear_empty_features!(tri)
+    return Adjacent2VertexMapAdjacentMapState(true, ∅, (∅, ∅), true)
 end
 
 struct GraphState <: AbstractTriangulationState
@@ -345,14 +358,14 @@ end
 function test_graph_contains_all_vertices(tri)
     all_vertices = Set{Int}()
     for T in each_triangle(tri) # need a method that doesn't use the graph
-        i, j, k = DT.triangle_vertices(T)
+        i, j, k = triangle_vertices(T)
         push!(all_vertices, i, j, k)
     end
     for i in all_vertices
-        flag = DT.has_vertex(tri, i)
+        flag = has_vertex(tri, i)
         !flag && return GraphState(flag, i)
     end
-    return GraphState(true, DT.∅)
+    return GraphState(true, ∅)
 end
 
 struct GraphAdjacentMapState <: AbstractTriangulationState
@@ -368,23 +381,23 @@ function Base.summary(state::GraphAdjacentMapState)
     end
 end
 function test_graph_matches_adjacent_map(tri)
-    E = DT.edge_type(tri)
+    E = edge_type(tri)
     adj_dict = get_adjacent(get_adjacent(tri))
     for e in each_edge(tri)
-        i, j = DT.edge_vertices(e)
-        if DT.has_ghost_triangles(tri) || !DT.is_ghost_edge(i, j)
-            eᵢⱼ = DT.construct_edge(E, i, j)
-            eⱼᵢ = DT.construct_edge(E, j, i)
-            flag = if !DT.has_multiple_sections(tri)
+        i, j = edge_vertices(e)
+        if has_ghost_triangles(tri) || !is_ghost_edge(i, j)
+            eᵢⱼ = construct_edge(E, i, j)
+            eⱼᵢ = construct_edge(E, j, i)
+            flag = if !has_multiple_sections(tri)
                 eᵢⱼ ∈ keys(adj_dict) && eⱼᵢ ∈ keys(adj_dict)
             else
-                DT.edge_exists(tri, i, j) && DT.edge_exists(tri, j, i)
+                edge_exists(tri, i, j) && edge_exists(tri, j, i)
             end
             !flag && return GraphAdjacentMapState(flag, Int.(edge_vertices(e)))
         end
     end
-    DT.clear_empty_features!(tri)
-    return GraphAdjacentMapState(true, (DT.∅, DT.∅))
+    clear_empty_features!(tri)
+    return GraphAdjacentMapState(true, (∅, ∅))
 end
 
 struct AdjacentMapGraphState <: AbstractTriangulationState
@@ -409,7 +422,7 @@ function Base.summary(state::AdjacentMapGraphState)
 end
 function test_adjacent_map_matches_graph(tri)
     for (e, k) in get_adjacent(get_adjacent(tri))
-        i, j = DT.edge_vertices(e)
+        i, j = edge_vertices(e)
         flag1 = all(∈(get_neighbours(tri, i)), (j, k))
         !flag1 && return AdjacentMapGraphState(flag1, Int(i), Int.((i, j)), Int.((j, k)), false, Int(k))
         flag2 = all(∈(get_neighbours(tri, j)), (k, i))
@@ -417,7 +430,7 @@ function test_adjacent_map_matches_graph(tri)
         flag3 = all(∈(get_neighbours(tri, k)), (i, j))
         !flag3 && return AdjacentMapGraphState(flag3, Int(k), Int.((i, j)), Int.((i, j)), true, Int(k))
     end
-    return AdjacentMapGraphState(true, DT.∅, (DT.∅, DT.∅), (DT.∅, DT.∅), true, DT.∅)
+    return AdjacentMapGraphState(true, ∅, (∅, ∅), (∅, ∅), true, ∅)
 end
 
 struct GraphTrianglesState <: AbstractTriangulationState
@@ -431,7 +444,7 @@ function Base.summary(state::GraphTrianglesState)
     if test_state(state)
         return "The graph correctly matches the triangle set."
     else
-        if state.bad_vertex == DT.∅
+        if state.bad_vertex == ∅
             return "The graph is inconsistent with the triangle set, as one of the vertices of $(state.bad_triangle) is not a vertex in the graph."
         else
             return "The graph is inconsistent with the triangle set. The triangle $(state.bad_triangle) is in the triangle set but either $(state.bad_edge[1]) or $(state.bad_edge[2]) are not in $(state.bad_vertex)'s neighbourhood."
@@ -441,7 +454,7 @@ end
 function test_graph_matches_triangles(tri)
     for T in each_triangle(tri)
         try
-            i, j, k = DT.triangle_vertices(T)
+            i, j, k = triangle_vertices(T)
             flag1 = all(∈(get_neighbours(tri, i)), (j, k))
             !flag1 && return GraphTrianglesState(flag1, Int.(triangle_vertices(T)), Int.((j, k)), Int(i))
             flag2 = all(∈(get_neighbours(tri, j)), (k, i))
@@ -449,11 +462,11 @@ function test_graph_matches_triangles(tri)
             flag3 = all(∈(get_neighbours(tri, k)), (i, j))
             !flag3 && return GraphTrianglesState(flag3, Int.(triangle_vertices(T)), Int.((i, j)), Int(k))
         catch e
-            e isa KeyError && return GraphTrianglesState(false, Int.(triangle_vertices(T)), (DT.∅, DT.∅), DT.∅)
+            e isa KeyError && return GraphTrianglesState(false, Int.(triangle_vertices(T)), (∅, ∅), ∅)
             rethrow(e)
         end
     end
-    return GraphTrianglesState(true, (DT.∅, DT.∅, DT.∅), (DT.∅, DT.∅), DT.∅)
+    return GraphTrianglesState(true, (∅, ∅, ∅), (∅, ∅), ∅)
 end
 
 struct SegmentState <: AbstractTriangulationState
@@ -477,31 +490,31 @@ function Base.summary(state::SegmentState)
 end
 function test_segments(tri)
     for e in each_segment(tri)
-        flag = DT.edge_exists(tri, e) || DT.edge_exists(tri, DT.reverse_edge(e))
+        flag = edge_exists(tri, e) || edge_exists(tri, reverse_edge(e))
         !flag && return SegmentState(flag, Int.(edge_vertices(e)), 0)
-        flag = DT.contains_edge(e, get_interior_segments(tri)) || DT.contains_edge(DT.reverse_edge(e), get_interior_segments(tri))
+        flag = contains_edge(e, get_interior_segments(tri)) || contains_edge(reverse_edge(e), get_interior_segments(tri))
         if !flag
-            flag = DT.contains_boundary_edge(tri, e) || DT.contains_boundary_edge(tri, DT.reverse_edge(e))
+            flag = contains_boundary_edge(tri, e) || contains_boundary_edge(tri, reverse_edge(e))
             !flag && return SegmentState(flag, Int.(edge_vertices(e)), 1)
         end
     end
     for e in get_interior_segments(tri)
-        flag = DT.edge_exists(tri, e) || DT.edge_exists(tri, DT.reverse_edge(e))
+        flag = edge_exists(tri, e) || edge_exists(tri, reverse_edge(e))
         !flag && return SegmentState(flag, Int.(edge_vertices(e)), 0)
     end
     int_segs = Set{NTuple{2,Int}}()
     all_segs = Set{NTuple{2,Int}}()
     for e in get_interior_segments(tri)
-        u, v = DT.edge_vertices(e)
+        u, v = edge_vertices(e)
         push!(int_segs, minmax(u, v))
     end
     for e in each_segment(tri)
-        u, v = DT.edge_vertices(e)
+        u, v = edge_vertices(e)
         push!(all_segs, minmax(u, v))
     end
     flag = int_segs ⊆ all_segs
-    !flag && return SegmentState(flag, (DT.∅, DT.∅), 2)
-    return SegmentState(true, (DT.∅, DT.∅), 0)
+    !flag && return SegmentState(flag, (∅, ∅), 2)
+    return SegmentState(true, (∅, ∅), 0)
 end
 
 struct BoundaryEdgeMapBoundaryNodesState <: AbstractTriangulationState
@@ -518,15 +531,15 @@ function Base.summary(state::BoundaryEdgeMapBoundaryNodesState)
         return "The boundary edge map is inconsistent with the boundary nodes. The edge $(state.bad_edge) maps to $(state.bad_pos) in the boundary edge map but $(state.bad_pos) corresponds to the edge $(state.bad_edge_2)."
     end
 end
-function test_boundary_edge_map_matches_boundary_nodes(tri::DT.Triangulation)
-    boundary_edge_map = DT.get_boundary_edge_map(tri)
+function test_boundary_edge_map_matches_boundary_nodes(tri::Triangulation)
+    boundary_edge_map = get_boundary_edge_map(tri)
     for (edge, pos) in boundary_edge_map
-        u = DT.get_boundary_nodes(DT.get_boundary_nodes(tri, pos[1]), pos[2])
-        v = DT.get_boundary_nodes(DT.get_boundary_nodes(tri, pos[1]), pos[2] + 1)
-        flag = u == DT.initial(edge)
+        u = get_boundary_nodes(get_boundary_nodes(tri, pos[1]), pos[2])
+        v = get_boundary_nodes(get_boundary_nodes(tri, pos[1]), pos[2] + 1)
+        flag = u == initial(edge)
         !flag && return BoundaryEdgeMapBoundaryNodesState(flag, Int.(edge_vertices(edge)), Int.((u, v)), pos)
     end
-    return BoundaryEdgeMapBoundaryNodesState(true, (DT.∅, DT.∅), (DT.∅, DT.∅), ())
+    return BoundaryEdgeMapBoundaryNodesState(true, (∅, ∅), (∅, ∅), ())
 end
 
 struct BoundaryNodesBoundaryEdgeMapState <: AbstractTriangulationState
@@ -543,53 +556,53 @@ function Base.summary(state::BoundaryNodesBoundaryEdgeMapState)
         return "The boundary nodes are inconsistent with the boundary edge map. The edge $(state.bad_edge) maps to $(state.bad_pos_2) but is at $(state.bad_pos)."
     end
 end
-function test_boundary_nodes_matches_boundary_edge_map(tri::DT.Triangulation)
-    boundary_nodes = DT.get_boundary_nodes(tri)
-    E = DT.edge_type(tri)
-    if DT.has_multiple_curves(tri)
-        nc = DT.num_curves(tri)
+function test_boundary_nodes_matches_boundary_edge_map(tri::Triangulation)
+    boundary_nodes = get_boundary_nodes(tri)
+    E = edge_type(tri)
+    if has_multiple_curves(tri)
+        nc = num_curves(tri)
         for i in 1:nc
             curve_nodes = get_boundary_nodes(tri, i)
-            ns = DT.num_sections(curve_nodes)
+            ns = num_sections(curve_nodes)
             for j in 1:ns
                 segment_nodes = get_boundary_nodes(curve_nodes, j)
-                ne = DT.num_boundary_edges(segment_nodes)
+                ne = num_boundary_edges(segment_nodes)
                 for k in 1:ne
                     left_node = get_boundary_nodes(segment_nodes, k)
                     right_node = get_boundary_nodes(segment_nodes, k + 1)
-                    edge = DT.construct_edge(E, left_node, right_node)
-                    pos = DT.get_boundary_edge_map(tri, edge)
+                    edge = construct_edge(E, left_node, right_node)
+                    pos = get_boundary_edge_map(tri, edge)
                     flag = pos == ((i, j), k)
                     !flag && return BoundaryNodesBoundaryEdgeMapState(flag, Int.(edge_vertices(edge)), pos, ((i, j), k))
                 end
             end
         end
-    elseif DT.has_multiple_sections(tri)
-        ns = DT.num_sections(tri)
+    elseif has_multiple_sections(tri)
+        ns = num_sections(tri)
         for i in 1:ns
             segment_nodes = get_boundary_nodes(tri, i)
-            ne = DT.num_boundary_edges(segment_nodes)
+            ne = num_boundary_edges(segment_nodes)
             for j in 1:ne
                 left_node = get_boundary_nodes(segment_nodes, j)
                 right_node = get_boundary_nodes(segment_nodes, j + 1)
-                edge = DT.construct_edge(E, left_node, right_node)
-                pos = DT.get_boundary_edge_map(tri, edge)
+                edge = construct_edge(E, left_node, right_node)
+                pos = get_boundary_edge_map(tri, edge)
                 flag = pos == (i, j)
                 !flag && return BoundaryNodesBoundaryEdgeMapState(flag, Int.(edge_vertices(edge)), pos, (i, j))
             end
         end
     else
-        ne = DT.num_boundary_edges(boundary_nodes)
+        ne = num_boundary_edges(boundary_nodes)
         for i in 1:ne
             left_node = get_boundary_nodes(boundary_nodes, i)
             right_node = get_boundary_nodes(boundary_nodes, i + 1)
-            edge = DT.construct_edge(E, left_node, right_node)
-            pos = DT.get_boundary_edge_map(tri, edge)
+            edge = construct_edge(E, left_node, right_node)
+            pos = get_boundary_edge_map(tri, edge)
             flag = pos == (get_boundary_nodes(tri), i)
             !flag && return BoundaryNodesBoundaryEdgeMapState(flag, Int.(edge_vertices(edge)), pos, (get_boundary_nodes(tri), i))
         end
     end
-    return BoundaryNodesBoundaryEdgeMapState(true, (DT.∅, DT.∅), (), ())
+    return BoundaryNodesBoundaryEdgeMapState(true, (∅, ∅), (), ())
 end
 
 abstract type AbstractTriangulationIteratorState <: AbstractTriangulationState end
@@ -687,8 +700,8 @@ end
 function GhostEdgeIteratorState(unique_flag, length_flag, output_flag, iterator_length, correct_length)
     return IteratorState(unique_flag, length_flag, output_flag, iterator_length, correct_length, :ghost_edge, :edge, :each_ghost_edge)
 end
-function test_iterators(tri::DT.Triangulation)
-    I = DT.integer_type(tri)
+function test_iterators(tri::Triangulation)
+    I = integer_type(tri)
     T = NTuple{3,I}
     E = NTuple{2,I}
     solid_triangles = T[]
@@ -701,16 +714,16 @@ function test_iterators(tri::DT.Triangulation)
     ghost_edges = E[]
     all_edges = E[]
     for T in each_triangle(tri)
-        i, j, k = DT.triangle_vertices(T)
+        i, j, k = triangle_vertices(T)
         push!(all_triangles, (i, j, k))
-        if DT.is_ghost_triangle(i, j, k)
+        if is_ghost_triangle(i, j, k)
             push!(ghost_triangles, (i, j, k))
         else
             push!(solid_triangles, (i, j, k))
         end
-        for (u, v) in DT.triangle_edges(i, j, k)
+        for (u, v) in triangle_edges(i, j, k)
             push!(all_edges, (u, v))
-            if DT.is_ghost_edge(u, v)
+            if is_ghost_edge(u, v)
                 push!(ghost_edges, (u, v))
             else
                 push!(solid_edges, (u, v))
@@ -718,7 +731,7 @@ function test_iterators(tri::DT.Triangulation)
         end
         for k in (i, j, k)
             push!(all_vertices, k)
-            if DT.is_ghost_vertex(k)
+            if is_ghost_vertex(k)
                 push!(ghost_vertices, k)
             else
                 push!(solid_vertices, k)
@@ -732,15 +745,15 @@ function test_iterators(tri::DT.Triangulation)
     unique_ghost_triangle_flag = allunique(ghost_triangles)
     unique_ghost_edge_flag = allunique(ghost_edges)
     for (i, e) in enumerate(all_edges)
-        u, v = DT.edge_vertices(e)
+        u, v = edge_vertices(e)
         all_edges[i] = (min(u, v), max(u, v))
     end
     for (i, e) in enumerate(solid_edges)
-        u, v = DT.edge_vertices(e)
+        u, v = edge_vertices(e)
         solid_edges[i] = (min(u, v), max(u, v))
     end
     for (i, e) in enumerate(ghost_edges)
-        u, v = DT.edge_vertices(e)
+        u, v = edge_vertices(e)
         ghost_edges[i] = (min(u, v), max(u, v))
     end
     unique!(all_edges)
@@ -770,9 +783,9 @@ function test_iterators(tri::DT.Triangulation)
     sort!(all_vertices)
     sort!(solid_vertices)
     sort!(ghost_vertices)
-    triangle_output_flag = DT.compare_triangle_collections(all_triangles, collect(each_triangle(tri)))
-    solid_triangle_output_flag = DT.compare_triangle_collections(solid_triangles, collect(each_solid_triangle(tri)))
-    ghost_triangle_output_flag = DT.compare_triangle_collections(ghost_triangles, collect(each_ghost_triangle(tri)))
+    triangle_output_flag = compare_triangle_collections(all_triangles, each_triangle(tri))
+    solid_triangle_output_flag = compare_triangle_collections(solid_triangles, each_solid_triangle(tri))
+    ghost_triangle_output_flag = compare_triangle_collections(ghost_triangles, each_ghost_triangle(tri))
     edge_output_flag = compare_edge_vectors(all_edges, each_edge(tri))
     solid_edge_output_flag = compare_edge_vectors(solid_edges, each_solid_edge(tri))
     ghost_edge_output_flag = compare_edge_vectors(ghost_edges, each_ghost_edge(tri))
@@ -813,14 +826,14 @@ function test_no_duplicate_segments(tri)
     interior_segments = get_interior_segments(tri)
     all_segments = get_all_segments(tri)
     for e in each_edge(interior_segments)
-        flag = !DT.contains_edge(DT.reverse_edge(e), interior_segments)
-        !flag && return DuplicateSegmentsState(flag, true, Int.(DT.edge_vertices(e)))
+        flag = !contains_edge(reverse_edge(e), interior_segments)
+        !flag && return DuplicateSegmentsState(flag, true, Int.(edge_vertices(e)))
     end
     for e in each_edge(all_segments)
-        flag = !DT.contains_edge(DT.reverse_edge(e), all_segments)
-        !flag && return DuplicateSegmentsState(flag, false, Int.(DT.edge_vertices(e)))
+        flag = !contains_edge(reverse_edge(e), all_segments)
+        !flag && return DuplicateSegmentsState(flag, false, Int.(edge_vertices(e)))
     end
-    return DuplicateSegmentsState(true, true, (DT.∅, DT.∅))
+    return DuplicateSegmentsState(true, true, (∅, ∅))
 end
 
 struct TriangulationState <: AbstractTriangulationState
@@ -850,11 +863,11 @@ struct TriangulationState <: AbstractTriangulationState
     ghost_edge_iterator_state::IteratorState
 end
 
-function TriangulationState(tri::DT.Triangulation)
-    has_ghosts = DT.has_ghost_triangles(tri)
-    DT.delete_ghost_triangles!(tri)
-    DT.add_ghost_triangles!(tri)
-    DT.clear_empty_features!(tri)
+function TriangulationState(tri::Triangulation)
+    has_ghosts = has_ghost_triangles(tri)
+    delete_ghost_triangles!(tri)
+    add_ghost_triangles!(tri)
+    clear_empty_features!(tri)
     state = TriangulationState(
         Adjacent2VertexMapAdjacentMapState(tri),
         Adjacent2VertexMapState(tri),
@@ -873,7 +886,7 @@ function TriangulationState(tri::DT.Triangulation)
         TriangleOrientationState(tri),
         test_iterators(tri)...
     )
-    !has_ghosts && DT.delete_ghost_triangles!(tri)
+    !has_ghosts && delete_ghost_triangles!(tri)
     return state
 end
 
@@ -898,7 +911,14 @@ function Base.show(io::IO, triangulation_state::TriangulationState)
     return
 end
 
-function validate_triangulation(tri; print_result=true)
+"""
+    validate_triangulation(tri::Triangulation; print_result=true) -> Bool 
+
+Tests if `tri` is a valid `Triangulation`. Returns `true` if so, 
+and `false` otherwise. If `print_result=true` and `tri` is not a 
+valid triangulation, all the issues with `tri` will be printed.
+"""
+function validate_triangulation(tri::Triangulation; print_result=true)
     state = TriangulationState(tri)
     print_result && !test_state(state) && println(state)
     return test_state(state)
diff --git a/test/Project.toml b/test/Project.toml
index a75aeaaec..676712c9d 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -17,6 +17,8 @@ GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
diff --git a/test/constrained_triangulation/segment_location.jl b/test/constrained_triangulation/segment_location.jl
index a6f2aa5fa..c46658514 100644
--- a/test/constrained_triangulation/segment_location.jl
+++ b/test/constrained_triangulation/segment_location.jl
@@ -1,6 +1,7 @@
 using ..DelaunayTriangulation
 const DT = DelaunayTriangulation
 using CairoMakie
+using Preferences
 
 @testset "Shewchuk Example: A small example with some collinearities" begin
     tri = shewchuk_example_constrained()
@@ -278,29 +279,31 @@ end
     end
 end
 
-@testset "A previously broken example" begin
-    for m in 1:100
-        a = -0.1
-        b = 0.1
-        c = -0.01
-        d = 0.01
-        nx = 25
-        ny = 25
-        tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
-        tri = triangulate(get_points(tri))
-        for i in 2:24
-            add_segment!(tri, i, 600 + i)
+if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+    @testset "A previously broken example" begin
+        for m in 1:100
+            a = -0.1
+            b = 0.1
+            c = -0.01
+            d = 0.01
+            nx = 25
+            ny = 25
+            tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
+            tri = triangulate(get_points(tri))
+            for i in 2:24
+                add_segment!(tri, i, 600 + i)
+            end
+            tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
+            tri = triangulate(get_points(tri))
+            e = (23, 71)
+            history = DT.PointLocationHistory{NTuple{3,Int},NTuple{2,Int},Int}()
+            jump_and_march(tri, get_point(tri, 71);
+                m=nothing, k=23, store_history=true, history=history)
+            collinear_segments = history.collinear_segments
+            DT.connect_segments!(collinear_segments)
+            DT.extend_segments!(collinear_segments, e)
+            @test collinear_segments == [(23, 47), (47, 71)]
         end
-        tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
-        tri = triangulate(get_points(tri))
-        e = (23, 71)
-        history = DT.PointLocationHistory{NTuple{3,Int},NTuple{2,Int},Int}()
-        jump_and_march(tri, get_point(tri, 71);
-            m=nothing, k=23,  store_history=true, history=history)
-        collinear_segments = history.collinear_segments
-        DT.connect_segments!(collinear_segments)
-        DT.extend_segments!(collinear_segments, e)
-        @test collinear_segments == [(23, 47), (47, 71)]
     end
 end
 
diff --git a/test/data_structures/curves.jl b/test/data_structures/curves.jl
index f2ee2ba2e..2f04225a3 100644
--- a/test/data_structures/curves.jl
+++ b/test/data_structures/curves.jl
@@ -5,6 +5,7 @@ using LinearAlgebra
 EllipticalArc = DT.EllipticalArc # shadow
 using CairoMakie
 using StableRNGs
+using Preferences
 using StructEquality
 using ForwardDiff
 using ReferenceTests
@@ -1978,36 +1979,38 @@ end
     end
 
     ## Closest point 
-    for spl in (spl, pspl)
-        for _ in 1:500
-            p = randn(2) * 30 |> Tuple
-            t′, q = DT.get_closest_point(spl, p)
-            @test q ⪧ spl(t′)
-            _t, _q = closest_point_on_curve(spl, p)
-            if !(spl == pspl)
-                @test _t ≈ t′ rtol = 1e-1 atol = 1e-1
-            elseif t ≠ 0 && t ≠ 1 && t′ ≠ 0 && t′ ≠ 1
-                @test _t ≈ t′ rtol = 1e-1 atol = 1e-1
+    if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+        for spl in (spl, pspl)
+            for _ in 1:500
+                p = randn(2) * 30 |> Tuple
+                t′, q = DT.get_closest_point(spl, p)
+                @test q ⪧ spl(t′)
+                _t, _q = closest_point_on_curve(spl, p)
+                if !(spl == pspl)
+                    @test _t ≈ t′ rtol = 1e-1 atol = 1e-1
+                elseif t ≠ 0 && t ≠ 1 && t′ ≠ 0 && t′ ≠ 1
+                    @test _t ≈ t′ rtol = 1e-1 atol = 1e-1
+                end
+                @test q ⪧ _q rtol = 1e-1 atol = 1e-1
+                @test DT.dist(p, _q) ≈ DT.dist(p, q) rtol = 1e-1 atol = 1e-1
             end
-            @test q ⪧ _q rtol = 1e-1 atol = 1e-1
-            @test DT.dist(p, _q) ≈ DT.dist(p, q) rtol = 1e-1 atol = 1e-1
-        end
-        for t in LinRange(0, 1, 150)
-            p = spl(t)
-            t′, q = DT.get_closest_point(spl, p)
-            @test q ⪧ spl(t′)
-            if !(spl == pspl)
-                @test t ≈ t′ rtol = 1e-1 atol = 1e-1
-            elseif t ≠ 0 && t ≠ 1
-                @test t ≈ t′ rtol = 1e-1 atol = 1e-1
+            for t in LinRange(0, 1, 150)
+                p = spl(t)
+                t′, q = DT.get_closest_point(spl, p)
+                @test q ⪧ spl(t′)
+                if !(spl == pspl)
+                    @test t ≈ t′ rtol = 1e-1 atol = 1e-1
+                elseif t ≠ 0 && t ≠ 1
+                    @test t ≈ t′ rtol = 1e-1 atol = 1e-1
+                end
+                @test p ⪧ q rtol = 1e-1 atol = 1e-1
+            end
+            @test DT.get_closest_point(spl, spl.control_points[1]) == (0.0, spl.control_points[1])
+            if spl ≠ pspl
+                @test DT.get_closest_point(spl, spl.control_points[end]) == (1.0, spl.control_points[end])
+            else
+                @test DT.get_closest_point(spl, spl.control_points[end]) == (0.0, spl.control_points[end])
             end
-            @test p ⪧ q rtol = 1e-1 atol = 1e-1
-        end
-        @test DT.get_closest_point(spl, spl.control_points[1]) == (0.0, spl.control_points[1])
-        if spl ≠ pspl
-            @test DT.get_closest_point(spl, spl.control_points[end]) == (1.0, spl.control_points[end])
-        else
-            @test DT.get_closest_point(spl, spl.control_points[end]) == (0.0, spl.control_points[end])
         end
     end
 
diff --git a/test/helper_functions.jl b/test/helper_functions.jl
index cb8d59565..2862d37f5 100644
--- a/test/helper_functions.jl
+++ b/test/helper_functions.jl
@@ -14,7 +14,7 @@ getxy((0.0, 0.0)) # avoid shadow
 
 const DT = DelaunayTriangulation
 
-include("triangulation_validation.jl")
+using DelaunayTriangulation: validate_triangulation
 
 function complicated_geometry()
     x1 = [collect(LinRange(0, 2, 4)),
@@ -2197,6 +2197,14 @@ function plot_small_angle_complexes(enricher)
     return fig
 end
 
+using DelaunayTriangulation:
+    test_adjacent2vertex_map_matches_adjacent_map,
+    test_state,
+    test_adjacent_map_matches_adjacent2vertex_map,
+    test_each_edge_has_two_incident_triangles,
+    test_triangle_orientation,
+    test_iterators
+
 export validate_triangulation
 export @_adj
 export _make_graph_from_adjacency
@@ -2210,11 +2218,11 @@ export complicated_geometry
 export validate_refinement
 export validate_statistics
 export validate_tessellation
-export compare_edge_vectors 
-export simple_geometry 
+export compare_edge_vectors
+export simple_geometry
 export fixed_shewchuk_example_constrained
 export ⊢
-export is_sink 
+export is_sink
 export is_conformal
 export _validate_offcenter
 export example_with_special_corners
@@ -2235,12 +2243,6 @@ export get_child_from_tree
 export traverse_tree
 export poor_triangulation_example
 export example_triangulation
-export test_adjacent2vertex_map_matches_adjacent_map
-export test_state 
-export test_adjacent_map_matches_adjacent2vertex_map
-export test_each_edge_has_two_incident_triangles
-export test_triangle_orientation
-export test_iterators
 export example_empty_triangulation
 export shewchuk_example_constrained
 export test_segment_triangle_intersections
@@ -2263,4 +2265,10 @@ export all_points_are_inside
 export all_diametral_circles_are_empty
 export compute_diametral_lens
 export get_points_in_diametral_lens
+export test_adjacent2vertex_map_matches_adjacent_map
+export test_state
+export test_adjacent_map_matches_adjacent2vertex_map
+export test_each_edge_has_two_incident_triangles
+export test_triangle_orientation
+export test_iterators
 end
\ No newline at end of file
diff --git a/test/interfaces/points.jl b/test/interfaces/points.jl
index 696f5580f..fc609e339 100644
--- a/test/interfaces/points.jl
+++ b/test/interfaces/points.jl
@@ -1,7 +1,7 @@
 using ..DelaunayTriangulation
 const DT = DelaunayTriangulation
 using Test
-using StaticArraysCore
+using StaticArrays
 using StatsBase
 import GeometryBasics: Point2f
 
@@ -235,3 +235,35 @@ end
     )
 end
 
+@testset "getz/_getz/getxyz/_getxyz" begin
+    p = (1.0, 2.0, 3.0)
+    q = (1.0f0, 2.0f0, 3.0f0)
+    r = [1.0f0, 4.0f0, 5.0f0]
+    s = @SVector [3, 7, 5]
+    t = [5.0, 13.0, -5.0]
+    
+    @test DT.getz(p) === 3.0 
+    @test DT._getz(p) === 3.0
+    @test DT.getxyz(p) === (1.0, 2.0, 3.0)
+    @test DT._getxyz(p) === (1.0, 2.0, 3.0)
+
+    @test DT.getz(q) === 3.0f0
+    @test DT._getz(q) === 3.0
+    @test DT.getxyz(q) === (1.0f0, 2.0f0, 3.0f0)
+    @test DT._getxyz(q) === (1.0, 2.0, 3.0)
+
+    @test DT.getz(r) === 5.0f0
+    @test DT._getz(r) === 5.0
+    @test DT.getxyz(r) === (1.0f0, 4.0f0, 5.0f0)
+    @test DT._getxyz(r) === (1.0, 4.0, 5.0)
+
+    @test DT.getz(s) === 5
+    @test DT._getz(s) === 5.0
+    @test DT.getxyz(s) === (3, 7, 5)
+    @test DT._getxyz(s) === (3.0, 7.0, 5.0)
+
+    @test DT.getz(t) === -5.0
+    @test DT._getz(t) === -5.0
+    @test DT.getxyz(t) === (5.0, 13.0, -5.0)
+    @test DT._getxyz(t) === (5.0, 13.0, -5.0)
+end
\ No newline at end of file
diff --git a/test/operations/delete_ghost_triangles.jl b/test/operations/delete_ghost_triangles.jl
index b57d6804c..0b4d963ff 100644
--- a/test/operations/delete_ghost_triangles.jl
+++ b/test/operations/delete_ghost_triangles.jl
@@ -3,7 +3,7 @@ const DT = DelaunayTriangulation
 using StructEquality
 using ..DelaunayTriangulation: Triangulation
 using StaticArrays
-
+using Preferences
 
 @testset "Deleting ghost triangles" begin
     tri, label_map, index_map = simple_geometry()
@@ -15,15 +15,19 @@ using StaticArrays
 end
 
 @testset "Making sure the ghost triangles are correctly adjusted on disjoint sets" begin
-    θ = LinRange(0, 2π, 500) |> collect
+    θ = LinRange(0, 2π, 100) |> collect
     θ[end] = 0
     xy = Vector{Vector{Vector{NTuple{2,Float64}}}}()
     cx = 0.0
     for i in 1:2
-        push!(xy, [[(cx + cos(θ), sin(θ)) for θ in θ]])
-        push!(xy, [[(cx + 0.5cos(θ), 0.5sin(θ)) for θ in reverse(θ)]])
+        push!(xy, [[(cx + cos(θ) + 1e-3rand(), sin(θ) + 1e-3rand()) for θ in θ]]) # use perturbations to also do this check without ExactPredicates loaded. iszero makes the boundary closed
+        push!(xy, [[(cx + 0.5cos(θ) + 1e-3rand(), 0.5sin(θ) + 1e-3rand()) for θ in reverse(θ)]])
         cx += 3.0
     end
+    xy[1][1][end] = xy[1][1][1]
+    xy[2][1][end] = xy[2][1][1] 
+    xy[3][1][end] = xy[3][1][1]
+    xy[4][1][end] = xy[4][1][1] 
     boundary_nodes, points = convert_boundary_points_to_indices(xy)
     tri = triangulate(points; boundary_nodes=boundary_nodes)
     @test validate_triangulation(tri)
diff --git a/test/operations/delete_holes.jl b/test/operations/delete_holes.jl
index 05bb77df2..8c965db82 100644
--- a/test/operations/delete_holes.jl
+++ b/test/operations/delete_holes.jl
@@ -3,7 +3,7 @@ const DT = DelaunayTriangulation
 using CairoMakie
 using StableRNGs
 using StaticArrays
-
+using Preferences
 
 @testset "Simple domain" begin
     p1 = (0.2, 0.2)
@@ -142,68 +142,70 @@ end
     validate_statistics(tri)
 end
 
-@testset "Interior holes that were already triangles" begin
-    p1 = (0.0, 0.0)
-    p2 = (1.0, 0.0)
-    p3 = (1.0, 1.0)
-    p4 = (0.0, 1.0)
-    p5 = (0.2, 0.2)
-    p6 = (0.8, 0.2)
-    p7 = (0.8, 0.8)
-    pts = [p1, p2, p3, p4, p5, p6, p7]
-    tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]], delete_holes=false, delete_ghosts=false)
-    @test num_ghost_triangles(tri) == 4
-    @test num_ghost_edges(tri) == 4 
-    points_to_process = DT.find_all_points_to_delete(tri)
-    @test points_to_process == Set((DT.𝒢,))
-    triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process)
-    @test DT.compare_triangle_collections(triangles_to_delete, Set((
-        (1, 4, DT.𝒢),
-        (4, 3, DT.𝒢),
-        (3, 2, DT.𝒢),
-        (2, 1, DT.𝒢),
-        (5, 6, 7)
-    )))
-    tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]])
-    @test validate_triangulation(tri)
-    validate_statistics(tri)
+if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+    @testset "Interior holes that were already triangles" begin
+        p1 = (0.0, 0.0)
+        p2 = (1.0, 0.0)
+        p3 = (1.0, 1.0)
+        p4 = (0.0, 1.0)
+        p5 = (0.2, 0.2)
+        p6 = (0.8, 0.2)
+        p7 = (0.8, 0.8)
+        pts = [p1, p2, p3, p4, p5, p6, p7]
+        tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]], delete_holes=false, delete_ghosts=false)
+        @test num_ghost_triangles(tri) == 4
+        @test num_ghost_edges(tri) == 4
+        points_to_process = DT.find_all_points_to_delete(tri)
+        @test points_to_process == Set((DT.𝒢,))
+        triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process)
+        @test DT.compare_triangle_collections(triangles_to_delete, Set((
+            (1, 4, DT.𝒢),
+            (4, 3, DT.𝒢),
+            (3, 2, DT.𝒢),
+            (2, 1, DT.𝒢),
+            (5, 6, 7)
+        )))
+        tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]])
+        @test validate_triangulation(tri)
+        validate_statistics(tri)
 
-    tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7], [7, 6, 5]]], delete_holes=false, delete_ghosts=false)
-    points_to_process = DT.find_all_points_to_delete(tri)
-    @test points_to_process == Set((DT.𝒢,))
-    triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process)
-    @test DT.compare_triangle_collections(triangles_to_delete, Set((
-        (1, 4, DT.𝒢),
-        (4, 3, DT.𝒢),
-        (3, 2, DT.𝒢),
-        (2, 1, DT.𝒢),
-        (5, 6, 7)
-    )))
-    tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7], [7, 6, 5]]])
-    @test validate_triangulation(tri)
-    validate_statistics(tri)
+        tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7], [7, 6, 5]]], delete_holes=false, delete_ghosts=false)
+        points_to_process = DT.find_all_points_to_delete(tri)
+        @test points_to_process == Set((DT.𝒢,))
+        triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process)
+        @test DT.compare_triangle_collections(triangles_to_delete, Set((
+            (1, 4, DT.𝒢),
+            (4, 3, DT.𝒢),
+            (3, 2, DT.𝒢),
+            (2, 1, DT.𝒢),
+            (5, 6, 7)
+        )))
+        tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7], [7, 6, 5]]])
+        @test validate_triangulation(tri)
+        validate_statistics(tri)
 
-    tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]], delete_holes=false, delete_ghosts=false)
-    points_to_process = DT.find_all_points_to_delete(tri)
-    @test points_to_process == Set((DT.𝒢,))
-    triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process)
-    @test DT.compare_triangle_collections(triangles_to_delete, Set((
-        (1, 4, DT.𝒢),
-        (4, 3, DT.𝒢),
-        (3, 2, DT.𝒢),
-        (2, 1, DT.𝒢),
-        (5, 6, 7)
-    )))
-    tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]])
-    @test validate_triangulation(tri)
-    validate_statistics(tri)
+        tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]], delete_holes=false, delete_ghosts=false)
+        points_to_process = DT.find_all_points_to_delete(tri)
+        @test points_to_process == Set((DT.𝒢,))
+        triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process)
+        @test DT.compare_triangle_collections(triangles_to_delete, Set((
+            (1, 4, DT.𝒢),
+            (4, 3, DT.𝒢),
+            (3, 2, DT.𝒢),
+            (2, 1, DT.𝒢),
+            (5, 6, 7)
+        )))
+        tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]])
+        @test validate_triangulation(tri)
+        validate_statistics(tri)
+    end
 end
 
 @testset "A previously broken example" begin
     a = 4 / 5
     t = LinRange(0, 2π, 6)
     x = @. a * (2cos(t) + cos(2t))
-    y = @. a * (2sin(t) - sin(2t)) # this used to be tri = generate_mesh(x, y, 5.0); points = get_points(tri); bn_nodes = get_boundary_nodes(tri)
+    y = @. a * (2sin(t) - sin(2t))
     points = [2.4 -0.152786404500042 -1.047213595499958 -1.047213595499958 -0.1527864045000426; 0.0 1.051462224238267 1.70130161670408 -1.70130161670408 -1.051462224238268]
     bn_nodes = [1, 2, 3, 4, 5, 1]
     tri = triangulate(points; boundary_nodes=bn_nodes, delete_ghosts=false, delete_holes=false)
@@ -226,7 +228,7 @@ end
     a = 4 / 5
     t = LinRange(0, 2π, 100)
     x = @. a * (2cos(t) + cos(2t))
-    y = @. a * (2sin(t) - sin(2t)) # this used to be tri = generate_mesh(x, y, 15.0); points = get_points(tri); bn_nodes = get_boundary_nodes(tri)
+    y = @. a * (2sin(t) - sin(2t)) 
     points = [2.4 2.390342532619651 2.361486661646358 2.313780262033189 2.247797423889998 2.16432999441036 2.064375923309772 1.949124594926212 1.819939377400058 1.678337662931399 1.525968712312314 1.364589651123909 1.196039993627106 1.022215093005981 0.8450389328830836 0.6664366846599798 0.4883074580921845 0.3124976685434042 0.1407754336467903 -0.02519360519360872 -0.1838686646186263 -0.333854062108855 -0.473917012361638 -0.6030027096667079 -0.7202464727404799 -0.824982776042241 -0.9167510419119772 -0.9952981201288541 -1.060577434849961 -1.112744832492037 -1.152151217102105 -1.179332111274063 -1.194994329879372 -1.2 -1.19534820274028 -1.182154550372295 -1.161629044931084 -1.13505259139796 -1.103752559560399 -1.069077803180918 -1.032373553014235 -0.994956601357817 -0.958091190190919 -0.9229660026456057 -0.8906726387622712 -0.8621859315182511 -0.8383464283873544 -0.8198453276894746 -0.8072121183067701 -0.8008051266355886 -0.8008051266355887 -0.8072121183067701 -0.8198453276894746 -0.8383464283873544 -0.8621859315182508 -0.890672638762271 -0.9229660026456055 -0.9580911901909194 -0.9949566013578169 -1.032373553014235 -1.069077803180918 -1.103752559560399 -1.135052591397961 -1.161629044931084 -1.182154550372295 -1.19534820274028 -1.2 -1.194994329879372 -1.179332111274063 -1.152151217102105 -1.112744832492038 -1.060577434849961 -0.9952981201288543 -0.9167510419119774 -0.8249827760422405 -0.7202464727404799 -0.6030027096667084 -0.473917012361639 -0.3338540621088567 -0.1838686646186263 -0.02519360519360934 0.1407754336467914 0.3124976685434047 0.4883074580921845 0.6664366846599791 0.8450389328830821 1.022215093005981 1.196039993627106 1.364589651123908 1.525968712312312 1.678337662931397 1.819939377400058 1.949124594926211 2.064375923309773 2.16432999441036 2.247797423889998 2.313780262033188 2.361486661646357 2.390342532619651 0.3170611563256605 0.3222593190661471 -0.6183519451946706 -0.6497375598908914 -0.03601883360826612 0.6813132651727021 -0.6442990922286358 0.677301056434819 -0.03212472989883586 1.049402030184789 -0.6932272745230642 -0.3561747556617257 1.050855112561903 -0.3586613134669976 -0.6921937990949052 -0.7738275096391242 1.420372082515759 -0.6465445728766354 -0.7401785386842585 -0.8130282799705698 1.553206818654827 -0.4939148325245979 -0.5486837455157717 -0.3610262116797264 -0.73633174592415 -0.5446205627063301 -0.5160523518705824 1.230246578448748 1.093304308222101 0.8980444010246713 1.226019949979472 -0.4997667391156656 -0.726253210863807 -0.6697924507551998 -0.4891744968539103 -0.4634488008923348 -0.301600826224346 -0.2783275986241442 -0.1177692287900377 -0.09393337228059367 -0.2546427830837441 -0.0703514373220399 0.08836743861011277 -0.2313237479039875 0.1137973151090641 0.2724811315356115 0.2996278178466897 0.4682688930880384 0.1380153936605681 0.406551246715369 0.6046552182453121 -0.426784127732563 -0.3923263038801806 -0.208465408072328 -0.3607313375439971 -0.1857892526190982 0.4960419968524275 -0.1415763437967015 0.03591030790050782 0.8685411342185969 -0.1926902902309472 -0.8592155035553608 1.716817787430798 -0.8576022838754369 -0.1921452564055306 -0.6683334786985783 0.8654953794762578 0.6952267131341282 -0.6287806122165458 0.4999252549098672 0.3218569037134379 -0.5986272435711641 0.1449039324769594 2.203335189185611 -1.101667594592808 -1.101667594592807 -0.0286998603756622 -0.3499716813694415 -0.4450559779573338 -0.04783795748406422 0.09480479796065325 2.114352958860066 -1.054192316972474 -1.054192316972475 -0.3317870997991825 1.39665732177026 -0.6206155967973439 -0.7760417249729167 0.4988340023743217 0.1461550962121355 0.3242397406248529 0.9036042833648704 -0.5300851339677813 -0.3652221764727946 -0.9373851776825108 1.875039225311675 -0.9376540476291644 -0.5066666554551625 0.2368123037070368 -0.02440461564986861 0.150646683235148 1.991415825682697 -0.9957149117466657 -0.9957009139360333 0.5336878756578392 0.7210970221965526 -0.1928809320498361 0.7689122343492278 -0.6467071096582683 1.285560087566004 -0.6388529779077368 -0.1969978949942105 0.1875732054682595 -0.4975220279896986 -0.7186754363850769 1.577166091039188 -0.8584906546541111; 0.0 0.000204308591503799 0.001629535973135443 0.005472026448394285 0.01287939060935179 0.02492716987386698 0.04259671987083449 0.06675468790926908 0.09813443380113843 0.1373197116453414 0.1847308933321656 0.2406139730886964 0.3050325470249727 0.3778629130904994 0.4587923858949879 0.5473208683088024 0.6427656684861544 0.7442694978083142 0.8508115330838226 0.9612213760111374 1.074195695220302 1.188317291935738 1.302076290158998 1.413893116908774 1.522142908049466 1.625180951076655 1.721368758301788 1.809100352482636 1.886828342268937 1.953089366954632 2.006528498919954 2.04592220767042 2.070199511290803 2.078460969082653 2.069995202699299 2.044292671697285 2.001056472471559 1.94020997634528 1.86190117239507 1.766503632611802 1.654614070392519 1.527046517275516 1.384823196404126 1.229162223576608 1.061462317070302 0.8832847449107668 0.6963327821298029 0.5024289901161495 0.3034906647750211 0.1015038293221599 -0.1015038293221605 -0.3034906647750203 -0.50242899011615 -0.6963327821298022 -0.8832847449107647 -1.061462317070301 -1.229162223576608 -1.384823196404126 -1.527046517275515 -1.654614070392518 -1.766503632611802 -1.86190117239507 -1.940209976345281 -2.001056472471559 -2.044292671697285 -2.069995202699299 -2.078460969082653 -2.070199511290803 -2.04592220767042 -2.006528498919954 -1.953089366954633 -1.886828342268938 -1.809100352482636 -1.721368758301788 -1.625180951076654 -1.522142908049467 -1.413893116908775 -1.302076290158999 -1.18831729193574 -1.074195695220302 -0.9612213760111378 -0.8508115330838222 -0.7442694978083139 -0.6427656684861544 -0.5473208683088027 -0.4587923858949887 -0.3778629130904991 -0.3050325470249729 -0.2406139730886968 -0.184730893332166 -0.137319711645342 -0.09813443380113851 -0.06675468790926926 -0.04259671987083431 -0.02492716987386698 -0.01287939060935179 -0.005472026448394374 -0.001629535973135488 -0.000204308591503799 -0.5473096449852685 0.5447620596126443 0.01428450391194042 -0.4309321583994598 0.7536074241825893 0.3536290172290472 0.4171992061645116 -0.354274248521617 -0.7697113089269275 0.1945973624956351 0.8115101356751746 -1.00610749817081 -0.1925650703607583 1.006349758355742 -0.8137846879949843 1.19333488738332 0.07348683780306886 -1.266821725186389 1.366146471157314 -1.324086653415268 -0.04174935346989145 1.135404524959176 0.9456563651429131 0.8229544359251958 -0.995445054751826 -0.9480022448317108 -0.741992356781049 -0.13995947020735 0.002345879688797254 -0.08111410706343182 0.130762025431598 -1.127145434944549 0.9963834095129508 0.6153059195085407 0.5329795234649529 0.3351323238863241 0.4528041803207281 0.2552428653861011 0.3751002341424148 0.1760741054849401 0.05622243244883154 -0.02310920543998908 0.0916097709503603 -0.1431517599572751 -0.1033611549989752 0.01112539166303821 -0.1796378655622734 -0.06423821903495724 -0.2933071885887684 0.1346998725151015 0.06389662809144629 -0.06697977586584479 -0.263147702667072 -0.3418086892373844 -0.4602230597841197 -0.5395791962871652 -0.260267079463506 0.5739282454666178 0.5007002518322551 0.2757413735876992 -0.8939205770864943 -1.486342121175829 0.0009313928164653013 1.487273513992295 0.8909822595409715 -0.6270118467477362 -0.2692763527002932 -0.1472140205276644 0.2179450297267495 0.4510351900656056 0.3329413377177828 -0.2182665505265999 -0.6631900778915942 0.0 1.908144246886934 -1.908144246886936 -0.4106426344765784 -0.6482825657476072 0.1353078965823503 -0.2192301054498285 0.3004465131589842 4.336808689942018e-17 1.825914653945082 -1.825914653945081 0.6415639095768699 -0.08968357289590935 1.254408379505216 -1.164673062563948 -0.4517076066327135 0.6541117112279523 -0.36426907951982 0.09997393244042915 0.731016147780223 -0.8325312305826258 1.623909218280243 0.0001552321360775167 -1.623753986144165 -0.5474230356714256 0.1741645772010534 -0.5950229828132197 -0.4758032489975526 8.081639737527319e-6 -1.724612653719712 1.724620735359448 0.2655807939077015 0.1867114129881211 -0.717357771233564 0.01594342553679897 -1.111060401383298 -0.004534585080420719 1.115594986463718 0.7308601982314702 0.4541306336569979 -0.3789760643386818 -1.406227077779554 0.08072235390443261 1.325504723875121]
     bn_nodes = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 1]
     tri = triangulate(points; boundary_nodes=bn_nodes)
diff --git a/test/operations/delete_point.jl b/test/operations/delete_point.jl
index 052a2989b..c9ac43f87 100644
--- a/test/operations/delete_point.jl
+++ b/test/operations/delete_point.jl
@@ -1,6 +1,7 @@
 using ..DelaunayTriangulation
 const DT = DelaunayTriangulation
 using StableRNGs
+using Preferences
 using StaticArrays
 
 @testset verbose = true "Deleting interior nodes" begin
@@ -46,51 +47,54 @@ using StaticArrays
             end
         end
     end
-    @testset "Lattice with a single boundary index" begin
-        for j in 1:20
-            rng1 = StableRNG(j)
-            a, b = sort(10randn(rng1, 2))
-            c, d = sort(15randn(rng1, 2))
-            nx = rand(rng1, 5:25)
-            ny = rand(rng1, 5:25)
-            tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
-            points = get_points(tri)
-            n = nx * ny
-            for k in 1:(n÷10)
-                rng2 = StableRNG(j + k * n)
-                i = rand(rng2, each_solid_vertex(tri) |> collect)
-                while DT.is_boundary_node(tri, i)[1]
+
+    if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+        @testset "Lattice with a single boundary index" begin
+            for j in 1:20
+                rng1 = StableRNG(j)
+                a, b = sort(10randn(rng1, 2))
+                c, d = sort(15randn(rng1, 2))
+                nx = rand(rng1, 5:25)
+                ny = rand(rng1, 5:25)
+                tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
+                points = get_points(tri)
+                n = nx * ny
+                for k in 1:(n÷10)
+                    rng2 = StableRNG(j + k * n)
                     i = rand(rng2, each_solid_vertex(tri) |> collect)
+                    while DT.is_boundary_node(tri, i)[1]
+                        i = rand(rng2, each_solid_vertex(tri) |> collect)
+                    end
+                    delete_point!(tri, i; rng=rng2)
+                    @test validate_triangulation(tri)
                 end
-                delete_point!(tri, i; rng=rng2)
-                @test validate_triangulation(tri)
             end
+            tri = triangulate_rectangle(0, 1, 0, 1, 25, 25; delete_ghosts=false, single_boundary=true)
+            add_segment!(tri, 7, 7 + 25)
+            @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 2)
+            @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 7)
+            @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 7 + 25)
+            @test_throws DT.InvalidVertexDeletionError delete_point!(tri, -1)
         end
-        tri = triangulate_rectangle(0, 1, 0, 1, 25, 25; delete_ghosts=false, single_boundary=true)
-        add_segment!(tri, 7, 7 + 25)
-        @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 2)
-        @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 7)
-        @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 7 + 25)
-        @test_throws DT.InvalidVertexDeletionError delete_point!(tri, -1)
-    end
-    @testset "Lattice with multiple boundary indices" begin
-        for j in 1:20
-            rng1 = StableRNG(j)
-            a, b = sort(10randn(rng1, 2))
-            c, d = sort(15randn(rng1, 2))
-            nx = rand(rng1, 5:25)
-            ny = rand(rng1, 5:25)
-            tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=false)
-            points = get_points(tri)
-            n = nx * ny
-            for k in 1:(n÷10)
-                rng2 = StableRNG(j + k * n)
-                i = rand(rng2, each_solid_vertex(tri) |> collect)
-                while DT.is_boundary_node(tri, i)[1]
+        @testset "Lattice with multiple boundary indices" begin
+            for j in 1:20
+                rng1 = StableRNG(j)
+                a, b = sort(10randn(rng1, 2))
+                c, d = sort(15randn(rng1, 2))
+                nx = rand(rng1, 5:25)
+                ny = rand(rng1, 5:25)
+                tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=false)
+                points = get_points(tri)
+                n = nx * ny
+                for k in 1:(n÷10)
+                    rng2 = StableRNG(j + k * n)
                     i = rand(rng2, each_solid_vertex(tri) |> collect)
+                    while DT.is_boundary_node(tri, i)[1]
+                        i = rand(rng2, each_solid_vertex(tri) |> collect)
+                    end
+                    delete_point!(tri, i; rng=rng2)
+                    @test validate_triangulation(tri)
                 end
-                delete_point!(tri, i; rng=rng2)
-                @test validate_triangulation(tri)
             end
         end
     end
diff --git a/test/operations/lock_convex_hull.jl b/test/operations/lock_convex_hull.jl
index 99aae9631..4aa922a43 100644
--- a/test/operations/lock_convex_hull.jl
+++ b/test/operations/lock_convex_hull.jl
@@ -4,7 +4,6 @@ using CairoMakie
 using DataStructures
 using StaticArrays
 
-
 @testset "lock_convex_hull" begin
     for _ in 1:150
         pts = rand(2, 500)
diff --git a/test/point_location/jump_and_march.jl b/test/point_location/jump_and_march.jl
index f9b7a90cd..d055e5f0a 100644
--- a/test/point_location/jump_and_march.jl
+++ b/test/point_location/jump_and_march.jl
@@ -4,7 +4,7 @@ using LinearAlgebra
 using Random
 using StableRNGs
 using StatsBase
-
+using Preferences
 
 
 tri, label_map, index_map = simple_geometry()
@@ -121,40 +121,42 @@ rep[3].y = mean([12.0, 6.0, 2.0, 4.0, 6.0, 10.0])
             rep[3].y = mean([12.0, 6.0, 2.0, 4.0, 6.0, 10.0])
         end
 
-        @testset "Test that we can find a point in every triangle" begin
-            for run in 1:6
-                for V in each_triangle(tri.triangles)
-                    rand() < 0.5 && continue # skip 50%
-                    if !DT.is_exterior_ghost_triangle(tri, triangle_vertices(V)...)
-                        i, j, k = triangle_vertices(V)
-                        p, q, r = get_point(tri, i, j, k)
-                        local c
-                        c = (p .+ q .+ r) ./ 3
-                        for k in DT.each_solid_vertex(tri)
-                            _V = jump_and_march(tri, c; k, concavity_protection=true)
-                            @test DT.is_positively_oriented(DT.triangle_orientation(tri, _V))
-                            if !DT.is_ghost_triangle(_V...)
-                                @test DT.compare_triangles(_V, V) &&
-                                      DT.is_inside(DT.point_position_relative_to_triangle(tri,
-                                    _V,
-                                    c))
-                            else
-                                local V1, V2
-                                V1 = DT.sort_triangle(V)
-                                V2 = DT.sort_triangle(_V)
-                                i1 = DT.geti(V1)
-                                i2 = DT.geti(V2)
-                                if i1 ≠ i2
-                                    i1 = i1 - 1
+        if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+            @testset "Test that we can find a point in every triangle" begin
+                for run in 1:6
+                    for V in each_triangle(tri.triangles)
+                        rand() < 0.5 && continue # skip 50%
+                        if !DT.is_exterior_ghost_triangle(tri, triangle_vertices(V)...)
+                            i, j, k = triangle_vertices(V)
+                            p, q, r = get_point(tri, i, j, k)
+                            local c
+                            c = (p .+ q .+ r) ./ 3
+                            for k in DT.each_solid_vertex(tri)
+                                _V = jump_and_march(tri, c; k, concavity_protection=true)
+                                @test DT.is_positively_oriented(DT.triangle_orientation(tri, _V))
+                                if !DT.is_ghost_triangle(_V...)
+                                    @test DT.compare_triangles(_V, V) &&
+                                          DT.is_inside(DT.point_position_relative_to_triangle(tri,
+                                        _V,
+                                        c))
+                                else
+                                    local V1, V2
+                                    V1 = DT.sort_triangle(V)
+                                    V2 = DT.sort_triangle(_V)
+                                    i1 = DT.geti(V1)
+                                    i2 = DT.geti(V2)
+                                    if i1 ≠ i2
+                                        i1 = i1 - 1
+                                    end
+                                    if i1 ≠ i2
+                                        i1 = i1 + 1
+                                    end
+                                    _V = DT.construct_triangle(typeof(V), i1, DT.getj(V1), DT.getk(V1))
+                                    @test DT.compare_triangles(_V, V) &&
+                                          DT.is_inside(DT.point_position_relative_to_triangle(tri,
+                                        _V,
+                                        c))
                                 end
-                                if i1 ≠ i2
-                                    i1 = i1 + 1
-                                end
-                                _V = DT.construct_triangle(typeof(V), i1, DT.getj(V1), DT.getk(V1))
-                                @test DT.compare_triangles(_V, V) &&
-                                      DT.is_inside(DT.point_position_relative_to_triangle(tri,
-                                    _V,
-                                    c))
                             end
                         end
                     end
@@ -162,14 +164,16 @@ rep[3].y = mean([12.0, 6.0, 2.0, 4.0, 6.0, 10.0])
             end
         end
 
-        @testset "Test that we don't break for points already in the triangulation" begin
-            for _ in 1:6
-                for k in DT.each_solid_vertex(tri)
-                    rand() < 1 / 2 && continue
-                    for j in DT.each_solid_vertex(tri)
-                        _V = jump_and_march(tri, get_point(tri, k); k=j)
-                        @test k ∈ triangle_vertices(_V)
-                        @test DT.is_positively_oriented(DT.triangle_orientation(tri, _V))
+        if !load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && tri_idx ≠ 3
+            @testset "Test that we don't break for points already in the triangulation" begin
+                for _ in 1:6
+                    for k in DT.each_solid_vertex(tri)
+                        rand() < 1 / 2 && continue
+                        for j in DT.each_solid_vertex(tri)
+                            _V = jump_and_march(tri, get_point(tri, k); k=j)
+                            @test k ∈ triangle_vertices(_V)
+                            @test DT.is_positively_oriented(DT.triangle_orientation(tri, _V))
+                        end
                     end
                 end
             end
diff --git a/test/predicates/general.jl b/test/predicates/general.jl
index 65da51d8e..e72a94879 100644
--- a/test/predicates/general.jl
+++ b/test/predicates/general.jl
@@ -4,7 +4,41 @@ using StaticArrays
 const DT = DelaunayTriangulation
 using ..DelaunayTriangulation: Certificate
 
+@testset "ext" begin
+    for _ in 1:50
+        p = @SVector rand(2)
+        q = @SVector rand(2)
+        @test ExactPredicates.ext(p, q) == DT.ext(p..., q...)
+    end
+end
+
+@testset "inp" begin
+    for _ in 1:50
+        p = @SVector rand(2)
+        q = @SVector rand(2)
+        @test ExactPredicates.inp(p, q) == DT.inp(p..., q...)
+    end
+end
+
+@testset "det" begin
+    for _ in 1:50
+        a, b, c, d = rand(4)
+        @test ExactPredicates.det(a, b, c, d) == DT.det(a, b, c, d)
 
+        a, b, c, d, e, f, g, h, i = rand(9)
+        @test ExactPredicates.det(a, b, c, d, e, f, g, h, i) == DT.det(a, b, c, d, e, f, g, h, i)
+    end
+end
+
+@testset "sgn" begin
+    @test DT.sgn(2.0) === 1
+    @test DT.sgn(2) === 1
+    @test DT.sgn(-1.0) === -1
+    @test DT.sgn(-1) === -1
+    @test DT.sgn(0.0) === 0
+    @test DT.sgn(-0.0) === 0
+    @test DT.sgn(0) === 0
+end
 
 @testset "Opposite signs" begin
     @test DT.opposite_signs(1, -1)
@@ -22,27 +56,35 @@ end
     @testset "Random points" begin
         for _ in 1:1500
             p, q, r = eachcol(rand(2, 3))
-            @test DT.orient_predicate(p, q, r) == orient(p, q, r)
+            @test DT.orient_predicate(p, q, r) === orient(p, q, r)
             @inferred DT.orient_predicate(p, q, r)
+            @test DT._orient_predicate(p, q, r) === orient(p, q, r)
+            @inferred DT._orient_predicate(p, q, r)
 
             p, q, r, s = eachcol(rand(3, 4))
-            @test DT.orient_predicate(p, q, r, s) == orient(p, q, r, s)
+            @test DT.orient_predicate(p, q, r, s) === orient(p, q, r, s)
             @inferred DT.orient_predicate(p, q, r, s)
+            @test DT._orient_predicate(p, q, r, s) === orient(p, q, r, s)
+            @inferred DT._orient_predicate(p, q, r, s)
 
             a, b, c, p = eachcol(rand(2, 4))
-            @test DT.incircle_predicate(a, b, c, p) == incircle(a, b, c, p)
+            @test DT.incircle_predicate(a, b, c, p) === incircle(a, b, c, p)
             @inferred DT.incircle_predicate(a, b, c, p)
+            @test DT._incircle_predicate(a, b, c, p) === incircle(a, b, c, p)
+            @inferred DT._incircle_predicate(a, b, c, p)
 
             a, b, p, q = eachcol(rand(2, 4))
-            @test DT.parallelorder_predicate(a, b, p, q) == parallelorder(a, b, p, q)
+            @test DT.parallelorder_predicate(a, b, p, q) === parallelorder(a, b, p, q)
             @inferred DT.parallelorder_predicate(a, b, p, q)
+            @test DT._parallelorder_predicate(a, b, p, q) === parallelorder(a, b, p, q)
+            @inferred DT._parallelorder_predicate(a, b, p, q)
 
             p, a, b = eachcol(rand(2, 3))
-            @test DT.sameside_predicate(a, b, p) == sameside(p, a, b)
+            @test DT.sameside_predicate(a, b, p) === sameside(p, a, b)
             @inferred DT.sameside_predicate(a, b, p)
 
             p, q, a, b = eachcol(rand(2, 4))
-            @test DT.meet_predicate(p, q, a, b) == meet(p, q, a, b)
+            @test DT.meet_predicate(p, q, a, b) === meet(p, q, a, b)
             @inferred DT.meet_predicate(p, q, a, b)
         end
     end
@@ -59,15 +101,18 @@ end
             a = (rand(x), rand(y))
             b = (rand(x), rand(y))
             c = (rand(x), rand(y))
-            @test DT.orient_predicate(p, q, r) == orient(p, q, r)
-            @test DT.incircle_predicate(a, b, c, p) == incircle(a, b, c, p)
+            @test DT.orient_predicate(p, q, r) === orient(p, q, r)
+            @test DT.incircle_predicate(a, b, c, p) === incircle(a, b, c, p)
             @test DT.sameside_predicate(a, b, p) == sameside(p, a, b)
-            @test DT.meet_predicate(p, q, a, b) == meet(p, q, a, b)
+            @test DT.meet_predicate(p, q, a, b) === meet(p, q, a, b)
+            @test DT._orient_predicate(p, q, r) === orient(p, q, r)
+            @test DT._incircle_predicate(a, b, c, p) === incircle(a, b, c, p)
             p = (rand(x), rand(y), rand(z))
             q = (rand(x), rand(y), rand(z))
             r = (rand(x), rand(y), rand(z))
             s = (rand(x), rand(y), rand(z))
             @test DT.orient_predicate(p, q, r, s) == orient(p, q, r, s)
+            @test DT._orient_predicate(p, q, r, s) == orient(p, q, r, s)
         end
     end
 
@@ -75,16 +120,23 @@ end
         p₁ = [5.7044025422189, 1.801603986463]
         p₂ = [8.3797127128527, 5.8924221838871]
         p₃ = [2.8875415689061, 6.2038339497809]
-        @test DT.orient_predicate(p₁, p₂, p₃) == 1
-        @test DT.orient_predicate(p₁, p₂, [10.0, 4.0]) == -1
+        @test DT.orient_predicate(p₁, p₂, p₃) === 1
+        @test DT.orient_predicate(p₁, p₂, [10.0, 4.0]) === -1
+        @test DT._orient_predicate(p₁, p₂, p₃) === 1
+        @test DT._orient_predicate(p₁, p₂, [10.0, 4.0]) === -1
         p₁ = [5.0, 1.0]
         p₂ = [5.0, 6.0]
-        @test DT.orient_predicate(p₁, p₂, [5.0, 5.0]) == 0
-        @test DT.orient_predicate(p₁, p₂, [5.0, 2.0]) == 0
-        @test DT.orient_predicate(p₂, p₁, [5.0, 2.0]) == 0
+        @test DT.orient_predicate(p₁, p₂, [5.0, 5.0]) === 0
+        @test DT.orient_predicate(p₁, p₂, [5.0, 2.0]) === 0
+        @test DT.orient_predicate(p₂, p₁, [5.0, 2.0]) === 0
+        @test DT._orient_predicate(p₁, p₂, [5.0, 5.0]) === 0
+        @test DT._orient_predicate(p₁, p₂, [5.0, 2.0]) === 0
+        @test DT._orient_predicate(p₂, p₁, [5.0, 2.0]) === 0
         p, q, r, s = (-4.16, -2.84, 0.0), (-4.65, -5.83, 0.0), (-1.12, -5.61, 0.0), (-7.83, 0.27, 4.6)
-        @test DT.orient_predicate(p, q, r, s) == -1
-        @test DT.orient_predicate(p, r, q, s) == 1
+        @test DT.orient_predicate(p, q, r, s) === -1
+        @test DT.orient_predicate(p, r, q, s) === 1
+        @test DT._orient_predicate(p, q, r, s) === -1
+        @test DT._orient_predicate(p, r, q, s) === 1
     end
 end
 
diff --git a/test/refinement/curve_bounded.jl b/test/refinement/curve_bounded.jl
index 927a1315f..9de320e17 100644
--- a/test/refinement/curve_bounded.jl
+++ b/test/refinement/curve_bounded.jl
@@ -10,7 +10,7 @@ using Test
 using StructEquality
 using DelimitedFiles
 using StableRNGs
-
+using Preferences
 const SACM = DT.SmallAngleComplexMember
 const SAC = DT.SmallAngleComplex
 
@@ -1450,6 +1450,9 @@ end
         point_sets = deepcopy.([points_I, points_II, points_III, points_IV, points_V, points_VI, points_VII, points_VIII, points_IX, points_X, points_XI, points_XII])
         curve_sets = deepcopy.([curve_I, curve_II, curve_III, curve_IV, curve_V, curve_VI, curve_VII, curve_VIII, curve_IX, curve_X, curve_XI, curve_XII])
         for i in eachindex(point_sets, curve_sets)
+            if !load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && i == 4 
+                continue 
+            end
             points, curve = deepcopy(point_sets[i]), deepcopy(curve_sets[i])
             tri = triangulate(points; boundary_nodes=curve, enrich=i ≤ 3)
             @test validate_triangulation(tri)
@@ -1517,6 +1520,9 @@ end
             (1e-1, 1e-2)
         ]
         for curve_idx in 4:lastindex(curve_sets)
+            if curve_idx ∈ (4, 6, 12) && !load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+                continue 
+            end
             for point_idx in 1:3
                 point_idx == 3 && curve_idx ≥ 5 && continue # no extra segments for curves ≥ 5
                 for (idx1, use_lens) in enumerate((false, true))
diff --git a/test/refinement/refine.jl b/test/refinement/refine.jl
index 989a6d22f..41f12e931 100644
--- a/test/refinement/refine.jl
+++ b/test/refinement/refine.jl
@@ -8,6 +8,7 @@ using ReferenceTests
 using Test
 using StructEquality
 using DelimitedFiles
+using Preferences
 using StableRNGs
 @struct_equal DT.InsertionEventHistory
 @struct_equal DT.RefinementConstraints
@@ -434,7 +435,9 @@ end
                         DT.undo_insertion!(tri, history)
                         @test tri.boundary_edge_map == orig_tri.boundary_edge_map
                         validate_statistics(tri)
-                        @test validate_triangulation(tri)
+                        if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+                            @test validate_triangulation(tri)
+                        end
                         @test tri == orig_tri
                         add_point!(tri, p, store_event_history=Val(true), event_history=history)
                         orig_tri = deepcopy(tri)
@@ -1186,7 +1189,7 @@ end
         tri_5 = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=false)
         tri_6 = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
         @static if VERSION ≥ v"1.10"
-        @inferred triangulate(rand(2, 250))
+            @inferred triangulate(rand(2, 250))
         end
         for (iii, tri) in enumerate((tri_1, tri_2, tri_3, tri_4, tri_5, tri_6))
             @info "Testing if encroached edges are detected. Run: $iii."
@@ -1707,92 +1710,94 @@ end
     end
 end
 
-@testset "enqueueing and splitting all encroached segments" begin
-    for iii in 1:100
-        for use_lens in (false, true)
-            for pass in 1:2
-                if pass == 1
-                    p1 = (0.0, 0.0)
-                    p2 = (1.0, 0.0)
-                    p3 = (0.0, 1.0)
-                    p4 = (1.0, 1.0)
-                    p5 = (0.5, 0.5)
-                    pts = [p1, p2, p3, p4, p5]
-                    C = Set{NTuple{2,Int}}()
-                    for i in 1:15
-                        θ = 2π * rand()
-                        r = 0.5sqrt(rand())
-                        x = 0.5 + r * cos(θ)
-                        y = 0.5 + r * sin(θ)
-                        push!(pts, (x, y))
-                        push!(C, (5, 5 + i))
-                    end
-                    tri = triangulate(pts; delete_ghosts=false, boundary_nodes=[1, 2, 4, 3, 1], segments=C)
-                    args = DT.RefinementArguments(tri, use_lens=use_lens, min_angle=25.0, seditious_angle=15.0, max_area=0.1)
-                    DT.enqueue_all_encroached_segments!(args, tri)
-                    @inferred DT.enqueue_all_encroached_segments!(args, tri)
-                    manual_enqueue = PriorityQueue{NTuple{2,Int},Float64}(Base.Order.Reverse)
-                    for e in each_edge(tri)
-                        if DT.contains_segment(tri, e...)
-                            flag = DT.is_encroached(tri, args, e)
-                            if flag
-                                manual_enqueue[e] = DT.edge_length_sqr(tri, e...)
+if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+    @testset "enqueueing and splitting all encroached segments" begin
+        for iii in 1:100
+            for use_lens in (false, true)
+                for pass in 1:2
+                    if pass == 1
+                        p1 = (0.0, 0.0)
+                        p2 = (1.0, 0.0)
+                        p3 = (0.0, 1.0)
+                        p4 = (1.0, 1.0)
+                        p5 = (0.5, 0.5)
+                        pts = [p1, p2, p3, p4, p5]
+                        C = Set{NTuple{2,Int}}()
+                        for i in 1:15
+                            θ = 2π * rand()
+                            r = 0.5sqrt(rand())
+                            x = 0.5 + r * cos(θ)
+                            y = 0.5 + r * sin(θ)
+                            push!(pts, (x, y))
+                            push!(C, (5, 5 + i))
+                        end
+                        tri = triangulate(pts; delete_ghosts=false, boundary_nodes=[1, 2, 4, 3, 1], segments=C)
+                        args = DT.RefinementArguments(tri, use_lens=use_lens, min_angle=25.0, seditious_angle=15.0, max_area=0.1)
+                        DT.enqueue_all_encroached_segments!(args, tri)
+                        @inferred DT.enqueue_all_encroached_segments!(args, tri)
+                        manual_enqueue = PriorityQueue{NTuple{2,Int},Float64}(Base.Order.Reverse)
+                        for e in each_edge(tri)
+                            if DT.contains_segment(tri, e...)
+                                flag = DT.is_encroached(tri, args, e)
+                                if flag
+                                    manual_enqueue[e] = DT.edge_length_sqr(tri, e...)
+                                end
                             end
                         end
-                    end
-                    compare_encroach_queues(args, manual_enqueue)
-                    empty!(args.events)
-                    DT.split_all_encroached_segments!(tri, args)
-                    @test isempty(args.queue.segments)
-                    @test !isempty(args.events.added_triangles)
-                    DT.enqueue_all_encroached_segments!(args, tri)
-                    @test isempty(args.queue.segments)
-                    if use_lens
-                        for e in each_segment(tri)
-                            @test !DT.is_encroached(tri, args, e)
+                        compare_encroach_queues(args, manual_enqueue)
+                        empty!(args.events)
+                        DT.split_all_encroached_segments!(tri, args)
+                        @test isempty(args.queue.segments)
+                        @test !isempty(args.events.added_triangles)
+                        DT.enqueue_all_encroached_segments!(args, tri)
+                        @test isempty(args.queue.segments)
+                        if use_lens
+                            for e in each_segment(tri)
+                                @test !DT.is_encroached(tri, args, e)
+                            end
+                        else
+                            @test is_conformal(tri)
                         end
+                        @test validate_triangulation(tri)
                     else
-                        @test is_conformal(tri)
-                    end
-                    @test validate_triangulation(tri)
-                else
-                    points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
-                    for r in LinRange(0.1, 0.9, 15)
-                        push!(points, (r, 0.01))
-                        push!(points, (r, 0.99))
-                        push!(points, (0.05, r))
-                        push!(points, (0.95, r))
-                    end
-                    for _ in 1:10
-                        push!(points, (rand(), rand()))
-                    end
-                    tri = triangulate(points; boundary_nodes=[1, 2, 3, 4, 1])
-                    args = DT.RefinementArguments(tri, use_lens=use_lens, min_angle=29.0, seditious_angle=19.0, max_area=0.05)
-                    DT.enqueue_all_encroached_segments!(args, tri)
-                    manual_enqueue = PriorityQueue{NTuple{2,Int},Float64}(Base.Order.Reverse)
-                    for e in each_edge(tri)
-                        if DT.contains_segment(tri, e...)
-                            flag = DT.is_encroached(tri, args, e)
-                            if flag
-                                manual_enqueue[e] = DT.edge_length_sqr(tri, e...)
+                        points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
+                        for r in LinRange(0.1, 0.9, 15)
+                            push!(points, (r, 0.01))
+                            push!(points, (r, 0.99))
+                            push!(points, (0.05, r))
+                            push!(points, (0.95, r))
+                        end
+                        for _ in 1:10
+                            push!(points, (rand(), rand()))
+                        end
+                        tri = triangulate(points; boundary_nodes=[1, 2, 3, 4, 1])
+                        args = DT.RefinementArguments(tri, use_lens=use_lens, min_angle=29.0, seditious_angle=19.0, max_area=0.05)
+                        DT.enqueue_all_encroached_segments!(args, tri)
+                        manual_enqueue = PriorityQueue{NTuple{2,Int},Float64}(Base.Order.Reverse)
+                        for e in each_edge(tri)
+                            if DT.contains_segment(tri, e...)
+                                flag = DT.is_encroached(tri, args, e)
+                                if flag
+                                    manual_enqueue[e] = DT.edge_length_sqr(tri, e...)
+                                end
                             end
                         end
-                    end
-                    compare_encroach_queues(args, manual_enqueue)
-                    empty!(args.events)
-                    DT.split_all_encroached_segments!(tri, args)
-                    @test isempty(args.queue.segments)
-                    @test !isempty(args.events.added_triangles)
-                    DT.enqueue_all_encroached_segments!(args, tri)
-                    @test isempty(args.queue.segments)
-                    if use_lens
-                        for e in each_segment(tri)
-                            @test !DT.is_encroached(tri, args, e)
+                        compare_encroach_queues(args, manual_enqueue)
+                        empty!(args.events)
+                        DT.split_all_encroached_segments!(tri, args)
+                        @test isempty(args.queue.segments)
+                        @test !isempty(args.events.added_triangles)
+                        DT.enqueue_all_encroached_segments!(args, tri)
+                        @test isempty(args.queue.segments)
+                        if use_lens
+                            for e in each_segment(tri)
+                                @test !DT.is_encroached(tri, args, e)
+                            end
+                        else
+                            @test is_conformal(tri)
                         end
-                    else
-                        @test is_conformal(tri)
+                        @test validate_triangulation(tri)
                     end
-                    @test validate_triangulation(tri)
                 end
             end
         end
@@ -2272,10 +2277,10 @@ end
                             @test validate_refinement(tri, args)
                             if _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(1, 3, 1, 2, 2)
                                 fig, ax, sc = triplot(tri)
-                                @test_reference "a_simple_convex_example_circle.png" fig by=psnr_equality(15)
+                                @test_reference "a_simple_convex_example_circle.png" fig by = psnr_equality(15)
                             elseif _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(2, 3, 1, 2, 2)
                                 fig, ax, sc = triplot(tri)
-                                @test_reference "a_simple_convex_example_lens.png" fig by=psnr_equality(15)
+                                @test_reference "a_simple_convex_example_lens.png" fig by = psnr_equality(15)
                             end
                         end
                     end
@@ -2300,10 +2305,10 @@ end
                             @test validate_refinement(tri, DT.RefinementArguments(tri; min_angle, min_area, seditious_angle, max_area, use_circumcenter=true, use_lens))
                             if _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(1, 3, 1, 3, 1)
                                 fig, ax, sc = triplot(tri)
-                                @test_reference "triangulation_with_a_hole_circle.png" fig by=psnr_equality(15)
+                                @test_reference "triangulation_with_a_hole_circle.png" fig by = psnr_equality(15)
                             elseif _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(2, 3, 1, 3, 1)
                                 fig, ax, sc = triplot(tri)
-                                @test_reference "triangulation_with_a_hole_lens.png" fig by=psnr_equality(15)
+                                @test_reference "triangulation_with_a_hole_lens.png" fig by = psnr_equality(15)
                             end
                         end
                     end
@@ -2390,10 +2395,10 @@ end
                             @test validate_refinement(tri; min_angle, min_area, max_area, use_circumcenter=true, seditious_angle, use_lens, check_conformal=false)
                             if _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(1, 3, 1, 3, 1)
                                 fig, ax, sc = triplot(tri)
-                                @test_reference "a_constrained_triangulation_with_multiple_holes_circle.png" fig by=psnr_equality(15)
+                                @test_reference "a_constrained_triangulation_with_multiple_holes_circle.png" fig by = psnr_equality(15)
                             elseif _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(2, 3, 1, 3, 1)
                                 fig, ax, sc = triplot(tri)
-                                @test_reference "a_constrained_triangulation_with_multiple_holes_lens.png" fig by=psnr_equality(15)
+                                @test_reference "a_constrained_triangulation_with_multiple_holes_lens.png" fig by = psnr_equality(15)
                             end
                         end
                     end
@@ -2428,10 +2433,10 @@ end
                             @test validate_refinement(tri; min_angle, min_area, max_area, seditious_angle, use_circumcenter=true, use_lens)
                             if _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(1, 4, 1, 3, 2)
                                 fig, ax, sc = triplot(tri)
-                                @test_reference "another_square_example_circle.png" fig by=psnr_equality(15)
+                                @test_reference "another_square_example_circle.png" fig by = psnr_equality(15)
                             elseif _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(2, 4, 1, 3, 2)
                                 fig, ax, sc = triplot(tri)
-                                @test_reference "another_square_example_lens.png" fig by=psnr_equality(15)
+                                @test_reference "another_square_example_lens.png" fig by = psnr_equality(15)
                             end
                         end
                     end
@@ -2475,10 +2480,10 @@ end
                             validate_statistics(tri)
                             if _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(1, 3, 1, 2, 2)
                                 fig, ax, sc = triplot(tri)
-                                @test_reference "avoid_bouncing_concentric_circular_shells_circle.png" fig by=psnr_equality(15)
+                                @test_reference "avoid_bouncing_concentric_circular_shells_circle.png" fig by = psnr_equality(15)
                             elseif _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(2, 3, 1, 2, 2)
                                 fig, ax, sc = triplot(tri)
-                                @test_reference "avoid_bouncing_concentric_circular_shells_lens.png" fig by=psnr_equality(15)
+                                @test_reference "avoid_bouncing_concentric_circular_shells_lens.png" fig by = psnr_equality(15)
                             end
                         end
                     end
@@ -2492,14 +2497,14 @@ end
             points = [(0.0, 0.0), (1.0, 0.0), (0.8, 0.2)]
             tri = triangulate(points)
             refine!(tri, use_circumcenter=true)
-            @test_reference "seditious_edge_testing_1.png" triplot(tri) by=psnr_equality(15)
+            @test_reference "seditious_edge_testing_1.png" triplot(tri) by = psnr_equality(15)
             @test validate_refinement(tri, use_circumcenter=true)
             validate_statistics(tri)
 
             points = [(0.0, 0.0), (1.0, 0.0), (0.8, 0.2)]
             tri = triangulate(points)
             refine!(tri, use_circumcenter=true, use_lens=false)
-            @test_reference "seditious_edge_testing_2.png" triplot(tri) by=psnr_equality(15)
+            @test_reference "seditious_edge_testing_2.png" triplot(tri) by = psnr_equality(15)
             @test validate_refinement(tri, use_circumcenter=true, use_lens=false)
             validate_statistics(tri)
 
@@ -2508,7 +2513,7 @@ end
             args = DT.RefinementArguments(tri, use_circumcenter=true, use_lens=false, max_area=1e-3get_area(tri))
             refine!(tri, args)
             @test validate_refinement(tri, args)
-            @test_reference "seditious_edge_testing_3.png" triplot(tri) by=psnr_equality(15)
+            @test_reference "seditious_edge_testing_3.png" triplot(tri) by = psnr_equality(15)
             validate_statistics(tri)
         end
     end
@@ -2544,7 +2549,7 @@ end
         @test validate_triangulation(tri)
         @test validate_refinement(tri; max_area=0.001, max_points=5000, use_circumcenter=true, use_lens=true)
         fig, ax, sc = triplot(tri)
-        @test_reference "triangulating_with_an_interior_hole.png" fig by=psnr_equality(15)
+        @test_reference "triangulating_with_an_interior_hole.png" fig by = psnr_equality(15)
     end
 
     @testset "Refining disjoint sets" begin
@@ -2568,50 +2573,52 @@ end
         validate_statistics(tri)
         @test validate_refinement(tri; min_area, max_area, use_circumcenter=true)
         fig, ax, sc = triplot(tri)
-        @test_reference "refining_disjoint_sets.png" fig by=psnr_equality(15)
+        @test_reference "refining_disjoint_sets.png" fig by = psnr_equality(15)
     end
 
-    @testset "Small angles" begin
-        ps = 0
-        fig = Figure(fontsize=52)
-        for i in 1:12
-            if i > 6
-                use_lens = true
-                i -= 6
-                ax = Axis(fig[2, i], title="Lens; $i", width=600, height=600)
-            else
-                use_lens = false
-                ax = Axis(fig[1, i], title="Circle; $i", width=600, height=600)
-            end
-            hidedecorations!(ax)
-            hidespines!(ax)
-            rng = StableRNG(i)
-            p1 = (0.0, 0.0)
-            p2 = (1.0, 0.0)
-            p3 = (0.0, 1.0)
-            p4 = (1.0, 1.0)
-            p5 = (0.5, 0.5)
-            pts = [p1, p2, p3, p4, p5]
-            C = Set{NTuple{2,Int}}()
-            for i in 1:20
-                θ = 2π * rand(rng)
-                r = 0.5sqrt(rand(rng))
-                x = 0.5 + r * cos(θ)
-                y = 0.5 + r * sin(θ)
-                push!(pts, (x, y))
-                push!(C, (5, 5 + i))
+    if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+        @testset "Small angles" begin
+            ps = 0
+            fig = Figure(fontsize=52)
+            for i in 1:12
+                if i > 6
+                    use_lens = true
+                    i -= 6
+                    ax = Axis(fig[2, i], title="Lens; $i", width=600, height=600)
+                else
+                    use_lens = false
+                    ax = Axis(fig[1, i], title="Circle; $i", width=600, height=600)
+                end
+                hidedecorations!(ax)
+                hidespines!(ax)
+                rng = StableRNG(i)
+                p1 = (0.0, 0.0)
+                p2 = (1.0, 0.0)
+                p3 = (0.0, 1.0)
+                p4 = (1.0, 1.0)
+                p5 = (0.5, 0.5)
+                pts = [p1, p2, p3, p4, p5]
+                C = Set{NTuple{2,Int}}()
+                for i in 1:20
+                    θ = 2π * rand(rng)
+                    r = 0.5sqrt(rand(rng))
+                    x = 0.5 + r * cos(θ)
+                    y = 0.5 + r * sin(θ)
+                    push!(pts, (x, y))
+                    push!(C, (5, 5 + i))
+                end
+                tri = triangulate(pts; rng, boundary_nodes=[1, 2, 4, 3, 1], segments=C)
+                refine!(tri; min_angle=27.3, min_area=0.0, use_circumcenter=true, rng, use_lens)
+                stats = statistics(tri)
+                ps += DT.get_largest_angle(stats) ≤ max(π - 2 * deg2rad(17.0), 2asin((sqrt(3) - 1) / 2)) # Corollary 8 of "When and Why Ruppert's Algorithm Works. In Twelfth International Meshing Roundtable, pp. 91–102, Santa Fe, NM, Sept 2003."
+                validate_statistics(tri)
+                @test validate_refinement(tri; min_angle=27.3, min_area=0.0, use_circumcenter=true, warn=false, use_lens, rng)
+                triplot!(ax, tri)
             end
-            tri = triangulate(pts; rng, boundary_nodes=[1, 2, 4, 3, 1], segments=C)
-            refine!(tri; min_angle=27.3, min_area=0.0, use_circumcenter=true, rng, use_lens)
-            stats = statistics(tri)
-            ps += DT.get_largest_angle(stats) ≤ max(π - 2 * deg2rad(17.0), 2asin((sqrt(3) - 1) / 2)) # Corollary 8 of "When and Why Ruppert's Algorithm Works. In Twelfth International Meshing Roundtable, pp. 91–102, Santa Fe, NM, Sept 2003."
-            validate_statistics(tri)
-            @test validate_refinement(tri; min_angle=27.3, min_area=0.0, use_circumcenter=true, warn=false, use_lens, rng)
-            triplot!(ax, tri)
+            resize_to_layout!(fig)
+            @test_reference "refinement_with_small_angles.png" fig
+            @test ps == 12
         end
-        resize_to_layout!(fig)
-        @test_reference "refinement_with_small_angles.png" fig
-        @test ps == 12
     end
 
     @testset "Another disjoint domain" begin
@@ -2643,7 +2650,7 @@ end
         validate_statistics(tri)
         @test validate_refinement(tri, max_area=1e-3A, min_area=0.0, use_circumcenter=true, use_lens=true, warn=false)
         fig, ax, sc = triplot(tri)
-        @test_reference "refining_disjoint_sets_2.png" fig by=psnr_equality(15)
+        @test_reference "refining_disjoint_sets_2.png" fig by = psnr_equality(15)
     end
 
     @testset "Tight example with a single triangle boundary interior" begin
@@ -2695,7 +2702,7 @@ end
         validate_statistics(tri)
         @test validate_refinement(tri; max_area=1e-3A, rng, use_circumcenter=true, warn=true)
         fig, ax, sc = triplot(tri)
-        @test_reference "tight_example_with_a_single_triangle_boundary_interior.png" fig by=psnr_equality(15)
+        @test_reference "tight_example_with_a_single_triangle_boundary_interior.png" fig by = psnr_equality(15)
     end
 
     @testset "A complicated example with tight walls and small angles" begin
@@ -2767,27 +2774,29 @@ end
         refine!(tri; max_area=1.0, rng, use_circumcenter=true)
         validate_statistics(tri)
         fig, ax, sc = triplot(tri)
-        @test_reference "complicated_example_with_tight_walls_and_small_angles.png" fig by=psnr_equality(15)
+        @test_reference "complicated_example_with_tight_walls_and_small_angles.png" fig by = psnr_equality(15)
     end
 
-    if !(get(ENV, "CI", "false") == "true")
-        @testset "Tasmania" begin
-            rng = StableRNG(123)
-            tassy_path = joinpath(dirname(dirname(pathof(DelaunayTriangulation))), "test", "tassy.txt")
-            tassy = readdlm(tassy_path)
-            ymax = @views maximum(tassy[:, 2])
-            tassy = [(x, ymax - y) for (x, y) in eachrow(tassy)]
-            reverse!(tassy)
-            unique!(tassy)
-            push!(tassy, tassy[begin])
-            boundary_nodes, points = convert_boundary_points_to_indices(tassy)
-            tri = triangulate(points; rng, boundary_nodes=boundary_nodes)
-            A = get_area(tri)
-            refine!(tri; max_area=1e-2A, use_circumcenter=true, rng)
-            validate_statistics(tri)
-            @test validate_refinement(tri; max_area=1e-2A, rng, use_circumcenter=true, warn=false)
-            fig, ax, sc = triplot(tri)
-            @test_reference "tasmania.png" fig by=psnr_equality(15)
+    if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+        if !(get(ENV, "CI", "false") == "true")
+            @testset "Tasmania" begin
+                rng = StableRNG(123)
+                tassy_path = joinpath(dirname(dirname(pathof(DelaunayTriangulation))), "test", "tassy.txt")
+                tassy = readdlm(tassy_path)
+                ymax = @views maximum(tassy[:, 2])
+                tassy = [(x, ymax - y) for (x, y) in eachrow(tassy)]
+                reverse!(tassy)
+                unique!(tassy)
+                push!(tassy, tassy[begin])
+                boundary_nodes, points = convert_boundary_points_to_indices(tassy)
+                tri = triangulate(points; rng, boundary_nodes=boundary_nodes)
+                A = get_area(tri)
+                refine!(tri; max_area=1e-2A, use_circumcenter=true, rng)
+                validate_statistics(tri)
+                @test validate_refinement(tri; max_area=1e-2A, rng, use_circumcenter=true, warn=false)
+                fig, ax, sc = triplot(tri)
+                @test_reference "tasmania.png" fig by = psnr_equality(15)
+            end
         end
     end
 
@@ -3051,6 +3060,6 @@ end
         validate_statistics(tri)
         @test validate_refinement(tri; min_angle=26.45, max_area=0.005A / 9, rng, use_circumcenter=true)
         fig, ax, sc = triplot(tri)
-        @test_reference "julia_logo.png" fig by=psnr_equality(15)
+        @test_reference "julia_logo.png" fig by = psnr_equality(15)
     end
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index bc52eb882..5e1c98102 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,3 +1,14 @@
+# setup LocalPreferences.toml 
+using Preferences 
+USE_EXACTPREDICATES = get(ENV, "USE_EXACTPREDICATES", "default")
+if USE_EXACTPREDICATES == "true"
+    set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => true)
+elseif USE_EXACTPREDICATES == "false" 
+    set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => false)
+end # if USE_EXACTPREDICATES == "default", do nothing
+
+@info "Testing with USE_EXACTPREDICATES = $USE_EXACTPREDICATES"
+
 # get all the compilation out of the way
 using BenchmarkTools
 using CairoMakie
@@ -41,7 +52,7 @@ function safe_include(filename; name=filename, push=true, verbose = true) # Work
     push && push!(ALL_TEST_SCRIPTS, normpath(filename))
     mod = @eval module $(gensym()) end
     @info "[$(ct())] Testing $name"
-    @testset verbose = verbose "$name" begin
+    @testset verbose = verbose "$(basename(name))" begin
         @eval mod using ..HelperFunctions
         @eval mod using ..Test
         Base.include(mod, filename)
@@ -50,7 +61,7 @@ end
 
 @testset verbose = true "DelaunayTriangulation.jl" begin
     @testset verbose = true "Aqua" begin
-        Aqua.test_all(DelaunayTriangulation; ambiguities=false, project_extras=false) # don't care about julia < 1.2
+        Aqua.test_all(DelaunayTriangulation; ambiguities=false, project_extras=false, stale_deps = load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)) # don't care about julia < 1.2
         Aqua.test_ambiguities(DelaunayTriangulation) # don't pick up Base and Core...
     end    
 
diff --git a/test/triangulation/bowyer_watson.jl b/test/triangulation/bowyer_watson.jl
index 888109f08..871dee261 100644
--- a/test/triangulation/bowyer_watson.jl
+++ b/test/triangulation/bowyer_watson.jl
@@ -4,8 +4,7 @@ using Random
 using StableRNGs
 using CairoMakie
 using DataStructures
-
-
+using Preferences
 
 @testset "Getting the correct order" begin
     points = rand(2, 50)
@@ -98,12 +97,14 @@ end
     @test validate_triangulation(tri)
 end
 
-@testset "Lots of collinearity" begin
-    _tri = triangulate_rectangle(-3.0, 2.0, 5.0, 17.3, 23, 57; single_boundary=true)
-    @test validate_triangulation(_tri)
-    for _ in 1:100
-        tri = DT.triangulate(_tri.points)
-        @test validate_triangulation(tri)
+if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+    @testset "Lots of collinearity" begin
+        _tri = triangulate_rectangle(-3.0, 2.0, 5.0, 17.3, 23, 57; single_boundary=true)
+        @test validate_triangulation(_tri)
+        for _ in 1:100
+            tri = DT.triangulate(_tri.points)
+            @test validate_triangulation(tri)
+        end
     end
 end
 
diff --git a/test/triangulation/constrained.jl b/test/triangulation/constrained.jl
index 6daee652e..53ed60f12 100644
--- a/test/triangulation/constrained.jl
+++ b/test/triangulation/constrained.jl
@@ -3,7 +3,7 @@ const DT = DelaunayTriangulation
 using CairoMakie
 using Random
 using StableRNGs
-
+using Preferences
 
 @testset "Random constrained Delaunay triangulations" begin
     for i in 1:25
@@ -91,81 +91,83 @@ end
     end
 end
 
-@testset "Lattice" begin
-    for m in 1:2
-        @info "Testing dense lattice. Run: $m"
-        rng = StableRNG(m)
-        a = 0.0
-        b = 5.0
-        c = -3.0
-        d = 7.0
-        nx = 13
-        ny = 20
-        tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
-        add_segment!(tri, 56, 162; rng)
-        for e in [(1, 249), (1, 250), (1, 251), (1, 26), (1, 39), (1, 52)]
-            add_segment!(tri, e; rng)
-        end
-        add_segment!(tri, 190, 99; rng)
-        for e in [(99, 113), (113, 101), (101, 115)]
-            add_segment!(tri, e; rng)
-        end
-        @test validate_triangulation(tri)
+if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+    @testset "Lattice" begin
+        for m in 1:2
+            @info "Testing dense lattice. Run: $m"
+            rng = StableRNG(m)
+            a = 0.0
+            b = 5.0
+            c = -3.0
+            d = 7.0
+            nx = 13
+            ny = 20
+            tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
+            add_segment!(tri, 56, 162; rng)
+            for e in [(1, 249), (1, 250), (1, 251), (1, 26), (1, 39), (1, 52)]
+                add_segment!(tri, e; rng)
+            end
+            add_segment!(tri, 190, 99; rng)
+            for e in [(99, 113), (113, 101), (101, 115)]
+                add_segment!(tri, e; rng)
+            end
+            @test validate_triangulation(tri)
 
-        a = -0.1
-        b = 0.1
-        c = -0.01
-        d = 0.01
-        nx = 25
-        ny = 25
-        tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
-        tri = triangulate(get_points(tri))
-        for i in 2:24
-            add_segment!(tri, i, 600 + i; rng)
-        end
-        @test validate_triangulation(tri)
-        tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
-        tri = triangulate(get_points(tri); rng)
-        for e in [(1, 28), (28, 103), (103, 180), (180, 625), (625, 523)]
-            add_segment!(tri, e; rng)
-        end
-        for e in [(437, 614), (527, 602), (528, 603), (555, 605)]
-            add_segment!(tri, e; rng)
+            a = -0.1
+            b = 0.1
+            c = -0.01
+            d = 0.01
+            nx = 25
+            ny = 25
+            tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
+            tri = triangulate(get_points(tri))
+            for i in 2:24
+                add_segment!(tri, i, 600 + i; rng)
+            end
+            @test validate_triangulation(tri)
+            tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
+            tri = triangulate(get_points(tri); rng)
+            for e in [(1, 28), (28, 103), (103, 180), (180, 625), (625, 523)]
+                add_segment!(tri, e; rng)
+            end
+            for e in [(437, 614), (527, 602), (528, 603), (555, 605)]
+                add_segment!(tri, e; rng)
+            end
+            @test validate_triangulation(tri)
         end
-        @test validate_triangulation(tri)
-    end
-    for m in 1:20
-        rng = StableRNG(m)
-        @info "Testing coarse lattice. Run: $m"
-        a = 0.0
-        b = 1.0
-        c = 0.0
-        d = 1.0
-        nx = 2
-        ny = 2
-        tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
-        add_segment!(tri, 1, 4; rng)
-        @test validate_triangulation(tri)
+        for m in 1:20
+            rng = StableRNG(m)
+            @info "Testing coarse lattice. Run: $m"
+            a = 0.0
+            b = 1.0
+            c = 0.0
+            d = 1.0
+            nx = 2
+            ny = 2
+            tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
+            add_segment!(tri, 1, 4; rng)
+            @test validate_triangulation(tri)
 
-        a = 0
-        b = 1
-        c = 0
-        d = 5
-        nx = 25
-        ny = 3
-        tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
-        tri = triangulate(get_points(tri); rng)
-        for i in 1:(nx-1)
-            u = i
-            v = 2nx
-            add_segment!(tri, u, v)
-        end
-        for i in 51:75
-            u = i
-            v = 26
-            add_segment!(tri, u, v; rng)
+            a = 0
+            b = 1
+            c = 0
+            d = 5
+            nx = 25
+            ny = 3
+            tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true)
+            tri = triangulate(get_points(tri); rng)
+            for i in 1:(nx-1)
+                u = i
+                v = 2nx
+                add_segment!(tri, u, v)
+            end
+            for i in 51:75
+                u = i
+                v = 26
+                add_segment!(tri, u, v; rng)
+            end
+            @test validate_triangulation(tri)
         end
-        @test validate_triangulation(tri)
     end
 end
 
@@ -185,34 +187,36 @@ end
     end
 end
 
-@testset "Triangulation with two curves" begin
-    for i in 1:10
-        @info "Testing triangulation of a domain with two curves: Run $i"
-        rng = StableRNG(i)
-        pts = [(rand(rng), rand(rng)) for _ in 1:50]
-        x1 = LinRange(0, 1, 100)
-        y1 = LinRange(0, 0, 100)
-        x2 = LinRange(1, 1, 100)
-        y2 = LinRange(0, 1, 100)
-        x3 = LinRange(1, 0, 100)
-        y3 = LinRange(1, 1, 100)
-        x4 = LinRange(0, 0, 100)
-        y4 = LinRange(1, 0, 100)
-        x = [x1..., x2..., x3..., x4...]
-        y = [y1..., y2..., y3..., y4...]
-        outer_square = [(x, y) for (x, y) in zip(x, y)] |> unique
-        push!(outer_square, outer_square[begin])
-        outer_square_x = [first.(outer_square)]
-        outer_square_y = [last.(outer_square)]
-        circ_pts = [(0.3cos(θ), 0.3sin(θ)) .+ 0.5 for θ in LinRange(2π, 0, 50)]
-        circ_pts[end] = circ_pts[1]
-        inner_circle_x = [first.(circ_pts)]
-        inner_circle_y = [last.(circ_pts)]
-        x = [outer_square_x, inner_circle_x]
-        y = [outer_square_y, inner_circle_y]
-        nodes, pts = convert_boundary_points_to_indices(x, y; existing_points=pts)
-        tri = triangulate(pts; boundary_nodes=nodes, rng)
-        @test validate_triangulation(tri)
+if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+    @testset "Triangulation with two curves" begin
+        for i in 1:10
+            @info "Testing triangulation of a domain with two curves: Run $i"
+            rng = StableRNG(i)
+            pts = [(rand(rng), rand(rng)) for _ in 1:50]
+            x1 = LinRange(0, 1, 100)
+            y1 = LinRange(0, 0, 100)
+            x2 = LinRange(1, 1, 100)
+            y2 = LinRange(0, 1, 100)
+            x3 = LinRange(1, 0, 100)
+            y3 = LinRange(1, 1, 100)
+            x4 = LinRange(0, 0, 100)
+            y4 = LinRange(1, 0, 100)
+            x = [x1..., x2..., x3..., x4...]
+            y = [y1..., y2..., y3..., y4...]
+            outer_square = [(x, y) for (x, y) in zip(x, y)] |> unique
+            push!(outer_square, outer_square[begin])
+            outer_square_x = [first.(outer_square)]
+            outer_square_y = [last.(outer_square)]
+            circ_pts = [(0.3cos(θ), 0.3sin(θ)) .+ 0.5 for θ in LinRange(2π, 0, 50)]
+            circ_pts[end] = circ_pts[1]
+            inner_circle_x = [first.(circ_pts)]
+            inner_circle_y = [last.(circ_pts)]
+            x = [outer_square_x, inner_circle_x]
+            y = [outer_square_y, inner_circle_y]
+            nodes, pts = convert_boundary_points_to_indices(x, y; existing_points=pts)
+            tri = triangulate(pts; boundary_nodes=nodes, rng)
+            @test validate_triangulation(tri)
+        end
     end
 end
 
diff --git a/test/utils.jl b/test/utils.jl
index b57e1cbd9..f837c929b 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -1224,17 +1224,19 @@ end
 end
 
 @testset "replace_(boundary/ghost)_triangle_with_(ghost/boundary)_triangle" begin
-      tri = triangulate(rand(2, 5000), delete_ghosts=false)
-      _tri = triangulate_rectangle(0.0, 1.0, 0.0, 1.0, 25, 25, delete_ghosts=false)
-      for tri in (tri, _tri)
-            for T in each_ghost_triangle(tri)
-                  T = DT.sort_triangle(T)
-                  i, j, k = triangle_vertices(T)
-                  w = get_adjacent(tri, j, i)
-                  V = (j, i, w)
-                  @test DT.replace_boundary_triangle_with_ghost_triangle(tri, V) == T
-                  @test DT.replace_ghost_triangle_with_boundary_triangle(tri, T) == V
-                  if (i, j) ∉ ((624, 625), (1, 26), (2, 1), (625, 600)) # corners of the lattice have two associated ghosts
+      for _ in 1:10
+            tri = triangulate(rand(2, 5000), delete_ghosts=false)
+            _tri = triangulate_rectangle(0.0, 1.0, 0.0, 1.0, 25, 25, delete_ghosts=false)
+            for tri in (tri, _tri)
+                  for T in each_ghost_triangle(tri)
+                        T = DT.sort_triangle(T)
+                        i, j, k = triangle_vertices(T)
+                        w = get_adjacent(tri, j, i)
+                        V = (j, i, w)
+                        kr = count(e -> DT.is_boundary_edge(tri, e...), DT.triangle_edges(V))
+                        kr == 2 && continue # two associated ghosts, test is not clear
+                        @test DT.replace_boundary_triangle_with_ghost_triangle(tri, V) == T
+                        @test DT.replace_ghost_triangle_with_boundary_triangle(tri, T) == V
                         V = (i, w, j)
                         @test DT.replace_boundary_triangle_with_ghost_triangle(tri, V) == T
                         @test DT.replace_ghost_triangle_with_boundary_triangle(tri, T) == (j, i, w)
diff --git a/test/voronoi/voronoi.jl b/test/voronoi/voronoi.jl
index a01956611..1639af91c 100644
--- a/test/voronoi/voronoi.jl
+++ b/test/voronoi/voronoi.jl
@@ -8,6 +8,7 @@ import GeometryBasics: Point2f
 using StaticArrays
 using LinearAlgebra
 using StructEquality
+using Preferences
 @struct_equal DT.Queue
 
 @testset "Unconstrained test" begin
@@ -194,13 +195,13 @@ end
 end
 
 @testset "Voronoi point location" begin
-    A = (-1.0, 7.0)
-    B = (4.0, 4.0)
-    C = (-2.0, -1.0)
-    D = (-1.0, 3.0)
-    E = (3.0, -1.0)
-    F = (1.0, 4.0)
-    G = (-3.0, 5.0)
+    A = (-1.0, 7.0) .+ (1e-6rand(), 1e-6rand()) # perturb to allow the tests to work even without ExactPredicates
+    B = (4.0, 4.0) .+ (1e-6rand(), 1e-6rand())
+    C = (-2.0, -1.0) .+ (1e-6rand(), 1e-6rand())
+    D = (-1.0, 3.0) .+ (1e-6rand(), 1e-6rand())
+    E = (3.0, -1.0) .+ (1e-6rand(), 1e-6rand())
+    F = (1.0, 4.0) .+ (1e-6rand(), 1e-6rand())
+    G = (-3.0, 5.0) .+ (1e-6rand(), 1e-6rand())
     pts = [A, B, C, D, E, F, G]
     tri = triangulate(pts; delete_ghosts=false, randomise=false)
     vor = voronoi(tri)
@@ -220,31 +221,33 @@ end
         @test findmin(all_dists)[2] == u
     end
 
-    points = [
-        (0.0, 0.0), (-1.0, 1.0), (-0.5, 1.0), (0.0, 1.0), (0.5, 1.0), (1.0, 1.0),
-        (1.0, 0.8), (1.0, 0.0), (1.0, -0.5), (1.0, -1.0),
-        (0.1, -1.0), (-0.8, -1.0), (-1.0, -1.0),
-        (-1.0, -0.7), (-1.0, -0.1), (-1.0, 0.6),
-        (-0.1, -0.8), (0.2, -0.8),
-        (-0.6, -0.4), (0.9, 0.0), (-0.5, 0.5), (-0.4, 0.6), (-0.1, 0.8)
-    ]
-    tri = triangulate(points, delete_ghosts=false)
-    vorn = voronoi(tri)
-    @test validate_tessellation(vorn)
-    xg = LinRange(-1, 1, 250)
-    yg = LinRange(-1, 1, 250)
-    x = vec([x for x in xg, _ in yg])
-    y = vec([y for _ in xg, y in yg])
-    for (ξ, η) in zip(x, y)
-        p = (ξ, η)
-        u = get_nearest_neighbour(vorn, p)
-        @inferred get_nearest_neighbour(vorn, p)
-        all_dists = [norm(p .- get_generator(vorn, i)) for i in sort(collect(each_generator(vorn)))]
-        k = findmin(all_dists)[2]
-        @test k == u
-        for m in DT.each_point_index(tri)
-            u = get_nearest_neighbour(vorn, p, try_points=m)
-            @test u == k
+    if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+        points = [
+            (0.0, 0.0), (-1.0, 1.0), (-0.5, 1.0), (0.0, 1.0), (0.5, 1.0), (1.0, 1.0),
+            (1.0, 0.8), (1.0, 0.0), (1.0, -0.5), (1.0, -1.0),
+            (0.1, -1.0), (-0.8, -1.0), (-1.0, -1.0),
+            (-1.0, -0.7), (-1.0, -0.1), (-1.0, 0.6),
+            (-0.1, -0.8), (0.2, -0.8),
+            (-0.6, -0.4), (0.9, 0.0), (-0.5, 0.5), (-0.4, 0.6), (-0.1, 0.8)
+        ]
+        tri = triangulate(points, delete_ghosts=false)
+        vorn = voronoi(tri)
+        @test validate_tessellation(vorn)
+        xg = LinRange(-1, 1, 250)
+        yg = LinRange(-1, 1, 250)
+        x = vec([x for x in xg, _ in yg])
+        y = vec([y for _ in xg, y in yg])
+        for (ξ, η) in zip(x, y)
+            p = (ξ, η)
+            u = get_nearest_neighbour(vorn, p)
+            @inferred get_nearest_neighbour(vorn, p)
+            all_dists = [norm(p .- get_generator(vorn, i)) for i in sort(collect(each_generator(vorn)))]
+            k = findmin(all_dists)[2]
+            @test k == u
+            for m in DT.each_point_index(tri)
+                u = get_nearest_neighbour(vorn, p, try_points=m)
+                @test u == k
+            end
         end
     end
 end
@@ -893,14 +896,16 @@ end
     @test flag / tot > 0.9
 end
 
-@testset "Lattice" begin
-    for _ in 1:100
-        tri = triangulate_rectangle(0, 1, 0, 1, 11, 11, delete_ghosts=false)
-        tri = triangulate(tri.points)
-        vorn = voronoi(tri)
-        @test validate_tessellation(vorn; check_convex=false)
-        vorn = voronoi(tri, clip=true)
-        @test validate_tessellation(vorn; check_convex=false, check_adjacent=false)
+if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)
+    @testset "Lattice" begin
+        for _ in 1:100
+            tri = triangulate_rectangle(0, 1, 0, 1, 11, 11, delete_ghosts=false)
+            tri = triangulate(tri.points)
+            vorn = voronoi(tri)
+            @test validate_tessellation(vorn; check_convex=false)
+            vorn = voronoi(tri, clip=true)
+            @test validate_tessellation(vorn; check_convex=false, check_adjacent=false)
+        end
     end
 end