@2cUniverse You can never make a valid comparison between two measurements both of which have been rounded down to 0. You must either use a finer measuring instrument, or increase the scale of the problem until the instrument returns positive values. The latter is usually easier in Maple.

iremFrac2a:=proc(n,d, b)
   local r1,r,i,li:
   li:=NULL:
   r:= irem(n*b,d):
   if NumberTheory:-AreCoprime(d,b) then
     while true do
        r1:= irem(b*r,d):
        li:=li,[[r,r],[r,r1]],[[r,r1],[r1,r1]]:
        if r1=n then break fi:
        r:=r1:
     end do:
   else print("Denom-Base not coprime "): return [] fi:
   [li]
end proc
:   
iremFrac2b:= proc(n, d, b) local N:= irem(n,d), r:= irem(N*b, d), s; 
    [if igcd(b,d) = 1 then 
        (do [[r,r], [r, (s:= irem(b*r, d))]], [[r,s], [s,s]] until (r:= s)=N)
    else
        print("Denom-Base not coprime")
    fi]
end proc
:
d:= nextprime(2^17):  b:= irem(rand(), d);
                          b := 116956

(L1, CT1, RT1):= CodeTools:-Usage(
    iremFrac2a(1, d, b), output= [output, cputime, realtime]
)[]:
memory used=14.24GiB, alloc change=-32.00MiB, 
cpu time=5.09m, real time=79.49s, gc time=4.93m

(L2, CT2, RT2):= CodeTools:-Usage(
    iremFrac2b(1, d, b), output= [output, cputime, realtime]
)[]:
memory used=14.89MiB, alloc change=0 bytes, 
cpu time=172.00ms, real time=165.00ms, gc time=0ns

#Verify that computed results are equal:
evalb(L1 = L2);
                              true

#time ratios (real and CPU):
RT1/RT2, CT1/CT2;
                    481.7333333, 1776.529070

To be fair, those ratios are not as extreme if I increase the scale of the problem "horizontally" (i.e., smaller test cases but more of them) rather than "vertically" (i.e., a single large test case), as above. I'll explain why in a moment.

d:= nextprime(2^10):  b:= irem~(['rand()'$64], d);
b := [173, 917, 905, 706, 476, 906, 759, 1001, 827, 720, 297, 35, 
  933, 582, 61, 289, 139, 245, 817, 496, 11, 477, 1001, 541, 824, 
  911, 787, 901, 990, 853, 562, 160, 248, 313, 631, 873, 630, 
  281, 565, 282, 85, 761, 311, 1021, 475, 123, 593, 632, 209, 
  683, 931, 296, 348, 381, 890, 772, 720, 565, 529, 277, 861, 
  286, 147, 912]

(L1, CT1, RT1):= CodeTools:-Usage(
    iremFrac2a~(1, d, b), output= [output, cputime, realtime]
)[]:
memory used=283.35MiB, alloc change=8.00KiB,
cpu time=6.34s, real time=1.70s, gc time=6.02s

(L2, CT2, RT2):= CodeTools:-Usage(
    iremFrac2b~(1, d, b), output= [output, cputime, realtime]
)[]:
memory used=13.95MiB, alloc change=0 bytes,
cpu time=172.00ms, real time=172.00ms, gc time=0ns

#Verify that computed results are equal:
evalb(L1 = L2);
                              true

#time ratios (real and CPU):
RT1/RT2, CT1/CT2;
                    9.877906977, 36.88372093

