diff --git a/Project.toml b/Project.toml
index 2abb41b..d941fb8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,20 +1,24 @@
 name = "ArrayDiff"
 uuid = "c45fa1ca-6901-44ac-ae5b-5513a4852d50"
-authors = ["Sophie Lequeu <slequeu@hotmail.com>", "Benoît Legat <benoit.legat@gmail.com>"]
 version = "0.1.0"
+authors = ["Sophie Lequeu <slequeu@hotmail.com>", "Benoît Legat <benoit.legat@gmail.com>"]
 
 [deps]
+Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
+Calculus = "0.5.2"
 DataStructures = "0.18, 0.19"
 ForwardDiff = "1"
 MathOptInterface = "1.40"
 NaNMath = "1"
 SparseArrays = "1.10"
+SpecialFunctions = "2.6.1"
 julia = "1.10"
diff --git a/src/ArrayDiff.jl b/src/ArrayDiff.jl
index faf460f..2480b3f 100644
--- a/src/ArrayDiff.jl
+++ b/src/ArrayDiff.jl
@@ -36,13 +36,14 @@ import NaNMath:
 include("Coloring/Coloring.jl")
 include("graph_tools.jl")
 include("sizes.jl")
+include("univariate_expressions.jl")
+include("operators.jl")
 include("types.jl")
 include("utils.jl")
 
 include("reverse_mode.jl")
 include("forward_over_reverse.jl")
 include("mathoptinterface_api.jl")
-include("operators.jl")
 include("model.jl")
 include("parse.jl")
 include("evaluator.jl")
diff --git a/src/mathoptinterface_api.jl b/src/mathoptinterface_api.jl
index 21c6804..9cc3b6e 100644
--- a/src/mathoptinterface_api.jl
+++ b/src/mathoptinterface_api.jl
@@ -4,7 +4,7 @@
 # Use of this source code is governed by an MIT-style license that can be found
 # in the LICENSE.md file or at https://opensource.org/licenses/MIT.
 
-_no_hessian(op::MOI.Nonlinear._UnivariateOperator) = op.f′′ === nothing
+_no_hessian(op::_UnivariateOperator) = op.f′′ === nothing
 _no_hessian(op::MOI.Nonlinear._MultivariateOperator) = op.∇²f === nothing
 
 function MOI.features_available(d::NLPEvaluator)
