@David Sycamore Just replace a(n) with a(n)/n in the previous plotting commands:

m:= 3000:
L:= [seq](a(n)/n, n= 1..m):
plots:-listplot(       
    L, style= point,  
    symbolsize= max(1, round(16/(1+ilog2(m))))
);

@Arastas I think that this simplified example may show you the crucial difference that makes it work. Compare the results of these two sum commands:

sum(t^2, t= k..k+13);
sum(t^2, t= k..k+n);

The 2nd uses symbolic summation; the 1st just crudely adds 14 terms together. The 1st could've used symbolic summation, but it chose not to. When the number of terms is specifically known to sum (as in the 1st case), it uses a heuristic to decide between symbolic summation and simple addition. For the case at hand, that heuristic turns out to not make the best choice.

The procedure pellsolve is missing from your posted code. Please provide it.

@stud_8013 Due to some bug in MaplePrimes, the number of Replies is shown as 0 for this Question (and it'll likely be shown as 1 after I post this), so I don't know if you saw that your Question has been Answered. And if you did see the Answer, I expected some sort of Reply (like whether it works for you).

@David Sycamore The sequence a(n)/n appears to have three distinct convergent subsequences, two converging to 3/2 (guess), one from above and one from below, and one converging to 3/4 (guess) from below.

@Christopher2222 You asked:

• But why did the author choose to show it as tan^(-1)?  Probably his convention of choice.

It is an unfortunate and ambiguous textbook convention to use -1 (with no parentheses) as a superscript on function names to indicate the functional inverse and to use positive integer superscripts with exactly the same typography to indicate algebraic exponentiation (rather than functional iteration). In Maple's 2D output, a (-1) superscript---in parentheses---on a function name indicates the functional inverse, but in cases such as tan where that inverse already has its own name (arctan), the tan@@(-1) gets replaced by arctan, unless you use unevaluation quotes or some form of inertness.

@David Sycamore I have no idea why those plot commands are not working for you.

Starting at about 500 points, one can see a remarkably regular pattern in the plot. Here's the plot for 1000 points:

At 10,000 points, there's no additional pattern apparent, just continuation of the 3 rays.

I have two quibbles with the wording (only) of your sequence's defining phrase, to wit:

• Lexicographicaly least sequence of nonnegative integers commencing 1,3,5,7 such that any four consecutive terms are mutually coprime.

I'd change it to

• Lexicographically least sequence of distinct nonnegative integers commencing 1, 3, 5, 7 such that any four consecutive terms are mutually pairwise coprime

My quibble is with the wording only, not with the underlying concepts.

@Thomas Dean 

Both my solution using ArrayTools:-AddAlongDimension and the one using ListTools:-LengthSplit implicitly rely on the fact that the pattern of indices in sumidx is extremely simple and thus sumidx isn't needed at all: We just need to add the 1st 10 elements of V, the 2nd 10 elements, etc. That's why they are so fast. I assume that in your actual application, there wouldn't be such a simple pattern, and you'd actually need the list sumidx. So, I don't think that it's valid to compare the method that uses LengthSplit with those using sumidx and for or seq for the summations.

Note that I made corrections to the Answer above; so, if you've tried that code already, please try it again.

Here's another method that I like much better than that. It also uses no indices; additionally, it uses no matrices or any other containers other than lists. The only difference is in the last 2 lines, which construct V and S. 

restart;
roll:= (n::posint)-> RandomTools:-Generate(list(posint(range= 9), n)):
nN:= (Jx:= 5)*(Ix:= 10): 
N:= sort(10*['seq'(3..2+Jx)$Ix] + roll(nN));
V:= 10*(roll(nN) + 10*iquo~(N, 10));
S:= add~([ListTools:-LengthSplit](V, Ix));

The output for S is the same as before.

@Carl Love The exact cause of your error is that you executed the code once, and then, without doing the restart, you executed it again. When a for loop terminates because its upper bound has been reached, its index variable is incremented past the upper bound. Thus, n=3 after your for loop. This n is then re-used as the index to sum. The index variable of sum can't have a numeric value.

If you execute your posted code straight through (including the restart), then you'll get a different error. This error is caused by your use of single quotes in the expression 'int'('R(x)' * ..., ...). No quotes are needed. Just use int(R(x) * ..., ...).

@rlopez Thanks for clarifying that. I was well aware of all the mathematical issues involved. I asked because I didn't know whether the OP's Question was simply ill-posed or whether "level curve f(x,y) = c" was some sort of "acceptable" or "known" textbook-problem abbreviation for "level curve f(x, y, 0) = c" that I simply had never seen before.

@Thomas Dean For each value of N, take the average of the 10 CPU times and 10 elapsed times. Set Tcpu and Telapsed to those averages. The raw data is not needed after that. For the larger values of N, just use whatever measurements you have; it needn't be 10.

@Thomas Dean I figured that you were doing integer factoring. Several points:

[Edit: In section 1 below, I originally incorrectly said that the n in the O formula represented the number of bits of the input number. Actually it's the input number itself. This has now been corrected.] 

1. There's no known polynomial-time algorithm for it (unless you include quantum computation on a scale that's not known to have been achieved yet). NFS is known to be "sub-exponential, super-polynomial". The best implementations are known to be specifically 

O(e^(4*(ln(n)*ln(ln(n))^2/9)^(1/3)))

where n is the number being factored (rather than the number of its bits or digits). To put this in terms of your computation, we let d = ceil(log[10](n)), so d is the number of base-10 digits. Using this in the formula above, simplifying, and discarding lower-order terms in the exponent, we get

O(A^((d*ln(d)^2)^(1/3)))

where A = exp(4*(ln(10)/9)^(1/3)), which is approximately 12.67.

2. For any algorithm, using parallelism (multi-processing, hyper-threading, distributed processing, etc.) will neither change the "big-O" nor help you to determine it. At best, parallelism gives you a constant factor A (0 < A < 1) reduction in elapsed time, and you "pay" for it with a constant factor B (B > 1) increase in total CPU time (and hence electricity usage), where necessarily A*B > 1 (often much greater and increasing with the number of processors used). Since "big-O" (formally known as part of Landau order notation) ignores constant factors, parallelism won't lower the order. (Although it's possible that poorly implemented parallelism could increase the order.) So, for this purpose, you should only use sequential (i.e., single-processor) computation.

