@David Sycamore You can iterate through the primes like this

N1:= 10^6:  N2:= 10^9:
P:= ithprime(N1-1):
for n from N1 to N2 do
   P:= nextprime(P);
   #blah blah blah
end do:

Throughout the loop, P will equal ithprime(n).

That N2 may not be an upper bound achievable in a reasonable amount of time. Once you find the correct P, you'll be done, right? Then make the code

N1:= 10^6: 
P:= ithprime(N1-1):
for n from N1 do
   P:= nextprime(P);
   #blah blah blah
   if ifoundit then break end if;
   #blah blah blah
end do:

Where ifoundit is some boolean indicator of whether the current P is the one that you want.

Regarding a stored list of primes in Maple: Reading (quickly) through the code of ithprime and its children `ithprime/small` and `ithprime/large`, it appears to me that Maple accumulates knowledge about primes during a given session; that this knowledge is stored in primary memory rather than files; and that this knowledge is lost when your session ends (through restart, for example). It's possible to save this knowledge across sessions by saving `ithprime/global/list` either in a file or a library. The technical details of how to do this are not complicated, but I don't feel like going into them right here right now. (Perhaps someone else will fill in those details.)

Regarding the speed of Maple's primality testing: I think that another post in this thread may have given the wrong impression, perhaps unintentionally, that it is slow. Maple's isprime is extremely fast. But the issue that's relevant to this thread is that knowing that a number is prime doesn't tell you its rank in the primes.

This shows that for numbers on the order of 10^12 (which is what rand() produces), isprime is essentially instantaneous:

restart:
gc():
Ps:= ['rand()' $ 2^12]: #4096 numbers to check
L:= CodeTools:-Usage(isprime~(Ps)):
memory used=8.39MiB, alloc change=0 bytes, cpu time=63.00ms, real time=53.00ms, gc time=0ns
nops(select(evalb, L)); #how many were prime?
                              163

Regarding counting and ranking primes: A very good approximation of the number of primes less than n (thus an approximation of the inverse function of ithprime) can obtained very quickly by

ApproxPi:= proc(n::posint)
   Digits:= max(trunc(evalhf(Digits)), ilog10(n)+3);
   round(Li(evalf(n)) - Li(2.)) #= Int(1/ln(x), x= 2..n)
end proc:

This approximation is much better than the often-used n/ln(n). There is a vast amount of literature available on how good this approximation is. I suggest starting with the Wikipedia article "Prime Number Theorem". And if you have a reliable way to approximate a continuous function's inverse, then the function itself can be approximated with fsolve. For example,

Digits:= 10:
round(fsolve(Li(n) - Li(2) = 10^7));
        179415211

So, ithprime(10^7) is approximately 179,415,211. I'll leave it to you to show that the relative error of this approximation is about 5x10^(-5). (Choosing a value to use for Digits so that fsolve actually returns something that'll round to an integer is a bit tricky.)

Text in plots (in legends, titles, captions, tickmarks, axis labels, etc.) can be neatly formatted by wrapping it with the keyword typeset (as if typeset were a procedure), like this

legend = [typeset(omega[0] = parse(cat(evalf(omega0), 0)))]

The purpose of the parse(cat(evalf[3](omega0), 0)) is simply to get the number to display as 10.0 in an upright font. If you'd be satisfied with simply 10., that part could be shortened to evalf(omega0).

NM, as shown by Tom, your confusion is simply caused by a misunderstanding of the meaning of the Maple reserved word global as well as the English adjective global. Both refer to the highest level, not a higher level; it's not what we call a relative or comparative adjective (i.e., one that can be used with -er or more-). The word lexical is used in Maple (not as a keyword or reserved word, but in documentation, error messages, etc.) to refer to variables that are local to or parameters of a procedure or module at a higher level than the current procedure or module.

My solution uses some simple calculus: If F(x,y) = f(x)*g(y) then

• f(x) = C1*exp(int(diff(F(x,y), x)/F(x,y), x))
• f(x) = C2*exp(int(diff(ln(F(x,y)), x), x))

