Carl Love

## 26827 Reputation

11 years, 281 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

## Trivial solution filtering...

I often need to filter a set/list of solution sets/lists to separate those that contain trivial equations. It can be done like this:

s:= [{f0(t, r) = 0, n(t) = n(t)}, {f0(t, r) = 0, n(t) = 0}, {f0(t, r) = 0, n(t) = 0}];
selectremove(t-> ormap(evalb, t), s);

[{f0(t, r) = 0, n(t) = n(t)}],  [{f0(t, r) = 0, n(t) = 0}, {f0(t, r) = 0, n(t) = 0}]

This code would be identical regardless of whether the solutions are sets or lists.

## Escaped local...

[This Answer was based on the original Question alone, before the OP posted any code. The Replies were written after the code was posted.]

I suspect that this is the correction for your problem: Let's suppose that you have a procedure that creates a table assigned to a variable T that is a local of P and you want that table to be the return value of P. Then due to the last name evaluation property of T, you should make the return value eval(T) rather than simply T. You should do this inside P. No second argument for the eval is necessary, but after you get this working with the simple eval(T), you might be able to make a slight efficiency improvement using an integer second argument for the eval depending on the what type of things are stored in T.

I suspect that this is the cause of your problem: Suppose that returns T without using eval. You make  your call P(...), and the displayed output makes it seem as if T has been returned. But that is still the T local to P. That's called an escaped local. (Although "escaped" might make them sound dangerous, they can be useful.) It is a distinct variable from the global T, which can use outside of P, such as in a worksheet. The two Ts have no connection whatsoever; they merely have names that are coincidentally spelled the same. I suspect that you've been trying to use something akin to T[k] from your worksheet level.

You can assign the returned table to any variable: If you use eval as I described in the first paragraph, and you make your call to P as T:= P(...), then you can index T at the worksheet level. But you can use any variable for this purpose; it needn't be T.

I am certain that these things have nothing to do with your problem: It makes absolutely no difference where the procedure P is defined (.mpl file, library archive, or directly in your worksheet). It makes no difference where you call P from (worksheet or from another procedure).

## Replace factorials with binomial...

When I read @dharr 's proposed GAMMA workaround this morning, before @acer looked into it, I knew that it couldn't possibly work in the desired way (i.e., under evalhf) because GAMMA(401.) is still well above the maximum value of a hardware float. But replacing the three factorials with a single call to binomial does work under evalhf:

F:= (n,k,p)-> 100*binomial(n,k)*p^k*(1-p)^(n-k):
evalhf(F(400, 5.1, 1/100));

15.1978021802672245

CodeTools:-Usage(plot(F(400, q, 1/100), q= 0..20, numpoints= 1000));
memory used=2.88MiB, alloc change=0 bytes, cpu time=31.00ms, real time=40.00ms, gc time=0ns

## Option continuous...

int(diff(y(x), x), x= 1..2, continuous)

## Option 'tryhard'...

If you use option tryhard, then it uses a completely different algorithm and that limitation doesn't apply.

## Iterators, no characters, any base...

I have a completely different approach that doesn't use any strings or character-based representations, works in any base, works for any desired "digit" count, and is faster than the fastest procedure from the original Question by a factor of 2 (although it doesn't quite make it to the desired "quarter of a minute").

Here's the procedure:

