diff --git a/README.md b/README.md
index f2e83baa..d3638edc 100644
--- a/README.md
+++ b/README.md
@@ -321,7 +321,7 @@ julia> ustrip(ua"degC", 295.15u"K")
 22.0
 ```
 
-### Arrays
+### [Arrays](@id home_arrays)
 
 For working with an array of quantities that have the same dimensions,
 you can use a `QuantityArray`:
@@ -364,6 +364,22 @@ julia> @btime $f.(qa) setup=(xa = randn(100000) .* u"km/s"; qa = QuantityArray(x
 
 So we can see the `QuantityArray` version saves on both time and memory.
 
+By default, DynamicQuantities will create a `QuantityArray` from an `AbstractArray`, similarly to how a `Quantity` is created from a scalar in the [Usage](@ref) examples:
+
+```@repl quantity-array
+using DynamicQuantities # hide
+
+x = [0.3, 0.4, 0.5]u"km/s"
+
+y = (42:45) * u"kg"
+```
+
+This can be overridden to produce a vector of `Quantity`s by explicitly broadcasting the unit:
+
+```@repl quantity-array
+z = [0.3, 0.4, 0.5] .* u"km/s"
+```
+
 ### Unitful
 
 DynamicQuantities allows you to convert back and forth from Unitful.jl:
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 47a056ee..619742f5 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -246,6 +246,7 @@ the same dimension) by passing an array and a single quantity:
 
 ```julia
 x = QuantityArray(randn(32), u"km/s")
+# or x = randn(32)u"km/s"
 ```
 
 or, by passing an array of individual quantities:
@@ -281,6 +282,8 @@ f_square(v) = v^2 * 1.5 - v^2
 println("Applying function to y_q: ", sum(f_square.(y_q)))
 ```
 
