diff --git a/src/AbstractTypes.jl b/src/AbstractTypes.jl
index ae25423a4c..bd036f22d3 100644
--- a/src/AbstractTypes.jl
+++ b/src/AbstractTypes.jl
@@ -76,6 +76,8 @@ abstract type IdentityMap <: SetMap end
 
 abstract type FPModuleHomomorphism <: FunctionalMap end
 
+abstract type SkewDerivation{D,S} <: Map{D,D,Map,SkewDerivation} where {S<:Map{D,D}} end
+
 # rings, fields etc, parameterised by an element type
 # these are the type classes of different kinds of
 # mathematical rings/fields/etc, which have a base ring,
@@ -106,6 +108,8 @@ abstract type MatRing{T<:NCRingElement} <: NCRing end
 
 abstract type FreeAssociativeAlgebra{T<:RingElement} <: NCRing end
 
+abstract type OreAlgebra{T<:RingElement} <: NCPolyRing{T} end
+
 # Abstract types for number fields, parmeterised by the element type of
 # the base field.
 abstract type NumField{T<:RingElement} <: Field end
@@ -151,6 +155,8 @@ abstract type MatRingElem{T<:NCRingElement} <: NCRingElem end
 
 abstract type FreeAssociativeAlgebraElem{T<:RingElement} <: NCRingElem end
 
+abstract type OreOperator{T<:RingElem} <: NCPolyRingElem{T} end
+
 abstract type NumFieldElem{T<:RingElement} <: FieldElem end
 
 abstract type SimpleNumFieldElem{T} <: NumFieldElem{T} end
diff --git a/src/Generic.jl b/src/Generic.jl
index e88ee316d5..121e0f99c2 100644
--- a/src/Generic.jl
+++ b/src/Generic.jl
@@ -94,6 +94,10 @@ include("generic/FreeAssociativeAlgebraGroebner.jl")
 
 include("generic/PolyRingHom.jl")
 
+include("generic/OreAlgebras.jl")
+include("generic/PolySkewDerivation.jl")
+include("generic/PolyFracFieldHom.jl")
+
 ###############################################################################
 #
 #   Temporary miscellaneous files being moved from Hecke.jl
diff --git a/src/exports.jl b/src/exports.jl
index fc1e5075cd..9761526437 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -620,3 +620,11 @@ export zero!
 export zero_matrix
 export zeros
 export zz
+
+export ore_extension
+export trivial_derivation
+export derivation
+export sigma_endomorphism
+export OreAlgebra
+export OreOperator
+export SkewDerivation
diff --git a/src/generic/GenericTypes.jl b/src/generic/GenericTypes.jl
index e0e4ab061b..f5a28ccae0 100644
--- a/src/generic/GenericTypes.jl
+++ b/src/generic/GenericTypes.jl
@@ -1607,3 +1607,64 @@ end
 function PolyRingAnyMap(d::D, c::C, cm::U, ig::V) where {D, C, U, V}
   return PolyRingAnyMap{D, C, U, V}(d, c, cm, ig)
 end
