inner_proc:=proc(M::integer,
                 A::Vector(datatype=float[8]),
                 B::Vector(datatype=float[8]))
local count::integer, i::integer;
count:=0;
for i from 1 to M do
  if sqrt(A[i]^2+B[i]^2) < 0.5 then
    count:=count+1;
  end if;
end do;
return 4.0*count/M;
end proc:

cinner_proc:=Compiler:-Compile(inner_proc):

p:=proc(N) # N number of points
local F,X,Y,n,inner_proc,P;
  n:=N;
  F:=Statistics:-RandomVariable(Uniform(-1/2, 1/2));
  X:=Statistics:-Sample(F,n);
  Y:=Statistics:-Sample(F,n);
  cinner_proc(n,X,Y);
end proc:

st:=time(): p(2*10^5); time()-st;
                      3.14218000000000020
                             0.063

I tried to resist the strong temptation to change how it works, and to be faithful to the given method.

acer

inner_proc:=proc(M::integer,
                 A::Vector(datatype=float[8]),
                 B::Vector(datatype=float[8]))
local count::integer, i::integer;
count:=0.0;
for i from 1 to M do
  if sqrt(A[i]^2+B[i]^2) < 0.5 then
    count:=count+1.0;
  end if;
end do;
return 4.0*count/M;
end proc:

cinner_proc:=Compiler:-Compile(inner_proc):

p:=proc(N) # N number of points
local F,X,Y,n,inner_proc,P;
  n:=N;
  F:=Statistics:-RandomVariable(Uniform(-1/2, 1/2));
  X:=Statistics:-Sample(F,n);
  Y:=Statistics:-Sample(F,n);
  cinner_proc(n,X,Y);
end proc:

st:=time(): p(2*10^5); time()-st;
                      3.14250000000000007
                             0.078

acer

What you are doing now is a bit like this (simplified to just 1 loop),

> restart:

> lc:=1:
> for ii from 1 to 3 do
>   T[ii]:=i->lc;
>   lc:=lc+1;
> end do:

> T[1](), T[2](), T[3]();
                                    4, 4, 4
 
> lc:=17:
> T[1](), T[2](), T[3]();
                                  17, 17, 17

But it seems like you want behaviour more like this,

> restart:

> lc:=1:
> for ii from 1 to 3 do
>   T[ii]:=unapply(lc,i);
>   lc:=lc+1;
> end do:

> T[1](), T[2](), T[3]();
                                    1, 2, 3
 
> lc:=17:
> T[1](), T[2](), T[3]();
                                    1, 2, 3

Is that closer to what you are after?

Notice also that there is no need to have the T indices be the same as the formal parameters of the operators (ie, i,j,k). Having it so may be just additionally confusing to you.

acer

Sorry, but 100*cos(20) is an exact quantity. (Recall, Maple is a computer algebra system, in which there is an important distinction since algebraic manipulation of exact quanitites is useful.) You can apply evalf or evalhf to it, to get a floating-point quantity.

Or, you could ensure that the argument to sin is itself a float, so that it returns a float. This is probably a better solution for you than explicit using explicit evalf/evalhf calls, if you intend to use Compile on the procedure.

Maple will not automatically cast an exact numeric quantity to a float when assigning into a float[8] Matrix, except in the cases that the exact quantity is an integer or a rational.

acer