+See [Home > Arrays](@ref home_arrays) for more.
+
 ### Fill
 
 We can also make `QuantityArray` using `fill`:
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 00000000..c650f3d4
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,482 @@
+
+
+![](https://github.com/SymbolicML/DynamicQuantities.jl/assets/7593028/10b1e6b8-f1c5-43bc-97e6-4ddb4c175293)
+
+
+[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://symbolicml.org/DynamicQuantities.jl/dev/)
+[![Build Status](https://github.com/SymbolicML/DynamicQuantities.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/SymbolicML/DynamicQuantities.jl/actions/workflows/CI.yml?query=branch%3Amain)
+[![Coverage](https://coveralls.io/repos/github/SymbolicML/DynamicQuantities.jl/badge.svg?branch=main)](https://coveralls.io/github/SymbolicML/DynamicQuantities.jl?branch=main)
+[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
+
+
+  
+DynamicQuantities defines a simple statically-typed `Quantity` type for Julia.
+Physical dimensions are stored as a *value*, as opposed to a parametric type, as in [Unitful.jl](https://github.com/PainterQubits/Unitful.jl).
+This can greatly improve both runtime performance, by avoiding type instabilities, and startup time, as it avoids overspecializing methods.
+
+- [Performance](#performance)
+- [Usage](#usage)
+  - [Constants](#constants)
+  - [Symbolic Units](#symbolic-units)
+    - [Custom Units](#custom-units)
+    - [Affine Units](#affine-units)
+  - [Arrays](#arrays)
+  - [Unitful](#unitful)
+- [Types](#types)
+
+## Performance
+
+DynamicQuantities can greatly outperform Unitful
+when the compiler cannot infer dimensions in a function:
+
+```julia
+julia> using BenchmarkTools, DynamicQuantities; import Unitful
+
+julia> dyn_uni = 0.2u"m/s"
+0.2 m s⁻¹
+
+julia> unitful = convert(Unitful.Quantity, dyn_uni)
+0.2 m s⁻¹
+
+julia> f(x, i) = x ^ i * 0.3;
+
+julia> @btime f($dyn_uni, 1);
+  2.708 ns (0 allocations: 0 bytes)
+
+julia> @btime f($unitful, 1);
+  2.597 μs (30 allocations: 1.33 KiB)
+```
+
+**Note the μ and n: this is a 1000x speedup!**
+Here, the DynamicQuantities quantity object allows the compiler to build a function that is type stable,
+while the Unitful quantity object, which stores its dimensions in the type, requires type inference at runtime.
+
+However, if the dimensions in your function *can* be inferred by the compiler,
+then you can get better speeds with Unitful:
+
+```julia
+julia> g(x) = x ^ 2 * 0.3;
+
+julia> @btime g($dyn_uni);
+  1.791 ns (0 allocations: 0 bytes)
+
+julia> @btime g($unitful);
+  1.500 ns (0 allocations: 0 bytes)
+```
+
+While both of these are type stable,
+because Unitful parametrizes the type on the dimensions, functions can specialize
+to units and the compiler can optimize away units from the code.
+
+## Usage
+
+You can create a `Quantity` object 
+by using the convenience macro `u"..."`:
+
+```julia
+julia> x = 0.3u"km/s"
+300.0 m s⁻¹
+
+julia> y = 42 * u"kg"
+42.0 kg
+```
+
+or by importing explicitly:
+
+```julia
+julia> using DynamicQuantities: kPa
+
+julia> room_temp = 100kPa
+100000.0 m⁻¹ kg s⁻²
+```
+
+Note that `Units` is an exported submodule, so you can
+also access this as `Units.kPa`. You may like to define
+
+```julia
+julia> const U = Units
+```
+
+so that you can simply write, say, `U.kPa` or `C.m_e`.
+
+This supports a wide range of SI base and derived units, with common
+prefixes.
+
+You can also construct values explicitly with the `Quantity` type,
+with a value and keyword arguments for the powers of the physical dimensions
+(`mass`, `length`, `time`, `current`, `temperature`, `luminosity`, `amount`):
+
+```julia
+julia> x = Quantity(300.0, length=1, time=-1)
+300.0 m s⁻¹
+```
+
+Elementary calculations with `+, -, *, /, ^, sqrt, cbrt, abs` are supported:
+
+```julia
+julia> x * y
+12600.0 m kg s⁻¹
+
+julia> x / y
+7.142857142857143 m kg⁻¹ s⁻¹
+
+julia> x ^ 3
+2.7e7 m³ s⁻³
+
+julia> x ^ -1
+0.0033333333333333335 m⁻¹ s
+
+julia> sqrt(x)
+17.320508075688775 m¹ᐟ² s⁻¹ᐟ²
+
+julia> x ^ 1.5
+5196.152422706632 m³ᐟ² s⁻³ᐟ²
+```
+
+Each of these values has the same type, which means we don't need to perform type inference at runtime.
+
+Furthermore, we can do dimensional analysis by detecting `DimensionError`:
+
+```julia
+julia> x + 3 * x
+1.2 m¹ᐟ² kg
+
+julia> x + y
+ERROR: DimensionError: 0.3 m¹ᐟ² kg and 10.2 kg² s⁻² have incompatible dimensions
+```
+
+The dimensions of a `Quantity` can be accessed either with `dimension(quantity)` for the entire `Dimensions` object:
+
+```julia
+julia> dimension(x)
+m¹ᐟ² kg
+```
+
+or with `umass`, `ulength`, etc., for the various dimensions:
+
+```julia
+julia> umass(x)
+1//1
+
+julia> ulength(x)
+1//2
+```
+
+Finally, you can strip units with `ustrip`:
+
+```julia
+julia> ustrip(x)
+0.2
+```
+
+You can also convert a quantity to a desired unit and *then* strip the units
+using a two-argument version of `ustrip`:
+
+```julia
+julia> ustrip(u"km", 1000u"m")
+1.0
+
+julia> ustrip(u"minute", 60u"s")
+1.0
+```
+
+This is equivalent to `ustrip(quantity / unit)` but performs dimension checks first.
+
+### Constants
+
+There are a variety of physical constants accessible
+via the `Constants` submodule:
+
+```julia
+julia> Constants.c
+2.99792458e8 m s⁻¹
+```
+
+which you may like to define as
+
+```julia
+julia> const C = Constants
+```
+
+These can also be used inside the `u"..."` macro:
+
+```julia
+julia> u"Constants.c * Hz"
+2.99792458e8 m s⁻²
+```
+
+Similarly, you can just import each individual constant:
+
+```julia
+julia> using DynamicQuantities.Constants: h
+```
+
+For the full list, see the [docs](https://symbolicml.org/DynamicQuantities.jl/dev/constants/).
+
+
+### Symbolic Units
+
+You can also choose to not eagerly convert to SI base units,
+instead leaving the units as the user had written them.
+For example:
+
+```julia
+julia> q = 100us"cm * kPa"
+100.0 cm kPa
+
+julia> q^2
+10000.0 cm² kPa²
+```
+
+You can convert to regular SI base units with
+`uexpand`:
+
+```julia
+julia> uexpand(q^2)
+1.0e6 kg² s⁻⁴
+```
+
+This also works with constants:
+
+```julia
+julia> x = us"Constants.c * Hz"
+1.0 Hz c
+
+julia> x^2
+1.0 Hz² c²
+
+julia> uexpand(x^2)
+8.987551787368176e16 m² s⁻⁴
+```
+
+You can also convert a quantity in regular base SI units to symbolic units with the `|>` infix operator
+```julia
+julia> 5e-9u"m" |> us"nm"
+5.0 nm
+```
+
+You can also convert between different symbolic units.
+(Note that you can write this more explicitly
+with `uconvert(us"nm", 5e-9u"m")`.)
+
+
+Finally, you can also import these directly:
+
+```julia
+julia> using DynamicQuantities.SymbolicUnits: cm
+```
+
+or constants:
+
+```julia
+julia> using DynamicQuantities.SymbolicConstants: h
+```
+
+Note that `SymbolicUnits` and `SymbolicConstants` are exported,
+so you can simply access these as `SymbolicUnits.cm` and `SymbolicConstants.h`,
+respectively.
+
+
+#### Custom Units
+
+You can create custom units with the `@register_unit` macro:
+
+```julia
+julia> @register_unit OneFiveV 1.5u"V"
+```
+
+and then use it in calculations normally:
+
+```julia
+julia> x = us"OneFiveV"
+1.0 OneFiveV
+
+julia> x * 10u"A" |> us"W"
+15.0 W
+
+julia> 3us"V" |> us"OneFiveV"
+2.0 OneFiveV
+```
+
+#### Affine Units
+
+You can also use "*affine*" units such as Celsius or Fahrenheit,
+using the `ua"..."` string macro:
+
+```julia
+julia> room_temp = 22ua"degC"
+295.15 K
+
+julia> freezing = 32ua"degF"
+273.15 K
+```
+
+These are regular `Quantity{Float64,Dimensions{...}}` objects, meaning that you
+can use them in the same way as regular quantities, including taking differences.
+
+To convert back, you can use the two-argument `ustrip` with the particular affine unit:
+
+```julia
+julia> ustrip(ua"degC", 295.15u"K")
+22.0
+```
+
+### [Arrays](@id home_arrays)
+
+For working with an array of quantities that have the same dimensions,
+you can use a `QuantityArray`:
+
+```julia
+julia> ar = QuantityArray(rand(3), u"m/s")
+3-element QuantityArray(::Vector{Float64}, ::Quantity{Float64, Dimensions{FixedRational{Int32, 25200}}}):
+ 0.2729202669351497 m s⁻¹
+ 0.992546340360901 m s⁻¹
+ 0.16863543422972482 m s⁻¹
+```
+
+This `QuantityArray` is a subtype `<:AbstractArray{Quantity{Float64,Dimensions{...}},1}`,
+meaning that indexing a specific element will return a `Quantity`:
+
+```julia
+julia> ar[2]
+0.992546340360901 m s⁻¹
+
+julia> ar[2] *= 2
+1.985092680721802 m s⁻¹
+
+julia> ar[2] += 0.5u"m/s"
+2.485092680721802 m s⁻¹
+```
+
+This also has a custom broadcasting interface which
+allows the compiler to avoid redundant dimension calculations,
+relative to if you had simply used an array of quantities:
+
+```julia
+julia> f(v) = v^2 * 1.5;
+
+julia> @btime $f.(xa) setup=(xa = randn(100000) .* u"km/s");
+  109.500 μs (2 allocations: 3.81 MiB)
+
+julia> @btime $f.(qa) setup=(xa = randn(100000) .* u"km/s"; qa = QuantityArray(xa));
+  50.917 μs (3 allocations: 781.34 KiB)
+```
+
+So we can see the `QuantityArray` version saves on both time and memory.
+
+By default, DynamicQuantities will create a `QuantityArray` from an `AbstractArray`, similarly to how a `Quantity` is created from a scalar in the [Usage](@ref) examples:
+
+```@repl quantity-array
+using DynamicQuantities # hide
+
+x = [0.3, 0.4, 0.5]u"km/s"
+
+y = (42:45) * u"kg"
+```
+
+This can be overridden to produce a vector of `Quantity`s by explicitly broadcasting the unit:
+
+```@repl quantity-array
+z = [0.3, 0.4, 0.5] .* u"km/s"
+```
+
+### Unitful
+
+DynamicQuantities allows you to convert back and forth from Unitful.jl:
+
+```julia
+julia> using Unitful: Unitful, @u_str; import DynamicQuantities
+
+julia> x = 0.5u"km/s"
+0.5 km s⁻¹
+
+julia> y = convert(DynamicQuantities.Quantity, x)
+500.0 m s⁻¹
+
+julia> y2 = y^2 * 0.3
+75000.0 m² s⁻²
+
+julia> x2 = convert(Unitful.Quantity, y2)
+75000.0 m² s⁻²
+
+julia> x^2*0.3 == x2
+true
+```
+
+## Types
+
+Both a `Quantity`'s values and dimensions are of arbitrary type. The default
+`Dimensions` (for the `u"..."` macro) performs exponent tracking for SI units,
+and `SymbolicDimensions` (for the `us"..."` macro) performs exponent tracking
+for all known unit and constant symbols, using a sparse array.
+
+You can create custom spaces dimension spaces by simply creating
+a Julia struct subtyped to `AbstractDimensions`:
+
+```julia
+julia> struct CookiesAndMilk{R} <: AbstractDimensions{R}
+           cookies::R
+           milk::R
+       end
+
+julia> cookie_rate = Quantity(0.9, CookiesAndMilk(cookies=1, milk=-1))
+0.9 cookies milk⁻¹
+
+julia> total_milk = Quantity(103, CookiesAndMilk(milk=1))
+103 milk
+
+julia> total_cookies = cookie_rate * total_milk
+92.7 cookies
+```
+
+Exponents are tracked by default with the type `R = FixedRational{Int32,C}`,
+which represents rational numbers with a fixed denominator `C`.
+This is much faster than `Rational`.
+
+```julia
+julia> typeof(0.5u"kg")
+Quantity{Float64, Dimensions{FixedRational{Int32, 25200}}}
+```
+
+You can change the type of the value field by initializing with a value
+explicitly of the desired type.
+
+```julia
+julia> typeof(Quantity(Float16(0.5), mass=1, length=1))
+Quantity{Float16, Dimensions{FixedRational{Int32, 25200}}}
+```
+
+or by conversion:
+
+```julia
+julia> typeof(convert(Quantity{Float16}, 0.5u"m/s"))
+Quantity{Float16, Dimensions{FixedRational{Int32, 25200}}}
+```
+
+For many applications, `FixedRational{Int8,6}` will suffice,
+and can be faster as it means the entire `Dimensions`
+struct will fit into 64 bits.
+You can change the type of the dimensions field by passing
+the type you wish to use as the second argument to `Quantity`:
+
+```julia
+julia> using DynamicQuantities
+
+julia> R8 = Dimensions{FixedRational{Int8,6}};
+
+julia> R32 = Dimensions{FixedRational{Int32,2^4 * 3^2 * 5^2 * 7}};  # Default
+
+julia> q8 = [Quantity{Float64,R8}(randn(), length=rand(-2:2)) for i in 1:1000];
+
+julia> q32 = [Quantity{Float64,R32}(randn(), length=rand(-2:2)) for i in 1:1000];
+
+julia> f(x) = @. x ^ 2 * 0.5;
+
+julia> @btime f($q8);
+  1.433 μs (3 allocations: 15.77 KiB)
+
+julia> @btime f($q32);
+  1.883 μs (4 allocations: 39.12 KiB)
+```
+
+
diff --git a/src/arrays.jl b/src/arrays.jl
index e0bf619c..f469ca70 100644
--- a/src/arrays.jl
+++ b/src/arrays.jl
@@ -14,6 +14,14 @@ and so can be used in most places where a normal array would be used, including
 
 # Constructors
 
+The most convenient way to create a `QuantityArray` is by multiplying your array-like object by the desired dimension(s), e.g.,
+
+```julia
+x = [3, 4, 5]u"km/s"
+```
+
+For more control, the following constructors are available:
+
 - `QuantityArray(v::AbstractArray, d::AbstractDimensions)`: Create a `QuantityArray` with value `v` and dimensions `d`,
   using `Quantity` if the eltype of `v` is numeric, and `GenericQuantity` otherwise.
 - `QuantityArray(v::AbstractArray{<:Number}, q::AbstractQuantity)`: Create a `QuantityArray` with value `v` and dimensions inferred
diff --git a/src/math.jl b/src/math.jl
index cfe80378..4dc8f83c 100644
--- a/src/math.jl
+++ b/src/math.jl
@@ -49,6 +49,11 @@ for (type, base_type, _) in ABSTRACT_QUANTITY_TYPES
         function Base.:/(l::AbstractDimensions, r::$type)
             new_quantity(typeof(r), inv(ustrip(r)), l / dimension(r))
         end
+
+        # Defining here instead of outside loop with UnionAbstractQuantity to avoid method ambiguities
+        Base.:*(A::AbstractArray, q::$type) = QuantityArray(A, q)
+        Base.:*(q::$type, A::AbstractArray) = A * q
+        Base.:/(A::AbstractArray, q::$type) = A * inv(q)
     end
 end
 
diff --git a/test/unittests.jl b/test/unittests.jl
index 189dab09..c467abe3 100644
--- a/test/unittests.jl
+++ b/test/unittests.jl
@@ -265,6 +265,7 @@ end
 end
 
 @testset "Arrays" begin
+    T_QA_GenericQuantity(T, N) = QuantityArray{T, N, D, Q, V} where {T, D<:Dimensions, Q<:UnionAbstractQuantity, V<:AbstractArray{T, N}}
     for T in [Float16, Float32, Float64], R in [Rational{Int16}, Rational{Int32}, SimpleRatio{Int}, SimpleRatio{SafeInt16}]
         D = Dimensions{R}
 
@@ -288,9 +289,9 @@ end
         @test ustrip(x + ones(T, 32))[32] == 2
         @test typeof(x + ones(T, 32)) <: GenericQuantity{Vector{T}}
         @test typeof(x - ones(T, 32)) <: GenericQuantity{Vector{T}}
-        @test typeof(ones(T, 32) * GenericQuantity(T(1), D, length=1)) <: GenericQuantity{Vector{T}}
-        @test typeof(ones(T, 32) / GenericQuantity(T(1), D, length=1)) <: GenericQuantity{Vector{T}}
-        @test ones(T, 32) / GenericQuantity(T(1), length=1) == GenericQuantity(ones(T, 32), length=-1)
+        @test typeof(ones(T, 32) * GenericQuantity(T(1), D, length=1)) <: T_QA_GenericQuantity(T, 1)
+        @test typeof(ones(T, 32) / GenericQuantity(T(1), D, length=1)) <: T_QA_GenericQuantity(T, 1)
+        @test ones(T, 32) / GenericQuantity(T(1), length=1) == QuantityArray(ones(T, 32), GenericQuantity(T(1), length=-1))
     end
 
     @testset "isapprox" begin
@@ -395,9 +396,11 @@ end
     end
 
     @testset "Multiplying ranges with units" begin
+        T_QA_StepRangeLen = QuantityArray{T, 1, D, Q, V} where {T, D<:Dimensions, Q<:UnionAbstractQuantity, V<:StepRangeLen}
+        T_QA_s_StepRangeLen = QuantityArray{T, 1, D, Q, V} where {T, D<:SymbolicDimensions, Q<:UnionAbstractQuantity, V<:StepRangeLen}
         # Test multiplying ranges with units
         x = (1:0.25:4)u"inch"
-        @test x isa StepRangeLen
+        @test x isa T_QA_StepRangeLen
         @test first(x) == 1u"inch"
         @test x[2] == 1.25u"inch"
         @test last(x) == 4u"inch"
@@ -405,7 +408,7 @@ end
 
         # Integer range (but real-valued unit)
         x = (1:4)u"inch"
-        @test x isa StepRangeLen
+        @test x isa T_QA_StepRangeLen
         @test first(x) == 1u"inch"
         @test x[2] == 2u"inch"
         @test last(x) == 4u"inch"
@@ -414,7 +417,7 @@ end
 
         # Test with floating point range
         x = (1.0:0.5:3.0)u"m"
-        @test x isa StepRangeLen
+        @test x isa T_QA_StepRangeLen
         @test first(x) == 1.0u"m"
         @test x[2] == 1.5u"m"
         @test last(x) == 3.0u"m"
@@ -427,7 +430,7 @@ end
 
         # Test with symbolic units
         x = (1:0.25:4)us"inch"
-        @test x isa StepRangeLen{<:Quantity{Float64,<:SymbolicDimensions}}
+        @test x isa T_QA_s_StepRangeLen
         @test first(x) == us"inch"
         @test x[2] == 1.25us"inch"
         @test last(x) == 4us"inch"
@@ -435,7 +438,7 @@ end
 
         # Test that symbolic units preserve their symbolic nature
         x = (0:0.1:1)us"km/h"
-        @test x isa AbstractRange
+        @test x isa QuantityArray
         @test first(x) == 0us"km/h"
         @test x[2] == 0.1us"km/h"
         @test last(x) == 1us"km/h"
@@ -443,14 +446,14 @@ end
 
         # Similarly, integers should stay integers:
         x = (1:4)us"inch"
-        @test_skip x isa StepRangeLen{<:Quantity{Int64,<:SymbolicDimensions}}
+        @test_skip x isa T_QA_s_StepRangeLen
         @test first(x) == us"inch"
         @test x[2] == 2us"inch"
         @test last(x) == 4us"inch"
         @test length(x) == 4
 
         # With RealQuantity:
-        @test_skip (1.0:4.0) * RealQuantity(u"inch") isa StepRangeLen{<:RealQuantity{Float64,<:SymbolicDimensions}}
+        @test_skip (1.0:4.0) * RealQuantity(u"inch") isa T_QA_s_StepRangeLen
         # TODO: This is not available as TwicePrecision interacts with Real in a way
         #       that demands many other functions to be defined.
     end