3. The "closeness" (as measured by the sum of the squared residuals or by any other metric) of a "fit" will always get better as the number of parameters being fitted approaches the number of data points. For example, a degree-d general polynomial (so, d+1 coefficients are being fitted) can always be found that's an exact fit (meaning the residuals are 0) for any set of d+1 data points with distinct independent-variable values. The statistical predictive value of such a fit is worthless. So, criteria in addition to small residuals are needed to assess the quality of a fit. When doing a big-O computation with a single independent variable (the n or N in this case), I don't think that you should ever use a model function with more than one term.

4. So, is there a method better than trial and error? Yes, start with the theoretical big-O known for the algorithm, which I gave above.

This Reply is "tangential"^[*1] to your main Question.

Your Question shows the plaintext input

series(tan^(-1)(x), x, 9)

and the output is the series for arctan(x), and I suspect that that was your intended function. In Maple's 1D input, tan^(-1)(x) is not arctan(x). (I'm not criticizing your post; I'm just curious how this expression came to be in this form and whether Maple, MaplePrimes, or your OS or web browser had anything to do with that.) So, I guess that you entered the command using 2D Input. My first question is How exactly did you create that 2D Input (my curiosity is only about the first argument)? My second question is How did you transcribe that to MaplePrimes (e.g., simple typing, copy/paste, something else)?

In 1D input

(tan@@(-1))(x) = arctan(x)  # i.e., the functional inverse of tangent
(tan^(-1))(x) = cot(x)  # i.e., cotangent, the multiplicative inverse of tangent
tan^(-1)(x) = 1/tan

In the last of those, tan is being treated as just a name, and (-1) is seen as a constant function whose argument is x. Since function application has higher operator precedence than ^, the expression is the same as tan^((-1)(x)) = tan^(-1) = 1/tan.

[*1] I realize of course that the actual function used is irrelevant to your Question, and all that matters is that the series have an O(...) term.