```(* a161786__2
This procedure returns the primes from ithprime(m) to ithprime(m+n-1), inclusive, that have
exactly D copies of any "digit" in their base-b representation.
The base b defaults to 10 but can be any nonunit positive integer..
The critical digit count defaults to 4.

My method is
1. Use simple combinatorial tools from the Iterator package to construct sequences of
base-b "digits" with the desired "digit" 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.

This code doesn't use strings or any cther charcater-based representations in any way,
It doesn't break down integers into digits; rather, it builds large integers from digits.

Comments of the form #*n refer to footnotes ar the end of the procedure.
*)
A161786__2:= proc(m::posint, n::posint, b::And(posint, Not(1)):= 10, D::nonnegint:= 4, \$)
option `Author: Carl Love <carl.j.love@gmail.com> 2023-Oct-13`;
uses It= Iterator, NT= NumberTheory;
local
(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), repunit:= add(B), rep, Rep, BC, k, d, C, p, mids:= bs\$D1-2,
d0:= select(x-> igcd(x,b)=1, bs), d1:= {\$map(iquo, p0..p1, B[ds])},
FactorSelector:=
if D1 = 1 then ()-> bs minus {d}
else ()-> map(`minus`, [`if`(C[1]=0, d0, bs), mids, `if`(C[D]=ds, d1, bs)], {d})[]
fi,
Isprime:= eval(`if`(p1 < 5*10^9, gmp_isprime, isprime))  #*1
;
# Short-circuit returns for corner cases:
if D1 < 0 then return {} fi;
if D1 = 0 then
return
if D>1 then `if`(p0 <= repunit and repunit <= p1 and Isprime(repunit), {repunit}, {})
else select(p-> p0 <= p and p <= p1 and Isprime(p), bs)
fi
fi;
if p0 < B[ds] then
return thisproc(m, (p:= NT:-pi(B[ds]))-m, b, D) union thisproc(p+1, n+m-p, b, D)
fi;

# Main code:
{
for C in It:-Combination(ds+1, D1) do
(
for d in bs minus `if`(C[D] = ds, {}, {0}) do
Rep:= d*rep;
(
for k in It:-CartesianProduct(FactorSelector()) do
p:= Rep + inner([seq](k), BC);  #*2
if p0 <= p and p <= p1 and Isprime(p) then p else fi
od
)
od
)
od
}
end proc

(* Footnotes:
*1: gmp_isprime is an undocumented builtin procedure that is faster
than isprime for most small cases.

*2: 'inner' is an undocumented builtin procedure that computes the
inner or "dot" product of 2 lists as if they were vectors.
*)
:```

And the worksheet:

 > restart:
(* a161786__2
 >
 > CodeTools:-Usage(A161786__2(10,10,10,1));

memory used=0.73MiB, alloc change=0 bytes, cpu time=0ns, real time=12.00ms, gc time=0ns

 > P2:= CodeTools:-Usage(A161786__2(10^6, 10^6)):

memory used=2.16GiB, alloc change=0 bytes, cpu time=21.25s, real time=19.99s, gc time=2.64s

 > #Compare with your procedure from the Question: A161786__1 := proc(m::nonnegint, n::posint := 1, ` \$`)::'Vector'(prime): uses ST= StringTools; options no_options: local p := ifelse(n = 1, 1, ithprime(n - 1)), vec := DEQueue();     to m do         if             member(                 4,                 rhs~({                     ST:-CharacterFrequencies(nprintf("%d", (p := nextprime(p))), 'digit')                 })             )         then             vec ,= p         fi     od;     Vector([vec[]], 'datatype' = integer[4]) end :
 > P1:= CodeTools:-Usage(A161786__1(10^6, 10^6)):

memory used=3.47GiB, alloc change=0 bytes, cpu time=42.97s, real time=39.27s, gc time=7.89s

 > evalb(P2 = {seq}(P1));

 >

I am vehemently opposed to any so-called "mathematics" based on the base-10 digits of numbers without generalizing it to all bases. To me, it's "junk math". It looks like OEIS will admit many such worthless sequences.

Your syntax is already perfectly correct, including your placement of derivatives in boundary conditions. To finish solving, replace degree with Pi/180, and set a specific numeric approximation for infty in the Parameters (I used 5). Then do

dsolve(eval({eq1, eq2, eq3, ics, bcs}, {Parameters}), numeric)

## Very simple system...

Your system of ODEs is very simple: four linear equations with constant coefficients for the homogenous parts and with the forcing functions all of the form C[k]*exp(-t) for 4 constants C[k]. The solutions are all very smooth, monotonic functions. Thus, there's very little error even for rk1 (a.k.a. Euler's method).

## Multi-assignment and multi-indexing...

