porting code from https://github.com/Craig-Macomber/N-dimensional-simplex-noise

simplex.cpp

#include <vector>
#include <cmath>
#include <iostream>
#include <numeric>
 
template<typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& obj) {
  os << "[";
  for(uint i = 0; i < obj.size(); ++i){
      os << obj.at(i);
      if(i+1 < obj.size()){
          os << ", ";
      }
  }
  os << "]";
  return os;
}

template<typename Ta, typename Tb>
std::vector<Ta> sort_with(const std::vector<Ta>& first, const std::vector<Tb>& second){
    std::vector<std::pair<Ta, Tb>> zip;
    for(uint i = 0; i < first.size() && i < second.size(); ++i){
        zip.emplace_back(first[i], second[i]);
    }
    
    struct {
        bool operator()(std::pair<Ta, Tb> a, std::pair<Ta, Tb> b)
        {   
            return a.second < b.second;
        }   
    } compare_second;
    
    std::sort(zip.begin(), zip.end(), compare_second);
    
    std::vector<Ta> r;
    for(auto pair : zip){
        r.push_back(pair.first);
    }
    return r;
} 
 
template<int dimensions>
struct simplex_noise{
    simplex_noise(int init_seed): seed(hash(init_seed)), d_range(dimensions){
        std::iota(d_range.begin(), d_range.end(), 0);
        
        // Compute skew constant. This is questionably correct
        skew_constant = (std::sqrt(dimensions + 1) - 1) / dimensions;
        
        // This is the proper relation between f and g in terms of d
        // that makes sure skewing and unskewing reverse eachother
        unskew_constant = skew_constant / (1 + (dimensions * skew_constant));
 
        // distace from corner to center of most distant face
        // this is the max influance distance for a vertex
        if (dimensions == 1) {
            corner_to_face_squared = 1.0;
//            corner_to_face_squared = std::pow(side_length, 2);
        } else {
            // simplex edge length
            float side_length = std::sqrt(dimensions) / ((dimensions * skew_constant) + 1);
            float a = std::pow(std::pow(side_length, 2) - std::pow(side_length / 2, 2), 0.5);
            
            if (dimensions == 2){
                corner_to_face_squared = std::pow(a, 2);
            } else{
                corner_to_face_squared = std::pow(a, 2) + std::pow((a/2), 2);
            }
        }
                
        // Precompute gradient vectors.
        // Make all possible vectors consisting of:
        // +1 and -1 components with a single 0 component:
        // Vecs from center to midddle of edges of hyper-cube
 
        if (dimensions > 1) {
            // inject 0s into vecs from vf to make needed vecs
            for(int z : d_range) {
                auto vecs_for_d = gen_gradients();
                for(auto gradient : vecs_for_d){
                    gradient.insert(gradient.begin() += z, 0);
                    gradients.push_back(gradient);
                }
            }
                        
            // All 1 or -1 version (hypercube corners)
            // Makes number of vecs higher, and a power of 2.
            //self.vecs=[v for z in xrange(d) for v in vf(d)]
            
            // Compensate for gradient vector lengths
            value_scaler = std::pow(dimensions - 1, -0.5);
            // use d instead of d-1 if using corners instead of edges
            
            // Rough estimated/expirmentally determined function
            // for scaling output to be -1 to 1
            value_scaler *= std::pow(dimensions - 1,-3.5) *100 + 13;
        } else {
//            skew_constant = 0; 
            gradients = {{1.0f}, {-1.0f}};
            value_scaler = 1.0;
//            corner_to_face_squared = 1.0;
//            unskew_constant = 1.0;
        }
        
        // shuffle the vectors using self.seed
       std::mt19937 rand(seed);
       std::shuffle(gradients.begin(), gradients.end(), rand);
    }
 
