Strangely, the faster gmp_isprime is a "undocumented protected name". I'm not sure why it wasn't officially documented.

I also try to use a faster primes generator; in theory, the following program should be quite efficient: (Assume that one has installed the most recent Primesieve under E:\Primesieve\bin\.) 

A161786__3 := proc(m::nonnegint, n::posint := 1, b(*ase*)::And(posint, Non(1)) := 10, o(*ccurrence*)::posint := 4, ` $`)::'Vector'(prime):
    description "Dependency: primesieve winget package.":
	options hfloat:
	local p::string := ifelse(n = 1, "1", ssystem(cat("\"E:/Primesieve/bin/primesieve.exe\" ", n - 1, " -n -p"))[2]), vec := rtable(1 .. 0, datatype = integer[4], subtype = Vector[column]), numberDecompose := (radix, uintlen) -> (thisproc(args) := convert(['tmp' = 'd', [''irem'('tmp', radix, ''tmp'')' $ uintlen]], procedure, 'locals' = ['tmp'::posint], 'params' = ['d'::posint])), v::integer[0 .. b - 1], C(*ounter*);
	if b = 10 then
		to m do
			p := ssystem(cat("\"E:/Primesieve/bin/primesieve.exe\" 1 ", p, " -n -p"))[-1];
			if member(o, rhs~({StringTools:-CharacterFrequencies(p, 'digit')})) then
				vec ,= parse(p)
			fi
		od
	elif b = 2 then
		to m do
			p := ssystem(cat("\"E:/Primesieve/bin/primesieve.exe\" 1 ", p, " -n -p"))[-1];
			if member(o, rhs~({StringTools:-CharacterFrequencies(Bits:-String(parse(p)), 'binary')})) then
				vec ,= parse(p)
			fi
		od
	elif b = 16 then
		to m do
			p := ssystem(cat("\"E:/Primesieve/bin/primesieve.exe\" 1 ", p, " -n -p"))[-1];
			if member(o, rhs~({StringTools:-CharacterFrequencies(convert(parse(p), 'hex'), 'hdigit')})) then
				vec ,= parse(p)
			fi
		od
	elif b = 8 then
		to m do
			p := ssystem(cat("\"E:/Primesieve/bin/primesieve.exe\" 1 ", p, " -n -p"))[-1];
			if member(o, rhs~({StringTools:-CharacterFrequencies(sprintf("%o", parse(p)), 'octal')})) then
				vec ,= parse(p)
			fi
		od
	elif b = 4 then
		to m do
			p := ssystem(cat("\"E:/Primesieve/bin/primesieve.exe\" 1 ", p, " -n -p"))[-1];
			if member(o, rhs~({StringTools:-CharacterFrequencies(String(op(Bits:-Split(parse(p), 2))), "0123")})) then
				vec ,= parse(p)
			fi
		od
	elif b = 2**(local k := ilog2(b)) then
		if k < 8 then
			to m do
				p := ssystem(cat("\"E:/Primesieve/bin/primesieve.exe\" 1 ", p, " -n -p"))[-1];
				if member(o, rhs~({StringTools:-CharacterFrequencies(convert(Bits:-Split(parse(p), k) +~ 1, 'bytes'))})) then
					vec ,= parse(p)
				fi
			od
		else
			to m do
				p := ssystem(cat("\"E:/Primesieve/bin/primesieve.exe\" 1 ", p, " -n -p"))[-1];
				C := table('sparse');
				for v in Bits:-Split(parse(p), k) do
					++C[v]
				od;
				if member(o, {entries(C, 'nolist')}) then
					vec ,= parse(p)
				fi
			od
		fi
	elif b < 2**8 then
		to m do
			p := ssystem(cat("\"E:/Primesieve/bin/primesieve.exe\" 1 ", p, " -n -p"))[-1];
			if member(o, rhs~({StringTools:-CharacterFrequencies(convert(numberDecompose(b, evalhf['hfloat'](trunc(log[b](p)) + 1))(p) +~ 1, 'bytes'))})) then
				vec ,= parse(p)
			fi
		od
	else
		to m do
			C := table('sparse');
			p := ssystem(cat("\"E:/Primesieve/bin/primesieve.exe\" 1 ", p, " -n -p"))[-1];
			for v in numberDecompose(b, ilog[b](p) + 1)(p) do
				C[v]++
			od;
			if member(o, entries(C, 'output = Array')) then
				vec ,= parse(p)
			fi
		od
	fi;
	vec
