diff --git a/pub/lib/expr.cpp b/pub/lib/expr.cpp
index 1548b73b75c26c88954887b3cf6d2c2e5bf65b54..1aa5316db17528afab2960360436d857194cb3d0 100644
--- a/pub/lib/expr.cpp
+++ b/pub/lib/expr.cpp
@@ -236,16 +236,16 @@ RObjVecInt CNode::to_index_list( int siz ) const
         PObjVecInt ret( new CObjVecInt );
         int ival = ri->to_i();
         if( ival<0 || ival>=siz )
-            return NULL;
+            return nullptr;
         ret->v.push_back( ival );
         return ret;
     } else if ( RObjVecInt riv = PCAST<const CObjVecInt>( val ) ) {
         for( int ival: riv->v )
             if( ival<0 || ival>=siz )
-                return NULL;
+                return nullptr;
         return riv;
     } else {
-        return NULL;
+        return nullptr;
     }
 }
 
diff --git a/pub/lib/fbase.cpp b/pub/lib/fbase.cpp
index 2c2d61d4d8e11c8bec17173fbbacf40614b8218c..fa9cbe39a1bc9fb60f65ba72726dd0a8fe42d786 100644
--- a/pub/lib/fbase.cpp
+++ b/pub/lib/fbase.cpp
@@ -28,7 +28,7 @@
 const double PI    = M_PI;    //4*atan(1.);
 const double TWOPI = 2*M_PI; //8*atan(1.);
 
-static triv::RandomNumberGenerator *rng = NULL;
+static triv::RandomNumberGenerator *rng = nullptr;
 
 
 //**************************************************************************************************