The reason is that when you create a sequence (or list or set) by appending to an existing sequence, each assignment reconstructs the entire sequence. Of course, the time for each such reconstruction is proportional to the length of the sequence. Thus, the time for the whole loop (just for this reconstruction minutiae, not for the actual desired computation) is proportional to the square of the number of loop iterations. Using the standard asymtotic-order ("big O") notation, we say that the time for this process is O(n^2) (where n is the number of iterations).  Assuming that the actual computation time for each sequence element is independent of its position in the sequence (which is very often the case, and it's certainly true in this case), that time is O(n) in total. It's mathematically certain that there's a value of n, say n_c, such that the O(n^2) process will dominate the O(n) process for all n > n_c. In more-technical mathematical language, for any A > 0, B > 0,

limit(A*n^2 / (B*n) , n= infinity) = infinity.

@2cUniverse Your new version can be coded like this:

iremFrac2:= proc(n, d, b) local N:= irem(n,d), r:= irem(N*b, d), s; 
    [if igcd(b,d) = 1 then 
        (do [[r,r], [r, (s:= irem(b*r, d))]], [[r,s], [s,s]] until (r:= s)=N)
    else
        print("Denom-Base not coprime")
    fi]
end proc
:

The N:= irem(n,d) is needed because if n >= d then the loop will never reach n.

Oh, I just realized that this may have not been obvious before: The number of loop iterations is always MutiplicativeOrder(b,d), but dividing that number by two isn't always relevant.

@sursumCorda GMP is open source and publically licensed. You can find numerous websites with documentation of the whole package or of invidual components such as "GMP isprime". The Maple help page ?GMP contains a link to the source code as it's been modified for use in Maple.

@Mike Mc Dermott Your example is too small for a difference to be shown. For any expression of practical size for which I've used optimize(..., tryhard) over the past 20 years, it has made a great improvement in the optimization. The extra time taken by the optimizer has always been trivial compared to the improvement in execution time. Here's an example:
 

Download OptimizeTryhard.mw

@sursumCorda That's a good find on your part. Here are the new timings using option compile on all Iterators, and running twice so that the compilation time isn't counted.

LL:= [seq]([$0..k], k= 1..8): #So, product has 9! 8-tuples.
for k to 6 do cp[k]:= CodeTools:-Usage(CP[k](LL)) od:
combinat:-cartprod:
memory used=189.79MiB, alloc change=0 bytes,
cpu time=2.66s, real time=2.15s, gc time=703.12ms

Iterator:-CartesianProduct:
memory used=173.28MiB, alloc change=5.64MiB, 
cpu time=2.66s, real time=2.09s, gc time=796.88ms

Iterator:-MixedRadixTuples:
memory used=106.80MiB, alloc change=2.87MiB, 
cpu time=2.56s, real time=2.29s, gc time=390.62ms

Iterator:-MixedRadixGrayCode:
memory used=106.80MiB, alloc change=0 bytes, 
cpu time=2.69s, real time=2.36s, gc time=421.88ms

Iterator:-MultiSeq:
memory used=386.46MiB, alloc change=0 bytes,
cpu time=4.70s, real time=3.46s, gc time=1.72s

Carl's own cartesian product iterator:
memory used=138.49MiB, alloc change=-8.50MiB,
cpu time=1.89s, real time=1.62s, gc time=359.38ms

@2cUniverse My ideas regarding remainder -1 weren't intended for all moduli. 97 is prime; 49136 is not. The idea will work for moduli of the forms p^k and 2*p^k for odd primes p. If the remainder -1 occurs, then it corresponds to d-1 and is obviously the maximum of the list. But for a highly composite modulus such as 49136, it's unlikely that -1 will occur as a remainder for some arbitrary base b.

Your procedure remainders is very inefficient because it builds a sequence by appending, and it also has a few other much-more-minor inefficiencies. Here's an improvement:

remainders:= proc(n, d, b) local n1:= irem(n,d), r:= n1; 
    `if`(igcd(b,d) = 1, [r, (do r:= irem(b*r, d) until r=n1)][..-2], FAIL)   
end proc:

The remainder 1 is usuallly put at the end of the list R and not the beginning, so that R[k] = b^k mod d. If you do it that way, the procedure can be simplified to

remainders:= proc(n, d, b) local n1:= irem(n,d), r:= n1; 
    `if`(igcd(b,d) = 1, [do r:= irem(b*r, d) until r=n1], FAIL)   
end proc:

@sursumCorda To my great surprise, combinat:-cartprod is the best of those 5, and Iterator:-CartesianProduct is by far the worst. (My timings below do not include the once-per-session compilation times for the Iterators.) To celebrate, I rewrote combinat:-cartprod and made a significant improvement in its time. Another difference among the methods is that they give the results in different orders.

Using my new CartProd iterator, I reduced the time for A161786__2(10^6, 10^6) to under 6 seconds and the time for A161786__2(1, 2*10^6) to under 7 seconds!

restart
:
CP[1]:= LL-> local nx:= combinat:-cartprod(LL)['nextvalue'];
    {to mul(nops~(LL)) do nx() od}
:
CP[2]:= LL-> local c;
    {for c in Iterator:-CartesianProduct(LL[]) do [seq](c) od}
:
CP[3]:= LL-> local c, k, j;
    {for c in Iterator:-MixedRadixTuples(nops~(LL)) do 
        [for k,j in c do LL[k][j+1] od]
    od}
:
CP[4]:= LL-> local c, k, j;
    {for c in Iterator:-MixedRadixGrayCode(nops~(LL)) do 
        [for k,j in c do LL[k][j+1] od]
    od}
:
CP[5]:= LL-> local c, k, j;
    {for c in Iterator:-MultiSeq(<[1$nops(LL)]>..<nops~(LL)>) do
        `?[]`~(LL, [seq](`[]`~(c)))
    od}
: 
#My improvement of combinat:-cartprod:
CartProd:= proc(LL::list({list, set})) 
option `Author: Carl Love <carl.j.love@gmail.com> 2023-Oct-24`;
local 
    n:= nops(LL), N:= nops~(LL), J:= rtable([1$n-1, 0], datatype= integer[4]),
    M:= rtable((1..n, 1..max(N)), [op]~(LL))
;
    proc() local i:= n, j;
        while J[i] = N[i] do J[i--]:= 1 od;
        J[i]++;
        [for i,j in J do M[i,j] od]
    end proc,
    mul(N)
end proc
:
CP[6]:= proc(LL) local (nx, N):= CartProd(LL);
    {to N do nx() od}
end proc
:
LL:= [seq]([$0..k], k= 1..8):

#I actually ran these tests twice after the restart 
#so that the Iterators would be compiled:
for k to 6 do cp[k]:= CodeTools:-Usage(CP[k](LL)) od:
memory used=189.79MiB, alloc change=2.77MiB,
cpu time=1.89s, real time=1.89s, gc time=0ns

memory used=407.77MiB, alloc change=5.64MiB,
cpu time=3.53s, real time=3.54s, gc time=0ns

memory used=106.80MiB, alloc change=2.87MiB, 
cpu time=3.25s, real time=2.37s, gc time=1.03s

memory used=106.80MiB, alloc change=5.64MiB,
cpu time=2.16s, real time=2.16s, gc time=0ns

memory used=386.45MiB, alloc change=5.64MiB, 
cpu time=2.92s, real time=2.93s, gc time=0ns

memory used=140.25MiB, alloc change=-19.78MiB,
cpu time=2.67s, real time=1.68s, gc time=1.17s

nops({entries}(cp, nolist));  #Show that all results are the same
                               1

Your worksheet shows two procedures for this---one named rho and one named InversePtMobius. Since those two procedures obviously give identical results, I don't understand what exactly you're asking.

@2cUniverse Since it's obvious that your list is the powers of 2 modulo 97, here are two more methods based on that:

NumberTheory:-ModularLog(-1, 2, 97);

NumberTheory:-MultiplicativeOrder(2, 97)/2;

Both of these return 24. You were getting 25 because you started your list at 2^0 rather than 2^1. The MultiplicativeOrder is likely more efficient because it makes use of the special status of -1 as the primitive square root of 1. The 2 in the denominator is to get the logarithm of -1, and it's independent of the base, 2, or the modulus, 97.  

I think that in your first line, you tried to enter a formula for h, but it got replaced by characters that look like rectangles with Xs in them.

@zenterix Using string constants as keys is robust and enhances code readabiliy. It seems unlikely to me that anyone with even a tiny knowledge of any standard programming language would think that "chevron2A" could have any value other than the string constant that it manifestly appears to be. And, unlike a variable, it's impossible for you to accidentally assign it a value, which would damage its use as a key.

Gradually over the years, more-and-more Maple library code has been using string constants as keys and keyword values.

In this case, your keys (aka indices) such as chevron2A are the escaped locals, and trying to use them at the worksheet level causes essentially the same problem as described in my 2nd paragraph above. As @acer said, making them module exports ameliorates the problem. (Indeed, module exports are essentially locals designed to be used in escaped form.) However, from the code that you've shown, I see no reason that they should be variables at all, whether local or export. Instead I recommend that they be strings, for example "chevron2A". This would correct all the problems that you've mentioned. 

I love Halloween and pumpkin carving, so I thank you for this, and I vote up. I just have a small numeric nit to pick; we are mathematicians afterall! Halloween/Samhain is the beginning of the darkest quarter of the year (north of the Equator), not the "darker half". The celebration of the end of that quarter is known as Groundhog Day / Imbolc / Candlemass, etc.

@Carl Love I streamlined the above code down to 26 lines and added over 50 lines of comments. Let me know if there's any question about how this code works or about why I chose any particular technique.

(* A161786__2
This procedure returns the primes from ithprime(m) to ithprime(m+n-1), inclusive, that hav exactly D
copies of any "digit" in their base-b representation. Henceforth "digit" always means base-b digit.

Defaults based on  _OEIS_ "A161786":
	The base b defaults to 10, but can be any nonunit positive integer.
	The critical digit count D defaults to 4, but can be any positive integer.

My method is
	1. Use simple combinatorial tools from the Iterator package to construct a sequence of
	digits with the desired count.
	2. Use simple arithmetic to compose the sequence into an integer.
	3. If the integer is in the desired range and prime, save it for return.

This code doesn't use strings or any other characater-based representations in any way.
It doesn't break down integers into digits; rather, it builds large integers from digits.
*)
A161786__2:= proc(m::posint, n::posint, b::And(posint, Not(1)):= 10, D::posint:= 4, $)
option `Author: Carl Love <carl.j.love@gmail.com> 2023-Oct-13`;
uses It= Iterator, NT= NumberTheory;
local  #*1 
	(p0,p1):= op(ithprime~([m, m+n-1])), ds:= ilog[b](p1), D1:= ds+1-D, bs:= {$0..b-1},
	B:= Array(0..ds, k-> b^k), rep1:= add(B), C, BC, rep, d, Rep, k, n1,
	d0:= select(x-> igcd(x,b)=1, bs), mids:= bs$D1-2, d1:= {$map(iquo, p0..p1, B[ds])},
	Prime?:= eval(`if`(p1 < 5*10^9, gmp_isprime, prime)),
	Accept:= p-> if p0 <= p and p <= p1 and Prime?(p) then p fi,
	DigitSets:= rcurry(op@eval, d= ()) @ 
		`if`(D1=1, ()-> [`if`(C[1]=0, d0, bs)], ()-> [`if`(C[1]=0, d0, bs), mids, `if`(C[D1]=ds, d1, bs)])
;
	if D1 <= 0 then `if`(D1 < 0, {}, Accept~(`if`(D>1, {rep1}, bs)))
	elif 0 in d1 then n1:= NT:-pi(B[ds]); thisproc~([m, n1+1], [n1-m, m+n-n1], b, D)
	else  # main loop:
		{for C in It:-Combination(ds+1, D1) do
			rep:= rep1 - add((BC:= [seq](B[[seq](C)])));
			(for d in `if`(C[1]=0, bs, d0) minus `if`(C[D1] = ds, {}, {0}) do
				Rep:= d*rep;
				(for k in It:-CartesianProduct(DigitSets()) do
					Accept(Rep + inner([seq](k), BC))
				od)
			od)
		od}
	fi[]
end proc

(* Code comments:

Locals:
	Static:
		p0:	smallest acceptable prime;
		p1:	largest acceptable prime;
		ds:	digit count of p1, counting the units digits as "zeroeth";
		D1:	number of digit positions NOT being specially counted;
		bs:	set of possible single digits;
		B:	array of powers of b such that B[k] = b^k;
		rep1: integer with desired number of digits, all of which are 1;
		d0:	set of allowed units digits for a multidigits prime;
		d1:	set of allowed leading digits (possibly including 0) in p1..p2;
		mids: middle digits: sequence of max(D1-2, 0) copies of bs;
		
	Variable:
		C:	iterates over all D1-combinations of digit positions (0-relative)
			(C is returned by Iterator:-Combination as a sorted vector);
		BC:	sorted list of powers of b corresponding to the postions in C;
		rep: integer with D 1-digits and D1 0-digits with the positions of the 0s given by C;
		d:	iterates over the allowed values for the specially counted digits, taking into
			account whether their positions are first and/or last;
		Rep:	d*rep, thus the digits become all D d's and D1 0s in some order;
		k:	iterates over all digit sequences that can replace the 0s in Rep
			(k is returned by Iterator:-CartesianProduct as a vector in digit-position order);
		n1:	new value of n potentially used for recursive call;

	Procedures:
		Prime?:
			Replaces isprime with the more-efficient undocumented builtin gmp_isprime if all
			its inputs will be < 5 billion;
		Accept:
			If its input is in the desired range and prime, return it;
		DigitSets:
			Constructs the sequence of allowed-digit sets to pass to Iterator:-CartesianProduct.
			C and d are unassigned at the time DigitSets is assigned, but they'll be assigned
			when it's called.
			
Body:
	The 1st branch of the if-statement is to handle cases where there are no digits other than those
	being specially counted. It's easier to handle these as separate cases.

	The 2nd branch handles the case where there are fewer digits in p0 than in p1. Since the main loop
	requires that their lengths be equal, this branch uses recursion to split the cases by length.

	Main loop:
		Consider all possible D1-combinations C of digit positions:
			BC is the powers of b corresponding to those positions;
			rep has all digits 1 or 0.
			Consider all possible digit values d for the specially counted value:
				Rep is rep with 1s replaced by d;
				Consider all possible digit combinations k to replace the 0s in Rep:
					`inner` is an undocumented builtin command to compute the inner- (aka dot-)
					product of lists as if they were vectors;
					so Rep + inner(..) is Rep with its 0s replaced by the digits in k. 
*)
:		

Download DigitCountPrimes.mw

@Jamie128 Yes, it can be animated by just making a few small changes to the code above:

R:= 2:  #max r of featured surface 
spokes:= 24:  #number of radial gridlines in base

plots:-animate(
    plot3d, #the feature
    [
        [r, theta, (1+k)*(csc(theta)/r)^2], r= 0..R, theta= -Pi..Pi,
        coords= cylindrical, style= surface, shading= zhue, transparency= 0.15
    ],
    k= -1..1,
    background= [
        plot3d(  #the base
            [r, theta, 0], r= 0..3/2*R, theta= -Pi..Pi, coords= cylindrical,
            grid= [7, 2*spokes+1], color= COLOR(RGB, 0.97$3), thickness= 2.5, glossiness= 0
        ),
        plots:-textplot3d(  #the spoke labels
            [seq](
                [7/4*R, i, 0, i],
                i= (2*Pi/spokes)*~([$0..spokes-1] -~ iquo(spokes,2))
            ), 
            coords= cylindrical, font= ["times", 10]
        )
    ],
    view= [map(`*`, -1..1, 2*R)$2, 0..20], axes= framed, 
    orientation= [40,40,0], labels= [r$2, z], labelfont= ["times", 16], 
    axis[1,2]= [tickmarks= [-R, R]]
);