It appears that the global bindings are present in what gets constructed by the 2D Math input inserted by the Matrix palette. This may be a hint that the problem is in parsing, or something related.

In other words, the global name is getting used, even when the name is supposed to have been rebound to the module export. That's a bug, but not a bug in your module package.

The following is just an example of how it might go wrong. It may well be that some other intermediate form is used, or some other mechanism. I'm just trying to illustrate that the "localness" can get lost.

restart:

m:=module() option package; export f;
     f:=proc(x) x^2; end proc:
   end module:

with(m):

f(3);

                               9

str := sprintf("%m",'f(3)'):

sscanf(str,"%m")[];

                              f(3)

eval(%); # this contains the global name :-f and not m:-f

                              f(3)

% - :-f(3);

                               0

You might want to consider just assigning your procedures to the global names :-Cos and :-Sin, and leaving them out of a module for now.

acer

For some examples you may wish to ensure that the objective depends upon the names, and remove constant names.

Below, `t` is a dummy variable, and likely not wanted.

> restart:

> expr := Pi*x*y+Int(sin(cos(x+t)),t=0..1);

                                       1
                                      /
                                     |
                   expr := Pi x y +  |   sin(cos(x + t)) dt
                                     |
                                    /
                                      0

> indets(expr,name);

                                 {Pi, t, x, y}

> depnames := indets(expr,And(And(name,Non(constant)),          
>                             'satisfies'(u->depends(expr,u))));

                              depnames := {x, y}

> seq(u=-1..1, u in depnames);

                           x = -1 .. 1, y = -1 .. 1

acer

restart:
with(LinearAlgebra):

A:=Matrix([[-2,-1],[9/5,4/5]]);

                                     [-2  -1]
                                     [      ]
                                A := [ 9   4]
                                     [ -   -]
                                     [ 5   5]

Apownodd:=simplify(MatrixPower(A,2*k+1)) assuming k::posint:

map(limit,Apownodd,k=infinity);

                                  [-9  -5]
                                  [--  --]
                                  [4   4 ]
                                  [      ]
                                  [ 9   5]
                                  [ -   -]
                                  [ 4   4]

map(identify,evalf(A^61));

                                  [-9  -5]
                                  [--  --]
                                  [4   4 ]
                                  [      ]
                                  [ 9   5]
                                  [ -   -]
                                  [ 4   4]

Apowneven:=simplify(MatrixPower(A,2*k)) assuming k::posint:

map(limit,Apowneven,k=infinity);

                                  [ 9   5]
                                  [ -   -]
                                  [ 4   4]
                                  [      ]
                                  [-9  -5]
                                  [--  --]
                                  [4   4 ]

map(identify,evalf(A^60));

                                  [ 9   5]
                                  [ -   -]
                                  [ 4   4]
                                  [      ]
                                  [-9  -5]
                                  [--  --]
                                  [4   4 ]

You may also be able to get there, with insight, by diagonalizing.

acer

Look at that 

acer

If rsolve succeeds then you can try and use `limit`.

eqs := { a(n)*(a(n-1)+3) = a(n-1)+2, a(1) = 2 }:

R:=rsolve(eqs, a(n) );

   /  /            n                      n       
 - |  |/     (1/2)\   (1/2)   /     (1/2)\   (1/2)
   \2 \\2 - 3     /  3      - \2 + 3     /  3     

                    n                 n\\//            n       
        /     (1/2)\      /     (1/2)\ || |/     (1/2)\   (1/2)
    + 2 \2 - 3     /  + 2 \2 + 3     / // \\2 - 3     /  3     

                  n                      n               n\
      /     (1/2)\   (1/2)   /     (1/2)\    /     (1/2)\ |
    - \2 + 3     /  3      + \2 - 3     /  + \2 + 3     / /

limit(R,n=infinity);

                                 (1/2)
                          4 - 2 3     
                          ------------
                            (1/2)     
                           3      - 1 

expand(rationalize(%));

                            (1/2)    
                           3      - 1

evalf(%);

                          0.7320508060

If rsolve fails (and if there is no easy fix-up to a workable polynomial form) then you might still be able to use `evalf/Limit` on the result of passing the `makeproc` option to `rsolve`. If that in turn succeeds then you might even get lucky with `identify`.

