Carl Love

Carl Love

28055 Reputation

25 Badges

12 years, 356 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

What is Q? Even if you know what Q is, it's clear from your output shown above that you haven't told Maple what Q is. Or were you hoping to get a solution for an arbitrary forcing function Q? [LMFAO] Maybe that's what Maple is trying to do....

I suspect that if a symbolic solution is possible at all, it'll be enormously complex with deeply nested integrals and RootOf expressions using your C__s as coefficients. And what's the point of having an enormously complex symbolic solution that's incomprehensible to the human mind? I'll admit that there is sometimes a point to that, but what's your reason? Chances are good that a parametric numeric solver, specific to this BVP system, would suit your purpose better than a symbolic solution. That numeric solver could be written in 10-15 minutes, including on-screen sliders to change the solution's plot as the parameters are varied.

 

The numerical solution of the original ODE system can be obtained easily, quickly, and straightforwardly:

restart:
F:= [S,E,J,R]:
IVP:= [
   diff(S(t),t) = mu*N - (beta*J(t)+mu)*S(t),
   diff(E(t),t) = beta*S(t)*J(t)-(mu+sigma)*E(t),
   diff(J(t),t) = sigma*E(t)-(mu+gamma1)*J(t),
   diff(R(t),t) = gamma1*J(t)-mu*R(t),
   S(0)=S0, E(0)=E0, J(0)=J0, R(0)=R0
]:
params:= [
   N= 5e7, beta=0.00001, mu=0.02, sigma=45.6, gamma1=73,
   S0=12500000, E0=50000, J0=30000, R0=37420000
]:
sol:= dsolve(IVP, numeric, parameters= lhs~(params), maxfun= -1):
sol(parameters= params):
plots:-display(
   <seq(
       plots:-odeplot(
          sol, [t,f(t)], t= 0..70, numpoints= 5000
       ),
       f= F
   )>, thickness= 3
);

But I don't know how to solve the system after you've applied "homotopy...".

In the code below, I used gamma1 for gamma and J for I because the symbols gamma and I are reserved in Maple for specific constants. If you really want to use gamma and I, it's easy to make a few adjustments to accomodate that.

restart:
Sys:= [
   S(k+1) = (mu*N - beta*Sum(S(l)*J(k-l), l= 0..k) - mu*S(k))/(k+1),
   E(k+1) = (beta*Sum(S(l)*J(k-l), l= 0..k) - (mu+sigma)*E(k))/(k+1),
   J(k+1) = (sigma*E(k) - (mu+gamma1)*J(k))/(k+1),
   R(k+1) = (gamma1*J(k) - mu*R(k))/(k+1)
]:
V:= [S,E,J,R]: #The unknown functions
ICs:= V(0) =~ [S0,E0,J0,R0]: #The symbolic initial values

K:= 100:        #highest value of k to compute and plot
params:= [
   N = 5e7,     #population
   mu = 0.02,     #?
   beta= 1e-5,     #?
   gamma1= 73, #?
   sigma= 45.6,  #?
   S0= 1.25e7,     #S(0)
   E0= 5e4,     #E(0)
   J0= 3e4,     #J(0)
   R0= 3.742e7     #R(0)
]:
assign(eval(ICs, params)); 
NumSys:= subs(Sum= add, unapply(eval(Sys, params), k)):
for _k from 0 to K-1 do assign(NumSys(_k)) od:
plot([seq([seq([k,f(k)], k= 0..K-1)], f= V)], legend= V, labels= [k, ``]);

Use this:

SP:= (n::posint)-> add(mul~(ifactors(n)[2])):
SP(12);

The [2] part remains the same regardless of n. It's an index into the list structure returned by ifactors. See ?ifactors, and note the (slight) distinction between it and the more commonly used ifactor. If you use ifactors  instead of ifactor, the format is more consistent, and you're less likely to encounter exceptional cases. The ifactor form is intended primarily for display, while the ifactors form is intended primarily for subsequent computation.

 

Your line

CC(i, j) := CC(i, j)+Testna(i, 1+(n-1)*3)*(Z(1, n)-Z(1, n+1))

should be

CC(i, j) := CC(i, j)+Testna(i, j+(n-1)*3)*(Z(1, n)-Z(1, n+1)),

the only difference being that I changed 1+... to j+....

Your whole thing could be made much simpler, and made to work for arbitrary K, like this:

