#
# Sieve of Eratosthenes from the Chapter 16 of the
# Maple Programming Guide where it is used as an example
# for using the Maple debugger - so Maple make no claim
# that it is partiularly "efficient"
#
  restart;
  sieve := proc(n::integer)
                local i, k, flags, count,twicei;
                count := 0;
                for i from 2 to n do
                    flags[i] := true
                end do;
                for i from 2 to n do
                    if flags[i] then
                       twicei := 2*i;
                       for k from twicei by i to n do
                           flags[k] := false;
                       end do;
                       count := count+1;
                    end if;
                end do;
                count;
           end proc:
   CodeTools:-Usage(sieve(1000000));

memory used=216.16MiB, alloc change=84.01MiB, cpu time=3.84s, real time=3.59s, gc time=670.80ms

#
# Sieve of Eratosthenes I wrote many years ago. It was
# originally developed to return all the primes up to n
# - not just how many primes there are.
#
# But still about 5X faster than the above
#
  restart;
  sieve := proc( Nval::integer)
                 local vec, k;
                 description "The Sieve of Eratosthenes";
                 uses ArrayTools, LinearAlgebra:
               #
               # Initialise a vector with all ones
               #
                 vec:= Vector( row, Nval, 1 ):
               #
               # set the first entry in the vector to zero
               #        
                 vec[1]:= 0:
               #
               # Set all even entries with index >=4, to zero
               # (eliminates all the even indices which can't
               # be prime)
               #
                 Fill(0, vec, 3, 2);
               #
               # For each odd number, say k,  eliminate all
               # of its odd multiples (other than 1*k)
               #
                 seq( Fill
                      ( 0, vec, k^2-1, k ),
                      k=3..floor( sqrt(Nval) ), 2
                    ):
               #
               # return the indices for all non-zero values
               #
                 return Transpose
                        ( SearchArray
                          ( vec )
                        );
           end proc:
  CodeTools:-Usage(numelems(sieve(1000000)));

memory used=89.29MiB, alloc change=47.64MiB, cpu time=842.00ms, real time=852.00ms, gc time=78.00ms

  restart;
  P:= [ 7.291666667*10^10*sin(2.154234962*10^(-11)*M*Unit('J')/(Unit('Pa')*Unit('mm')^3))*Unit('Pa')*Unit('m')*Unit('mm')^3/(M*Unit(J)),
       -7.291666667*10^10*(-1 + cos(2.154234962*10^(-11)*M*Unit('J')/(Unit('Pa')*Unit('mm')^3)))*Unit('Pa')*Unit('m')*Unit('mm')^3/(M*Unit('J')),
        M = 0..72
      ];
  plot([map(convert, P[1..2], unit_free)[], P[3]]);

 restart;
 with(plots):
 A := (x,r) -> 2*sqrt(-x^2 + r^2)*x:
 Explore
 ( display
   ( [ plot
       ( sqrt(r^2 - p^2),
         p=-r..r,
         color=blue,
         scaling=constrained
       ),
       plot
       ( [-x, -x, x, x, -x],
         [0 , sqrt(-x^2 + r^2), sqrt(-x^2 + r^2), 0, 0],
         title=typeset( "Opt x= %1, r=%2, Area = %2", x, r, A(x,r)),
         scaling=constrained
       )
     ]
   ),
   parameters=[   x= 0.0..4.0,
                  r=0.0..4.0
              ],
   initialvalues=[ x=1.0,
                   r=2.0                    
                 ]    
 );

  a:=2:
  convert(a, string);
  convert('a', string);

restart;
with(PDEtools):
with(DEtools):
declare(u(y, z));

declare(v(y, z));

declare(w(y, z));

declare(p(y, z));

declare(R(z));