eqs := { a(n) = (a(n-1)+2)/(a(n-1)+3), a(1) = 2 }:

h:=rsolve(eqs, a, makeproc ):

evalf(Limit('h'(x),x=infinity));

                          0.7320508076

identify(%);

                            (1/2)    
                           3      - 1

acer

Did you mean that you tried it on the 2D Math output from evaluating expand((x+y)^19);

I tried that with exponent 18, and then 19. Ie. I selected the whole output (without right-clicking) and waited to see if the Smart Popup menu appeared. It did, for exponent 18, but not for 19.

Then I issued stopat(SubexpressionMenu); and tried again. The debugger popped up for 18, but not for 19. This indication that Library code is not getting called on the larger example seems like a hint that there may be some kind of size cutoff in the GUI. Your problematic example doesn't look overly large to me.

In trying to discover whether there was a size cutoff I looked at the length of dotm (sprintf to "%m" format) of a call to Typesetting:-Typeset of the output expression, and tried a few other expressions, but that didn't seem to be the key metric. Maybe something else? This might be the wrong tack.

acer

Below I use i, i+2, and i^2 for your "point" P(1,2,3), so that it depends on the value of `i`.

Another possibility for referring to such a "point" might be to use a Vector, V, in which case its entries might be referenced instead as V[1], V[2], V[3].

Using expressions,

restart:
F:=6*x^2 + 3*x*y^2 - 1.5*z^3;
H:=VectorCalculus:-Hessian(F,[x,y,z]);
for i from 1 to 3 do
   eval(F,[x=i,y=i+2,z=i^2]);
   eval(H,[x=i,y=i+2,z=i^2]);
end do;

Using operators (procedures),

restart:
F:=(x,y,z)->6*x^2 + 3*x*y^2 - 1.5*z^3:
H:=unapply(VectorCalculus:-Hessian(F(x,y,z),[x,y,z]),[x,y,z]):
for i from 1 to 3 do
   F(i,i+2,i^2);
   H(i,i+2,i^2);
end do;

acer

restart:

T:=sqrt(2/3)*Matrix([[cos(omega*t),cos(omega*t-(2/3)*Pi),cos(omega*t+(2/3)*Pi)],
                     [-sin(omega*t),-sin(omega*t-(2/3)*Pi),-sin(omega*t+(2/3)*Pi)],
                     [1/sqrt(2),1/sqrt(2),1/sqrt(2)]]):

q:=T . map(diff,T^%T,t):

map(expand,map(convert,q,expln));

                             [  0    -omega  0]
                             [                ]
                             [omega    0     0]
                             [                ]
                             [  0      0     0]

acer

The problem is numeric (floating-point) evaluation error. Under the fast (hardware float) `evalhf` interpreter your `L` evaluates to Float(undefined) for values of `k` closer to 100.

Even under regular (software float) `evalf` there are numeric issues when done at lower working precision (ie. Digits).

There are several ways around this. The attached worksheet shows three related and relatively simple ways of getting around it by avoiding evalhf mode computation.

Download error1.mw

acer

How about having `f` return unevaluated if its arguments are not of type numeric.

restart:                                          

simp:=plot(sin(x), x=0..2*Pi):                    

f:=(xx,yy)->`if`(yy::numeric and xx::numeric,     
                 `if`(yy > 0, [xx,4*yy], [xx,yy]),
                 'procname'(xx,yy)):              

trns:=plottools:-transform(f):                                               

trns(simp);    

acer

Here are two approaches, the second of which gets that precise form.

(There are more ways to do this particular example, but from just this one stated example and no more description it is not possible to completely characterize some general task that you may or may not have in mind.)

Download denom.mw

acer

It's not clear what you want, exactly. But you may be able to adjust this to suit.

restart:

