What is the cause of artifacts on the Z axis in my 3D perlin noise?

I am currently working on a 3D Perlin noise implementation in C++, and there are strange block-like artifacts when I visualize it as a video of 2D slices through the 3D volume as defined by this code with just one octave:

#include <png++/png.hpp>
#include <memory>
#include <string>

#include "Octave.h"

constexpr unsigned IMAGE_SIZE = 512;

std::string numberToString(unsigned n, unsigned digits);

int main() {
    std::mt19937_64 rnd(0);
    auto octave = std::make_unique<Octave>(&rnd, 32, 1); //make_unique because Octave objects are too big to fit on the stack

    for(unsigned z=0;z<625;++z){
        std::cout << z << "/625" << std::endl;
        png::image<png::rgb_pixel> image(IMAGE_SIZE, IMAGE_SIZE);
        for(unsigned x=0;x<IMAGE_SIZE;++x){
            for(unsigned y=0;y<IMAGE_SIZE;++y){
                unsigned brightness = (octave->noise(x, z*(64.0/300.0), y)*.5+.5)*255;
                image[y][x] = png::rgb_pixel(brightness, brightness, brightness);
            }
        }
        image.write("output/perlin-" + numberToString(z, 4) + ".png");
    }


    return 0;
}

std::string numberToString(unsigned n, unsigned digits){
    std::string string = std::to_string(n);
    while(string.length() < digits){
        string = "0" + string;
    }
    return string;
}

Each image from the above program is a single frame in a video. A video showing the output of this problem (converted to a video via ffmpeg) is available here: https://youtu.be/f7NxYo8U7TQ. Youtube added some compression artifacts in that video that are not present in the program's output.

As is clear in the above video, as the slice moves along the Z axis, square artifacts are present: . Depending on which plane slices are made in (eg. XZ or ZY), the artifacts slightly change in nature. It appears that they are worse just after a roll over in the frequency of the one rendered octave.

What is the cause of this?

Octave.h:

#pragma once

#include <array>
#include <random>

class Octave{
public:
    Octave(std::mt19937_64 *rnd, double frequency, double amplitude);
    double noise(double x, double y, double z);

private:
    static constexpr int PERMUTATION_TABLE_PART_SIZE = 1000000;
    static constexpr int PERMUTATION_TABLE_PART_COUNT = 3;
    static constexpr int PERMUTATION_TABLE_SIZE = PERMUTATION_TABLE_PART_SIZE*PERMUTATION_TABLE_PART_COUNT;
    std::array<int, PERMUTATION_TABLE_SIZE> m_permutationTable;
    double m_frequency, m_amplitude;

    double influence(int x, int y, int z, double distanceX, double distanceY, double distanceZ);

    inline static double square(double d) { return d*d; }
    inline static double vectorLength(double x, double y, double z);
    inline static double interpolate(double a, double b, double x);

};

Octave.cpp:

#include "Octave.h"

#include <utility>
#include <cmath>
#include <iostream>

Octave::Octave(std::mt19937_64 *rnd, double frequency, double amplitude) : m_frequency(frequency), m_amplitude(amplitude) {
    //fill in basic array
    for(int i=0;i<PERMUTATION_TABLE_PART_SIZE;++i){
        for(int j=0;j<PERMUTATION_TABLE_PART_COUNT;++j){
            m_permutationTable[i+PERMUTATION_TABLE_PART_SIZE*j] = i;
        }
    }

    //shuffle array
    for(int i=0;i<PERMUTATION_TABLE_SIZE;++i){
        int swapWith = ((*rnd)() % (PERMUTATION_TABLE_SIZE-i))+i;
        std::swap(m_permutationTable[i], m_permutationTable[swapWith]);
    }
}

double Octave::noise(double x, double y, double z) {
    x /= m_frequency;
    y /= m_frequency;
    z /= m_frequency;

    int intX = std::floor(x);
    int intY = std::floor(y);
    int intZ = std::floor(z);

    double floatX = x - intX;
    double floatY = y - intY;
    double floatZ = z - intZ;

    double influence1 = influence(intX, intY, intZ, floatX, floatY, floatZ);
    double influence2 = influence(intX+1, intY, intZ, floatX-1, floatY, floatZ);
    double influence3 = influence(intX+1, intY+1, intZ, floatX-1, floatY-1, floatZ);
    double influence4 = influence(intX, intY+1, intZ, floatX, floatY-1, floatZ);
    double influence5 = influence(intX, intY, intZ+1, floatX, floatY, floatZ-1);
    double influence6 = influence(intX+1, intY, intZ+1, floatX-1, floatY, floatZ);
    double influence7 = influence(intX+1, intY+1, intZ+1, floatX-1, floatY-1, floatZ-1);
    double influence8 = influence(intX, intY+1, intZ+1, floatX, floatY-1, floatZ-1);

    double frontUpperInterpolatedValue = interpolate(influence4, influence3, floatX);
    double backUpperInterpolatedValue = interpolate(influence8, influence7, floatX);
    double frontLowerInterpolatedValue = interpolate(influence1, influence2, floatX);
    double backLowerInterpolatedValue = interpolate(influence5, influence6, floatX);

    double upperInterpolatedValue = interpolate(frontUpperInterpolatedValue, backUpperInterpolatedValue, floatZ);
    double lowerInterpolatedValue = interpolate(frontLowerInterpolatedValue, backLowerInterpolatedValue, floatZ);

    return interpolate(lowerInterpolatedValue, upperInterpolatedValue, floatY)*m_amplitude;
}

double Octave::influence(int x, int y, int z, double distanceX, double distanceY, double distanceZ) {
    //create un-normalized gradient vector
    //the ordering of x, y, and z is arbitrary but different to produce different x y and z
    double gradientX = (m_permutationTable[m_permutationTable[m_permutationTable[x]+y]+z]/static_cast<double>(PERMUTATION_TABLE_PART_SIZE))*2-1;
    double gradientY = (m_permutationTable[m_permutationTable[m_permutationTable[y]+x]+z]/static_cast<double>(PERMUTATION_TABLE_PART_SIZE))*2-1;
    double gradientZ = (m_permutationTable[m_permutationTable[m_permutationTable[y]+x]+z]/static_cast<double>(PERMUTATION_TABLE_PART_SIZE))*2-1;

    //normalize gradient vector
    double gradientVectorInverseLength = 1/vectorLength(gradientX, gradientY, gradientZ);
    gradientX *= gradientVectorInverseLength;
    gradientY *= gradientVectorInverseLength;
    gradientZ *= gradientVectorInverseLength;

    //compute dot product
    double dot = gradientX*distanceX+gradientY*distanceY+gradientZ*distanceZ;

    return dot;
}

double Octave::vectorLength(double x, double y, double z) {
    return std::sqrt(square(x)+square(y)+square(z));
}

double Octave::interpolate(double a, double b, double x) {
    return (b-a)*(6*x*x*x*x*x-15*x*x*x*x+10*x*x*x)+a;
}

Solution

I determined the issue. I forgot to do floatZ-1 when calculating influence6.