K:= 4:
Z:= Vector[row](K+1, n-> z-2*(n-1));
DZ:= Z[..-2] - Z[2..]; # minus Delta Z
Testna:= LinearAlgebra:-RandomMatrix(K-1, K*(K-1)):
CC:= Matrix((K-1)$2, (i,j)-> add(Testna[i, j+(K-1)*(n-1)]*DZ[n], n= 1..K)); 

[The third line of the above code has been corrected due to an error discovered by Tom, below.]

So, there's actually no explicit loops required because the looping is built-in to the Matrix and add commands, and it's usually more efficient as well as much more convenient to rely on these built-in features.

The preferred form of Vector and Matrix indexing uses square brackets rather than parentheses. So, it's Z[n] and Testna[i,j] rather than your Z(1,n) and Testna(i, j). Your parentheses also work, in this case, but I think that it'd be safer for a beginner to stick to square brackets, as the parentheses form has a special meaning that I'd rather not go into here.

Using Maple's Standard GUI (which you're probably already using anyway), try this:

(xrange,yrange):= (-4*Pi..4*Pi, -2..2):
Explore(plot(A*sin(B*x), x= xrange, view= [xrange,yrange]), A= yrange, B= -2..2);

and play with the A and B sliders that'll appear on the screen beneath the plot.

This specific command doesn't animate, but I wanted you to "explore" it first before we animate it. And I think that any pedagogic purpose served by an animation is better served by the above command, assuming that your purpose is pedagogic. Explore has a vast number of options, including many for animation. See ?Explore.

 

You can use Diff as an inert form instead of diff. Replacing diff with %diff should've worked though.

When you call it an "operator", I think that you mean that you intend to "operate" on an unknown function f, making a new function. Is that correct?

It is done by using the output option of Usage, as in

Time:= CodeTools:-Usage(..., output= cputime);

A global variable and a local variable with the same apparent name are still completely different variables to Maple. Surely you know that already (at least on some level) because it's the same in virtually all computer languages. What's your reason for making C1 local to foo? I wouldn't try doing pattern matching with local variables at all. I could probably get something to work, but I think that it'd get rather convoluted.

The command that you need is subs (short for substitute).

eq:= x(t) = _C1*sin(sqrt(k)*t/sqrt(m)); #You don't need to type this.

subs([ _C1= A,  sqrt(k)= W*sqrt(m)], eq);

Note that the first substitution is rather obvious, but the second is more subtle. It has to be that way because the sqrt(k) and sqrt(m) are split in the original. The left sides of substitutions are rather limited (they need to be vaguely "atomic" ), but the right sides can be anything.

To other readers: Yeah, I know that I'm not using "atomic" in the Maple sense of that word defined at ?type,atomic.

The attached worksheet shows you how you can carefully use the printlevel environment variable to begin to research the dsolve code that does the parsing and syntax error checking.

The reason that I use such extreme care in doing this is that the amount of information that you get with using printlevel can easily become overwhelming---hundreds or thousands of screenfuls---which may leave you with no other option than killing your whole Maple session.
 

Step 1: Make up an ODE that you think that dsolve will reject, and make sure that that is the case.

#End all top-level lines with : during this process:
#Always put restart in its own execution group.
restart:

dsolve(diff(y(x),x) = y(x,x)):

Error, (in dsolve) found the indeterminate function y with different arguments {y(x), y(x, x)}

 

Okay, now I'm satisfied that dsolve rejects that input for the reason that I thought that it would.

 

Step 2: Make sure that that rejection happens quickly:

restart:

CodeTools:-Usage(
   proc()
      try dsolve(diff(y(x),x) = y(x,x)) catch: end try
   end proc
   ()
):

memory used=0.58MiB, alloc change=0 bytes, cpu time=16.00ms, real time=15.00ms, gc time=0ns

 

Okay, now I'm satisfied that it happened quickly.

 

Step 3: Redo it with a high setting of printlevel. This will give you a bunch of procedure names

that you can look at with showstat. Just ignore the procedures named unknown for now. The green

really useful information is at the bottom, starting with the first report of ERROR.

 

Note that the last green line before the pink error message tells you the line number in dsolve

itself (67 in this case) where the error was last invoked. You can look up the line numbers with

showstat.

restart:

#Note that the next 3-lines are a single execution group.
#The 40 in the 1st line is somewhat arbitrary. Be careful,
#the amount of information that you get increases literally
#exponentially as you increase printlevel. It quickly becomes
#overwhelming.
printlevel:= 40:
dsolve(diff(y(x),x) = y(x,x)):
printlevel:= 1: #Make sure that's a 1.