diff --git a/pub/lib/fit.cpp b/pub/lib/fit.cpp
index c93f4491579579a66fa88690b80f78409c7c18c1..e5f3939a9b43e1f46a66a413778f4643be80db4a 100644
--- a/pub/lib/fit.cpp
+++ b/pub/lib/fit.cpp
@@ -148,7 +148,7 @@ static void fit_evaluate( const double* par, int m_dat, const void *data,
     try {
         FitDatTyp *mydata;
         mydata = (FitDatTyp*) data;
-        if( mydata->timeout && time(NULL) > mydata->timeout )
+        if( mydata->timeout && time(nullptr) > mydata->timeout )
             throw S("reached fit time out");
         POlc fc = mydata->fc;
         ROld fd = mydata->fd;
@@ -211,7 +211,7 @@ string NCurveFit::fit_one_spec( POlc fc, ROld fd, int k, int j, const lm_control
 
         // Levenberg-Marquardt routine from library lmfit:
 
-        data.timeout = time(NULL) + omp_get_num_threads()*maxtime;
+        data.timeout = time(nullptr) + omp_get_num_threads()*maxtime;
 
         lm_status_struct status;
         lmmin( npfree, &(Par[0]), nd, &data, fit_evaluate, &control, &status );
@@ -281,7 +281,7 @@ static void fit_evaluate_glo( const double* par, int m_dat, const void *data,
 {
     try {
         FitGloDatTyp* mydata = (FitGloDatTyp*) data;
-        if( mydata->timeout && time(NULL) > mydata->timeout )
+        if( mydata->timeout && time(nullptr) > mydata->timeout )
             throw S("reached fit time out");
         POlc fc = mydata->fc;
         ROld fd = mydata->fd;
@@ -372,7 +372,7 @@ void NCurveFit::fit_global( POlc fc, ROld fd, int k, const lm_control_struct& co
 
     // Levenberg-Marquardt routine from library lmfit:
 
-    data.timeout = time(NULL) + fc->nJ()*maxtime;
+    data.timeout = time(nullptr) + fc->nJ()*maxtime;
 
     lm_status_struct status;
     lmmin( npfree, &(Par[0]), nd, &data, fit_evaluate_glo, &control, &status );
diff --git a/pub/lib/fregistry.cpp b/pub/lib/fregistry.cpp
index d1e835c91a8bdf6069143800acc4a0c81bee43c8..5ca2b0d6c3c8b5589d5e786a5332c3ef090015bd 100644
--- a/pub/lib/fregistry.cpp
+++ b/pub/lib/fregistry.cpp
@@ -17,12 +17,12 @@
 #include "coord.hpp"
 
 
-//! Returns a pointer to the CFunc of given name, or NULL.
+//! Returns a pointer to the CFunc of given name, or nullptr.
 
 const CFunc* SFuncRegistry::find( string tag ) const
 {
     auto pos = FMap.find(tag);
-    return pos == FMap.end() ? NULL : pos->second;
+    return pos == FMap.end() ? nullptr : pos->second;
 }
 
 
diff --git a/pub/lib/geni.cpp b/pub/lib/geni.cpp
index acdf9a0d332db6d5aff5c8daae6eb4f9b44ef950..1146584fa8ac8f871db4f31fba2151e8d151c4c3 100644
--- a/pub/lib/geni.cpp
+++ b/pub/lib/geni.cpp
@@ -271,7 +271,7 @@ CGeni::CGeni( string _name, int _narg, string _default_arg_str, geni_eval _eval,
 const CGeni* NGeni::find( string key )
 {
     auto pos = gmap.find(key);
-    return pos == gmap.end() ? NULL : pos->second;
+    return pos == gmap.end() ? nullptr : pos->second;
 }
 
 void CGeni::register_me() const
diff --git a/pub/lib/integrate.cpp b/pub/lib/integrate.cpp
index bac1cf624eafe920520511ae2f56051d998e4151..12cf9be18fd123a651b085059535f0ee022209aa 100644
--- a/pub/lib/integrate.cpp
+++ b/pub/lib/integrate.cpp
@@ -135,7 +135,7 @@ CCvin::CCvin( string _txt, int _narg, string _def, cvin_eval _eval, string _com
 const CCvin* NCvin::find( string key )
 {
     auto pos = hmap.find(key);
-    return pos == hmap.end() ? NULL : pos->second;
+    return pos == hmap.end() ? nullptr : pos->second;
 }
 
 void CCvin::register_me() const
diff --git a/pub/lib/node.cpp b/pub/lib/node.cpp
index 88e1703e93de29a936b5690c7ba6603f1665cf01..97474ec3d1ee15548ae6edbcb2e49f242c6057cc 100644
--- a/pub/lib/node.cpp
+++ b/pub/lib/node.cpp
@@ -91,7 +91,7 @@ CNodeFun3::CNodeFun3( const CFunc* _fu, RNode _a0, RNode _a1, RNode _a2 )
 RObj CNodeFun::tree_val( const CContext& ctx ) const
 {
     string basetypes;
-    const CTypedFunc* tf = NULL;
+    const CTypedFunc* tf = nullptr;
     try {
          // evaluate arguments, determine if all results are scalar.
         vector<RObj> pa(narg);
@@ -119,7 +119,7 @@ RObj CNodeFun::tree_val( const CContext& ctx ) const
                     (*(func_0_s)(f))( pa[0]->to_s() );
                 else
                     throw S("BUG: unexpected intypes");
-                return RObjStr( new CObjStr( "__NULL__" ) );
+                return RObjStr( new CObjStr( "__nullptr__" ) );
             } else if ( tf->outtype=="i" ) {
                 int val;
                 if       ( tf->intypes=="i" )
@@ -535,7 +535,7 @@ RObj CNodeGeni::tree_val( const CContext& ctx ) const
         }
         // Now evaluate the function.
         double val;
-        geni->eval( &val, NULL, n, pav ); // TODO: compute errors
+        geni->eval( &val, nullptr, n, pav ); // TODO: compute errors
         return RObjDbl( new CObjDbl( val ) );
 
     } catch( string& ex ) {
@@ -625,7 +625,7 @@ RObj CNodeCvin::tree_val( const CContext& ctx ) const
         int j = PCAST<const CRef2>(ref)->get_j( ctx, f->nJ() );
         // now evaluate the function
         double val;
-        cvin->eval( &val, NULL, fc, k, j, a ); // TODO: compute errors?
+        cvin->eval( &val, nullptr, fc, k, j, a ); // TODO: compute errors?
         return PObjDbl( new CObjDbl( val ) );
     } catch( string& ex ) {
         throw ex+" in NodeCvin "+tree_info()+" {context: "+ctx.context_info()+"}";
@@ -667,7 +667,7 @@ CNodeVal::CNodeVal( int _val ) : val( PObjInt( new CObjInt( _val ) ) ) {}
 
 bool CNodeVal::looks_like_indices() const
 {
-    return PCAST<const CObjInt>(val)!=NULL;
+    return PCAST<const CObjInt>(val)!=nullptr;
 }
 
 RObj CNodeVal::tree_val( const CContext& ctx ) const
@@ -1087,7 +1087,7 @@ RObjVecDbl CNodeConv::convolve( const CContext& ctx, double theshift, const RSpe
        ret = dret;
     } else {
         ret = PObjVecDbl( new CObjVecDbl(n) );
-        dret = NULL;
+        dret = nullptr;
     }
     // Convolute theory with resolution.
     if( conv_step ){ // accelerated version with equidistant grids
@@ -1180,7 +1180,7 @@ RObjVecDbl CNodePConv::convolve( const CContext& ctx, double theshift, const RSp
        }
     } else {
         ret = PObjVecDbl( new CObjVecDbl(n) );
-        dret = NULL;
+        dret = nullptr;
     }
     // Convolute theory with resolution.
     if( conv_step ){ // accelerated version with equidistant grids
@@ -1294,7 +1294,7 @@ RObjVecDbl CNodeDirac::convolve( const CContext& ctx, double theshift, const RSp
         ret = dret;
     } else {
         ret = PObjVecDbl( new CObjVecDbl(n) );
-        dret = NULL;
+        dret = nullptr;
     }
     // Delta function copies resolution.
     vector<double> tt(n), yy(n), dd(n);
@@ -1303,7 +1303,7 @@ RObjVecDbl CNodeDirac::convolve( const CContext& ctx, double theshift, const RSp
         tt[i] = (*(ctx.vt))[i]-theshift;
     // interpolation
     try {
-        sv->intpol( tt, &yy, dret ? &dd : NULL );
+        sv->intpol( tt, &yy, dret ? &dd : nullptr );
     } catch( string& s ) {
         throw "error while evaluating delta function with shift "+S(theshift)+": "+s;
     }
diff --git a/pub/lib/olf.cpp b/pub/lib/olf.cpp
index be1ff3457ac4e78d2a375849192e3577f53262a3..eabf3a6cdf07fdb9b37feca7c6664454020e49a5 100644
--- a/pub/lib/olf.cpp
+++ b/pub/lib/olf.cpp
@@ -555,7 +555,7 @@ RObjVec COlc::eval_curve_syscall( const vector<double>& vt, int k, int j,
         if ( S.size()<2 )
             throw errmsg + " returns only " + S(S.size()) + " valid lines";
         //cout << "DEB cvv: intpol "<<vt.front()<<" .. "<<vt.back()<<"\n";
-        S.intpol( vt, &v, NULL );
+        S.intpol( vt, &v, nullptr );
     } else
         throw S("BUG: invalid evaMode");
     FF.close();
diff --git a/pub/lib/reg.cpp b/pub/lib/reg.cpp
index 7a8192ba708d63757ab4be4d79f46237a99e2f68..0244c6987480983624aa1a83d1f70689265957ea 100644
--- a/pub/lib/reg.cpp
+++ b/pub/lib/reg.cpp
@@ -31,7 +31,7 @@ double NReg::lastresult;
 PObj SRegistry::find( string key ) const
 {
     auto pos = Map.find(key);
-    return pos == Map.end() ? NULL : pos->second;
+    return pos == Map.end() ? nullptr : pos->second;
 }
 
 PObj SRegistry::find_or_fail( string key ) const
diff --git a/pub/lib/toplevel.cpp b/pub/lib/toplevel.cpp
index 068ff8a3e1c7c81545e470ed45c60fb87b81e50f..9a9f34fa196d6250c168b83b0122143de5c6faaf 100644
--- a/pub/lib/toplevel.cpp
+++ b/pub/lib/toplevel.cpp
@@ -103,7 +103,7 @@ void CFrida::execute_cmd( const string cmdline )
             NOlm::IterateO fiter( NOlm::IterTOL );
             if( fiter.size()==0 || !T->k_dependent() ) {
                 RObj ret = T->tree_point_val();
-                if( !CNode::eval( "defined (\"$silent\")" )->to_b() && ret->to_s() != "__NULL__" )
+                if( !CNode::eval( "defined (\"$silent\")" )->to_b() && ret->to_s() != "__nullptr__" )
                     cout << ret->to_s() << "\n";
             } else {
                 while ( (fiter()) ) {