diff --git a/docs/source/gradient_solver.rst b/docs/source/gradient_solver.rst
index 1356e74..dde9d7e 100644
--- a/docs/source/gradient_solver.rst
+++ b/docs/source/gradient_solver.rst
@@ -33,10 +33,10 @@
 .. function:: bool FirstOrderFunction::Evaluate(const double* const parameters, double* cost, double* gradient) const
 
    Evaluate the cost/value of the function. If ``gradient`` is not
-   ``NULL`` then evaluate the gradient too. If evaluation is
+   ``nullptr`` then evaluate the gradient too. If evaluation is
    successful return, ``true`` else return ``false``.
 
-   ``cost`` guaranteed to be never ``NULL``, ``gradient`` can be ``NULL``.
+   ``cost`` guaranteed to be never ``nullptr``, ``gradient`` can be ``nullptr``.
 
 .. function:: int FirstOrderFunction::NumParameters() const
 
diff --git a/docs/source/gradient_tutorial.rst b/docs/source/gradient_tutorial.rst
index 0bbdee4..3fef6b6 100644
--- a/docs/source/gradient_tutorial.rst
+++ b/docs/source/gradient_tutorial.rst
@@ -40,7 +40,7 @@
       const double y = parameters[1];
 
       cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
-      if (gradient != NULL) {
+      if (gradient != nullptr) {
         gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x;
         gradient[1] = 200.0 * (y - x * x);
       }
diff --git a/docs/source/interfacing_with_autodiff.rst b/docs/source/interfacing_with_autodiff.rst
index b79ed45..02f58b2 100644
--- a/docs/source/interfacing_with_autodiff.rst
+++ b/docs/source/interfacing_with_autodiff.rst
@@ -181,7 +181,7 @@
                            double* residuals,
                            double** jacobians) const {
        if (!jacobians) {
-         ComputeDistortionValueAndJacobian(parameters[0][0], residuals, NULL);
+         ComputeDistortionValueAndJacobian(parameters[0][0], residuals, nullptr);
        } else {
          ComputeDistortionValueAndJacobian(parameters[0][0], residuals, jacobians[0]);
        }
diff --git a/docs/source/nnls_modeling.rst b/docs/source/nnls_modeling.rst
index 257c5eb..593fcab 100644
--- a/docs/source/nnls_modeling.rst
+++ b/docs/source/nnls_modeling.rst
@@ -108,29 +108,29 @@
    that contains the :math:`i^{\text{th}}` parameter block that the
    ``CostFunction`` depends on.
 
-   ``parameters`` is never ``NULL``.
+   ``parameters`` is never ``nullptr``.
 
    ``residuals`` is an array of size ``num_residuals_``.
 
-   ``residuals`` is never ``NULL``.
+   ``residuals`` is never ``nullptr``.
 
    ``jacobians`` is an array of arrays of size
    ``CostFunction::parameter_block_sizes_.size()``.
 
-   If ``jacobians`` is ``NULL``, the user is only expected to compute
+   If ``jacobians`` is ``nullptr``, the user is only expected to compute
    the residuals.
 
    ``jacobians[i]`` is a row-major array of size ``num_residuals x
    parameter_block_sizes_[i]``.
 
-   If ``jacobians[i]`` is **not** ``NULL``, the user is required to
+   If ``jacobians[i]`` is **not** ``nullptr``, the user is required to
    compute the Jacobian of the residual vector with respect to
    ``parameters[i]`` and store it in this array, i.e.
 
    ``jacobians[i][r * parameter_block_sizes_[i] + c]`` =
    :math:`\frac{\displaystyle \partial \text{residual}[r]}{\displaystyle \partial \text{parameters}[i][c]}`
 
-   If ``jacobians[i]`` is ``NULL``, then this computation can be
+   If ``jacobians[i]`` is ``nullptr``, then this computation can be
    skipped. This is the case when the corresponding parameter block is
    marked constant.
 
@@ -914,7 +914,7 @@
 
        std::vector<LocalParameterization*> local_parameterizations;
        local_parameterizations.push_back(my_parameterization);
-       local_parameterizations.push_back(NULL);
+       local_parameterizations.push_back(nullptr);
 
        std::vector parameter1;
        std::vector parameter2;
@@ -1109,8 +1109,8 @@
    Given a loss function :math:`\rho(s)` and a scalar :math:`a`, :class:`ScaledLoss`
    implements the function :math:`a \rho(s)`.
 
-   Since we treat a ``NULL`` Loss function as the Identity loss
-   function, :math:`rho` = ``NULL``: is a valid input and will result
+   Since we treat a ``nullptr`` Loss function as the Identity loss
+   function, :math:`rho` = ``nullptr``: is a valid input and will result
    in the input being scaled by :math:`a`. This provides a simple way
    of implementing a scaled ResidualBlock.
 
@@ -1587,7 +1587,7 @@
    the parameter blocks it expects. The function checks that these
    match the sizes of the parameter blocks listed in
    ``parameter_blocks``. The program aborts if a mismatch is
-   detected. ``loss_function`` can be ``NULL``, in which case the cost
+   detected. ``loss_function`` can be ``nullptr``, in which case the cost
    of the term is just the squared norm of the residuals.
 
    The user has the option of explicitly adding the parameter blocks
@@ -1752,7 +1752,7 @@
    parameter blocks it expects. The function checks that these match
    the sizes of the parameter blocks listed in parameter_blocks. The
    program aborts if a mismatch is detected. loss_function can be
-   NULL, in which case the cost of the term is just the squared norm
+   nullptr, in which case the cost of the term is just the squared norm
    of the residuals.
 
    The parameter blocks may be passed together as a
@@ -1791,10 +1791,10 @@
 
       Problem problem;
 
-      problem.AddResidualBlock(new MyUnaryCostFunction(...), NULL, x1);
-      problem.AddResidualBlock(new MyBinaryCostFunction(...), NULL, x2, x1);
-      problem.AddResidualBlock(new MyUnaryCostFunction(...), NULL, v1);
-      problem.AddResidualBlock(new MyBinaryCostFunction(...), NULL, v2);
+      problem.AddResidualBlock(new MyUnaryCostFunction(...), nullptr, x1);
+      problem.AddResidualBlock(new MyBinaryCostFunction(...), nullptr, x2, x1);
+      problem.AddResidualBlock(new MyUnaryCostFunction(...), nullptr, v1);
+      problem.AddResidualBlock(new MyBinaryCostFunction(...), nullptr, v2);
 
 .. function:: void Problem::AddParameterBlock(double* values, int size, LocalParameterization* local_parameterization)
 
@@ -1871,7 +1871,7 @@
 
    Get the local parameterization object associated with this
    parameter block. If there is no parameterization object associated
-   then `NULL` is returned
+   then `nullptr` is returned
 
 .. function:: void Problem::SetParameterLowerBound(double* values, int index, double lower_bound)
 
@@ -2018,7 +2018,7 @@
 .. function:: bool Problem::Evaluate(const Problem::EvaluateOptions& options, double* cost, vector<double>* residuals, vector<double>* gradient, CRSMatrix* jacobian)
 
    Evaluate a :class:`Problem`. Any of the output pointers can be
-   `NULL`. Which residual blocks and parameter blocks are used is
+   `nullptr`. Which residual blocks and parameter blocks are used is
    controlled by the :class:`Problem::EvaluateOptions` struct below.
 
    .. NOTE::
@@ -2032,10 +2032,10 @@
 
         Problem problem;
         double x = 1;
-        problem.Add(new MyCostFunction, NULL, &x);
+        problem.Add(new MyCostFunction, nullptr, &x);
 
         double cost = 0.0;
-        problem.Evaluate(Problem::EvaluateOptions(), &cost, NULL, NULL, NULL);
+        problem.Evaluate(Problem::EvaluateOptions(), &cost, nullptr, nullptr, nullptr);
 
       The cost is evaluated at `x = 1`. If you wish to evaluate the
       problem at `x = 2`, then
@@ -2043,7 +2043,7 @@
       .. code-block:: c++
 
          x = 2;
-         problem.Evaluate(Problem::EvaluateOptions(), &cost, NULL, NULL, NULL);
+         problem.Evaluate(Problem::EvaluateOptions(), &cost, nullptr, nullptr, nullptr);
 
       is the way to do so.
 
diff --git a/docs/source/nnls_tutorial.rst b/docs/source/nnls_tutorial.rst
index fb7059c..6c89032 100644
--- a/docs/source/nnls_tutorial.rst
+++ b/docs/source/nnls_tutorial.rst
@@ -111,7 +111,7 @@
      // auto-differentiation to obtain the derivative (jacobian).
      CostFunction* cost_function =
          new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
-     problem.AddResidualBlock(cost_function, NULL, &x);
+     problem.AddResidualBlock(cost_function, nullptr, &x);
 
      // Run the solver!
      Solver::Options options;
@@ -212,7 +212,7 @@
   CostFunction* cost_function =
     new NumericDiffCostFunction<NumericDiffCostFunctor, ceres::CENTRAL, 1, 1>(
         new NumericDiffCostFunctor);
-  problem.AddResidualBlock(cost_function, NULL, &x);
+  problem.AddResidualBlock(cost_function, nullptr, &x);
 
 Notice the parallel from when we were using automatic differentiation
 
@@ -220,7 +220,7 @@
 
   CostFunction* cost_function =
       new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
-  problem.AddResidualBlock(cost_function, NULL, &x);
+  problem.AddResidualBlock(cost_function, nullptr, &x);
 
 The construction looks almost identical to the one used for automatic
 differentiation, except for an extra template parameter that indicates
@@ -261,7 +261,7 @@
       residuals[0] = 10 - x;
 
       // Compute the Jacobian if asked for.
-      if (jacobians != NULL && jacobians[0] != NULL) {
+      if (jacobians != nullptr && jacobians[0] != nullptr) {
         jacobians[0][0] = -1;
       }
       return true;
@@ -358,13 +358,13 @@
   // Add residual terms to the problem using the using the autodiff
   // wrapper to get the derivatives automatically.
   problem.AddResidualBlock(
-    new AutoDiffCostFunction<F1, 1, 1, 1>(new F1), NULL, &x1, &x2);
+    new AutoDiffCostFunction<F1, 1, 1, 1>(new F1), nullptr, &x1, &x2);
   problem.AddResidualBlock(
-    new AutoDiffCostFunction<F2, 1, 1, 1>(new F2), NULL, &x3, &x4);
+    new AutoDiffCostFunction<F2, 1, 1, 1>(new F2), nullptr, &x3, &x4);
   problem.AddResidualBlock(
-    new AutoDiffCostFunction<F3, 1, 1, 1>(new F3), NULL, &x2, &x3)
+    new AutoDiffCostFunction<F3, 1, 1, 1>(new F3), nullptr, &x2, &x3)
   problem.AddResidualBlock(
-    new AutoDiffCostFunction<F4, 1, 1, 1>(new F4), NULL, &x1, &x4);
+    new AutoDiffCostFunction<F4, 1, 1, 1>(new F4), nullptr, &x1, &x4);
 
 
 Note that each ``ResidualBlock`` only depends on the two parameters
@@ -496,7 +496,7 @@
    CostFunction* cost_function =
         new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
             new ExponentialResidual(data[2 * i], data[2 * i + 1]));
-   problem.AddResidualBlock(cost_function, NULL, &m, &c);
+   problem.AddResidualBlock(cost_function, nullptr, &m, &c);
  }
 
 Compiling and running `examples/curve_fitting.cc
@@ -568,7 +568,7 @@
 
 .. code-block:: c++
 
-   problem.AddResidualBlock(cost_function, NULL , &m, &c);
+   problem.AddResidualBlock(cost_function, nullptr , &m, &c);
 
 to
 
@@ -697,7 +697,7 @@
             bal_problem.observations()[2 * i + 0],
             bal_problem.observations()[2 * i + 1]);
    problem.AddResidualBlock(cost_function,
-                            NULL /* squared loss */,
+                            nullptr /* squared loss */,
                             bal_problem.mutable_camera_for_observation(i),
                             bal_problem.mutable_point_for_observation(i));
  }