Conti := diff(u(y, z), y) + u(y, z)/(y - R(z)) + diff(w(y, z), z) + diff(R(z), z)*diff(w(y, z), y) = 0;
NSU := v(y, z)^2/(y - R(z)) - diff(p(y, z), y) = 0;
NSV := ((u(y, z)*diff(v(y, z), y) + u(y, z)*v(y, z)/(y - R(z)) + w(y, z)*(diff(v(y, z), z) + diff(R(z), z)*diff(v(y, z), y))) + (-diff(v(y, z), y)/(y - R(z)) - diff(v(y, z), y, y) + v(y, z)/(y - R(z))^2)) + diff(R(z), z)*diff(w(y, z), y) = 0;
NSW := u(y, z)*diff(w(y, z), z) + w(y, z)*(diff(w(y, z), z) + diff(R(z), z)*diff(w(y, z), y)) + diff(p(y, z), z) + diff(R(z), z)*diff(p(y, z), y) - diff(w(y, z), y)/(y - R(z)) - diff(w(y, z), y, y) - diff(R(z), z)*diff(w(y, z), y) = 0;
sys := [Conti, NSU, NSV, NSW];
deteqs := Infinitesimals(sys);

  restart:
  Explore
  ( plot
    ( [0, x, x+w, x+w+x],
      [l*cos(theta), 0, 0, l*cos(theta)],
      title=typeset( "Area=", (w+x)*l*cos(theta)),
      scaling=constrained
    ),
    parameters=[ x=1..10,
                 w=1..10,
                 l=1..10,
                 theta=0..Pi/2
               ],
    initialvalues=[ x=5,
                    w=5,
                    l=5,
                    theta=Pi/4
                  ]
  );

  restart:
  with(plots):
  with(plottools):
  with(Student[Calculus1]):
  display( arc( [0,0], 2, 0..Pi), scaling=constrained, color=red);
  VolumeOfRevolution( sqrt(4-x^2), x=-2..2);

  restart;
  fig:= proc( alpha )
              uses geometry, plots:
              local R:=2, xv:
              _EnvHorizontalName:= x:
              _EnvVerticalName:= y:
            #
            # Origin and centre circle
            #
              point(C, R*cos(alpha), R*sin(alpha));
              circle(c, [ point(orig, 0,0), R ]  ):
              segment(s0, [orig, C]):
            #
            # Construction lines
            #
              segment(s1, [ point(ps1, R, -R), point(ps2, -5*R, 5*R) ]):
              segment(s2, [ point(ps3, -R, -R), point(ps4, 5*R, 5*R) ]):
            #
            # Connecting rods
            #
              xv:= rhs
                   ( solve
                     ( [ x<0,
                         sqrt((x-coordinates(C)[1])^2+(-x-coordinates(C)[2])^2)=4*R
                       ],
                       x
                     )[]
                   ):
              point(A, xv, -xv):
              xv:= rhs
                   ( solve
                     ( [ x>0,
                         sqrt((x-coordinates(C)[1])^2+(x-coordinates(C)[2])^2)=4*R
                       ],
                       x
                     )[]
                   ):
              point(B, xv, xv):
              segment(s3, [A, C]):
              segment(s4, [B, C]):
            #
            # Cylinder walls
            #
              segment(s5, [point(ps5, 2*R, R), point(ps6, 5*R, 4*R)] ):
              reflection(s6, s5, line(l1, y=x)):
              reflection(s7, s5, line(l2, x=0)):
              reflection(s8, s6, l2):
            #
            # Draw everything
            #
              display
              ( [ draw( [ A(color=orange, symbol=solidcircle, symbolsize=24),
                          B(color=orange, symbol=solidcircle, symbolsize=24),
                          C(color=orange, symbol=solidcircle, symbolsize=24),
                          orig(color=orange, symbol=solidcircle, symbolsize=24),
                          c(color=red, linestyle=dash),
                          s0(color=green, thickness=4),
                          s1(color=blue, linestyle=dash),
                          s2(color=blue, linestyle=dash),
                          s3(color=green, thickness=4),
                          s4(color=green, thickness=4),
                          s5(color=black, thickness=8),
                          s6(color=black, thickness=8),
                          s7(color=black, thickness=8),
                          s8(color=black, thickness=8)
                        ]
                      ),
                 textplot
                 ( [ coordinates(A)[], "A"],
                   align=[above, left],
                   font=[times, bold, 16]
                 ),
                 textplot
                 ( [ coordinates(B)[], "B"],
                   align=[above,right],
                   font=[times, bold, 16]
                 ),
                 textplot
                 ( [ coordinates(C)[], "C"],
                   align=above,
                   font=[times, bold, 16]
                 )
                ],
                axes=none
             ):
      end proc:
#
# Produce animation
#
  plots:-display
         ( [ seq
             ( fig(j),
               j=0..evalf(2*Pi), evalf(Pi/20)
             )
           ],
           insequence=true
         );