end:

But unfortunately, the overhead of communication between Maple and Cmd (i.e., ssystem) appears to seriously weaken its advantage … I hope there will be a faster generator/iterator for prime numbers in a furture Maple release.

@Carl Love There is no problem in recent versions, but since `inttrans` is a "table-based package" (and the use of “inttrans:-invlaplace” is invalid in legacy Maple), for backward compatibility, I think that using “inttrans[':-invlaplace']” is somewhat better. 

Besides, the documentation of colondash is outdated. 
Member Selection Operator says that 

Using Packages also says that

Actually, those old packages are now compatible with the use of binary `:-` operator.
A few years ago, there used to be an option (Was this information helpful?—No: Tell us what we can do better.) at the bottom of the help page, yet today it disappears. 

@Carl Love Oops, I just found it and edited my reply. Thanks.

@Carl Love I believe that this is related. 

restart;
libname := subs(FileTools:-JoinPath([kernelopts(':-toolboxdir' = "Physics Updates"), "lib"]) = NULL, [(tmp := libname)])[]:
inttrans[':-invlaplace'](1/(s + a), s, t);
                           exp(-a t)

libname := tmp:
forget(inttrans[':-invlaplace']);
inttrans[':-invlaplace'](1/(s + a), s, t);
                            exp(a t)

In addition, 

> inttrans['invlaplace'](inttrans['laplace'](ln(t^2+1), t, s), s, t); # Maple 2023.1 
 = 
                              / 2    \
                          2 ln\t  + 1/

So Maple thinks that ㏑(t²＋1)＝2㏑(t² + 1)…? 

Additional examples: 
 

Download ilaplace.mw

(1) is correct, but (2), (3), (4), and (5) are wrong.
And for uknown reason, Maple cannot calculate (6), (7), (9), and (10).

@Thomas Richard Right. But consequently, if one has installed Python before, one will have to add a multiple Python on the machine (when installing Maple). Will such a Python installation be optional in future Maple releases?

I think that you mean: 

Python:-ImportModule("import math, numpy as np");
Python:-DefineFunction("def Concentration_calculation(C_0, Q, V_r, m_b, rho, R, Gamma_i, delta_t=1):\n\tt = np.arange(0, 360*60, delta_t)\n\tC_i = [C_0]\n\tfor i in range(len(t)-1):\n\t\tdC = -(Q/V_r)*(1-math.exp(-3*m_b*math.sqrt(Gamma_i/(rho*t[i+1]))/(math.sqrt(math.pi)*Q*R)))*C_i[i]\n\t\tC_i.append(C_i[i] + dC*delta_t)\n\treturn t, C_i"):

@Carl Love Many thanks. This is an impressive procedure. (At least the computation can be done in one-third of a minute!) 

The fourth parameter `D` being 0 appears to mean that not all digits can be used, but there exists a bug in your program:

A161786__2(1, 1, 10, 0); # though it is not hard to fix it 
Error, (in A161786__2) Array index out of range

.

@dharr Thanks. In practice, I find that most use of “convert(p, 'base', 10)” can be simply replaced by faster “[`convert/base`:-MakeSplit(length(p), 1, 10)(p)]”, but to use it, one has to set "kernelopts('opaquemodules' = false):" in advance (which is somewhat dangerous). Besides, after converting integers to strings, I also try to call StringTools:-RegMatch directly, but Maple's regular expression engine seems to be not optimized enough ("Regular expressions" are extremely arcane, but they are known for being very powerful as well.). 
As for the `ithprime`, I believe that its inefficiency is because it's "implemented in the high-level Maple programming language". But I don't know why such a fundamental function (for a math software) is not built into the core of the Maple computation engine …. 

@mmcdara Thanks. I think that your `B1` and `B2` can be obtained more directly (using NumberTheory:-PrimeCounting(N1) + 1 and NumberTheory:-PrimeCounting(10*N2 - 1)).
But in this question, the desired interval is p_1000000～p_1999999 (corresponding to B1 = 10^6 and B2 = 2*10^6 - 1). (I consider testing this range because it may be suited for a benchmark problem—neither too large nor too small.) Anyway, though I ask for a large-scale problem, I'm happy to find that your code performs better on small numbers.

