     int total_interpolation_points = n_total_boxes * n_interpolation_points * n_interpolation_points;
     int n_fft_coeffs_half = n_interpolation_points * n_boxes_per_dim;
     int n_fft_coeffs = 2 * n_interpolation_points * n_boxes_per_dim;
+    int n_interpolation_points_1d = n_interpolation_points * n_boxes_per_dim;
 
     // FIT-TSNE Device Vectors
     thrust::device_vector<int> point_box_idx_device(N);

            potentialsQij_device, num_points, num_nodes, n_terms);

	time(_time_norm)
        // Copy data back to the GPU

        delete[] box_lower_bounds;
        delete[] box_upper_bounds;
        delete[] y_tilde_spacings;
        delete[] y_tilde;
        delete[] x_tilde;
        delete[] fft_kernel_tilde;
        
        // Calculate Attractive Forces
        tsnecuda::bh::ComputeAttractiveForces(gpu_opt,

    thrust::device_vector<int> output_indices(
        n_terms * n_interpolation_points * n_interpolation_points * N);
    thrust::device_vector<float> chargesQij_device(N * n_terms);
    thrust::device_vector<float> box_lower_bounds_device(2 * n_total_boxes);
    thrust::device_vector<float> box_upper_bounds_device(2 * n_total_boxes);
    thrust::device_vector<thrust::complex<float>> fft_kernel_tilde_device(2 * n_interpolation_points_1d * 2 * n_interpolation_points_1d);
       
    // Easier to compute denominator on CPU, so we should just calculate y_tilde_spacing on CPU also
    float h = 1 / (float) n_interpolation_points;
    float y_tilde_spacings[n_interpolation_points];
    y_tilde_spacings[0] = h / 2;
    for (int i = 1; i < n_interpolation_points; i++) {
        y_tilde_spacings[i] = y_tilde_spacings[i - 1] + h;
    }
    float denominator[n_interpolation_points];
    for (int i = 0; i < n_interpolation_points; i++) {
        denominator[i] = 1;
        for (int j = 0; j < n_interpolation_points; j++) {
            if (i != j) {
                denominator[i] *= y_tilde_spacings[i] - y_tilde_spacings[j];
            }
        }
    }
    thrust::device_vector<float> y_tilde_spacings_device(y_tilde_spacings, y_tilde_spacings + n_interpolation_points);
    thrust::device_vector<float> denominator_device(denominator, denominator + n_interpolation_points);
    
    // Dump file
    float *host_ys = nullptr;
    std::ofstream dump_file;