M:=Matrix(10, 10,
 {(1, 1) = 1, (1, 2) = 1, (1, 3) = 1, (1, 4) = 1, (1, 5) = 1, (1, 6) = 1,
 (1, 7) = 1, (1, 8) = 1, (1, 9) = 1, (1, 10) = 1, (2, 1) = 1, (2, 2) = 1,
 (2, 3) = 2, (2, 4) = 2, (2, 5) = 2, (2, 6) = 2, (2, 7) = 1, (2, 8) = 1,
 (2, 9) = 1, (2, 10) = 1, (3, 1) = 1, (3, 2) = 2, (3, 3) = 2, (3, 4) = 2,
 (3, 5) = 2, (3, 6) = 2, (3, 7) = 2, (3, 8) = 2, (3, 9) = 1, (3, 10) = 1,
 (4, 1) = 1, (4, 2) = 2, (4, 3) = 2, (4, 4) = 2, (4, 5) = 2, (4, 6) = 2,
 (4, 7) = 2, (4, 8) = 2, (4, 9) = 1, (4, 10) = 1, (5, 1) = 2, (5, 2) = 2,
 (5, 3) = 2, (5, 4) = 2, (5, 5) = 2, (5, 6) = 2, (5, 7) = 2, (5, 8) = 2,
 (5, 9) = 2, (5, 10) = 1, (6, 1) = 2, (6, 2) = 2, (6, 3) = 2, (6, 4) = 2,
 (6, 5) = 2, (6, 6) = 2, (6, 7) = 2, (6, 8) = 2, (6, 9) = 2, (6, 10) = 1,
 (7, 1) = 2, (7, 2) = 2, (7, 3) = 2, (7, 4) = 2, (7, 5) = 2, (7, 6) = 2,
 (7, 7) = 2, (7, 8) = 2, (7, 9) = 2, (7, 10) = 1, (8, 1) = 2, (8, 2) = 2,
 (8, 3) = 2, (8, 4) = 2, (8, 5) = 2, (8, 6) = 2, (8, 7) = 2, (8, 8) = 2,
 (8, 9) = 2, (8, 10) = 1, (9, 1) = 2, (9, 2) = 2, (9, 3) = 2, (9, 4) = 2,
 (9, 5) = 2, (9, 6) = 2, (9, 7) = 2, (9, 8) = 2, (9, 9) = 2, (9, 10) = 1,
 (10, 1) = 2, (10, 2) = 2, (10, 3) = 2, (10, 4) = 2, (10, 5) = 2, (10, 6) = 2,
 (10, 7) = 2, (10, 8) = 2, (10, 9) = 2, (10, 10) = 1}):

m,n:=LinearAlgebra:-Dimensions(M):

Plines:=plots:-display(
  seq(seq(`if`(M[m+1-j,i-1]<>M[m+1-j,i],
               plottools:-line([i-0.5,j-1+0.5],[i-0.5,j+0.5]),NULL),
          j=1..m),i=2..n),
  seq(seq(`if`(M[m+1-j,i]<>M[m+2-j,i],
               plottools:-line([i-0.5,j-1+0.5],[i-0.5+1,j-1+0.5]),NULL),
          j=2..m),i=1..n),
              view=[1..m,1..n]):

Ppoints:=plots:-display(
  plots:-sparsematrixplot(Matrix(m,n,(i,j)->M[m+1-i,j]-1),color=red),
  plots:-sparsematrixplot(Matrix(m,n,(i,j)->M[m+1-i,j]-2),color=green)):

plots:-display( Plines, Ppoints );

acer

