It wasn't clear to me what you were asking, but I figured that you want to divide each entry of the first Array A by the corresponding entry of the second Array B, and then raise those results to the corresponding entries of Array C.

A := Array(1..2,1..2,[[1,2],[3,4]]):
B := Array(1..2,1..2,[[-1,-2],[7,9]]):
C := Array(1..2,1..2,[[1,2],[-3,-4]]):

zip(`^`,zip(`/`,A,B),C);

(A/~B)^~C;

But maybe you meant something else, eg. mapping some separate actions to only similar subranges of the Arrays. Just let us know. Maybe with a simple 2x2 example of what you want.

That's a good question. But answering it, more generally, would depend on what you mean by "correct". It looks like you are not content with the principal branch of the natural logarithm being used for fractional powering of negative real numbers.

The 5th bullet point of the Description of ?arithop mentions that the principal branch of ln is used when computing x^y as exp(y*ln(x)) for complex x. But it is also used for negative x (and that help page should be more clear). For example

> eval( x^y, [x=-1.1,y=4/3] );
                 -0.5677540645 - 0.9833788836 I

> eval( exp(y*ln(x)), [x=-1.1,y=4/3] );
                 -0.5677540643 - 0.9833788837 I

Suppose that you had terms like (k-1)^(4/3) instead. Would you want the corrective code to figure out that the substitution should automatically be for k-1 instead of for just k? What if you had a mix of rational powers of both k and k-1?

