How do you account for the curvature at either cone's apex?

@AOdrzywolek @vv Both of you (VV and the OP) seem to be unaware of this paragraph from help page ?evalhf (where I've added the emphasis):

• The evalhf function converts all its arguments to hardware floats, computes the answer and converts the answer to a Maple float result. Hence the user never deals with hardware floating-point numbers....

In short, evalhf does not return hfloats, which is exactly what you seem to be trying to get it to do; rather, it returns sfloats.

If you want hfloats to write to your file, then put the code that creates and writes those numbers in a procedure with option hfloat. Do not use evalhf within this procedure, because even within such a procedure, evalhf still returns sfloat​​​​​​​s.

@Kitonum Indeed, one could very easily extract that numeric (x,y,z) double-loop data directly from the plot.

@lcz Like I said, there will be  no memory savings from using iterator if you save the graphs. Nonetheless, here are two ways to do it:

restart:
#
#Use iterator and save graphs by diameter--Array method
#
CountDiamterClasses4:= proc(order::posint)
uses GT= GraphTheory;
local  
    prc:= GT:-NonIsomorphicGraphs(
        order, output= iterator, outputform= graph, restrictto= connected
    ),
    d, R:= table([seq(d= Array(1..0), d= 1..order-1)]), graphmark
;
    while (graphmark:= prc()) <> FAIL do R[GT:-Diameter(graphmark)],= graphmark od;
    {seq}~(R)
end proc
:
CodeTools:-Usage(CountDiamterClasses4(9)):
memory used=11.99GiB, alloc change=0.91GiB, cpu time=2.90m, 
real time=2.57m, gc time=26.03s

nops~(%);
   TABLE([1 = 1, 2 = 91518, 3 = 148229, 4 = 19320, 5 = 1818, 
     6 = 180, 7 = 13, 8 = 1])

restart:
#
#Use iterator, and save graphs---table method.
#
CountDiameterClasses5:= proc(order::posint)
uses GT= GraphTheory;
local  
    prc:= GT:-NonIsomorphicGraphs(
        order, output= iterator, outputform= graph, restrictto= connected
    ),
    R, graphmark
;
    while (graphmark:= prc()) <> FAIL do R[GT:-Diameter(graphmark)][graphmark]:= () od;
    {indices}~(R, 'nolist')
end proc
:
CodeTools:-Usage(CountDiameterClasses5(9)):
memory used=12.00GiB, alloc change=0.91GiB, cpu time=2.90m, 
real time=2.56m, gc time=26.69s

nops~(%);
   TABLE([1 = 1, 2 = 91518, 3 = 148229, 4 = 19320, 5 = 1818, 
     6 = 180, 7 = 13, 8 = 1])

@lcz Here's a comparison of three methods for counting the diameter classes of graphs of order 9. The first uses an iterator, the second is what I used in my Answer, and the third uses an iterator and the naive list-building method, without saving the graphs. You'll see, as expected, the major benefit of iterator is a reduction in memory allocation (reported as "alloc change" by CodeTools:-Usage).

restart:
#
#Use iterator and count accumulators:
#
CountDiameterClasses1:= proc(order::posint)
uses GT= GraphTheory;
local  
    prc:= GT:-NonIsomorphicGraphs(
        order, output= iterator, outputform= graph, restrictto= connected
    ),
    R:= table(sparse), graphmark
;
    while (graphmark:= prc()) <> FAIL do ++R[GT:-Diameter(graphmark)] od;
    eval(R)
end proc
:
CodeTools:-Usage(CountDiameterClasses1(9));
memory used=11.85GiB, alloc change=36.00MiB, cpu time=2.41m, 
real time=2.41m, gc time=4.41s

TABLE(sparse, [1 = 1, 2 = 91518, 3 = 148229, 4 = 19320, 5 = 1818, 
  6 = 180, 7 = 13, 8 = 1])

restart:
#
#Use direct generation of graphs:
#
CountDiamterClasses2:= proc(order::posint)
uses GT= GraphTheory, LT= ListTools;
    nops~(
        LT:-Classify(
            GT:-Diameter, 
            [GT:-NonIsomorphicGraphs](
                order, output= graphs, outputform= graph, restrictto= connected
            )
        )
    )
end proc
:  
CodeTools:-Usage(CountDiamterClasses2(9));
memory used=11.95GiB, alloc change=1.00GiB, cpu time=2.90m, 
real time=2.47m, gc time=32.62s

   TABLE([1 = 1, 2 = 91518, 3 = 148229, 4 = 19320, 5 = 1818, 
     6 = 180, 7 = 13, 8 = 1])

restart:
#
#Use iterator, but also use naive method of building lists inside loop.
#
CountDiameterClasses3:= proc(order::posint)
uses GT= GraphTheory;
local  
    prc:= GT:-NonIsomorphicGraphs(
        order, output= iterator, outputform= graph, restrictto= connected
    ),
    d, R:= table([seq(k=(), k= 1..order-1)]), graphmark
;
    while (graphmark:= prc()) <> FAIL do R[GT:-Diameter(graphmark)],= 0 od;
    (nops@`[]`)~(R)
end proc
:
CodeTools:-Usage(CountDiameterClasses3(9));
memory used=126.32GiB, alloc change=66.64MiB, cpu time=5.79m, 
real time=5.92m, gc time=30.47s

   TABLE([1 = 1, 2 = 91518, 3 = 148229, 4 = 19320, 5 = 1818, 
     6 = 180, 7 = 13, 8 = 1])

@AOdrzywolek Why do you want to "export" results in "binary"? "Export" implies that the file is meant to be read by some other program. What is this primitive program that can't read normal ASCII data?

@lcz The code that you posted in your Question was saving every graph generated into one of the lists. Since all the graphs are saved, there is no memory benefit from using the iterator. If instead you simply wanted to count the graphs of each diameter, or do anything else other than save them, then you can benefit from the iterator.

The reason that your code took more time has nothing to do with the GraphTheory package. It took more time because you were building lists inside a loop. In Maple, that's extremely inefficient, and you should never do it. Instead, you should build an Array or table inside the loop, and if needed, convert it to a list outside the loop.

@emendes Putting [] after a list or set (only those!) returns the underlying sequence, equivalent to op(...). The syntax allows [] to be put after almost anything, but unless it's a list or a set, it mostly does nothing other than create a mess. 

Well, clearly this one wasn't removed, so it's not automatic.

The most common reason for removal is duplication, i.e., the same question or a too-similar question asked by the same user.

It is often I who removes them, but I don't recall removing any of yours. Anyone with a reputation over 500 can moderate this forum.

@vv What you said about the 18 digits is true; however, not every possible sequence of 18 base-10 digits has a distinct representation as an hfloat, as your example posted immediately above shows. (In other words, your input and output both map to the same hfloat.) But every 15-digit number does map to a distinct hfloat, and 15 is the greatest integer for which this is true, which is why 15 is the standard returned by trunc(evalhf(Digits)).

You wrote:

• I don't see any contradiction using evalhf when Digits is 30. It could be e.g. a partial result, or a first approximation.

The OP's very first example clearly shows the nonsense that can result. Note that that first result is returned as an sfloat at 20 digits, not 18! Yes, it may be appropriate to use evalhf to get a first approximation, but (if using the results at the top level) you should change the value of Digits as you do it.

Don't repeat Questions as new Questions. They'll be deleted. If you want to raise your Question up in the "recents" stack, just post a followup Reply to it.

@emendes Your code is correct.

Your B could be made more efficient and more robust (easily extended to other situations) by 

op1:= {op~(1, parms)[]}:
B:= s-> {op~(1,s)[]} = op1:

There's no need to use S; you could use parms throughout.

The ,= is the new appending assignment operator. It's documented on help page ?assignment (or you can do simply ?,=).

@fatemeh1090 To try solve(..., parametric, real) with the new inequalities, they must be in polynomial form. The only one which isn't is u > 0, because of the sqrt in your specification of u. It can easily (by hand, for example) be shown to be equivalent to c*q+m > 0, which is polynomial​​​​​. If I put all those inequalties into the same solve that I had before, the command is

pwsol:= solve(
    {
       pp, 
       c*q+m > 0, 1+m-q > 0, q > 0, m > -1, m < 1, s > 0, c > 0,
       (1+m-q)^2 - 4*(m+c*q) > 0, q^2 + (-4*c-2*m-2)*q + (m-1)^2 >= 0
    }, 
    q, parametric, real
);

I let this command run (in Maple 2020) for over 3 hours of cputime and use 5 Gig of memory. Eventually it stopped due to some small bug. So, this is too complicated for the current capabilities of the command.

However, as I said, there are essentially only two possible solutions. The fact that you don't get 0 when you back-substitute them and simplify is largely irrelevant.

In conclusion, the title of this Question is incorrect in that there is no error at all in your results. Rather, it's the normal result of an equation with both a radical and parameters.

Other than the initial values, all the terms are of the form 3*2^m-1, not just the prime terms (modulo a small calculation error that you made).

@Jean-Claude Arbaut I was just about to post an Answer, but I see that you've already figured out what I was going to say.

What you called "type promotion" is often called floating-point contagion in Maple circles. In other words, once a symbolic expression has any float in it, symbolic expressions derived from it tend to become "infected" with even more floats. So, generally it's a bad idea to use evalf on expressions with free variables. But when you use seq, the numeric values of n are passed to u due to seq's special evaluation rule.