restart;
with(plots):
with(plottools):
AB := 39;
BC := 140;
BD := 140;
local D; #not local(D)!
Vdot := proc(U, V) local i; add(U[i]*V[i], i = 1 .. 2); end proc;
dist := proc(M, N) sqrt(Vdot(expand(M - N), expand(M - N))); end proc;

Fig := proc(alpha)
local cir, R, BC, BD, AC, AD, lAC, A, lBC, lAB, lBD, beta, B, C, Cc, Dd, D, Aa, Bb, F1, F2, d, k, h, i, Pb, Ph, pb1, ph1, Qb, Qh, qh1, qb1, p1, P1, p2, P2, p3, P3, p4, P4, q1, Q1, q2, Q2, q3, Q3, q4, Q4, cy1, cy2, cy3, cy4, tA, tB, tC, tD;
A := [0, 0]; R := 39; d := 83; BC := 140; BD := 140;
B := [R*cos(alpha), R*sin(alpha)];
k := BC/R; h := 1/2*sqrt(2);
Ph := [h*(R + BC), h*(R + BC)];
Pb := [h*(-R + BC), h*(-R + BC)];
Qh := [-h*(R + BC), h*(R + BC)];
Qb := [-h*(-R + BC), h*(-R + BC)];
P1 := [Ph[1] - 1/2*d*h, Ph[2] + 1/2*d*h];
P2 := [Ph[1] + 1/2*d*h, Ph[2] - 1/2*d*h];
P3 := [Pb[1] - 1/2*d*h, Pb[2] + 1/2*d*h];
P4 := [Pb[1] + 1/2*d*h, Pb[2] - 1/2*d*h];
Q1 := [Qh[1] + 1/2*d*h, Qh[2] + 1/2*d*h];
Q2 := [Qh[1] - 1/2*d*h, Qh[2] - 1/2*d*h];
Q3 := [Qb[1] + 1/2*d*h, Qb[2] + 1/2*d*h];
Q4 := [Qb[1] - 1/2*d*h, Qb[2] - 1/2*d*h];
cir := circle(A, R, color = black, linestyle = longdash);
F1 := plot(x, x = -R .. R + BC, color = black, linestyle = longdash);
F2 := plot(-x, x = -R - BC .. R, color = black, linestyle = longdash);
AC := R*(cos(alpha) + sqrt(k^2 - sin(alpha)^2));
C := [h . AC, h . AC];
AD := R*(cos(Pi - alpha) + sqrt(k^2 - sin(Pi - alpha)^2));
D := [-h*AD, h*AD]; lBC := plot([B, C], color = red, thickness = 4);
lAB := plot([A, B], color = red, thickness = 4); print(evalf(dist(B, C)), evalf(dist(B, D)));
lBD := plot([B, D], color = red, thickness = 4);
pb1 := pointplot(Pb, symbol = solidcircle, symbolsize = 5, color = black);
ph1 := pointplot(Ph, symbol = solidcircle, symbolsize = 5, color = black);
qb1 := pointplot(Qb, symbol = solidcircle, symbolsize = 5, color = black);
qh1 := pointplot(Qh, symbol = solidcircle, symbolsize = 5, color = black);
p1 := pointplot(P1, symbol = solidcircle, symbolsize = 10, color = black);
p2 := pointplot(P2, symbol = solidcircle, symbolsize = 10, color = black);
p3 := pointplot(P3, symbol = solidcircle, symbolsize = 10, color = black);
p4 := pointplot(P4, symbol = solidcircle, symbolsize = 10, color = black);
q1 := pointplot(Q1, symbol = solidcircle, symbolsize = 10, color = black);
q2 := pointplot(Q2, symbol = solidcircle, symbolsize = 10, color = black);
q3 := pointplot(Q3, symbol = solidcircle, symbolsize = 10, color = black);
q4 := pointplot(Q4, symbol = solidcircle, symbolsize = 10, color = black);
Aa := pointplot(A, symbol = solidcircle, symbolsize = 12, color = blue);
Bb := pointplot(B, symbol = solidcircle, symbolsize = 12, color = blue);
Cc := pointplot(C, symbol = solidcircle, symbolsize = 12, color = blue);
Dd := pointplot(D, symbol = solidcircle, symbolsize = 12, color = blue);
cy1 := plot([P1, P3], color = black, thickness = 8); cy2 := plot([P2, P4], color = black, thickness = 8);
cy3 := plot([Q1, Q3], color = black, thickness = 8); cy4 := plot([Q2, Q4], color = black, thickness = 8);
tA := textplot([0, 0, "A"], 'align' = {'above', 'right'});
tB := textplot([B[1], B[2], "B"], 'align' = {'above', 'right'});
tC := textplot([C[1], C[2], "C"], 'align' = {'above', 'right'});
tD := textplot([D[1], D[2], "D"], 'align' = {'above', 'right'});
display([cir, F1, F2, pb1, ph1, qb1, qh1, p1, p2, p3, p4, q1, q2, q3, q4, Aa, Bb, Cc, Dd, lAB, lBC, lBD, cy1, cy2, cy3, cy4, tA, tB, tC, tD], scaling = constrained);
end proc;