    float noise(std::array<float, dimensions> coord, bool getDerivative){
        //""" loc is a list of coordinates for the sample position """
        // Perform the skew operation on input space will convert the

        // regular simplexs to right simplexes making up huper-cubes
        auto s = std::accumulate(coord.begin(), coord.end(), 0) * skew_constant;
        
        // Skew and round loc to get origin of containing hypercube in skewed space
        std::vector<int> skewed_coord; // intSkewLoc
        for(auto v : coord){
            skewed_coord.push_back((int) std::floor(v + s));
        }

        // Unskewing factor for intSkewLoc to get to input space
        auto t = std::accumulate(skewed_coord.begin(), skewed_coord.end(), 0) * unskew_constant;
        // skewed simplex origin unskewed to input space would be:
        // cellOrigin=[v-t for v in intSkewLoc]
        
        // Distance from unskewed simplex origin (intSkewLoc[i]-t) to loc,
        // all in input space
        std::vector<float> cell_dist; // cellDist
        for(int i : d_range){
            cell_dist.push_back(coord[i] - skewed_coord[i] + t);
        }

        // Indexs of items in cellDist, largest to smallest
        // To find correct vertexes of containing simplex, the containing hypercube
        // is traversed one step of +1 on each axis, in the order given
        // by greatest displacement from origin of hyper cube first.
        // This order is stored in distOrder: The order to traverse the axies
        auto dist_order = sort_with(d_range, cell_dist);
        dist_order.push_back(-1);
        std::reverse(dist_order.begin(), dist_order.end());
                
        // Copy intSkewLoc to work through verts of simplex
        // intLoc will hold the current vertex,
        // and intSkewLoc will stay the containing origin's vertex
        // these are still the skewed right simplexes/ hypercube space vertex indexes
        auto skewed_coord_copy(skewed_coord); // intLoc

        // our accumulator of noise
        float n(0.0);
        // skewOffset holds addational skew that needs to get added.
        // It will be self.g * the how many +1s have been added to all the axies total
        // Thus, adding it to the current vertex skewes it to be in input space
        // relative to the simplex's origin vertex:
        // intSkewLoc in skewed space.
        float skew_offset(0.0); // skewOffset

        for(int v : dist_order){
            // Move to the next corner of simplex of not on first corner
            if (v != -1){
                skewed_coord_copy[v] += 1;
            }

            // get u: loc's position relative to the current vertex, in input space
            std::vector<float> u;
            for(int i : d_range){
                u.push_back(cell_dist[i] - (skewed_coord_copy[i] - skewed_coord[i]) + skew_offset);
            }
            // t is the factor for attenuating the effect from the current vertexes gradient
            // its based on distance squared
            // if distance squared exceeds self.cornerToFaceSquared, there is no contribution
            auto t = corner_to_face_squared;
            for(float a : u){
                // Accumulate negative distance squared into t
                t -= (a * a);
            }
            
            if(t > 0){
                // fech a pseudorandom vec from self.vecs using intLoc
                auto index = seed;
                for (int i : d_range){
                    index += skewed_coord_copy[i] << ((5*i)%16);
                }
                auto gradient = gradients[(hash(index)) % gradients.size()];
                                
                // dot product of vertex to loc vector and vector for gradient
                float gr(0);
                for (int i : d_range){
                    gr += gradient[i]*u[i];
                }
                
                // add current vertexes contribution
                float t4= std::pow(t,4);
                n += gr*t4;
            }
            skew_offset += unskew_constant;
        }
        return n * value_scaler;
    }
 
    private:
    long hash(uint value){
        long hash = 255 & value;
        hash += (hash << 10);
        hash ^= (hash >> 6);

        hash += 255 & (value>>(8));
        hash += (hash << 10);
        hash ^= (hash >> 6);
        
        hash += 255 & (value>>(16));
        hash += (hash << 10);
        hash ^= (hash >> 6);
        
        hash += (hash << 3);
        hash ^= (hash >> 11);
        return hash+(hash << 15);
    }
 
    std::vector<std::vector<float>> gen_gradients (){
        // Little helper generator function for making gradient vecs
        // Makes vecs of all -1 or 1, of d dimensions 
        std::vector<std::vector<float>> vectors;
        for(int i = 0; i < std::pow(2, dimensions - 1); ++i){
            std::vector<float> vector;
            for(int d = 0; d < dimensions - 1; ++d){
                vector.push_back((i >> d) % 2 * 2 -1);
            }
            vectors.push_back(vector);
        }
        return vectors;
    };
    int seed;
    float skew_constant, unskew_constant, corner_to_face_squared, value_scaler;
    std::vector<std::vector<float>> gradients;
    std::vector<int> d_range;
};