{--> enter dsolve, args = diff(y(x), x) = y(x, x)

{--> enter type/ODEtools/ODE, args = diff(y(x), x) = y(x, x)
<-- exit type/ODEtools/ODE (now in dsolve) = true}
{--> enter tools/atomize_names, args = {diff(y(x), x) = y(x, x)}
{--> enter tools/get_names, args = {diff(y(x), x) = y(x, x)}
<-- exit tools/get_names (now in tools/atomize_names) = {diff, x, y}}
<-- exit tools/atomize_names (now in dsolve) = {}}
{--> enter ModuleLoad, args =
{--> enter unknown, args = _c
{--> enter unprotect, args = _c
<-- exit unprotect (now in unknown) = }
{--> enter unassign, args = _c
<-- exit unassign (now in unknown) = }
{--> enter protect, args = _c
<-- exit protect (now in unknown) = }
<-- exit unknown (now in ModuleLoad) = }
{--> enter unknown, args = _s
{--> enter unprotect, args = _s
<-- exit unprotect (now in unknown) = }
{--> enter unassign, args = _s
<-- exit unassign (now in unknown) = }
{--> enter protect, args = _s
<-- exit protect (now in unknown) = }
<-- exit unknown (now in ModuleLoad) = }
{--> enter unknown, args = _p
{--> enter unprotect, args = _p
<-- exit unprotect (now in unknown) = }
{--> enter unassign, args = _p
<-- exit unassign (now in unknown) = }
{--> enter protect, args = _p
<-- exit protect (now in unknown) = }
<-- exit unknown (now in ModuleLoad) = }
{--> enter unknown, args = _xi
{--> enter unprotect, args = _xi
<-- exit unprotect (now in unknown) = }
{--> enter unassign, args = _xi
<-- exit unassign (now in unknown) = }
{--> enter protect, args = _xi
<-- exit protect (now in unknown) = }
<-- exit unknown (now in ModuleLoad) = }
{--> enter unknown, args = PDESolStruc
{--> enter unprotect, args = PDESolStruc
<-- exit unprotect (now in unknown) = }
{--> enter unassign, args = PDESolStruc
<-- exit unassign (now in unknown) = }
{--> enter protect, args = PDESolStruc
<-- exit protect (now in unknown) = }
<-- exit unknown (now in ModuleLoad) = }
{--> enter unknown, args = _u
{--> enter unprotect, args = _u
<-- exit unprotect (now in unknown) = }
{--> enter unassign, args = _u
<-- exit unassign (now in unknown) = }
{--> enter protect, args = _u
<-- exit protect (now in unknown) = }
<-- exit unknown (now in ModuleLoad) = }
{--> enter unknown, args = _l
{--> enter unprotect, args = _l
<-- exit unprotect (now in unknown) = }
{--> enter unassign, args = _l
<-- exit unassign (now in unknown) = }
{--> enter protect, args = _l
<-- exit protect (now in unknown) = }
<-- exit unknown (now in ModuleLoad) = }
{--> enter unknown, args = _alpha
{--> enter unprotect, args = _alpha
<-- exit unprotect (now in unknown) = }
{--> enter unassign, args = _alpha
<-- exit unassign (now in unknown) = }
{--> enter protect, args = _alpha
<-- exit protect (now in unknown) = }
<-- exit unknown (now in ModuleLoad) = }
{--> enter unknown, args = _beta
{--> enter unprotect, args = _beta
<-- exit unprotect (now in unknown) = }
{--> enter unassign, args = _beta
<-- exit unassign (now in unknown) = }
{--> enter protect, args = _beta
<-- exit protect (now in unknown) = }
<-- exit unknown (now in ModuleLoad) = }
{--> enter ODEtools/init, args =
{--> enter unknown, args = _xi
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _eta
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _chi
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _mu
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y1
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y2
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y3
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y4
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y5
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y6
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y7
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y8
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y9
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y10
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y11
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y12
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y13
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y14
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y15
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y16
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y17
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y18
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y19
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _y20
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = _s
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = symgen/y
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter unknown, args = symgen/x
<-- exit unknown (now in ODEtools/init) = 0}
{--> enter protect, args = _xi
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _eta
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _chi
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _mu
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y1
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y2
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y3
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y4
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y5
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y6
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y7
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y8
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y9
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y10
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y11
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y12
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y13
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y14
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y15
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y16
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y17
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y18
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y19
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _y20
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = _s
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = symgen/y
<-- exit protect (now in ODEtools/init) = }
{--> enter protect, args = symgen/x
<-- exit protect (now in ODEtools/init) = }
{--> enter unknown, args =
<-- exit unknown (now in ODEtools/init) = y}
{--> enter unknown, args =
<-- exit unknown (now in ODEtools/init) = x}
{--> enter unknown, args =
<-- exit unknown (now in ODEtools/init) = Y}
{--> enter unknown, args =
<-- exit unknown (now in ODEtools/init) = X}
{--> enter protect, args = y, x, Y, X
<-- exit protect (now in ODEtools/init) = }
<-- exit ODEtools/init (now in ModuleLoad) = }
<-- exit ModuleLoad (now in dsolve) = }
{--> enter PDEtools:-Library:-HasAnticommutativeVariables, args = diff(y(x), x) = y(x, x)
{--> enter UsePhysics, args = diff(y(x), x) = y(x, x)
<-- exit UsePhysics (now in PDEtools:-Library:-HasAnticommutativeVariables) = false}
<-- exit PDEtools:-Library:-HasAnticommutativeVariables (now in dsolve) = false}
{--> enter convert/global, args = []
{--> enter ProcessOptions, args = 1, [[]], convert/global/Defaults, Options
{--> enter unknown, args = []
{--> enter type/local, args = []
<-- exit type/local (now in unknown) = false}
<-- exit unknown (now in ProcessOptions) = []}
<-- exit ProcessOptions (now in convert/global) = }
<-- exit convert/global (now in dsolve) = []}
{--> enter ODEtools/odsolve, args = diff(y(x), x) = y(x, x), atomizenames = false
{--> enter ODEtools/integrator, args =
<-- exit ODEtools/integrator (now in ODEtools/odsolve) = PDEtools/int}
{--> enter ODEtools/get_yn, args = diff(y(x), x) = y(x, x)
{--> enter type/suffixed, args = x, _y, nonnegint
<-- exit type/suffixed (now in ODEtools/get_yn) = false}
{--> enter ODEtools/get_yn, args = {diff, x, y, y(x)}
{--> enter type/suffixed, args = diff, _y, nonnegint
<-- exit type/suffixed (now in ODEtools/get_yn) = false}
value remembered (in ODEtools/get_yn): type/suffixed(x, _y, nonnegint) -> false
{--> enter type/suffixed, args = y, _y, nonnegint
<-- exit type/suffixed (now in ODEtools/get_yn) = false}
{--> enter ODEtools/get_yn, args = {x, y}
<-- exit ODEtools/get_yn (now in ODEtools/get_yn) = {}}
<-- exit ODEtools/get_yn (now in ODEtools/get_yn) = {}}
<-- exit ODEtools/get_yn (now in ODEtools/odsolve) = {}}
{--> enter ODEtools/info, args = [diff(y(x), x) = y(x, x)], ode, FAIL, y, x, diff_ord, hint, WAY, extra, int_scheme, methods, singsol
{--> enter PDEtools/usediff, args = diff(y(x), x) = y(x, x)
{--> enter type/PDEtools/D, args = diff(y(x), x)
<-- exit type/PDEtools/D (now in PDEtools/usediff) = false}
{--> enter type/PDEtools/D, args = y(x)
<-- exit type/PDEtools/D (now in PDEtools/usediff) = false}
{--> enter type/PDEtools/D, args = y(x, x)
<-- exit type/PDEtools/D (now in PDEtools/usediff) = false}
<-- exit PDEtools/usediff (now in ODEtools/info) = diff(y(x), x) = y(x, x)}
{--> enter type/ODEtools/ODE, args = diff(y(x), x) = y(x, x)
<-- exit type/ODEtools/ODE (now in ODEtools/info) = true}
{--> enter Physics:-ModuleLoad, args =
{--> enter unprotect, args = Physics/ModuleIsLoaded
<-- exit unprotect (now in Physics:-ModuleLoad) = }
{--> enter TypeTools:-AddType, args = name, Or(symbol, indexed)
{--> enter known, args = name
<-- exit known (now in TypeTools:-AddType) = false}
{--> enter forgetAll, args =
<-- exit forgetAll (now in TypeTools:-AddType) = }
<-- exit TypeTools:-AddType (now in Physics:-ModuleLoad) = }
<-- exit Physics:-ModuleLoad (now in ODEtools/info) = }
{--> enter type/ODEtools/F(x), args = diff(y(x), x)
<-- exit type/ODEtools/F(x) (now in ODEtools/info) = false}
{--> enter type/ODEtools/F(x), args = y(x)
{--> enter type/PDEtools/known_function, args = y(x)
<-- exit type/PDEtools/known_function (now in type/ODEtools/F(x)) = false}
<-- exit type/ODEtools/F(x) (now in ODEtools/info) = true}
value remembered (in ODEtools/info): type/PDEtools/D(y(x)) -> false
value remembered (in ODEtools/info): type/PDEtools/D(y(x, x)) -> false
{--> enter unknown, args = y(x), y
<-- exit unknown (now in ODEtools/info) = y(x)}
{--> enter unknown, args = y(x, x), y
<-- exit unknown (now in ODEtools/info) = y(x, x)}
value remembered (in ODEtools/info): type/PDEtools/D(y(x, x)) -> false
<-- ERROR in ODEtools/info (now in ODEtools/odsolve) = found the indeterminate function %1 with different arguments %2, y, {y(x), y(x, x)}}
<-- ERROR in ODEtools/odsolve (now in dsolve) = found the indeterminate function %1 with different arguments %2, y, {y(x), y(x, x)}}
<-- ERROR in dsolve (now at top level) = found the indeterminate function %1 with different arguments %2, y, {y(x), y(x, x)}}
 #(dsolve,67): error