In your example, `c` has been assigned the value that you see. The fact that the value came from evaluating a+b is lost (except mathematically, but you don't seem to be wanting that), by which I mean that in your example the literal 'a+b' is not a part of the value assigned to name `c`.

If you reverse the order of operations, then yes you can evaluate `c` to one level, and get the unevaluated 'a+b'.

> restart:

> c:=a+b;

                                  c := a + b

> a:=5*x^2-3*y^3:
 
> b:=3*x^2+10*y^3:

> c;

                                     2      3
                                  8 x  + 7 y

> eval(c,1);

                                     a + b

In my example above it is full evaluation of `c` that produces and displays that expression in x and y.

Another variation is,

> restart:

> a:=5*x^2-3*y^3:      
> b:=3*x^2+10*y^3:     

> c:='a+b';

                                  c := a + b

> c;

                                     2      3
                                  8 x  + 7 y

> eval(c,1);

                                     a + b

You could also study the evaluation rules within procedures, which differ from those for the "top level." The manuals should help here.

You might also consider the somewhat related topic of command history.

acer

One can programmatically access the body of a procedure using the `eval` command. Eg. to subs into it, or get its inert form, or assign a copy of it to something else, etc.

By default one can also use just the `print` command to display the body of a procedure. This differs from `eval` in that it only prints. The return value from this operation is NULL, so it doesn't give programmatic access. (See also ?verboseproc)

To display line numbers, one can use the `showstat` command.

> restart:

> u:=x->x^2;
                      
                            2
                      x -> x 

> f:=eval(u);

                            2
                      x -> x 

> f(3);

                         9

> print(u);

                            2
                      x -> x

> showstat(u);

u := proc(x)
   1   x^2
end proc

Note that all of these will display without showing any comments present in the original source.

acer

I don't really understand why you wouldn't post your full example in the original question. All these followups with additional examples just makes this thread hard to read (or organize).

This is as far as I got with your 3rd example. You haven't stated what assumptions you might be prepared to accept on the variables.

Trying to combine or simplify under simultaneous pair of assumptions that numerator and denominator of the final results radicals were both positive just seemed to consume time and memory, and I interrupted it.

restart:
ex:=-(1/6)*(4*e+a^2*g*2^(1/2))^(1/2)*2^(1/4)*(-64*EllipticF((g*a^2

+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2),

 I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))*e^4

+2*EllipticF((g*a^2+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2

-g^2*a^4)^(1/2), I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2

-g^2*a^4)^(1/2))*g^3*a^6*2^(1/2)*e+g^4*a^8*EllipticE((g*a^2

+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2),

 I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))

+8*EllipticF((g*a^2+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2

-g^2*a^4)^(1/2), I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2

-g^2*a^4)^(1/2))*e^2*g^2*a^4+2*g^3*a^6*2^(1/2)*e*EllipticE((g*a^2

+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2)

, I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))

-16*g*a^2*2^(1/2)*e^3*EllipticE((g*a^2+2*2^(1/2)*e)^(1/2)*(-g*a^2

+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2), I*(8*e^2+g^2*a^4

+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))-16*EllipticF((g*a^2

+2*2^(1/2)*e)^(1/2)*(-g*a^2+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2),

 I*(8*e^2+g^2*a^4+4*g*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))*g

*a^2*2^(1/2)*e^3-8*g^2*a^4*EllipticE((g*a^2+2*2^(1/2)*e)^(1/2)*(-g*a^2

+2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2), I*(8*e^2+g^2*a^4+4*g

*a^2*2^(1/2)*e)^(1/2)/(8*e^2-g^2*a^4)^(1/2))*e^2)/((g*a^2+2*2^(1/2)

*e)^2*g^(1/2)*(8*e^2-g^2*a^4)^(1/2)):

Q:=combine(ex,radical) assuming -g*a^2+2*2^(1/2)*e > 0:
simplify(Q,size) assuming 8*e^2-g^2*a^2 > 0:
simplify(%) assuming 8*e^2-g^2*a^2 > 0:
simplify(%,size);

                              /       /                         
                              |       |                         
              1               | (1/4) |/   2      (1/2)  \    2 
----------------------------- |2      |\g a  + 2 2      e/ g a  
                            2 |       |                         
   (1/2) /   2      (1/2)  \  |       |                         
6 g      \g a  + 2 2      e/  \       \                         

           /                                  (1/2)\     
           |  /   2    2  4        2  (1/2)  \     |     
           |I \8 e  + g  a  + 4 g a  2      e/     |     
  EllipticE|---------------------------------------| + 2 
           |                        (1/2)          |     
           |          /   2    2  4\               |     
           \          \8 e  - g  a /               /     

           /                                  (1/2)\       
           |  /   2    2  4        2  (1/2)  \     |       
           |I \8 e  + g  a  + 4 g a  2      e/     |   /   
  EllipticK|---------------------------------------| e \4 e
           |                        (1/2)          |       
           |          /   2    2  4\               |       
           \          \8 e  - g  a /               /       

                 \                                             \
                 |                    (1/2)               (1/2)|
      2    (1/2)\| /       2    (1/2)\      /   2    2  4\     |
   + a  g 2     /| \4 e + a  g 2     /      \8 e  - g  a /     |
                 |                                             |
                 |                                             |
                 /                                             /

acer