@Teep You should add this procedure to my code, right after AssignToClusters:

#Returns the partition of the point indices.
PartitionIntoClusters:= proc(U::Matrix)
local CL:= AssignToClusters(U), k;
    [seq([ListTools:-SearchAll(k,CL)], k= 1..max(CL))]
end proc
:

Running this on the same example as before, I get

CP:= PartitionIntoClusters(U);
[[4, 5, 6, 10, 14, 15, 19, 20], [1, 2, 3, 7, 8, 9, 11, 12, 13, 16, 17, 18],
[21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40]]

You may get a different partition due to randomization of the initial configuration done in my algorithm.

The partitioning of the original matrix of points can then be done by

map2(index, XM, CP);

I should've mentioned this at the time: The continuous red curve in my plot (along the top of the histogram) is not meant to be a proper probability function. It's just a simple connecting of the dots of the discrete function. Sorry if this caused any confusion. 

My guess is that there was a temporary slowdown in the Maplesoft license server. I recall a few times in the past week where startup seemed slow, but the problem seems to have passed. Does the problem persist for you? Here's an experiment: Compare startup times with and without an Internet connection. 

That's a rather banal Post. Were you just testing the MaplePrimes posting system?

@JAMET A related curiousity is that it was recently (2015) proved (by Samir Siksek) that "Every integer greater 454 is the sum of at most seven positive cubes." But, they're not necessarily distinct.

@Teep Let's consider the example of dividing the 40 points in matrix XM into 3 clusters. Then the following commands will partition XM into 3 matrices, one for each cluster:

(C,U):= FuzzyCMeans(XM, clusters= 3):
CL:= AssignToClusters(U):
(C1,C2,C3):= seq(XM[[ListTools:-SearchAll(k,CL)]], k= 1..3); 

In other words, if cluster label k appears in the i^th position in list CL, then the i^th row of XM is a point in cluster number k.

Hmm, I don't like this behavior. It shows simplify and expand ignoring the inertness! The only command that should be allowed to ignore inertness is value.

Also try `%%sum` to see another curiosity.

@JAMET Change root~(3,c) to iroot~(c,3). But note that most cases have no solution. If you include negative cubes, most will have a solution, but they'll take longer to find.

The commands option package and with(TestLibrary) are superfluous to your example. They don't cause any trouble here, but they are heavy-handed commands that (IMO) should be avoided when they're not needed.

All of the rounding modes mentioned by Tom, and other modes, are available in Maple through the environment variable Rounding. See help page ?Rounding.

So, with this, I believe that Maple can easily emulate the rounding behavior of any other software.

But the whole mathematical theory of rounding errors is built on the concept of "significant digits" which is akin to "relative error" but is less well related to "digits after the decimal place" (which is akin to "absolute error".)

I am eager to see the results of applying this idea to your graph, which has far more degree-2 nodes than any randomly generated graph of the same order and size.

@JAMET In order for me to diagnose the problem, you'll need to upload as an attached file an executed worksheet showing the error message. Use the green uparrow on the toolbar.

The most likely culprit is that you haven't correctly transcribed the line found:= false;.

@JAMET If you're truly using Maple 2018 like your header says, this shouldn't make a difference, but try changing add(c) to add(x, x= c).

@JAMET Put your Replies after the Answer that they pertain to.

@JAMET The code above will only run in 1D input.