diff --git a/src/model.jl b/src/model.jl
index c1f9480..acc1234 100644
--- a/src/model.jl
+++ b/src/model.jl
@@ -53,7 +53,7 @@ function register_operator(
         elseif haskey(registry.multivariate_operator_to_id, op)
             error("Operator $op is already registered.")
         end
-        operator = Nonlinear._UnivariateOperator(op, f...)
+        operator = _UnivariateOperator(op, f...)
         push!(registry.univariate_operators, op)
         push!(registry.registered_univariate_operators, operator)
         registry.univariate_operator_to_id[op] =
diff --git a/src/operators.jl b/src/operators.jl
index f293928..a6f89eb 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -20,6 +20,158 @@ const DEFAULT_MULTIVARIATE_OPERATORS = [
     :row,
 ]
 
+function _validate_register_assumptions(
+    f::Function,
+    name::Symbol,
+    nb_args::Integer,
+)
+    # Assumption 1: check that `f` can be called with `Float64` arguments.
+    arg = nb_args == 1 ? 0.0 : zeros(nb_args)
+    if hasmethod(f, Tuple{typeof(arg)})
+        y = f(arg)
+    else
+        error(
+            "Unable to register the function :$name.\n\n" *
+            "The function must be able to be called with $nb_args Float64 " *
+            "arguments, but no method was found for this.",
+        )
+    end
+    if !(y isa Real)
+        error(
+            "Expected return type of `Float64` from the user-defined " *
+            "function :$(name), but got `$(typeof(y))`.",
+        )
+    end
+    # Assumption 2: check that `f` can be differentiated using `ForwardDiff`.
+    try
+        if nb_args == 1
+            ForwardDiff.derivative(f, 0.0)
+        else
+            ForwardDiff.gradient(x -> f(x...), zeros(nb_args))
+        end
+    catch err
+        if err isa MethodError
+            error(
+                "Unable to register the function :$name.\n\n" *
+                _FORWARD_DIFF_METHOD_ERROR_HELPER,
+            )
+        end
+        # We hit some other error, perhaps we called a function like log(-1).
+        # Ignore for now, and hope that a useful error is shown to the user
+        # during the solve.
+    end
+    return
+end
+
+function _checked_derivative(f::F, op::Symbol) where {F}
+    return function (x)
+        try
+            return ForwardDiff.derivative(f, x)
+        catch err
+            _intercept_ForwardDiff_MethodError(err, op)
+        end
+    end
+end
+
+"""
+    check_return_type(::Type{T}, ret::S) where {T,S}
+
+Overload this method for new types `S` to throw an informative error if a
+user-defined function returns the type `S` instead of `T`.
+"""
+check_return_type(::Type{T}, ret::T) where {T} = nothing
+
+function check_return_type(::Type{T}, ret) where {T}
+    return error(
+        "Expected return type of $T from a user-defined function, but got " *
+        "$(typeof(ret)).",
+    )
+end
+
+struct _UnivariateOperator{F,F′,F′′}
+    f::F
+    f′::F′
+    f′′::F′′
+    function _UnivariateOperator(
+        f::Function,
+        f′::Function,
+        f′′::Union{Nothing,Function} = nothing,
+    )
+        return new{typeof(f),typeof(f′),typeof(f′′)}(f, f′, f′′)
+    end
+end
+
+function _UnivariateOperator(op::Symbol, f::Function)
+    _validate_register_assumptions(f, op, 1)
+    f′ = _checked_derivative(f, op)
+    return _UnivariateOperator(op, f, f′)
+end
+
+function _UnivariateOperator(op::Symbol, f::Function, f′::Function)
+    try
+        _validate_register_assumptions(f′, op, 1)
+        f′′ = _checked_derivative(f′, op)
+        return _UnivariateOperator(f, f′, f′′)
+    catch
+        return _UnivariateOperator(f, f′, nothing)
+    end
+end
+
+function _UnivariateOperator(::Symbol, f::Function, f′::Function, f′′::Function)
+    return _UnivariateOperator(f, f′, f′′)
+end
+
+struct OperatorRegistry
+    # NODE_CALL_UNIVARIATE
+    univariate_operators::Vector{Symbol}
+    univariate_operator_to_id::Dict{Symbol,Int}
+    univariate_user_operator_start::Int
+    registered_univariate_operators::Vector{_UnivariateOperator}
+    # NODE_CALL_MULTIVARIATE
+    multivariate_operators::Vector{Symbol}
+    multivariate_operator_to_id::Dict{Symbol,Int}
+    multivariate_user_operator_start::Int
+    registered_multivariate_operators::Vector{
+        MOI.Nonlinear._MultivariateOperator,
+    }
+    # NODE_LOGIC
+    logic_operators::Vector{Symbol}
+    logic_operator_to_id::Dict{Symbol,Int}
+    # NODE_COMPARISON
+    comparison_operators::Vector{Symbol}
+    comparison_operator_to_id::Dict{Symbol,Int}
+    function OperatorRegistry()
+        univariate_operators = copy(MOI.Nonlinear.DEFAULT_UNIVARIATE_OPERATORS)
+        multivariate_operators = copy(DEFAULT_MULTIVARIATE_OPERATORS)
+        logic_operators = [:&&, :||]
+        comparison_operators = [:<=, :(==), :>=, :<, :>]
+        return new(
+            # NODE_CALL_UNIVARIATE
+            univariate_operators,
+            Dict{Symbol,Int}(
+                op => i for (i, op) in enumerate(univariate_operators)
+            ),
+            length(univariate_operators),
+            _UnivariateOperator[],
+            # NODE_CALL
+            multivariate_operators,
+            Dict{Symbol,Int}(
+                op => i for (i, op) in enumerate(multivariate_operators)
+            ),
+            length(multivariate_operators),
+            MOI.Nonlinear._MultivariateOperator[],
+            # NODE_LOGIC
+            logic_operators,
+            Dict{Symbol,Int}(op => i for (i, op) in enumerate(logic_operators)),
+            # NODE_COMPARISON
+            comparison_operators,
+            Dict{Symbol,Int}(
+                op => i for (i, op) in enumerate(comparison_operators)
+            ),
+        )
+    end
+end
+
 function eval_logic_function(
     ::OperatorRegistry,
     op::Symbol,
@@ -34,6 +186,23 @@ function eval_logic_function(
     end
 end
 
+function _generate_eval_univariate()
+    exprs = map(Nonlinear.DEFAULT_UNIVARIATE_OPERATORS) do op
+        return :(
+            return (
+                value_deriv_and_second($op, x)[1],
+                value_deriv_and_second($op, x)[2],
+            )
+        )
+    end
+    return Nonlinear._create_binary_switch(1:length(exprs), exprs)
+end
+
+@eval @inline function _eval_univariate(id, x::T) where {T}
+    $(_generate_eval_univariate())
+    return error("Invalid id for univariate operator: $id")
+end
+
 function eval_multivariate_function(
     registry::OperatorRegistry,
     op::Symbol,
@@ -165,17 +334,44 @@ function eval_multivariate_hessian(
     return true
 end
 
+function eval_univariate_function(operator::_UnivariateOperator, x::T) where {T}
+    ret = operator.f(x)
+    check_return_type(T, ret)
+    return ret::T
+end
+
+function eval_univariate_gradient(operator::_UnivariateOperator, x::T) where {T}
+    ret = operator.f′(x)
+    check_return_type(T, ret)
+    return ret::T
+end
+
+function eval_univariate_hessian(operator::_UnivariateOperator, x::T) where {T}
+    ret = operator.f′′(x)
+    check_return_type(T, ret)
+    return ret::T
+end
+
+function eval_univariate_function_and_gradient(
+    operator::_UnivariateOperator,
+    x::T,
+) where {T}
+    ret_f = eval_univariate_function(operator, x)
+    ret_f′ = eval_univariate_gradient(operator, x)
+    return ret_f, ret_f′
+end
+
 function eval_univariate_function_and_gradient(
     registry::OperatorRegistry,
     id::Integer,
     x::T,
 ) where {T}
     if id <= registry.univariate_user_operator_start
-        return Nonlinear._eval_univariate(id, x)::Tuple{T,T}
+        return _eval_univariate(id, x)::Tuple{T,T}
     end
     offset = id - registry.univariate_user_operator_start
     operator = registry.registered_univariate_operators[offset]
-    return Nonlinear.eval_univariate_function_and_gradient(operator, x)
+    return eval_univariate_function_and_gradient(operator, x)
 end
 
 function eval_multivariate_gradient(
diff --git a/src/types.jl b/src/types.jl
index 9f5bd2c..bd735a1 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -133,57 +133,6 @@ struct _FunctionStorage
     end
 end
 
-struct OperatorRegistry
-    # NODE_CALL_UNIVARIATE
-    univariate_operators::Vector{Symbol}
-    univariate_operator_to_id::Dict{Symbol,Int}
-    univariate_user_operator_start::Int
-    registered_univariate_operators::Vector{MOI.Nonlinear._UnivariateOperator}
-    # NODE_CALL_MULTIVARIATE
-    multivariate_operators::Vector{Symbol}
-    multivariate_operator_to_id::Dict{Symbol,Int}
-    multivariate_user_operator_start::Int
-    registered_multivariate_operators::Vector{
-        MOI.Nonlinear._MultivariateOperator,
-    }
-    # NODE_LOGIC
-    logic_operators::Vector{Symbol}
-    logic_operator_to_id::Dict{Symbol,Int}
-    # NODE_COMPARISON
-    comparison_operators::Vector{Symbol}
-    comparison_operator_to_id::Dict{Symbol,Int}
-    function OperatorRegistry()
-        univariate_operators = copy(MOI.Nonlinear.DEFAULT_UNIVARIATE_OPERATORS)
-        multivariate_operators = copy(DEFAULT_MULTIVARIATE_OPERATORS)
-        logic_operators = [:&&, :||]
-        comparison_operators = [:<=, :(==), :>=, :<, :>]
-        return new(
-            # NODE_CALL_UNIVARIATE
-            univariate_operators,
-            Dict{Symbol,Int}(
-                op => i for (i, op) in enumerate(univariate_operators)
-            ),
-            length(univariate_operators),
-            MOI.Nonlinear._UnivariateOperator[],
-            # NODE_CALL
-            multivariate_operators,
-            Dict{Symbol,Int}(
-                op => i for (i, op) in enumerate(multivariate_operators)
-            ),
-            length(multivariate_operators),
-            MOI.Nonlinear._MultivariateOperator[],
-            # NODE_LOGIC
-            logic_operators,
-            Dict{Symbol,Int}(op => i for (i, op) in enumerate(logic_operators)),
-            # NODE_COMPARISON
-            comparison_operators,
-            Dict{Symbol,Int}(
-                op => i for (i, op) in enumerate(comparison_operators)
-            ),
-        )
-    end
-end
-
 """
     Model()
 
diff --git a/src/univariate_expressions.jl b/src/univariate_expressions.jl
new file mode 100644
index 0000000..5ac28c2
--- /dev/null
+++ b/src/univariate_expressions.jl
@@ -0,0 +1,79 @@
+#! format:off
+#
+# This file is machine generated by running univariate_expressions_generator.jl
+
+using SpecialFunctions
+
+value_deriv_and_second(::typeof(+), x) = (+x, one(x), zero(x))
+value_deriv_and_second(::typeof(-), x) = (-x, -one(x), zero(x))
+value_deriv_and_second(::typeof(abs), x) = (abs(x), ifelse(x >= 0, one(x), -one(x)), zero(x))
+value_deriv_and_second(::typeof(sign), x) = (sign(x), zero(x), zero(x))
+value_deriv_and_second(::typeof(sqrt), x) = (sqrt(x), 0.5 / sqrt(x), (0.5 * -(0.5 / sqrt(x))) / sqrt(x) ^ 2)
+value_deriv_and_second(::typeof(cbrt), x) = (cbrt(x), 0.3333333333333333 / cbrt(x) ^ 2, (0.3333333333333333 * -(2 * (0.3333333333333333 / cbrt(x) ^ 2) * cbrt(x))) / (cbrt(x) ^ 2) ^ 2)
+value_deriv_and_second(::typeof(abs2), x) = (abs2(x), 2x, 2 * one(x))
+value_deriv_and_second(::typeof(inv), x) = (inv(x), -(abs2(inv(x))), -(-(abs2(inv(x))) * (2 * inv(x))))
+value_deriv_and_second(::typeof(log), x) = (log(x), 1 / x, -1 / x ^ 2)
+value_deriv_and_second(::typeof(log10), x) = (log10(x), (1 / x) / 2.302585092994046, (-1 / x ^ 2) / 2.302585092994046)
+value_deriv_and_second(::typeof(log2), x) = (log2(x), (1 / x) / 0.6931471805599453, (-1 / x ^ 2) / 0.6931471805599453)
+value_deriv_and_second(::typeof(log1p), x) = (log1p(x), 1 / (1 + x), -1 / (1 + x) ^ 2)
+value_deriv_and_second(::typeof(exp), x) = (exp(x), exp(x), exp(x))
+value_deriv_and_second(::typeof(exp2), x) = (exp2(x), 0.6931471805599453 * exp2(x), 0.6931471805599453 * (0.6931471805599453 * exp2(x)))
+value_deriv_and_second(::typeof(expm1), x) = (expm1(x), exp(x), exp(x))
+value_deriv_and_second(::typeof(sin), x) = (sin(x), cos(x), -(sin(x)))
+value_deriv_and_second(::typeof(cos), x) = (cos(x), -(sin(x)), -(cos(x)))
+value_deriv_and_second(::typeof(tan), x) = (tan(x), 1 + tan(x) ^ 2, 2 * (1 + tan(x) ^ 2) * tan(x))
+value_deriv_and_second(::typeof(sec), x) = (sec(x), sec(x) * tan(x), (sec(x) * tan(x)) * tan(x) + sec(x) * (1 + tan(x) ^ 2))
+value_deriv_and_second(::typeof(csc), x) = (csc(x), -(csc(x)) * cot(x), -(-(csc(x)) * cot(x)) * cot(x) + -(csc(x)) * -((1 + cot(x) ^ 2)))
+value_deriv_and_second(::typeof(cot), x) = (cot(x), -((1 + cot(x) ^ 2)), -(2 * -((1 + cot(x) ^ 2)) * cot(x)))
+value_deriv_and_second(::typeof(sind), x) = (sind(x), 0.017453292519943295 * cosd(x), 0.017453292519943295 * (-0.017453292519943295 * sind(x)))
+value_deriv_and_second(::typeof(cosd), x) = (cosd(x), -0.017453292519943295 * sind(x), -0.017453292519943295 * (0.017453292519943295 * cosd(x)))
+value_deriv_and_second(::typeof(tand), x) = (tand(x), 0.017453292519943295 * (1 + tand(x) ^ 2), 0.017453292519943295 * (2 * (0.017453292519943295 * (1 + tand(x) ^ 2)) * tand(x)))
+value_deriv_and_second(::typeof(secd), x) = (secd(x), 0.017453292519943295 * secd(x) * tand(x), 0.017453292519943295 * (0.017453292519943295 * secd(x) * tand(x)) * tand(x) + 0.017453292519943295 * secd(x) * (0.017453292519943295 * (1 + tand(x) ^ 2)))
+value_deriv_and_second(::typeof(cscd), x) = (cscd(x), -0.017453292519943295 * cscd(x) * cotd(x), -0.017453292519943295 * (-0.017453292519943295 * cscd(x) * cotd(x)) * cotd(x) + -0.017453292519943295 * cscd(x) * (-0.017453292519943295 * (1 + cotd(x) ^ 2)))
+value_deriv_and_second(::typeof(cotd), x) = (cotd(x), -0.017453292519943295 * (1 + cotd(x) ^ 2), -0.017453292519943295 * (2 * (-0.017453292519943295 * (1 + cotd(x) ^ 2)) * cotd(x)))
+value_deriv_and_second(::typeof(asin), x) = (asin(x), 1 / sqrt(1 - x ^ 2), -(-(2x) * (0.5 / sqrt(1 - x ^ 2))) / sqrt(1 - x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(acos), x) = (acos(x), -1 / sqrt(1 - x ^ 2), -(-(-(2x) * (0.5 / sqrt(1 - x ^ 2)))) / sqrt(1 - x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(atan), x) = (atan(x), 1 / (1 + x ^ 2), -(2x) / (1 + x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(asec), x) = (asec(x), (1 / abs(x)) / sqrt(x ^ 2 - 1), nothing)
+value_deriv_and_second(::typeof(acsc), x) = (acsc(x), (-1 / abs(x)) / sqrt(x ^ 2 - 1), nothing)
+value_deriv_and_second(::typeof(acot), x) = (acot(x), -1 / (1 + x ^ 2), -(-(2x)) / (1 + x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(asind), x) = (asind(x), 57.29577951308232 / sqrt(1 - x ^ 2), (57.29577951308232 * -(-(2x) * (0.5 / sqrt(1 - x ^ 2)))) / sqrt(1 - x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(acosd), x) = (acosd(x), -57.29577951308232 / sqrt(1 - x ^ 2), (-57.29577951308232 * -(-(2x) * (0.5 / sqrt(1 - x ^ 2)))) / sqrt(1 - x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(atand), x) = (atand(x), 57.29577951308232 / (1 + x ^ 2), (57.29577951308232 * -(2x)) / (1 + x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(asecd), x) = (asecd(x), (57.29577951308232 / abs(x)) / sqrt(x ^ 2 - 1), nothing)
+value_deriv_and_second(::typeof(acscd), x) = (acscd(x), (-57.29577951308232 / abs(x)) / sqrt(x ^ 2 - 1), nothing)
+value_deriv_and_second(::typeof(acotd), x) = (acotd(x), -57.29577951308232 / (1 + x ^ 2), (-57.29577951308232 * -(2x)) / (1 + x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(sinh), x) = (sinh(x), cosh(x), sinh(x))
+value_deriv_and_second(::typeof(cosh), x) = (cosh(x), sinh(x), cosh(x))
+value_deriv_and_second(::typeof(tanh), x) = (tanh(x), sech(x) ^ 2, 2 * (-(tanh(x)) * sech(x)) * sech(x))
+value_deriv_and_second(::typeof(sech), x) = (sech(x), -(tanh(x)) * sech(x), -(sech(x) ^ 2) * sech(x) + -(tanh(x)) * (-(tanh(x)) * sech(x)))
+value_deriv_and_second(::typeof(csch), x) = (csch(x), -(coth(x)) * csch(x), -(-(csch(x) ^ 2)) * csch(x) + -(coth(x)) * (-(coth(x)) * csch(x)))
+value_deriv_and_second(::typeof(coth), x) = (coth(x), -(csch(x) ^ 2), -(2 * (-(coth(x)) * csch(x)) * csch(x)))
+value_deriv_and_second(::typeof(asinh), x) = (asinh(x), 1 / sqrt(1 + x ^ 2), -((2x) * (0.5 / sqrt(1 + x ^ 2))) / sqrt(1 + x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(acosh), x) = (acosh(x), 1 / sqrt(x ^ 2 - 1), -((2x) * (0.5 / sqrt(x ^ 2 - 1))) / sqrt(x ^ 2 - 1) ^ 2)
+value_deriv_and_second(::typeof(atanh), x) = (atanh(x), 1 / (1 - x ^ 2), -(-(2x)) / (1 - x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(asech), x) = (asech(x), (-1 / x) / sqrt(1 - x ^ 2), ((1 / x ^ 2) * sqrt(1 - x ^ 2) - (-1 / x) * (-(2x) * (0.5 / sqrt(1 - x ^ 2)))) / sqrt(1 - x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(acsch), x) = (acsch(x), (-1 / abs(x)) / sqrt(1 + x ^ 2), nothing)
+value_deriv_and_second(::typeof(acoth), x) = (acoth(x), 1 / (1 - x ^ 2), -(-(2x)) / (1 - x ^ 2) ^ 2)
+value_deriv_and_second(::typeof(deg2rad), x) = (deg2rad(x), 0.017453292519943295 * one(x), 0 * one(x))
+value_deriv_and_second(::typeof(rad2deg), x) = (rad2deg(x), 57.29577951308232 * one(x), 0 * one(x))
+value_deriv_and_second(::typeof(erf), x) = (erf(x), (2 * exp(-x * x)) / 1.7724538509055159, (2 * ((-x + -x) * exp(-x * x))) / 1.7724538509055159)
+value_deriv_and_second(::typeof(erfinv), x) = (erfinv(x), 0.8862269254527579 * exp(erfinv(x) * erfinv(x)), 0.8862269254527579 * (((0.8862269254527579 * exp(erfinv(x) * erfinv(x))) * erfinv(x) + erfinv(x) * (0.8862269254527579 * exp(erfinv(x) * erfinv(x)))) * exp(erfinv(x) * erfinv(x))))
+value_deriv_and_second(::typeof(erfc), x) = (erfc(x), (-2 * exp(-x * x)) / 1.7724538509055159, (-2 * ((-x + -x) * exp(-x * x))) / 1.7724538509055159)
+value_deriv_and_second(::typeof(erfcinv), x) = (erfcinv(x), -0.8862269254527579 * exp(erfcinv(x) * erfcinv(x)), -0.8862269254527579 * (((-0.8862269254527579 * exp(erfcinv(x) * erfcinv(x))) * erfcinv(x) + erfcinv(x) * (-0.8862269254527579 * exp(erfcinv(x) * erfcinv(x)))) * exp(erfcinv(x) * erfcinv(x))))
+value_deriv_and_second(::typeof(erfi), x) = (erfi(x), (2 * exp(x * x)) / 1.7724538509055159, (2 * ((x + x) * exp(x * x))) / 1.7724538509055159)
+value_deriv_and_second(::typeof(gamma), x) = (gamma(x), digamma(x) * gamma(x), trigamma(x) * gamma(x) + digamma(x) * (digamma(x) * gamma(x)))
+value_deriv_and_second(::typeof(lgamma), x) = (lgamma(x), digamma(x), trigamma(x))
+value_deriv_and_second(::typeof(digamma), x) = (digamma(x), trigamma(x), polygamma(2, x))
+value_deriv_and_second(::typeof(invdigamma), x) = (invdigamma(x), inv(trigamma(invdigamma(x))), (inv(trigamma(invdigamma(x))) * polygamma(2, invdigamma(x))) * -(abs2(inv(trigamma(invdigamma(x))))))
+value_deriv_and_second(::typeof(trigamma), x) = (trigamma(x), polygamma(2, x), nothing)
+value_deriv_and_second(::typeof(airyai), x) = (airyai(x), airyaiprime(x), x * airyai(x))
+value_deriv_and_second(::typeof(airybi), x) = (airybi(x), airybiprime(x), x * airybi(x))
+value_deriv_and_second(::typeof(airyaiprime), x) = (airyaiprime(x), x * airyai(x), airyai(x) + x * airyaiprime(x))
+value_deriv_and_second(::typeof(airybiprime), x) = (airybiprime(x), x * airybi(x), airybi(x) + x * airybiprime(x))
+value_deriv_and_second(::typeof(besselj0), x) = (besselj0(x), -(besselj1(x)), -((besselj0(x) - besselj(2, x)) / 2))
+value_deriv_and_second(::typeof(besselj1), x) = (besselj1(x), (besselj0(x) - besselj(2, x)) / 2, (-(besselj1(x)) - (besselj(1, x) - besselj(3, x)) / 2) / 2)
+value_deriv_and_second(::typeof(bessely0), x) = (bessely0(x), -(bessely1(x)), -((bessely0(x) - bessely(2, x)) / 2))
+value_deriv_and_second(::typeof(bessely1), x) = (bessely1(x), (bessely0(x) - bessely(2, x)) / 2, (-(bessely1(x)) - (bessely(1, x) - bessely(3, x)) / 2) / 2)
+value_deriv_and_second(::typeof(erfcx), x) = (erfcx(x), 2 * x * erfcx(x) - 1.1283791670955126, 2 * erfcx(x) + 2 * x * (2 * x * erfcx(x) - 1.1283791670955126))
+value_deriv_and_second(::typeof(dawson), x) = (dawson(x), 1 - (2x) * dawson(x), -((2 * dawson(x) + (2x) * (1 - (2x) * dawson(x)))))
diff --git a/src/univariate_expressions_generator.jl b/src/univariate_expressions_generator.jl
new file mode 100644
index 0000000..36f6454
--- /dev/null
+++ b/src/univariate_expressions_generator.jl
@@ -0,0 +1,76 @@
+# Copyright (c) 2017: Miles Lubin and contributors
+# Copyright (c) 2017: Google Inc.
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+# JuMP requires first and second-order symbolic derivatives for univariate
+# functions. Calculus.jl is a package for computing them, and JuMP 1.0.0 used to
+# dynamically compute them during precompilation. The benefit of this approach
+# is that if the rules in Calculus are updated, then hey will automatically flow
+# into JuMP. However, the Calculus package has not seen any development for two
+# years, and JuMP uses only a very small part of the total package. Therefore,
+# this script statically builds the list of univariate functions supported by
+# JuMP and writes them to univariate_expressions.jl. This has the benefit of
+# removing a dependency, and saves the recalculation of these derivatives by
+# every user.
+#
+# If you want to rebuild this list of derivatives in future, you will need v0.5
+# of Calculus: `import Pkg; Pkg.pkg"add Calculus@0.5"`.
+
+import Calculus
+
+function _differentiate(f)
+    try
+        return Calculus.simplify(Calculus.differentiate(f, :x))
+    catch
+        return nothing
+    end
+end
+
+_to_expr(f) = f === nothing ? :(nothing) : f isa Number ? :($f * one(x)) : f
+
+open("src/univariate_expressions.jl", "w") do io
+    print(
+        io,
+        """
+#! format:off
+#
+# This file is machine generated by running univariate_expressions_generator.jl
+
+""",
+    )
+
+    println(
+        io,
+        "value_deriv_and_second(::typeof(+), x) = (+x, one(x), zero(x))",
+    )
+    println(
+        io,
+        "value_deriv_and_second(::typeof(-), x) = (-x, -one(x), zero(x))",
+    )
+    println(
+        io,
+        "value_deriv_and_second(::typeof(abs), x) = (abs(x), ifelse(x >= 0, one(x), -one(x)), zero(x))",
+    )
+    println(
+        io,
+        "value_deriv_and_second(::typeof(sign), x) = (sign(x), zero(x), zero(x))",
+    )
+
+    for (op, _) in Calculus.symbolic_derivatives_1arg()
+        f = Expr(:call, op, :x)
+        df = _differentiate(f)
+        d2f = _differentiate(df)
+        println(
+            io,
+            "value_deriv_and_second(::typeof($op), x) = (",
+            _to_expr(f),
+            ", ",
+            _to_expr(df),
+            ", ",
+            _to_expr(d2f),
+            ")",
+        )
+    end
+end
diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 1af6554..86993c1 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -529,6 +529,25 @@ function test_objective_norm_of_mtx_vector_product()
     return
 end
 
+function test_objective_univariate_operator()
+    model = ArrayDiff.Model()
+    x = MOI.VariableIndex(1)
+    ArrayDiff.set_objective(model, :(sin($x)))
+    evaluator = ArrayDiff.Evaluator(model, ArrayDiff.Mode(), [x])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 0]
+    @test sizes.size_offset == [0, 0]
+    @test sizes.size == []
+    @test sizes.storage_offset == [0, 1, 2]
+    x = [pi / 4]
+    @test MOI.eval_objective(evaluator, x) ≈ sqrt(2) / 2
+    g = ones(1)
+    MOI.eval_objective_gradient(evaluator, g, x)
+    @test g[1] ≈ cos(pi / 4)
+    return
+end
+
 end  # module
 
 TestArrayDiff.runtests()
diff --git a/test/Project.toml b/test/Project.toml
index d66f113..d768ff2 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,7 +1,10 @@
 [deps]
 ArrayDiff = "c45fa1ca-6901-44ac-ae5b-5513a4852d50"
+Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
\ No newline at end of file