Your can simplify your procedure substantially while correcting the error that you asked about. In the following, I use the optional argument to list the indices that you want to perform this operation on (because it's not clear to me where your index-pattern 1, 2, 4, 5, 7 comes from). The index -1 refers to the last vector position regardless of the actual number of entries.

```oneStep_egcd2:= proc(prev::Vector, curr::Vector, J::list(posint):= [1,2,4,5,7])
if curr[-1] <> 0 then
(prev[], curr[J]):= (curr[], prev[J] - curr[J]*iquo(prev[-1], curr[-1]))
end if
end proc
:
oneStep_egcd2(prev, curr);
```

Maple's multi-assignments happen essentially simultaneously without needing a temporary intermediary variable. For example, compare these:

```(a,b):= (1,2):  (a,b):= (b,a):  a, b;  #multi-assignment
2, 1

(a,b):= (1,2):  a:= b: b:= a:  a, b;   #sequential assignments
2, 2

```

## A better algorithm...

I think that I have a better algorithm for this. I think it has a much better average-case time. It guarantees return of the minimal possible (if there's any y larger than its corresponding x). It does this by stepping sequentially through y, stopping at the first that's greater than its corresponding x. And, it's a much simpler algorithm theoretically.

```MinY:= proc(Coeffs::list(rational), m::list(posint))
local c, M:= ilcm(m), g, X:= 0, Y:= 0;
`mod`:= modp; #Use nonnegative remainders only
M/= (g:= igcd((c:= chrem(Coeffs, m)), M)); #chrem is Chinese Remainder Algorithm
c:= g/c mod M;
to M do until (X:= X+c mod M) < (Y+= g);
`if`(X<Y, [X,Y], FAIL)
end proc
:
# -x_coeff/y_coeff for each equation:
Coef:= [-154/69, -13/716, -23/3059, -2295/4522, -6479/5396]:

# Moduli: The 3rd and 4th have their power reduced because the coeffs are reducible.
M:= [7^3, 13^3, 23^2, 17^2, 29^3]
:
seq(CodeTools:-Usage(MinY(Coef[..k], M[..k])), k= 1..nops(M));
memory used=2.80KiB, alloc change=0 bytes,
cpu time=0ns, real time=2.00ms, gc time=0ns

memory used=3.09KiB, alloc change=0 bytes,
cpu time=0ns, real time=0ns, gc time=0ns

memory used=12.67KiB, alloc change=0 bytes,
cpu time=0ns, real time=0ns, gc time=0ns

memory used=16.86KiB, alloc change=0 bytes,
cpu time=0ns, real time=0ns, gc time=0ns

memory used=95.96MiB, alloc change=0 bytes,
+cpu time=2.55s, real time=2.03s, gc time=796.88ms

[7, 49], [691, 819], [15893, 18837], [1103, 26390], [186957261, 190687133]

```

## Array indices must be integers...

Array indices must be integers. You're trying to use Pi/3 as an index.

## The power series...

A fairly simple power series for the incomplete upper Gamma function is

GAMMA(a,z) = GAMMA(a) - z^a*Sum((-1)^n*z^n/n!/(n+a), n= 0..infinity)

provided that is not a nonpositive integer. Although this series involves a possibly noninteger power of z, that can be factored out of the power series proper, as shown.

## MemoryInUse alternative...

Here's my alternative version of MemoryInUse that runs 68 times faster on my computer:

MemUse:= proc() local r; add(r[3], r= kernelopts('memusage')) end proc:

This measures exactly the same thing as MmaTranslator:-Mma:-MemoryInUse. There is something akin to a "Heisenberg Uncertainity Principle" involved here: It's impossible to measure memory usage without having some effect on it. And of course that's compounded by the fact that the device doing the measuring is the same as the device being measured. So don't expect identical numbers.

It's not incorrect for the memory usage number to sometimes decrease. This is due to "garbage collection"---the recycling of memory addresses whose contents the kernel determines will no longer be needed. Every time that you see a change in the "Memory: " or "Time: " number on the "status bar" (bottom border of your Maple window), a garbage collection has occurred.

Although there is some possibility that "remember tables" (which are the things that the forget command forgets) are causing your memory trouble, I think that that chance is only slight. So you should focus on other causes for now.

I think that you've gotten the mistaken idea that remember tables are primarily used with trig functions. Actually, any piece of Maple code could involve hundreds of remember tables, most of which will be associated with procedures that the user didn't directly call. But, most such remember tables are cleared at every garbage collection (this is called option system).

## No equivalence operator...

The Maple language doesn't have a builtin operator for logical equivalence. I guess that the 2D Input parser translates a variety arrow-type logical symbols into the language's implies operator (which prettyprints as a character like =>).

But, any character that has an HTML code (and thus begins with `&...in ASCII) can be used as an inert infix operator, so what you want prettyprinted can be achieved by

`&iff;`(not (P and Q), ``(not P) or ``(not Q));

The quotes on ``(not P) and ``(not Q) are needed to prevent automatic simplification.

You may want to get rid of some parentheses from the above (at which point the expression will only be good for display, not computation). That can be done like this

use T=Typesetting in
T:-mrow(
T:-Typeset(not (P and Q)),
T:-mo(`&iff;`),
T:-Typeset(``(not P) or ``(not Q))
)
end use;

## Linear regression...

The problem initially seems nonlinear, but it can be easily converted (completely) to a linear-regression problem. By "completely", I mean that the following is the complete original problem, not just some linear approximation of it. The linear form gives you much more control.

 > restart: Digits:= 15: St:= Statistics :
 > Zscore:= proc(V::Vector) uses St=Statistics; local mu:= St:-Mean(V), sigma:= St:-StandardDeviation(V);     mu, sigma, (V -~ mu)/sigma end proc :
 > #Dependent variable: Hvals:= <210.84, 218.19, 219.17, 222.27, 226.35, 228.45, 230.44, 234.65, 239.75>: H__avg, H__sd, H__z:= Zscore(Hvals) :
 > #Independent variable: Qvals:= <1800, 1750, 1730, 1700, 1650, 1635, 1600, 1570, 1500>: Q__avg, Q__sd, Q__z:= Zscore(Qvals) :
 > n:= numelems(Qvals) :
 > #Proposed model, with (apparently) A known beforehand and Q__0, k__f, k__1, and Q__s to be #determined. model:=     subs(         Q= Q__N*Q__sigma + Q__mu,         H__N = (A*(1 - Q/Q__0) - k__f*Q^2 - k__1*(Q - Q__s)^2 - H__mu)/H__sigma     ) :
 > pv:= St:-LinearFit(     [1,Q,Q^2], Q__z, H__z, Q, summarize,     output= [parametervector, residualsumofsquares] );

Summary:
----------------
Model: .52025574e-1-1.0044218*Q-.58528770e-1*Q^2
----------------
Coefficients:
Estimate  Std. Error  t-value   P(>|t|)
Parameter 1    0.0520    0.0414      1.2565    0.2556
Parameter 2   -1.0044    0.0319     -31.5038   0.0000
Parameter 3   -0.0585    0.0325     -1.8025    0.1215
----------------

Solve for the symbolic paarameters in the model by equating coefficients of the powers of Q. We have 3 coefficients and 4 parameters, so I'll just set the value of k__1 (arbitrary choice) to the MathCAD value, 1e-5.

 > Known_params:=     [A, k__1, Q__mu, Q__sigma, H__mu, H__sigma]     =~ [350.156, 1e-5, Q__avg, Q__sd, H__avg, H__sd] ; PV:= solve(     eval(         {seq(pv[1]=~ coeff~(expand(rhs(model)), Q__N, [0,1,2]))},         Known_params     ) );

Compute the sum of the squared residuals. It should be the same as returned by the LinearFit command, allowing for some round off error.

 > add(     eval(         eval(H__N - rhs(model), [PV[2][], Known_params[]])^2,         [H__N= H__z[k], Q__N= Q__z[k]]     ),     k= 1..n );

Compute the sum of the squared residuals using MathCAD's parameter values.

 > add(     eval(         eval(             H__N - rhs(model),             [                 ([Q__0, k__f, Q__s]=~ [6962.32773, 2e-5, 2198.86103])[],                 Known_params[]             ]         )^2,         [H__N= H__z[k], Q__N= Q__z[k]]     ),     k= 1..n );