Error, (in dsolve) found the indeterminate function y with different arguments {y(x), y(x, x)}

 

 locals defined as: n_ODEs = 1, ARGS = [], METHOD = METHOD, args2 = (diff(y(x), x) = y(x, x)), zz = [], n = n, funcs = funcs, funcs_in_args = funcs_in_args, ans = ans, type_of_conversion = type_of_conversion, use_physics = false, j = j, vars = vars

 

 


 

Download printlevel.mw

The ONLY documented purpose of the quit command is to terminate (or "crash") your entire command-line-interface Maple session. It's purpose is not to end any individual procedure. It's purpose is not to display information about resource usage.

The optional number (the 12 in your case) is returned to the operating system. It's purpose is so that you can send some type of message to the operating system. What the messages mean is entirely up to you; the actual number is meaningless to Maple. There are standard or conventional predefined meanings for 0 through 5, but whether you adhere to those is up to you. The 12 was simply a number chosen at random to make a help page example.

This command has no documented purpose in worksheet Maple, and I suspect that nothing good can come from using it there.

If you want to end a Maple procedure via a command in that procedure, use return or error.

I seriously doubt that there's been any change to this command since Maple 7.

A solution essentially the same as Fabio's is available directly from menu commands. Go to Tools -> Options -> General. Answer the first question, "How do you want Maple to handle the creation of new Math Engines?". I use the last answer: "Ask each time", and I have that set globally and never need to bother with going to the menu. But you may find it more convenient to just "Share" and save that as "Apply to session."

