Optimal solution to finding all points within a certain distance

I have to find the nearest points within a particular distance for every point. For every 3d-point provided (labeled 1...n) I have to show all the nearest neighbors within a certain distance r.

I solved it, but my algorithm is O(n^2) and I'm not sure how to improve it.

I tried thinking about sorting all the points but I would still need to iterate over all x elements that are within distance r, making it a O(nx log n) algorithm, which isn't any better for large x.

#include <iostream>
#include <vector>
#include <cassert>
#include <cmath>
#include <fstream>
using namespace std;

struct Point {
    double x, y, z;
    double square(double v) {
        return v * v;
    }
    double operator-(Point const& rhs) const {
        return sqrt( square(rhs.x - x) + square(rhs.y - y) + square(rhs.z - z) );
    }
};

int main() {
    float r;
    int n;
    std::vector<Point> points;
    std::ifstream s("file2.txt");
    std::cin.rdbuf(s.rdbuf());
    assert(s.is_open());
    assert(std::cin >> r >> n);
    
    for (int i = 0; i < n; i++) {
        double x, y, z;
        assert(std::cin >> x >> y >> z);
        points.push_back({x, y, z});
    }

    std::vector<std::vector<int>> adj(n, vector<int>());
    
    for (int i = 0; i < n; i++) {
        std::cout << i << ": ";
        for (int j = i + 1; j < n; j++) {
            if (points[i] - points[j] <= r) {
                adj[i].push_back(j);
                adj[j].push_back(i);
            }
        }
        for (int j = 0; j < adj[i].size(); j++) {
            cout << adj[i][j] << " ";
        }
        cout << endl;
    }
}

This is the output, which is correct:

0: 1 2 3 8 
1: 0 4 6 8 
2: 0 4 5 8 
3: 0 5 6 8 
4: 1 2 7 8 
5: 2 3 7 8 
6: 1 3 7 8 
7: 4 5 6 8 
8: 0 1 2 3 4 5 6 7

I don't want the answer but a nudge in the right direction. Or is this the best I can do given I have to print all the points out?