In reply to Bruno Mario Cesana.

Thank you for sharing your opinion. I was very explicit about the details that I implemented. Users who do not want to exclude slopes that are exactly -1 can replace the statements

idx = loc(G[,3] ^= -1);
S = G[idx,3];

with

S = G[,3];

> which is the greatest number of X and Y values that your macro can manage?

There are no enforced bounds on the size, but since this implementation uses a brute-force O(n^2) computation for the pairwise differences, it begins to get slow as the sample size increases. Here are some timing on my PC:

N Time (s)
====== ========
1000 0.13
2000 0.48
4000 2.29
5000 3.40
8000 10.48
10000 17.20

By fitting a quadratic to those points, you can estimate times for other sizes.
You probably know that Raymaekers and Dufey (2022) published a faster method.