> with(LinearAlgebra):
>
> A:=Matrix([[A11,A12],[A21,A22]]):
>
> InvA := Matrix([[(A11-A12&*A22^(-1)&*A21)^(-1),
>           -A11^(-1)&*A12&*(A22-A21&*A11^(-1)&*A12)^(-1)],
>         [-A22^(-1)&*A21&*(A11-A12&*A22^(-1)&*A21)^(-1),
>          (A22-A21&*A11^(-1)&*A12)^(-1)]]);
InvA :=
 
    [             1
    [--------------------------- ,
    [      //        1 \       \
    [A11 - ||A12 &* ---| &* A21|
    [      \\       A22/       /
 
     // 1        \                 1             \]
    -||--- &* A12| &* ---------------------------|]
     |\A11       /          //        1 \       \|]
     |                A22 - ||A21 &* ---| &* A12||]
     \                      \\       A11/       //]
 
    [ // 1        \                 1             \
    [-||--- &* A21| &* ---------------------------| ,
    [ |\A22       /          //        1 \       \|
    [ |                A11 - ||A12 &* ---| &* A21||
    [ \                      \\       A22/       //
 
                 1             ]
    ---------------------------]
          //        1 \       \]
    A22 - ||A21 &* ---| &* A12|]
          \\       A11/       /]
 
>
> map(normal, A . subs(`&*`=`*`,InvA) );
                                   [1    0]
                                   [      ]
                                   [0    1]
 
>
> data:={A11=RandomMatrix(2,2), A12=RandomMatrix(2,2),
>        A21=RandomMatrix(2,2), A22=RandomMatrix(2,2)}:
>
> A:=Matrix(eval([[A11,A12],[A21,A22]],data)):
>
> InvA:=Matrix(eval(subs(`&*`=`.`,
>           [[(A11-A12&*A22^(-1)&*A21)^(-1),
>             -A11^(-1)&*A12&*(A22-A21&*A11^(-1)&*A12)^(-1)],
>            [-A22^(-1)&*A21&*(A11-A12&*A22^(-1)&*A21)^(-1),
>             (A22-A21&*A11^(-1)&*A12)^(-1)]] ),data) ):
>
> A . InvA;
                              [1    0    0    0]
                              [                ]
                              [0    1    0    0]
                              [                ]
                              [0    0    1    0]
                              [                ]
                              [0    0    0    1]

acer

@Hamidreza Is there a timing difference if you compute

add(add(evalf(f[i, j](i, j, x, y)), i = 1 .. 50), j = 1 .. 10);

instead of

evalf(add(add(f[i, j](i, j, x, y), i = 1 .. 50), j = 1 .. 10));

Also, could you use evalhf instead of evalf, in the former of those two?

And what is D1()? Is it always the same, or does it change for each i and j (using globals)? If unchanging, then can you pull it out of the summation as a multiplicative constant, or evalf(D()) up front? That too could make a difference.

And, taking another tack, how about first finding Q:=sum(...,i=1..N) to get a formula, finding Q50:=eval(Q,N=50), and then computing add(evalf(Q50),j=1..100) or similar?

If you raise Digits, and then recompute, you can get some evidence that those relatively small imaginary components are merely artecfacts of doing the computation with floating-point numbers (ie. as an approximation).

For this example, you might be able to convince yourself of this with an assignment of only Digits:=16.

acer

Your procedure `inertia` expects to find the global xyzc to be an indexable object with three scalar values. That's how  `inertia` uses xyzc. And it works fine, the first time you call it.

But then you assign to xyzc (at the top level) something which is a list of three equations. This occurs when you compute Equate(cog(...),[xc,yc,zc]), the result of which you assign to global xyzc. So then xyzc is no longer in a form for which some of the commands internal to `inertia` are valid or permissible. The three entries of xyzc are now equations rather than scalar values. So the next invocation to `inertia` fails with an error, going wrong when it tries to square the entries of xyzc.

For example,

> (4.99=xc)^2;
Error, (in simpl/reloprod) invalid terms in product

acer

Must a local optimization (least squares) method find the global optimum for any nonlinear fit?

You could try a few things, such as setting it up yourself as a least squares problem, and then using different starting points or using a global optimizer.

I'd like to see Maple automatically detect the presence of GlobalOptimization, and use it within NonLinearFit if found. (Using yet another global optimizer might be possible but could get tricky, as it would have to accept the same syntax as Optimization and also somehow be registered for detection.)

acer

Thanks Alec! That's better for Primes, but worse for me, as the inability-to-post problem looks to be at my end.

A related question: is there anyone for whom the 2DMath displayed in Alec's response looks good? For me, it doesn't even look as good as the 2DMath in the Online Help using the exact same browser and session. It's small, and more jagged than even lack of antialising usually gives.

acer

Give it the "bug" tag, Alec.

It's an outright bug that assumptions are put on x merely by calling `series`, present in 13.01 (and probably 13.02 too). Which makes it present in the version most people are using at this date.

"Du musst Amboss oder Hammer sein." - Goethe

acer

Does LinearAlgebra:-SylvesterSolve do this?

You started with this,

> Z.C.X + C = D;

                              (Z . C . X) + C = D

So, now if you right-multiply by the inverse of X, you get,

> Z.C.X.X^(-1) + C.X^(`-1`) = D.X^(`-1`);

                                        -1           -1

                       (Z . C) + (C . (X  )) = D . (X  )

Now, if you read the help-page for LinearAlgebra:-SylvesterSolve, it says that it solves Matrix equations of the form A . X + isgn * X . B = scale*C.

And you seem to have that form, except that your C is their X, your Z is their A, your X^(-1) is their B, and your D.X^(-1) is their C. If that's correct, then how about calling,

LinearAlgebra:-SylvesterSolve(Z, X^(-1), D.X^(-1) );

You'd want to check that I didn't make a mistake in all that. If it doesn't work, maybe you could post or upload your actual Matrices.

acer

Who thinks that the original student's numeric task cannot be implemented in Maple within a factor of four or better of the Excel timing? How about a factor of 2?

Nick, please upload the Maple program if you are so permitted. If it cannot be achieved, then at least people will know more about what needs improving.

acer

Good question. When you pass an rtable (ie. any of Matrix, Array, or Vector) as an argument to a procedure call it gets passed similarly to what you called "by reference".

The procedure gets the same rtable object, not a copy, of what was passed. Hence changes made to the entries of that rtable parameter, while inside that procedure, also change the entries in the object at the higher level.

acer

@Nick lists are not ordered, in the sense that Maple retains their entries in the original order in which you input them or create the object. This is in contrast with sets.

Also, the syntax for the Vector constructor is with round-brackets or angle-brackets, like so,

Vector([1,2,3]);

<1,2,3>;

You have several choices about how you create and use your Array of Vectors. If you know the number of Vectors in advance, then you could simply use M:=Matrix(3,120) say, and code it so that references to the jth Vector was like M[1..3,j] or M[i,j] etc. But if you do that, then any time you had to pull out a column (ie, a 3-Vector) you might be creating a new object which has some performance issues if done many times. Of course, you might be able write your code to always reference the M[i,j] entries, and so avoid "pulling out" the Vectors; it depends on whether you have to pass them separately to another routine. (I don't want to barrage you with too much info, but there is also a round-bracket mechanism to help deal with grabbing subselections of Matrices.)

You could also do exactly what you wrote. Ie, first create A:=Array(1..120) then then assign to A[i] any Vector(3) that you created. This might avoid repeated or unnecessary creation/collection of objects. In this scenario, you would access the ith component of the jth Vector as A[j][i].

You might also want to look at the Portal page on such objects, for clarification of their differences. Note especially that Maple lists are not really mutable structures as are 1D Arrays

Finally, if you do use a 1D structure that contains Vector objects instead of a 2D Matrix then you have a choice between plain ol' Vector and the version of Vector available under the VectorCalculus package. If you load that package using `with` then you automatically get its version when you call Vector(). If you use the former you can easily do Curl and all that stuff with quite simple constructions of commands. If you use the latter then you get such vector calculus stuff for free as commands in the VectorCalculus package. But note that not everyone finds the coordinate system stuff that is built into VectorCalculus Vectors to fit their expectations.

Hope that helps some.

acer