Fig(Pi/3);

Warning, A new binding for the name `D` has been created. The global instance of this name is still accessible using the :- prefix, :-`D`.  See ?protect for details.

display([seq(Fig((2*alpha*Pi)/50), alpha = 0 .. 50)], insequence = true, axes = none);

 restart;
 ans:=6 *a[0]+2 *k^2 *b[1]*x^2 *y+k^2 *b[1] *y *lambda+6* k^2 *b[2]* x^3 *y-c *b[1] *x* y-2* c* b[2] *x^2 *y-c *b[2] *y *lambda+(6 *b[2]^2 *x^4 *lambda)/(lambda^2 *sigma+mu^2)+(6 *b[2]^2 *x^2 *lambda^2)/(lambda^2 *sigma+mu^2)+(6 *b[1]^2 *lambda* x^2)/(lambda^2 *sigma+mu^2)-12 *a[0] *b[2] *x *y-12* (a[1] * x)* b[2] *x* y-12* (a[2] * (x^2)) *b[2] *x* y+(6 *b[1]^2 *lambda^2)/(lambda^2 *sigma+mu^2)+k^2 *(a[1] * (-2 *x* (mu* y-x^2-lambda)-mu *x* y))+2 *k^2* (a[2] * ((mu* y-x^2-lambda)^2+x* (-2 *x* (mu* y-x^2-lambda)-mu *x* y)))+c *(a[1] * (mu *y-x^2-lambda))+2 *c* (a[2] * (x *(mu *y-x^2-lambda)))+5* k^2 *b[2] *x *y *lambda-(c *b[2] *lambda^2 *mu)/(lambda^2 *sigma+mu^2)+(12 *b[1] *lambda *b[2] *x^3)/(lambda^2 *sigma+mu^2)+(12 *b[1] *lambda^2 *b[2] *x)/(lambda^2 *sigma+mu^2)-(12 *b[1]^2 *lambda *mu* y)/(lambda^2 *sigma+mu^2)+(k^2 *b[1] *lambda^2 *mu)/(lambda^2* sigma+mu^2)-12 *a[0] *b[1] *y-12* (a[1] * x) *b[1] *y-12* (a[2] * (x^2))* b[1] *y-6 *a[0]^2-12 *a[0] *(a[1] * x)-12 *a[0]* (a[2] * (x^2))-6* (a[1] * x)^2-12* (a[1] * x) *(a[2] * (x^2))-6* (a[2] * (x^2))^2-(2 *k^(2)* b[1] *lambda* mu^2 *y)/(lambda^2 *sigma+mu^2)+(k^2 *b[1] *lambda *mu *x^2)/(lambda^2 *sigma+mu^2)+(6* k^2 *b[2] *lambda *mu *x^3)/(lambda^2 *sigma+mu^2)+(6* k^2 *b[2] *lambda^2 *mu *x)/(lambda^2 *sigma+mu^2)+(2* c *b[2] *lambda *mu^2* y)/(lambda^2 *sigma+mu^2)-(c *b[2] *lambda* mu* x^2)/(lambda^2 *sigma+mu^2)-(12 *b[2]^2 *x^2 *lambda* mu* y)/(lambda^2 *sigma+mu^2)-(12* k^2 *b[2] *lambda* mu^2 *x *y)/(lambda^2 *sigma+mu^2)-(24 *b[1] *lambda *b[2] *x *mu* y)/(lambda^2 *sigma+mu^2)+6* (a[1] * x)+6* (a[2] * (x^2))+6 *b[2] *x *y+6 *b[1] *y;

P1 := coeff(coeff(ans, x, 4), y, 0);