How about 

RealDomain:-solve({x + abs(y + 1) = 1, y + abs(z + 2) = 1, z + abs(x - 2) = 1});
       {x = -abs(y+1)+1, z = -abs(y+1), -1 <= y, y <= 1}

@vv For example, 

But one can easily check that (x, y)＝(78917749904961, 78919591651475) is better. Regrettably, I have no time to wait for a larger R. (Note that minimizing x＋y will lead to a bigger solution (x, y) = (158339291, 112245157230352).)

@Carl Love Here is the Maple code: 

mfoo := proc(` $`)::posint:
	local bound::posint, s::object := Session('logic' = "QF_LIA");
	uses SMTLIB;
	s:-Assert(And((154*x + 69*y) - 0 = 7**3*m__0, (13*x + 716*y) - 0 = 13**3*m__1, (23*x + 3059*y) - 0 = 23**3*m__2, x > 0, x < y));
	for bound do
		s:-Push();
		s:-Assert(y = bound);
		ifelse(s:-Satisfiable(), RETURN(bound), s:-Pop())
	od
end:

CodeTools:-Usage(mfoo(), iterations = 1);
 = 
memory used=0.97GiB, alloc change=0 bytes, cpu time=88.03s, real time=88.82s, gc time=2.95s

                             18837

And a practically word-for-word translation: 

Python:-Verbose((*not *)false);
ssystem(sprintf("PowerShell -Command \"& '%spython.exe' -m pip install -U --force-reinstall z3-solver --no-cache-dir\"", StringTools:-Subs("bin" = "Python", kernelopts(bindir)))): # Remember to uninstall `z3-solver` later! 
Python:-Import("from z3 import *"):
Python:-DefineFunction("def pfoo() -> Int:\n\tx, y = Ints('x y'); m = IntVector('m', 3)\n\ts = SolverFor(\"QF_LIA\")\n\ts.insert((154*x + 69*y) - 0 == 7**3*m[0], (13*x + 716*y) - 0 == 13**3*m[1], (23*x + 3059*y) - 0 == 23**3*m[2], x > 0, x < y)\n\tbound = 0\n\twhile True:\n\t\ts.push()\n\t\ts.insert(y == (bound := bound + 1))\n\t\tif s.check() != sat:\n\t\t\ts.pop()\n\t\telse:\n\t\t\treturn bound"):
CodeTools:-Usage(Python:-EvalFunction("pfoo"), 'iterations' = 1);
memory used=1.20KiB, alloc change=0 bytes, cpu time=31.00ms, real time=8.60s, gc time=0ns

                             18837

There exist some other state-of-the-art SMT solvers (e.g., Yices 2, Boolector, Bitwuzla, and cvc5); nevertheless, the first three are not convenient to use on Windows (whereas the last one provides pre-built binaries for Windows). If I have time, I shall test cvc5 (whose input can be obtained by SMTLIB:-ToString) in Maple. 

@Carl Love Thanks for your advice.
As I see it, the program should not be too "special" since if the OP plans to additionally minimize x＋y instead of the current y (Note that their results are different. One is (x, y) = (15893, 18837); the other is (x, y) = (1103, 26390).), the objective expression will no longer be a multiple of 7×13 (though some other tricks will remain available). So I attempt to concentrate on defining OP's central task and write code that captures concepts rather than burning in a particular generation of specific algorithms. 

Anyway, I believe that there is something wrong with Maple's SMTLIB package. According to 

showstat(SMTLIB:-Execute, 1);

SMTLIB:-Execute := proc(smt::string, {getassignment::truefalse := false, getmodel::truefalse := false, getproof::truefalse := false, getunsatcore::truefalse := false, getvalue::sequential(name) := [], solver::Or(identical(maple,z3),string) := ':-z3', z3options::sequential(equation) := []})
   1   if solver = (':-z3') then
           ...
       elif solver = (':-maple') then
           ...
       else
           ...
       end if
end proc

Maple by default makes use of the external Z3 prover. However, even the naïve Python analogue runs considerably faster than the raw Maple implementation: 

What happened in Maple?! Maybe certain overheads of linking and calling external routines become critical here. (Recall that Python's native loops are known for their poor performance!)

@Axel Vogt I'm not sure what happened. On my computer: 
 

Download question237110.mw

This is far slower than @vv's one. (Why?)