where C1 and C2 are constants (i.e., expressions depending on neither x nor y). Only one of the two formulas above needs to be used. My procedure (depending on how it's called) will try the first method and, if that doesn't work, try the second method.

My procedure works for all the examples that I could find in this thread, although I'm sure that you could construct some input that it won't work on. It (theoretically) works for any number of factors.

Download Separable.mw

MaplePrimes refuses to display the worksheet above (although you can still download it). It shows the execution of my procedure on all the examples from this thread. Here's the procedure code:

Separate:= proc(
   F::algebraic, V::list(name):= [indets(F, And(name, Not(constant)))[]], 
   {
      nosymbolic::truefalse:= false,
      method::identical(ln, divide, tryhard):= 'divide',
      _pass::{1,2,3,4}:= 1 #for internal use; only for recursive calls
   }
)
local
   Sym:= `if`(nosymbolic, [][], 'symbolic'),
   FF:= `if`(_pass=1, factor(F), F),
   setV:= {V[]},
   K:= [$1..nops(V)],
   ReDo:= subs(
      _P= procname, 
     ()-> 
        if method = 'tryhard' or _pass in {2,3} then
           _P(
              FF, V, 
              ':-method'= `if`(_pass=1, 'divide', 'ln'),
              ':-_pass'= _pass+1,
              ':-nosymbolic'= not nosymbolic
           )
        else
           FAIL
        fi
   ),
   lnF, DlnF, Fs, MF, C, r
;
   if method = 'tryhard' or _pass > 1 then
      userinfo(1, procname, "Trying ", _options)
   fi;  
   if method in {'divide', 'tryhard'} then
      DlnF:= map(v-> diff(FF,v)/FF, V)
   else
      lnF:= simplify(ln(FF), Sym);
      DlnF:= map(v-> diff(lnF,v), V)
   fi;
   Fs:= zip((df,v)-> simplify(exp(int(simplify(df, Sym), v)), Sym), DlnF, V);
   MF:= mul(Fs);
   C:= simplify(FF/MF, Sym);
   if 
      depends(C, setV) and ormap(v-> simplify(diff(C,v), Sym) <> 0, setV) or 
      ormap(k-> depends(Fs[k], setV minus {V[k]}), K) 
   then 
      ReDo()          
   else
      if simplify(C*MF - FF, Sym) <> 0 then
         r:= ReDo();
         if r = FAIL then
            userinfo(1, procname, "Warning: Can't verify equality.")
         fi
      fi;  
      `if`(nops(Fs)=0, [C], [C*Fs[1], Fs[2..][]])
   fi
end proc:

In your simple example, direct substitution of ln@abs for ln via subs or eval (as shown by VV) works. But the idiom for more general cases (even slightly more general) is subsindets or evalindets. (The difference between the two is analogous to the difference between subs and two-argument eval, so they usually produce the same results.) I find patmatch and applyrule unreliable, very difficult to use, and inefficient.

Suppose that I change your requirement so that you want to replace ln with ln@abs only if the argument of ln depends on certain variables. Try these examples:

e1:= 7*ln(arcsin(x)) - ln(x-1)*sin(x)/2 - ln(x+1)/2 + ln(f):
e2:= 7*ln(arcsin(x)) - ln(x-ln(3*ln(y)))*sin(x)/2 - ln(x+1)/2 + ln(f):
Sln:= (e, x::{name, set(name)})-> subsindets(e, 'ln'(dependent(x)), ln@abs@op):
Sln(e1, x); Sln(e2, x); Sln(e2, {x,y});

The second argument to subsindets or evalindets selects the type of subexpression that you want to transform. It's identical to what you'd use as the second argument to indets. The third argument is a procedure or operator that transforms the items selected by the second argument. In this example, I've used ln@abs@op, which is equivalent to e-> ln(abs(op(e))). The op is to extract the argument of the ln expression.

In Maple, square brackets cannot be used instead of parentheses for algebraic groupings. Thus, you need to edit the objective function. Unfortunately, since square brackets do have several other legitimate meanings in Maple, your error was not caught at the time that you entered p.

(This paragraph applies if you want the curve to go exactly through the given points (which is somewhat implied by the word interpolation).) The situation that you present is an extremely well-known drawback of polynomial interpolation. The usual alternative is to use a spline. This is a piecewise curve of low-degree polynomials, usually cubic, and usually made so that the resulting curve is differentiable at the split points (called knots). The drawback of a spline is that unlike a polynomial, it is not infinitely differentiable. So, use Maple command CurveFitting:-Spline.

(This paragraph applies if you only have a finite set of points.) If you just want a curve that gets as close as possible to the points but doesn't necessarily go through them exactly, then many more options are available, but you need to first decide on a template or parameterized formula for the curve, and then start with Statistics:-Fit. For example, you can find the best cubic polynomial that approximates the data. Or if you know that you want an increasing function without inflection points, the model may be y = C*x^P where C and P are parameters to be determined (Statistics:-PowerFit handles this specific case).

(This paragraph applies if you have an essentially infinite number of points; that is, you have some means (perhaps very complicated means) of numerically computing the function values to arbitrary accuracy over some interval.) So, you want some nice, relatively simple, algebraic function that approximates your not-so-nice numeric procedure over some closed and bounded interval. Solutions of this sort are in the numapprox package.

I can confirm that the situation is as you described. I think that this is bug, because the underscore-prefix convention is supposed to apply to automatically generated global symbols, not global symbols specified by the user. Here is a workaround:

OdeTest:= proc(Sol,ODE)
local 
   Names:= 
      [(indets(Sol, And(name, Not(constant))) 
         minus indets(ODE, And(name, Not(constant)))
       )[]
      ],
   _Cn:= [_C||(1..nops(Names))]
;
   subs(_Cn=~ Names, odetest(subs(Names=~ _Cn, Sol), ODE, args[3..]))
end proc:

Now just use OdeTest any place where you use odetest.

If you simply want it to pause indefinitely (rather than for a pre-specified length of time) so that you can do something else and then get back to it, then just use any command that expects user input, such as

readline(terminal);

If you trace the code that ensues from evaluating Ei(1, 1.), you'll see that the relevant code is in

showstat(`evalf/Ei/taylor`);

A simple alternating series approach is used. First, the value of Digits is increased slightly in line 2. The formula for the amount of increase is complicated (but only 1 line), but the end result is that if Digits starts with its default value, 10, then it's increased to 13. Second, there are some initializations in lines 3-6. Third, a very simple for loop in lines 14-17 accumulates the sum of an alternating series until reaching the first term that doesn't change the sum (at Digits=13). That's 15 terms in this case. Finally, the sum is added to Psi(1.) - ln(x), where x=1. but Psi(1.) is constant, independent of x.

The relevant alternating series can be obtained by

convert(Ei(1,x), FormalPowerSeries);

(That's equivalent to the series that Mariusz showed in his Answer above.)

And the code that I detailed above is equivalent to

S:= -Sum((-1)^k*x^k/(k*k!), k= 1..infinity) - gamma - ln(x);
evalf[13](eval(S, [x=1., Sum= add, infinity= 15])):
evalf[10](%);

except that the terms of the series are generated iteratively, so that no actual powering or factorial computations are done, and no foreknowledge of the number of terms is needed.

In the declarations section of foo, include the line:

uses Shorter_name= B;

Then call boo as Shorter_name:-boo.

Never use with inside a procedure or module. Use uses instead.

If you want to actually remove the assumptions on alpha, do

alpha:= 'alpha';

If you simply want alpha to be displayed as alpha rather than alpha~ (without actually removing the assumptions), do

interface(showassumed= 0):

Put the appropriate legend option in each individual odeplot. Then display will merge the legends correctly, even if there's more than one legend in each individual plot. So,

p1:= odeplot(lp1, ..., legend= ["step"]);
...
p2:= odeplot(lp2, ..., legend= ["impulse"]);
display([p1, p2]);

For the code Test that you show, the call must be Test(f) and not Test(f(x)). If you use Text(f(x)), you'll get weird results, but not an error message. It's worthwhile to force an error message like this:

Test:= proc(f::procedure) f(1) + f(2) end proc:

Now the call Test(f(x)) will give you an error message. The key is the ::procedure. The procedure is what Maple calls a type, and f::procedure will force f to be that type.^[1] There's no practical limit on how complicated or intricate a type specification can be.

It's easier to pass a procedure f than an algebraic expression such as f(x). But if you have a compelling need to pass f(x) instead, let me know; it's not that complicated. The complexity is that the procedure Test will either need to be told or need to figure out for itself what the variable x in the expression is. It would be a very bad practice (although Maple would allow it) to code Test to just assume that the variable is global x.

[1] Both f and Test are considered by Maple to be of type procedure. In ordinary mathematical English, we may refer to f as a function. That is fine, but f is not of type function in the way that that word would be used in actual Maple code.

You're confusing QR decomposition, a matrix factorization algorithm, with the QR (eigenvector) algorithm, which may use QR decomposition as an iterative step (but in practice usually doesn't). See the Wikipedia article "QR algorithm". If, for educational purposes, you want to implement the QR eigenvector algorithm in a way that actually uses QR decomposition, the steps for doing that are given in that Wikipedia article. But, note that that is not now considered a practical algorithm.

Please don't direct any followup questions at me if they go on and on with ridiculously verbose and ultimately meaningless code, as is your long-established style.