I like Acer's way; indeed, it's the way that I first thought of. But here's another way which you may like better, especially if you like to see your polynomials instead of seeing H.

`diff/`:= ()-> ``(diff(args)):
`eval/`:= (f,e)-> eval(expand(f), e):
H:= ``@orthopoly:-H:

#Usage examples:
P:= H(1,x)*H(2,x);
P1:= diff(P,x);
eval(P1, x= 7);

You can use expand to expand these polynomials in the normal way. Two-argument eval will also force exapnsion, but subs will not. So, I believe that that satisfies both of your requests, right?

I thought that I'd handle both the base-3 case and the base-n case. The base-3 case can be handled in a special way because there's never a carry needed for the addition. (The same is true for bases 2 and 4, of course.) The procedure ZPad below pads one of the two lists with the number of zeros needed to match the length of the other list.

B10:= (N::list(nonnegint), b::posint)-> add(N*~b^~[$0..nops(N)-1]):

#The general case:
Nim:= (a::nonnegint, b::nonnegint, B::posint)->
   B10(irem~(convert(add(B10~(convert~([a,b], 'base', B), 10)), 'base', 10), B), B)
:
#The base-3 case:
ZPad:= (a::list, b::list)-> 
   (n-> [[a[], 0$max(-n,0)], [b[], 0$max(n,0)]])(nops(a)-nops(b))
:
`&Nim3`:= (a::nonnegint, b::nonnegint)-> 
   B10(irem~(add(ZPad(convert~([a,b], 'base', 3)[])), 3), 3)
: 

Usage:

Nim(11, 21, 3);
or
11 &Nim3 21;

Does that do what you want?

I can make ZPad a little more efficient if that's desired.

First 175 176 177 178 179 180 181 Last Page 177 of 395