@fff97 The paper you linked gives some information on how to implement this in purely hardware-float arithmetic, which is what Matlab uses. It'll work for reasonably sized d (the starting digit position). You'd need to implement the modular exponentiation for the expression

16 &^ (d-k) mod (8*k+j),

which'll only take a few lines of code. But this'll only work if d is small enough so that the computation will happen in exact-integer arithmetic. I think that that means that (8*d+6)^2 will need to fit into the 53-bit mantissa of a hardware float, which means (I think) that (8*d+6)^2 < 2^52.

@vv Thank you for the compliment.

Regarding 0 &^ 0: For purely integer arithmetic, I believe that 0^0 should be 1. What are some arguments to make it undefined? We have 0! = 1 and binomial(0,0) = 1. And for the program above, making 0 &^ 0 equal 1 produces the correct result, which allows me to extend the BBP algorithm from the paper (where it's stated that d > 0) to the case d=0.

For real-number and complex-number arithmetic, I have no problem with keeping 0^0 undefined. Unfortunately, Maple's automatic simplification makes it 1.

@Deank1905 Yes, my technique can be directly applied to your situation:

LCM:= 1;
for i from 1 to 10 do
    # Your code that computes an individual r goes here:
    # ...
   LCM:= ilcm(LCM, r)
end do; 
LCM;

This technique doesn't "get all of the values of r out of the loop" because that's not necessary to compute their lowest common multiple. But if you still do want to get the values of r out of the loop for some other reason, then I can show you how to do that.

@Bosco Emmanuel The weird thing about what you say is that the error message that you received is precisely what it would've been had you not included the ~s in your input. You could do this instead:

M2:= map(rhs, Matrix(map(convert, [solution], list)));

The ~ is (essentially) just an abbreviation for map. (There are a few technical differences; it's not a purely syntactic abbreviation.) What Maple version are you using? This use of ~ was introduced in Maple 13 I believe.

@Bosco Emmanuel You didn't transcribe Tom's solution correctly. You wrote

M2 := rhs(Matrix(convert([solution], list)));

whereas he wrote

M2 := rhs~(Matrix(convert~([solution], list)));

@Bosco Emmanuel The number of rows is nops([solution]).

@Muhammad Usman You never say what the contents (or entries) of the matrices named A are supposed to be. I don't see how your specific question

A[n]:= ?

can be answered without that information.

Nowhere do you say what goes into the matrices A. You only say what the As go into.

@mmcdara Yeah, I think that PartialSums is misplaced in ListTools. It's more of a mathematical command than a list-manipulation command.

@ravenHound You are confusing the extremal values of N with the values of t at which those extrema are attained. If you plug those t back into N, you'll get the other two values.

To get both the values of N and t for the extrema directly from the plot, do

P:= plot(N, 0..20):
A:= op([1,1], P): #The plot's point matrix
A[[(min[index],max[index])(A[.., 2])], ..];

Of course, these values that are taken directly from the plot are just approximations because a plot is just a finite set of points (213 points in this case). However, the plot command does usually increase the point density a little bit near the extrema.

@ravenHound Your analyses are totally correct, and I'm convinced that you understand the underlying math.

Regarding the Maple: Yes, the graph will theoretically be a straight line with slope 2/3. But in reality, using actual data and doing the log-log transformation, you'll get a not-quite-straight (at best) line with slope approximately 2/3. To find the line that best fits that data, you use Maple command Statistics:-PowerFit.

The professor's instruction to "experiment with the graphing calculator to see how the value of a affects the graph" is silly to my mind: It's probably obvious to you how a multiplied constant affects a graph.

@mmcdara I think that your interpretation is the correct one, but I don't understand why/how dsolve solved it. There's more than one independent variable with respect to which derivatives are taken, so I'd expect only pdsolve to work.

@LustForLife You probably realize that your procedure redoes the sums and products of the initial segments many times. This can be easily avoided by using more elementwise operations:

BigProc:= (H::list, T::list)-> 
   (P-> (1 +~ ListTools:-PartialSums(P)) /~ (1 +~ ListTools:-PartialSums(1/~P)))
      (H*~T)
:

@Carl Love Here's an update that detects limit cycles of any length. Unlike Tom's, this one doesn't require a linear search of all of the previous values to detect a repeat. However, it does require that each value be stored twice.

KaprekarSeq:= proc(x::posint, b::posint:= 10)
local j, r:= Vector(1, [x]), kx:= Kaprekar(x), J:= table([x=1]);
   for j from 2 while not assigned(J[kx]) do
      J[kx]:= j;
      r(j):= kx; 
      kx:= Kaprekar(kx, b)
   end do;
   convert(r, list), convert(r[J[kx]..], list)
end proc:

Kaprekar:= proc(x::posint, b::posint:= 10)
local X:= convert(x, base, b), N:= b^nops(X), S:= sort(X);
   convert(S, base, b, N)[] - convert(ListTools:-Reverse(S), base, b, N)[]
end proc:

Example usage:

KaprekarSeq(12345);

      [12345, 41976, 82962, 75933, 63954, 61974], [82962, 75933, 63954, 61974]

The first returned list is the entire sequence; the second is the limit cycle.

@tomleslie Yes, that was intentional; I mentioned it in my introduction. However, inspired by you, I've improved my code to detect any cycle, and I'll post that update in a few minutes.