There are several choices of approach, for your simple example. (I add a different constant to each, just so that they are each easily visible and don't overlay each other when plotted together.)

y := 3/8*2^(1/3)*k^(4/3):
w := eval(y, k = sqrt(k^2)) + 1:
z := eval(y, k = abs(k) ) + 2:
v := subsindets(y,`^`,t->RealDomain:-`^`(op(t))) + 3:

plot([y,w,z,v],k=-3..3);

Now let's make it slightly tougher, with a power of k-1 instead. Note that the plot might not be what you want unless you substitute for k-1. But what if you had to write code which wouldn't know what to substitute for?

y := (k-1)^(8/5):
w := eval(y, k = sqrt(k^2)) + 1: # how to know to sub for k-1 instead?
w2 := algsubs(k-1 = sqrt((k-1)^2), y) + 1.2:
v := subsindets(y,:-`^`,t->RealDomain:-`^`(op(t))) + 2:

plot([y,w,w2,v],k=-3..3);

Harder still is a mix of powers. But now no single substitution works if just applied to the whole expression.

y := k^(4/3) + (k-1)^(4/3):
w := eval(y, k = sqrt(k^2)) + 1: # not right?
w2 := algsubs(k-1 = sqrt((k-1)^2), y) + 1.2: # not right?
v := subsindets(y,:-`^`,t->RealDomain:-`^`(op(t))) + 2:

plot([y,w,w2,v],k=-3..3);

Notice that by using RealDomain:-`^` the corrective code looks the same each time above, and works, and doesn't depend on what's in the expression.

I did it by substituting with RealDomain:-`^` after the fact of constructing expression y, because you might encounter other issues if you loaded all/part of the RealDomain package up front. An alternative is to do something like with(RealDomain) before the unseen computation that produced expression y. But that can introduce other difficulties I won't get into. Hence I offer that post-computation substitution as a suggestion. Compare this with the example at top

> with(RealDomain):

> eval( x^y, [x=-1.1,y=4/3] );
                          1.135508125

One usual way to do this would be to apply prep2trans to `a`, and then to apply CodeGeneration[C] both to that result and to `g` separately,  to obtain C source for a pair of routines which could be compiled and linked together.

But you appear to have your heart set on a single, inlined routine. With that in mind, here is one way to get it. I'll throw in another procedure `b` similar to `a` just for fun. (This works for your simple example, but may have difficulties for more complicate examples.)

> a := x->piecewise(x<0,x^2,2*x);

                                        /        2     \
                     a := x -> piecewise\x < 0, x , 2 x/

> b := x->piecewise(x<0,x^2,2*x);

                                        /        2     \
                     b := x -> piecewise\x < 0, x , 2 x/

> g := proc(x) local ans1, ans2;
>        ans1 := a(x);
>        ans2 := b(x);
>        ans1+ans2;
>      end proc;:

> codegen[prep2trans](
>   codegen[prep2trans](
>     subs([a=subsop(3=op([inline,op(3,eval(a))]),eval(a)),
>           b=subsop(3=op([inline,op(3,eval(b))]),eval(b))],
>                  eval(g))));

               proc(x)
               local ans1, ans2;
                 ans1 := if x < 0 then x^2 else 2*x end if; ;
                 ans2 := if x < 0 then x^2 else 2*x end if; ;
                 ans1 + ans2;
               end proc;

> CodeGeneration[C](%);

void cg (int x)
{
  void ans1;
  void ans2;
  ans1 = (  if (x < 0)
    x * x;
  else
    2 * x;
);
  ans2 = (  if (x < 0)
    x * x;
  else
    2 * x;
);
  return(ans1 + ans2);

Do you just mean like this?

> WhatDay:=proc(Year,Month,Day)
>   StringTools:-ParseTime(
>      "%Y-%m-%d",
>      cat(Year,"-",Month,"-",Day) )[weekDayName];
> end proc:

> WhatDay(2010,9,8);
                          "Wednesday"

ps. Never use || unless you have to.

Maple makes a strong distinction between exact quantities and floating-point approximations.

The suggestion to wrap f(a)*f(m) with a call to evalf will help with that error message. But notice how it produces an exact rational set of results for `a` and `b` (which you may or may not want).

So an alternative is to start with `a` (and/or `b`) as floats, eg. 0.0 and 1.0 instead of 1 and 0 which are exact rationals. Of you could assign evalf((a+b)/2) or (a+b)/2.0 to `m` instead. For your simple `f` any float value for a, b, or m would be contagious and cause `f` to also return a float value.

Since nothing else in your original code introduces `m` (or evaluates later values for `a` or `b`) as a float approximation then if you start with exact values then you end up with exact values.

The comparator `<` under evalb (eg called by `if` or `while`) will not use sqrt(2) as a float. That's why you got that error message, because it could not tell whether the relation held without its doing some approximation or more sophisticated analysis. But under `is` it could tell. So yet another approach, while retaining purely use of exact quantities, could be

if is(f(a)*f(m) > 0)=true then ...

notice that I compared the result from `is` to `true`, since `is` can return stuff other than true/false.

piecewise(seq(op([1-q<=t and t<q, f[q](t+1-q)[i][j]]),q=1..3),0);

S:=seq(op([1-q<=t and t<q, f[q](t+1-q)[i][j]]),q=1..3);
piecewise(S,0);

This might have run out of memory.

One way to try and keep down the memory cost is to perform the factorization of the Matrix inplace as well as the linear solving step. Note that this overwrites the Matrix.

Something like this:

M := Matrix(.....): # or however you build it
B := Vector(.....): # or however built. The RHS of M.X=B

(p,lu) := LinearAlgebra:-LUDecomposition(M,output=NAG,inplace=true):

LinearAlgebra:-LinearSolve([p,lu], B, inplace):

The above does the linear solving inplace on B (overwriting Vector B with solution X). I think that you had that part of it already. But I believe that the above example also overwrites M with the superimposed L and U of the LU factorization, also doing that step inplace. So it should save memory allocation on the O(N^2) order of the size of high precision float M.

Lastly, if you can, also check this with an old version like Maple 10. Some recent versions had quite severe problems with high precision float LinearAlgebra, by somehow being unable to do garbage collection during the external calls. An O(N^3) operation can produce very many intermediate float values which, while not themselves stored in any input/output Matrix, altogether take a LOT of memory. If the system is not able to collect any of those intermediate floats produced during the arithmetic then at N=2000 they may be too much to handle. If this turns out to be your problem, then I don't know of an immediate fix for a recent release.

On the face of it, your example didn't look to me like it was going wrong because X[1,1] was a Matrix element. It looked more like a confusion between properties (assumptions) and types.

You could try your code with this.

test:=proc(x::satisfies(t->is(t,RealRange(0,1))))
 print(x);
end;

A simpler example:

> restart:
> assume(x,RealRange(0,1));

> is(x,RealRange(0,1));
                              true

> type(x,RealRange(0,1));
                             false

> type(x,satisfies(t->is(t,RealRange(0,1))));
                              true

Changing topic a little, I now try to pass the entire Matrix (which costs no more than passing a particular entry) with a similar type check. I find this puzzling behaviour:

> restart:
> type(Matrix([[2,3],[5,7]]),'Matrix(prime)'); # ok
                              true

> for i to 2 do: for j to 2 do:
>    assume(x[i,j],RealRange(0,1));
> od; od:

> X:=Matrix(2,2,(i,j)->x[i,j]):

> seq(type(Z,satisfies(t->is(t,RealRange(0,1)))), Z in X);
                     true, true, true, true

> type(X,'Matrix(satisfies(t->is(t,RealRange(0,1))))'); #??
                             false

Thanks for uploading your worksheet.

It looks to me like the problematic unevaluated fsolve call contains arguments with HFloats in them. (Ie, a special kind of scalar "hardware" double precision float, technically of Maple type `hfloat`, which is relatively new.) And it appears that is messing up fsolve.

Those HFloats printed just like regular software floats in the message you posted, and so copying and pasting that as subsequent input led to Maple's doing something quite different from what you yourself had done. Using actual HFloats, your example fails in Maple 14. But using regular floats (sfloats) it succeeds.

So one workaround is to convert any hfloats in the arguments to fsolve into regular floats. You procedure could be edited to do something like this:

BV := fsolve(op(subsindets([
        {cndtn1, cndtn2}, {bv1L = bv1Lappr .. Mi, bv1R = Mi .. bv1Rappr}
                ],hfloat,SFloat@op)));

Now, it's still possible that inside some loop in your code, a call to fsolve might fail to find a root for some other reason. One way to treat that is to put in a type-check, to see if the result of an fsolve call is returned as an unevaluated fsolve call (ie, it failed) instead of a set of equations of name=numeric, etc.

if type(BV, 'specfunc(anything,fsolve)' then
   #lprint(BV);
   #error "failed fsolve attempt %1", BV;
   #whatever you else you'd fall back to...
end if;

> V:=Vector([1,5,3,10,2,1,0,4]):

> r,c,v:=ArrayTools:-SearchArray(V):

> seq(`if`(v[i]<3,r[i],NULL),i=rtable_dims(v));
                            1, 5, 6

The command SearchArray is fascinating, as the number of objects in the expression sequence that it returns depends on the number of names in the expression sequence on the left-hand-side of the assignment statement.

I guess that you don't want any ugly single-left (name) quotes to appear in what gets displayed in the MathContainer. If that's the case then some usual solutions involving ``() or wrapped single-left quotes might not satisfy you.

One way to bring about the desired effect might be to convert some of the names in the expression to escaped local names that merely appear the same as the global names. So, for your example you could convert all the names in the numerator to locals. Here's example code, inserting both the "unsimplified" expression, and then following that up with a call to extract, "simplify", and re-insert it.

P:=module() option package; export `&/`;
 `&/`:=proc(x,y) local t;
     :-`/`(subs(seq(t=convert(t,`local`),
                    t in indets(x,name)),x),y);
  end proc:
end module:
with(P):

DocumentTools:-SetProperty('MathContainer0', 'expression', (a^2*b)&/a);

DocumentTools:-SetProperty('MathContainer0', 'expression',
   convert(DocumentTools:-GetProperty('MathContainer0', 'expression'),`global`));

I had package P export `&/` instead of `/`, so that the latter still works normally.

You could try this, although you might not like the letf-single (name) quotes that display (wrapping it).

`Δy`

Or, if you don't mind the space between the symbols you could just use multiplication.

Delta*y

Another way, using MathML and the 'value' option of SetProperty, could be this (use the whole thing, including quotation marks).

"<math><mi>&Delta;y</mi></math>"

So the command would then be

DocumentTools:-SetProperty('MathContainer0', 'value', "<math><mi>&Delta;y</mi></math>");

@John May He may be using a Document insteda of a Worksheet. Inside a Document, you can use the Insert -> Execution Group from the top menubar to get the red > command cursor. When using that the results are as described. It's probably by design that it "works" that way.

For example.

Download printlevel_proble.mw

We should note Erik's own comment that one of the benefits of having the `with` inside the worksheet is that it makes the worksheet portable amongst other users. It makes sense: if the worksheet needs the packages loaded in order to work as intended, then it is the worksheet which should somehow bring about the `with` loading.

Having said that, of all the approaches mentioned I like Joe's most, to use an option of the GUI's launcher. At some universities the IT support is somehow able to install customized application launchers on students desktops. It could be done for a class, too, even if by running an insall script just once at the start of the semester.

And repeated -c options might not be the only way, using the launcher. isn't there also a -i option, to force use of a particular file as the initialization file instead of the usual maple.ini? Or is that only on Unix?

Trying to somehow clear the entire environment, by removing all assignments, etc, seems hard to do exhaustively. Erik only mentioned inadvertant assignments to variables. But there are also remember tables and caches of system routines to consider. That includes both remember tables of evalf, limit, etc, but also mistakes like sin(x):=2.3 which don't show up in anames(user).

Use `depends` and not `has` for your test. You probably want to avoid `has` spuriously recognizing terms where the name is buried irrelevantly inside a structure. For example, when the name in question is the dummy variable of integration or summation in an Int or Sum inert (or int or sum unevaluated( function call.

You want to map over products and sums, and not into individual terms. So use a type-check like type(...,{`+`,`*`}) similarly to how Doug did above.