+
+################################################################################
+#
+#   R-algebra morphisms of R(x)
+#
+################################################################################
+
+mutable struct PolyFracFieldAnyMap{
+    D <: Union{FracField{<:PolyRingElem},RationalFunctionField{<:RingElement,<:PolyRingElem}},
+    C <: NCRing,
+    V <: Map{<:PolyRing,C}} <: Map{D,C,Map,PolyFracFieldAnyMap}
+
+  domain::D
+  morphism::V
+
+  function PolyFracFieldAnyMap{D,C,V}(d::D,phi::V) where {D,C,V}
+    return new{D,C,V}(d,phi)
+  end
+end
+
+function PolyFracFieldAnyMap(d::D,i::V) where {D,C,V<:Map{<:PolyRing,C}}
+  return PolyFracFieldAnyMap{D,C,V}(d,i)
+end
+
+
+################################################################################
+#
+#   Univariate Ore algebra, its element type skew derivations
+#
+################################################################################
+@attributes mutable struct OreAlgebra{T<:RingElem} <: AbstractAlgebra.OreAlgebra{T}
+  base_ring::Ring
+  D::Symbol
+  δ::Map(SkewDerivation)
+
+  function OreAlgebra{T}(R::Ring, D::Symbol, δ; cached=true) where T<:RingElem
+    return get_cached!(OreID, (R,D,δ), cached) do
+      new{T}(R,D,δ)
+    end
+  end
+end
+
+const OreID = CacheDictType{Tuple{Ring,Symbol,Map(SkewDerivation)},NCRing}()
+
+mutable struct OreOperator{T<:RingElem} <: AbstractAlgebra.OreOperator{T}
+  parent::OreAlgebra{T}
+  coeffs::Vector{T}
+  length::Int
+end
+
+mutable struct PolySkewDerivation{D,S} <: SkewDerivation{D,S}
+  domain::D
+  σ::S
+  intermediate_cache::Vector{<:NCRingElem}
+
+  function PolySkewDerivation{D,S}(dom::D,σ::S,coeff::T) where {D,S,T<:NCRingElem}
+    return new{D,S}(dom,σ,[coeff])
+  end
+end
+
+Univariateish = Union{PolyRing,FracField{<:PolyRingElem},RationalFunctionField}
diff --git a/src/generic/OreAlgebras.jl b/src/generic/OreAlgebras.jl
new file mode 100644
index 0000000000..53e34cd7a7
--- /dev/null
+++ b/src/generic/OreAlgebras.jl
@@ -0,0 +1,307 @@
+function ore_extension(R::DT, D::Symbol, δ) where {DT<:Ring}
+  σ = sigma_endomorphism(δ)
+  OO = ore_extension_only(R,D,δ)
+  return OO,gen(OO)
+end
+ore_extension_only(R::T, D::Symbol, δ) where T<:Ring = OreAlgebra{elem_type(R)}(R, D, δ)
+
+function (R::OreAlgebra{T})() where T
+  return zero(R)
+end
+
+function (R::OreAlgebra{T})(c::T) where T
+  iszero(c) && return zero(R)
+  return OreOperator{T}(R,[c],1)
+end
+
+function (R::OreAlgebra{T})(c::Union{Integer,Rational,AbstractFloat}) where T
+  C = base_ring(R)
+  return R(C(c))
+end
+
+function (R::OreAlgebra{T})(a::OreOperator{T}) where T
+ parent(a) != R && error("Unable to coerce polynomial")
+ return a
+end
+
+function (R::OreAlgebra{T})(c::Vector{T}) where T<:RingElem
+  return OreOperator{T}(R,c,length(c))
+end
+
+function (R::OreAlgebra{T})(c::Vector{U}) where {T<:RingElem,U<:RingElem}
+  return OreOperator{T}(R,c,length(c))
+end
+
+function (R::OreAlgebra{T})(c::Vector{U}) where {T<:RingElem,U<:Integer}
+  C = base_ring(R)
+  return OreOperator{T}(R,C.(c),length(c))
+end
+
+###############################################################################
+#
+#   OreAlgebra
+#
+###############################################################################
+
+elem_type(::Type{OreAlgebra{T}}) where T = OreOperator{T}
+base_ring_type(::Type{OreAlgebra{T}}) where T = parent_type(T)
+
+base_ring(R::OreAlgebra) = R.base_ring
+
+var(R::OreAlgebra) = R.D
+symbols(R::OreAlgebra) = [var(R)]
+
+function gen(R::OreAlgebra{T}) where T
+  C = base_ring(R)
+  return OreOperator{T}(R,[zero(C),one(C)],2)
+end
+
+ngens(R::OreAlgebra) = 1
+gens(R::OreAlgebra) = [gen(R)]
+
+function zero(R::OreAlgebra{T}) where T
+  return OreOperator{T}(R,T[],0)
+end
+
+function one(R::OreAlgebra{T}) where T
+  return OreOperator{T}(R,T[one(base_ring(R))],1)
+end
+
+derivation(R::OreAlgebra) = R.δ
+sigma_endomorphism(R::OreAlgebra) = sigma_endomorphism(derivation(R))
+
+###############################################################################
+#
+#   OreOperator
+#
+###############################################################################
+
+parent_type(::Type{OreOperator{T}}) where T = OreAlgebra{T}
+parent(a::OreOperator{T}) where T = a.parent
+
+is_domain_type(::Type{OreOperator{T}}) where T = false
+is_exact_type(::Type{OreOperator{T}}) where T = is_exact_type(T)
+
+promote_rule(::Type{OreOperator{T}},::Type{OreOperator{T}}) where T<:RingElement = OreOperator{T}
+function promote_rule(::Type{OreOperator{T}},::Type{U}) where {T<:RingElement,U<:RingElement}
+  promote_rule(T, U) == T ? OreOperator{T} : Union{}
+end
+
+function deepcopy_internal(a::OreOperator{T}, dict::IdDict) where T
+  C = [deepcopy_internal(c,dict) for c in coefficients(a)]
+
+  return OreOperator{T}(parent(a), C, length(C))
+end
+
+coefficients(a::OreOperator{T}) where T = a.coeffs[1:length(a)]
+length(a::OreOperator{T}) where T = a.length
+
+iszero(a::OreOperator{T}) where T = a.length == 0
+function isone(a::OreOperator{T}) where T
+  return length(a) == 1 && first(coefficients(a)) |> isone
+end
+
+function canonical_unit(a::OreOperator{T}) where T
+  return one(parent(a))
+end
+
+function coeff(a::OreOperator{T},i) where T
+  if 0 <= i && i <= order(a)
+    return a.coeffs[i+1]
+  else
+    return zero(base_ring(a))
+  end
+end
+
+setcoeff!(a::OreOperator{T}, n::Int, c::Integer) where T = setcoeff!(a,n,base_ring(a)(c))
+function setcoeff!(a::OreOperator{T}, n::Int, c::T) where T
+  i = n+1
+  if i > length(a.coeffs)
+    resize!(a.coeffs,i)
+    a.coeffs[a.length+1:i] .= zero(base_ring(a))
+  end
+  a.coeffs[i] = c
+  a.length = normalise(a,length(a.coeffs))
+  return a
+end
+
+function leading_coefficient(a::OreOperator{T}) where T
+  iszero(a) && return base_ring(a)()
+  return coefficients(a) |> last
+end
+function leading_term(a::OreOperator{T}) where T
+  iszero(a) && return parent(a)()
+  A = parent(a)
+  C = base_ring(a)
+  c = leading_coefficient(a)
+  k = order(a)
+
+  return A([fill(zero(C),k); c])
+end
+function leading_monomial(a::OreOperator{T}) where T
+  iszero(a) && return parent(a)()
+  A = parent(a)
+  C = base_ring(a)
+  k = order(a)
+
+  return A(T[(zero(C) for _ in 1:k)..., one(C)])
+end
+
+order(a::OreOperator{T}) where T = iszero(a) ? -1 : length(a)-1
+
+function +(a::OreOperator{T},b::OreOperator{T}) where T<:RingElem
+  check_parent(a,b)
+
+  la = length(a)
+  lb = length(b)
+  l = min(la,lb)
+  L = max(la,lb)
+
+  new_c = Vector{T}(undef,L)
+
+  new_c[1:l] .= coefficients(a)[1:l] + coefficients(b)[1:l]
+  c = la > lb ? coefficients(a) : coefficients(b)
+  new_c[(l+1):L] = c[(l+1):L]
+
+  if all(iszero,new_c)
+    return zero(parent(a))
+  else
+    i = findlast(!iszero,new_c)
+    return OreOperator{T}(parent(a), new_c, i)
+  end
+end
+
+function -(a::OreOperator{T}) where T
+  R = parent(a)
+  c = -a.coeffs
+
+  return OreOperator{T}(R,c,length(c))
+end
+
+function -(a::OreOperator{T},b::OreOperator{T}) where T
+  return a + -b
+end
+
+function *(a::OreOperator{T},b::T) where T<:RingElem
+  return a*parent(a)(b)
+end
+
+function *(a::T,b::OreOperator{T}) where T<:RingElem
+  R = parent(b)
+  res = a.*coefficients(b)
+  return OreOperator{T}(R,res,length(b))
+end
+
+function *(a::OreOperator{T},b::OreOperator{T}) where T
+  check_parent(a,b)
+  (iszero(a) || iszero(b)) && return zero(a)
+
+  R = parent(a)
+  k = base_ring(R)
+  δ = derivation(R)
+  σ = sigma_endomorphism(R)
+
+  res = zeros(parent(a)|>base_ring,length(b))
+  q = coefficients(b)
+  l = length(q)
+  for p in coefficients(a)
+    res += p.*q
+    push!(res,zero(R|>base_ring))
+    q = [δ.(q); zero(k)] .+ [zero(k); σ.(q)]
+  end
+
+  i = findlast(!iszero,res)
+  if isnothing(i)
+    return zero(R)
+  end
+
+  return OreOperator{T}(R,res[1:i],i)
+end
+
+function ^(a::OreOperator{T},i::Int) where T
+  if i < 0
+    throw(DomainError(i, "exponent must be >= 0"))
+  elseif i == 0
+    return one(parent(a))
+  elseif i == 1
+    return a
+  else
+    return reduce(*,Iterators.repeated(a,i-1);init=a)
+  end
+end
+
+
+function divexact_right(a::OreOperator{T},b::OreOperator{T}; check=true) where T
+  iszero(b) && throw(DivideError())
+  q,r = rdivrem(a,b)
+
+  !check && return q
+  @req iszero(r) "not an exact division"
+  return q
+end
+
+function ==(a::OreOperator{T},b::OreOperator{T}) where T
+  fl = check_parent(a,b,false)
+  !fl && return false
+  if a.length != b.length
+    return false
+  end
+  return all(splat(==),zip(a.coeffs,b.coeffs))
+end
+
+function ==(a::OreOperator{T}, c::T) where T
+  if iszero(c)
+    return a.length == 0
+  elseif a.length == 1
+    return leading_coefficient(a) == c
+  else
+    return false
+  end
+end
+==(c::T,a::OreOperator{T}) where T = a == c
+
+function is_unit(a::OreOperator{T}) where T
+  iszero(a) && return false
+  isone(a) && return true
+  order(a) == 0 && leading_coefficient(a)|>is_unit && return true
+end
+
+function zero!(a::OreOperator{T}) where T
+  a.coeffs = elem_type(base_ring(a))[]
+  a.length = 0
+
+  return a
+end
+
+function one!(a::OreOperator{T}) where T
+  a.coeffs = [one(base_ring(a))]
+  a.length = 1
+
+  return a
+end
+
+function set_length!(a::OreOperator{T},n::Int) where T
+  resize!(a.coeffs, n)
+  a.length = n
+
+  return a
+end
+
+function normalise(a::OreOperator{T}, n::Int) where T
+  n = min(n,length(a.coeffs))
+  i = findlast(!iszero, a.coeffs[1:n])
+  isnothing(i) && return 0
+  return i
+end
+
+function fit!(a::OreOperator{T}, n::Int) where T
+  if n > length(a)
+    old_n = length(a)+1
+    resize!(a.coeffs, n)
+    a.coeffs[old_n:n] .= base_ring(a)|>zero
+  end
+
+  return
+end
+
+
diff --git a/src/generic/PolyFracFieldHom.jl b/src/generic/PolyFracFieldHom.jl
new file mode 100644
index 0000000000..9d20d985ef
--- /dev/null
+++ b/src/generic/PolyFracFieldHom.jl
@@ -0,0 +1,58 @@
+################################################################################
+#
+#  Field access
+#
+################################################################################
+domain(f::PolyFracFieldAnyMap{D,C,V}) where {D,C,V} = f.domain
+codomain(f::PolyFracFieldAnyMap{D,C,V}) where {D,C,V} = f.morphism.codomain
+underlying_morphism(f) = f.morphism
+
+################################################################################
+#
+#  String I/O
+#
+################################################################################
+function Base.show(io::IO, f::PolyFracFieldAnyMap)
+  io = pretty(io)
+  if is_terse(io)
+    print(io, "Ring homomorphism")
+  else
+    print(io, "Hom: ")
+    print(terse(io), Lowercase(), domain(f), " -> ")
+    print(terse(io), Lowercase(), codomain(f))
+  end
+end
+
+function AbstractAlgebra.show_map_data(io::IO, f::PolyFracFieldAnyMap)
+  println(io)
+  println(io, "defined by", Indent())
+  R = domain(f)
+  g = gen(R)
+  print(io, g, " -> ", f(g))
+  print(io, Dedent())
+end
+
+################################################################################
+#
+#  Constructor
+#
+################################################################################
+function hom(R::D, S::NCRing, image) where {D<:FracField{<:PolyRingElem}}
+  r = base_ring(R)
+  return PolyFracFieldAnyMap(R,hom(r,S,image))
+end
+
+function hom(R::D, S::NCRing, image) where {D<:RationalFunctionField{<:RingElement,<:PolyRingElem}}
+  r = R() |> numerator |> parent
+  return PolyFracFieldAnyMap(R,hom(r,S,image))
+end
+
+################################################################################
+#
+#  Evaluation functions
+#
+################################################################################
+function (f::PolyFracFieldAnyMap{D,C,V})(a) where {D,C,V<:Map{<:PolyRing,C}}
+  ϕ = underlying_morphism(f)
+  return ϕ.(numerator(a))/ϕ.(denominator(a))
+end
diff --git a/src/generic/PolySkewDerivation.jl b/src/generic/PolySkewDerivation.jl
new file mode 100644
index 0000000000..5af0385ce9
--- /dev/null
+++ b/src/generic/PolySkewDerivation.jl
@@ -0,0 +1,112 @@
+
+function trivial_derivation(R::D,σ::S) where {D<:Univariateish,S}
+  return PolySkewDerivation{D,S}(R,σ,zero(R))
+end
+
+function derivation(R::D,σ::S,c::T) where {D<:Univariateish,S<:Map{D,D},T}
+  @req elem_type(R) == T "incompatible coefficient type"
+  return PolySkewDerivation{D,S}(R,σ,c)
+end
+derivation(R::D) where {D<:Univariateish} = derivation(R,identity_map(R),one(R))
+derivation(R::D,σ::S) where {D<:Univariateish,S<:Map{D,D}} = derivation(R,σ,one(R))
+
+function show(io::IO, d::PolySkewDerivation)
+  io = pretty(io)
+  if is_terse(io)
+    print(io, "Skew-derivation")
+  else
+    print(io, "σ-Der: ")
+    print(terse(io), Lowercase(), domain(d), " -> ")
+    print(terse(io), Lowercase(), codomain(d))
+  end
+end
+
+function show(io::IO, d::PolySkewDerivation{D,<:Map(IdentityMap)}) where D
+  io = pretty(io)
+  if is_terse(io)
+    print(io, "Derivation")
+  else
+    print(io, "Der: ")
+    print(terse(io), Lowercase(), domain(d), " -> ")
+    print(terse(io), Lowercase(), codomain(d))
+  end
+end
+
+function show_map_data(io::IO, d::PolySkewDerivation)
+  println(io)
+  if coefficient(d) |> iszero
+    print(io, "which is identically zero ")
+    println(io, "and with commutation rule", Indent())
+  elseif !(coefficient(d) |> isone)
+    println(io, "with coefficient ", coefficient(d))
+    println(io, "and commutation rule", Indent())
+  else
+    println(io, "with commutation rule", Indent())
+  end
+  print(io, sigma_endomorphism(d))
+  show_map_data(io,sigma_endomorphism(d))
+end
+
+function show_map_data(io::IO, d::PolySkewDerivation{D,<:Map(IdentityMap)}) where D
+  if !(coefficient(d) |> isone)
+    println(io)
+    print(io, "with coefficient ", coefficient(d))
+  end
+end
+
+domain(d::PolySkewDerivation) = d.domain
+codomain(d::PolySkewDerivation) = d.domain
+coefficient(d::PolySkewDerivation) = first(d.intermediate_cache)
+
+sigma_endomorphism(δ::PolySkewDerivation{D,S}) where {D,S} = δ.σ
+
+function (d::PolySkewDerivation{D,S})(a) where {D,S<:Map(IdentityMap)}
+  c = coefficient(d)
+  iszero(c) && return domain(d)()
+  return c*derivative(a)
+end
+
+function (d::PolySkewDerivation{D,S})(a::T) where {D,S<:Map{D,D,<:Map,<:Any},T<:PolyRingElem}
+  c = coefficient(d)
+  iszero(c) && return domain(d)()
+
+  x = parent(a) |> gen
+  σ = sigma_endomorphism(d)
+  cached_degree = length(d.intermediate_cache+1)
+  for i in cached_degree:degree(a)
+    res = last(d.intermediate_cache)*x + σ(x^(i-1))*c
+    push!(d.intermediate_cache, res)
+  end
+
+  return sum(σ.(Iterators.drop(coefficients(a),1)) .* d.intermediate_cache; init=parent(a)())
+end
+
+function (d::PolySkewDerivation{D,S})(a::T) where {
+  D<:FracField{<:PolyRingElem},
+  S<:Map{D,D,<:Map,<:Any},
+  T<:FracElem{<:PolyRingElem}}
+  
+  c = coefficient(d)
+  iszero(c) && return domain(d)()
+
+  R = codomain(d)
+  σ = sigma_endomorphism(d)
+  p = numerator(a)
+  q = denominator(a)
+  return R((d(p) - σ(q)*d(p))//(σ(q)*q))::elem_type(R)
+end
+
+function (d::PolySkewDerivation{D,S})(a::T) where {
+  D<:RationalFunctionField{<:RingElement,<:PolyRingElem},
+  S<:Map{D,D,<:Map,<:Any},
+  T<:RationalFunctionFieldElem{<:RingElement,<:PolyRingElem}}
+
+  c = coefficient(d)
+  iszero(c) && return domain(d)()
+
+  R = codomain(d)
+  σ = sigma_endomorphism(d)
+  p = numerator(a)
+  q = denominator(a)
+  return R((d(p) - σ(q)*d(p))//(σ(q)*q))::elem_type(codomain(d))
+end
diff --git a/src/generic/exports.jl b/src/generic/exports.jl
index 8379a4c06c..b9459c723c 100644
--- a/src/generic/exports.jl
+++ b/src/generic/exports.jl
@@ -133,3 +133,8 @@ export unit
 export upscale
 export weights
 export zero
+
+export ore_extension
+export trivial_derivation
+export derivation
+export sigma_endomorphism
diff --git a/test/NCRings-test.jl b/test/NCRings-test.jl
index c89f8e5bf5..f6ff9316c0 100644
--- a/test/NCRings-test.jl
+++ b/test/NCRings-test.jl
@@ -3,6 +3,8 @@ include("generic/NCPoly-test.jl")
 include("generic/FreeAssociativeAlgebra-test.jl")
 include("generic/FreeAssociativeAlgebraGroebner-test.jl")
 
+#include("generic/OreOperator-test.jl")
+
 @testset "NCRings.oftype" begin
    F = GF(3)
    Fx, x = polynomial_ring(F, "x")
diff --git a/test/generic/OreOperator-test.jl b/test/generic/OreOperator-test.jl
new file mode 100644
index 0000000000..8a938daefa
--- /dev/null
+++ b/test/generic/OreOperator-test.jl
@@ -0,0 +1,187 @@
+@testset "Generic.OreOperator" begin
+  R,x = polynomial_ring(QQ,:x)
+  @testset "Differential operators" verbose=true begin
+    δ = derivation(R)
+    AD,D = ore_extension(R,:D,δ)
+
+    @test order(D) == 1
+    @test length(D) == 2
+
+    a0 = zero(AD)
+    @test order(a0) == -1
+    @test length(a0) == 0
+    @test a0.coeffs == elem_type(R)[]
+    @test iszero(a0)
+    @test !is_unit(a0)
+
+    @test leading_coefficient(a0) |> iszero
+    @test leading_term(a0)  |> iszero
+    @test leading_monomial(a0) |> iszero
+
+    a1 = one(AD)
+    @test order(a1) == 0
+    @test length(a1) == 1
+    @test is_unit(a1)
+
+    @test leading_coefficient(a1) |> isone
+    @test leading_term(a1) |> isone
+    @test leading_monomial(a1) |> isone
+
+    @test a0 * a1 == a0
+    @test a1 * a1 == a1
+
+    a2 = D^2
+    @test order(a2) == 2
+    @test coefficients(a2) == R.([0,0,1])
+
+    a3 = 3*a2 + D
+    @test leading_coefficient(a3) == R(3)
+    @test leading_monomial(a3) == D^2
+    @test leading_term(a3) == leading_coefficient(a3)*leading_monomial(a3)
+
+#    @testset "Division with remainder" verbose=true begin
+#      @test rdivrem(D,D) == (one(AD),zero(AD))
+#      @test rdivrem(D^2+D,D) == (D+1,zero(AD))
+#    end
+
+    @testset "Arithmetic" verbose=true begin
+      c4 = [x, x^3, x^2, x^5]
+      a4 = (c4[1]*D + c4[2])*(c4[3]*D + c4[4])
+      @test coefficients(a4) == [c4[1]*δ(c4[4]) + c4[2]*c4[4], c4[1]*δ(c4[3]) + c4[2]*c4[3] + c4[1]*c4[4], c4[1]*c4[3]]
+      @test length(a4) == 3
+      @test order(a4) == 2
+    end
+  end
+
+  @testset "Euler operators" verbose=true begin
+    σ = identity_map(R)
+    δ = derivation(R,σ,x)
+    AE,θ = ore_extension(R,:θ,δ)
+
+    a1 = θ*x
+    @test order(a1) == 1
+    @test coefficients(a1) == [x,x]
+
+    a2 = θ*x^2
+    @test order(a2) == 1
+    @test coefficients(a2) == [2*x^2, x^2]
+  end
+
+  @testset "Forward differences" verbose=true begin
+    σ = hom(R,R,x+1)
+    δ = derivation(R,σ)
+    AF,F = ore_extension(R,:F,δ)
+
+    a1 = F*x
+    @test order(a1) == 1
+    @test coefficients(a1) == [one(R),x+1]
+
+    a2 = F*x^2
+    @test order(a2) == 1
+    @test coefficients(a2) == [2*x + 1, x^2 + 2*x + 1]
+  end
+
+  @testset "Recurrence operators" verbose=true begin
+    σ = hom(R,R,x+1)
+    δ = trivial_derivation(R,σ)
+    AS,S = ore_extension(R,:S,δ)
+
+#    @testset "Division with remainder" verbose=true begin
+#      @test rdivrem(S,S) == (one(AS),zero(AS))
+#      @test rdivrem(S^2+S,S) == (S+1,zero(AS))
+#    end
+
+    @testset "Arithmetic" verbose=true begin
+      c1 = [x, x^3, x^2, x^5]
+      a1 = (c1[1]*S + c1[2])*(c1[3]*S + c1[4])
+      @test coefficients(a1) == [c1[2]*c1[4],c1[1]*σ(c1[4]) + c1[2]*c1[3],c1[1]*σ(c1[3])]
+      @test length(a1) == 3
+      @test order(a1) == 2
+    end
+  end
+end
+
+#@testset "Univariate Ore algebra over $name" verbose=true for (name,(R,x)) in [
+#    ("univariate polynomial fraction field", begin
+#      R,x = polynomial_ring(QQ,:x)
+#      K = fraction_field(R)
+#      K,K(x)
+#     end),
+#    ("univariate polynomial function field", rational_function_field(QQ)),
+#  ]
+#
+#  @testset "Differential operators" verbose=true begin
+#    δ = derivation(R)
+#    σ = sigma_endomorphism(δ)
+#    AD,D = ore_extension(R,:D,δ)
+#
+#    U = AD([3x+1,2x-3,-(2x+1),3x+5])
+#    V = AD([-(3x+5),x+2])
+#    q,r = rdivrem(U,V)
+#    @test q*V + r == U
+#
+#    U = (x-1)*D^5 + 5*D^4
+#    V = (x-1)*D^3 + (6 - 3x)*D^2 + (27x^2 + 9x -42)*D + 117 - 54x
+#    g = grcd(U,V)
+#    @test g == D - (2x^2-x-7)//((x+2)*(x+1)*(x-1))
+#
+#    U = (2x+3)*x*D^2 + 2*(4x^2 + 3x - 3)*D + 2*(4x^2 - 3)
+#    V = (x+1)*(x-1)*D^2 + 2*(2x^2 - x - 2)*D + 2*(2x^2 - 2x - 1)
+#    f = lclm(U,V)
+#    @test iszero(rrem(f,U))
+#    @test iszero(rrem(f,V))
+#
+#    @testset "Hermite normal form" verbose=true begin
+#      H1 = AD[D-1 D; 0 D^2]
+#      @test hnf(H1) == H1
+#
+#      H2 = AD[D-(3x^2+4)//(x*(3x+2)) 6*(x-2)//(x*(3x+2)); -(3-x)//(3x+2) D-11//(3x+2)]
+#      @test hnf(H2) == AD[1 (3x+2)//(x-3)*D-11/(x-3); 0 D^2+(6-x^2)//(x*(x-3))*D+3*(x-2)//(x*(x-3))]
+#    end
+#
+#    @testset "Uncoupling" verbose=true begin
+#      H3 = R[(3x^2+4)//(x*(3x+2)) -6*(x-2)//(x*(3x+2)); (3-x)//(3x+2) 11//(3x+2)]
+#      P3 = R[1 0; (3x^2+4)//(x*(3x+2)) -6*(x-2)//(x*(3x+2))]
+#      @test uncouple(H3,δ) == P3
+#      @test gauge_transform(H3,δ) == (σ.(P3)*H3 + δ.(P3))*inv(P3)
+#    end
+#  end
+#  @testset "Recurrence operators" verbose=true begin
+#    σ = hom(R,R,x+1)
+#    δ = trivial_derivation(R,σ)
+#    AS,S = ore_extension(R,:S,δ)
+#
+#    U = AS([3x+1,2x-3,-(2x+1),3x+5])
+#    V = AS([-(3x+5),x+2])
+#    q,r = rdivrem(U,V)
+#    @test q*V + r == U
+#
+#    a1 = (x*S + x) * (2*x*S + x)
+#    @test a1 == x*2*(x+1)*S^2 + (x*(x+1) + x*2*x)*S + x*x
+#
+#    U = S^7 + 5S^6 + 9S^5 +5S^4 - 5S^3 - 9S^2 - 5S - 1
+#    V = (x+5)*S^5 + 6S^4 - 6S^3 - 2*(x-1)*S^2 - (x+19)*S + 2*(x+6)
+#    g = grcd(U,V)
+#    @test g == S^2 + 10*(x+4)//((x+2)*(2x+7))*S - (2x+9)*(x+6)//((2x+7)*(x+2))
+#
+#    U = (x+2)*(2x-1)*S^2 - 8*(x^2+3x-1)*S + 4*(x+4)*(2x+1)
+#    V = (x+4)*(x+3)*S^2 - 2*(x+5)*(2x+5)*S + 4*(x+4)*(x+5)
+#    f = lclm(U,V)
+#    @test iszero(rrem(f,U))
+#    @test iszero(rrem(f,V))
+#
+#    @testset "Division with remainder" verbose=true begin
+#      @test rdivrem(S,S) == (one(AS),zero(AS))
+#      @test rdivrem(S^2+S,S) == (S+1,zero(AS))
+#    end
+#
+#    @testset "Arithmetic" verbose=true begin
+#      c1 = [x, x^3, x^2, x^5]
+#      a1 = (c1[1]*S + c1[2])*(c1[3]*S + c1[4])
+#      @test coefficients(a1) == [c1[2]*c1[4],c1[1]*σ(c1[4]) + c1[2]*c1[3],c1[1]*σ(c1[3])]
+#      @test length(a1) == 3
+#      @test order(a1) == 2
+#    end
+#  end
+#end
+#