introMaple.mws

>    restart;

  A Short Introduction to Maple

>   

 Arithmetic Operations

 Variable Assignement

>    a := 100;

a := 100

>    b := 200;

b := 200

>    a+b;

300

>    a mod 3;

1

>    12*23;

276

>    12^23;

6624737266949237011120128

>    30!;

265252859812191058636308480000000

>    100/20;

5

>    100/30;

10/3

>    trunc(100/30);

3

>    round(100/30);

3

>    trunc(-4.3);

-4

>    round(-4.3);

-4

>    floor(2.3); floor(-2.3);

2

-3

>    ceil(2.3); ceil(-2.3);

3

-2

>    abs(-4);

4

Quotes, Concatenation operator

>    a:=3;

a := 3

>    a:='a';

a := 'a'

>    a;

a

>    sqrt(s^2);

(s^2)^(1/2)

>    assume(s>0);   sqrt(s^2);

s

>    assume(s<0);   sqrt(s^2);

-s

>    s:='s';

s := 's'

>    a:=(x+y)^3;

a := (x+y)^3

  In Maple 5, the concatenation operator was the dot.

 In Maple 6, the concatenation operator is the two vertical bars.

>    a||1;

a1

>    a||1||2;

a12

  Definition of strings

>    c:="abcd";

c :=

  Interface, Code Generation

  Customize the user interface

>    interface(version);

`Maple Worksheet Interface, Maple 8.00, IBM INTEL LINUX, Apr 22 2002 Build ID 110847`

>    interface(screenwidth=100);

>    interface( quiet=true);

>    interface( quiet=false);

>    interface( verboseproc = 2 );

>    eval(isprime);

proc (n) local btor, nr, p, r; options remember, system, `Copyright (c) 1993 Gaston Gonnet, Wissenschaftliches Rechnen, ETH Zurich. All rights reserved.`; if not type(n,'integer') then if type(n,'numer...

  LaTeX interface

>    expr:=expand((x+y)^10);

expr := x^10+10*x^9*y+45*x^8*y^2+120*x^7*y^3+210*x^6*y^4+252*x^5*y^5+210*x^4*y^6+120*x^3*y^7+45*x^2*y^8+10*x*y^9+y^10

>    latex(expr);

{x}^{10}+10\,{x}^{9}y+45\,{x}^{8}{y}^{2}+120\,{x}^{7}{y}^{3}+210\,{x}^{6}{y}^{4}+252\,{x}^{

5}{y}^{5}+210\,{x}^{4}{y}^{6}+120\,{x}^{3}{y}^{7}+45\,{x}^{2}{y}^{8}+10\,x{y}^{9}+{y}^{10}

  C and Fortran code generation

>    with(codegen);

Warning, the protected name MathML has been redefined and unprotected

[C, GRAD, GRADIENT, HESSIAN, JACOBIAN, MathML, WebEQ, cost, declare, dontreturn, eqn, fortran, horner, intrep2maple, joinprocs, makeglobal, makeparam, makeproc, makevoid, maple2intrep, optimize, packar...
[C, GRAD, GRADIENT, HESSIAN, JACOBIAN, MathML, WebEQ, cost, declare, dontreturn, eqn, fortran, horner, intrep2maple, joinprocs, makeglobal, makeparam, makeproc, makevoid, maple2intrep, optimize, packar...

>    C(1-x/2+3*x^2-x^3+x^4);

      t0 = 1.0-x/2.0+3.0*x*x-x*x*x+x*x*x*x;

>    fortran(1-x/2+3*x^2-x^3+x^4);

      t0 = 1-x/2+3*x**2-x**3+x**4

>    C([a=10.3*x-y,b=z/t*log(z)]);

      a = 0.103E2*x-y;

      b = z/t*log(z);

>    fortran([a=10.3*x-y,b=z/t*log(z)]);

      a = 0.103D2*x-y

      b = z/t*log(z)

>    cost(expand((x+y)^5));

5*additions+28*multiplications

 Types  (whattype, type, hastype)

  Maple has a type system. Using it is not obligatory though.

>    a:=1; whattype(a);

a := 1

integer

>    b:=1.2; whattype(b);

b := 1.2

float

>    c:=2-3*I; whattype(c);

c := 2-3*I

complex

>    c:=cos(x); whattype(c);

c := cos(x)

function

>    whattype(1/2);

fraction

>    whattype(x);

symbol

   For lazy typists, there is an alias mechanism in Maple.

>    alias(w=whattype);

w

>    w(a+b);

float

>    w(x+y); w(x*y); w(x^3);

`+`

`*`

`^`

>    w(w);

symbol

Besides whattype, there are two more commands that handle types

>    type(a,integer); type(b,integer);

true

false

>    hastype(a,integer); hastype(b,integer);

true

false

The hastype(expr,t) returns true if expr has a subexpression of the type t

>    pol:=(x+x*y)^2;

pol := (x+x*y)^2

>    w(pol);

`^`

>    type(pol,`*`);

false

>    hastype(pol,`*`);

true

>    type(pol,`+`);

false

>    hastype(pol,`+`);

true

>    pol:=expand((x+x*y)^2);

pol := x^2+2*x^2*y+x^2*y^2

>    w(pol);

`+`

>    type(pol,`*`);

false

>    hastype(pol,`^`);

true

 Access parts of a Maple expression (op, nops, indets)

   Tree structure of an expression in Maple.

 Internal representation of mathematical expressions.

>    a:=x+y;

a := x+y

>    nops(a);

2

>    op(1,a);

x

>    op(2,a);

y

>    b:=sqrt(a);

b := (x+y)^(1/2)

>    nops(b);

2

>    op(1,b);

x+y

>    op(2,b);

1/2

>    w(b);

w((x+y)^(1/2))

>    nops(op(2,b));

2

>    op(1,op(2,b)); op(2,op(2,b));

1

2

 Because fractions are quite common,

 there are two special commands to get

 the numerator and the denominator.

>    numer(1/2); denom(1/2);

1

2

   There is a more compact way to  write this

>    op([2,1],b); op([2,2],b);

1

2

>    op(0,b);

`^`

>    op(1..3,[5,6,7,8]);

5, 6, 7

>    op(b);

x+y, 1/2

>    indets(x+y^2);

{x, y}

>    a:=x+z^3:  indets(a);

{x, z}

>    indets(sqrt(x)); nops(sqrt(x));

{x, x^(1/2)}

2

>    indets(sqrt(x)+y^(1/3));

{x, y, x^(1/2), y^(1/3)}

>    indets(exp(x+y));  indets(log(x+y));

{x, y, exp(x+y)}

{x, y, ln(x+y)}

 Sequences, Lists and Sets

  A sequence is defined as a succession of elements separated by commas

>    1,2,3;

1, 2, 3

>    whattype(%);

exprseq

We can create fairly complicated sequences

>    seq(i,i=1..10);

1, 2, 3, 4, 5, 6, 7, 8, 9, 10

>    seq(i^2,i=1..10);

1, 4, 9, 16, 25, 36, 49, 64, 81, 100

>    seq(x^i-1,i=1..7);

x-1, x^2-1, x^3-1, x^4-1, x^5-1, x^6-1, x^7-1

  A list  is defined as a sequence enclosed in (square) brackets

>    l:=[1,2,3];

l := [1, 2, 3]

>    whattype(%);

list

  You can use seq to define new lists

>    l:=[seq(ifactor(i),i=1..6)];

l := [1, ``(2), ``(3), ``(2)^2, ``(5), ``(2)*``(3)]

   To access the i-th element of a list

>    l[4];

``(2)^2

  A set is defined as a sequence enclosed in (curly) braces

>    s:={1,2,3};

s := {1, 2, 3}

>    whattype(%);

set

  Sets cannot contain multiple elements as opposed to lists

>    l:=[1,1,2,3];

l := [1, 1, 2, 3]

>    s:={1,1,2,3};

s := {1, 2, 3}

 Apply a function to every element of a list (map)

>    l1:=[1,9,64];

l1 := [1, 9, 64]

>    map(sqrt,l1);

[1, 3, 8]

>    l2:=[-1,4,-6];

l2 := [-1, 4, -6]

>    map(abs,l2);

[1, 4, 6]

>    l3:=[10/3,2,3/7];

l3 := [10/3, 2, 3/7]

>    map(evalf,l3);

[3.333333333, 2., .4285714286]

>    map(evalf[15],l3);

[3.33333333333333, 2., .428571428571429]

>    f:=x->x mod 3;

f := proc (x) options operator, arrow; `mod`(x,3) end proc

>    map(f,[1,10,12]);

[1, 1, 0]

>    map(x->x^2+1,[1,10,12]);

[2, 101, 145]

>    map(x->x^3-1,[x,y+1,z+2]);

[x^3-1, (y+1)^3-1, (z+2)^3-1]

>    map(factor,%);

[(x-1)*(x^2+x+1), y*(y^2+3*y+3), (z+1)*(z^2+5*z+7)]

When the function has arguments,

  these appear at the end.

>    l:=[1+x^2+x^2*y+x-x^2*y^3,x^3-x^3*y+x^2];

l := [1+x^2+x^2*y+x-x^2*y^3, x^3-x^3*y+x^2]

>    map(collect,l,x);

[1+(1+y-y^3)*x^2+x, (1-y)*x^3+x^2]

 Operations on sets (union, intersect, member)

>    s1 := {1,2,3};   s2 := {1,4,5};

s1 := {1, 2, 3}

s2 := {1, 4, 5}

>    s1 union s2;

{1, 2, 3, 4, 5}

>    s1 intersect s2;

{1}

>    member(2,s1); member(2,s2);

true

false

>    s1 minus s2;

{2, 3}

>    s2 minus s1;

{4, 5}

>    s1 minus {1};

{2, 3}

>    s2 union {6};

{1, 4, 5, 6}

 Operations on lists

 Add an element to a list

>    m1:=[1,2,3];

m1 := [1, 2, 3]

>    m1:=[op(m1),4];

m1 := [1, 2, 3, 4]

>    m1:=[5,op(m1),6];

m1 := [5, 1, 2, 3, 4, 6]

>    m1:=[op(1..3,m1),100,op(4..nops(m1))];

m1 := [5, 1, 2, 100, 4, 6]

  This technique is not efficient when

   we are dealing with big lists.

 Delete an element from a list

>    m2:= subsop(2=NULL,m1);

m2 := [5, 2, 100, 4, 6]

  Note that the initial list has been left unaltered.

 m2 is a simply copy of m1, with the second element removed.

>    m1;

[5, 1, 2, 100, 4, 6]

 Modify an element of a list

>    m1[2]:=20;

m1[2] := 20

>    m1;

[5, 20, 2, 100, 4, 6]

>    subsop(2=25,m1);

[5, 25, 2, 100, 4, 6]

>    m1;

[5, 20, 2, 100, 4, 6]

 Sums and Products (sum, prod)

>    sum(i,i=1..10^6);

500000500000

>    add(i,i=1..10^6);

500000500000

   sum looks for a closed form and then substitutes the

 numerical values of the bounds of summations

  add performs a loop from the lower bound to the upper bound

>    sum(i^2,i=1..1000000000000);

333333333333833333333333500000000000

>    sum(i,i=1..n);

1/2*(n+1)^2-1/2*n-1/2

>    factor(%);

1/2*n*(n+1)

>    factor(sum(i^3,i=1..n));

1/4*n^2*(n+1)^2

>    factor(sum(i^3,i=a..b));

-1/4*(a-1-b)*(a+b)*(a^2-a+b^2+b)

>    product(i,i=1..10);   10!;

3628800

3628800

   Infinite sums and products are handled

>    sum('1/i^2','i=1..infinity');

1/6*Pi^2

>    sum('1/k!', 'k'=0..infinity);

exp(1)

>    product('1-z^2/n^2','n=1..infinity');

product:   "Cannot show that 1-z^2/n^2 has no zeros on [1,infinity]"

1/z/Pi*sin(Pi*z)

  Sometimes the results are not immediately readable

>    product('1/i^2','i=1..n');

1/(GAMMA(n+1)^2)

>    i:='i':k:='k':
sum(combinat[binomial](i+k-1,k-1),i=0..n);

(n+1)/k*binomial(n+k,k-1)

>    sum('(2*z)/(z^2-n^2)','n=1..infinity');

Psi(1-z)-Psi(z+1)

  use single quotes to avoid premature evaluation

>    i:=2;

i := 2

>    sum(i,i=1..1000);

Error, (in sum) summation variable previously assigned, second argument evaluates to 2 = 1 .. 1000

>    sum('i','i'=1..1000);

500500

in the following sum we access the command binomial  from the

package combinat in the long form, avoiding the with(combinat)

>    n:='n';
sum(combinat[binomial](n,k),k=0..n);

n := 'n'

2^n

  the sum is actually equal to  (1+1)^n (binomial theorem)

   we can have nested sums and products

>    factor(sum(sum('i*j','i=1..n'),'j=1..n'));

1/4*n^2*(n+1)^2

>    i:='i': j:='j': n:='2':

>    sum(sum(n!/(i!*j!*(n-i-j)!),i=0..n),j=0..n-i);

-288/(5+7^(1/2)*I)/(-5+7^(1/2)*I)

>    evalc(%);

9

the sum is actually equal to  (1+1+1)^2 = 9, by the trinomial theorem

  now let us try to verify  that (1+1+1)^3 = 27:

>    i:='i': j:='j': n:='3':

>    sum(sum(n!/(i!*j!*(n-i-j)!),i=0..n),j=0..n-i);

(2*(i-2)/i*(i-3)*(-4+i)*((((-4+i)!+1/6*(-1+i)!)*(i-2)!+(-4+i)!*(-1+i)!)*(i-3)!+1/2*(i-2)!*(-4+i)!*(-1+i)!)*(i-5/2)*(-1+i)*(i^2+1/2*i+8)/(2+i)/(i+1)/(4-i)!/(i-2)!/(i-3)!/(-4+i)!/(-1+i)!*hypergeom([2, -1...
(2*(i-2)/i*(i-3)*(-4+i)*((((-4+i)!+1/6*(-1+i)!)*(i-2)!+(-4+i)!*(-1+i)!)*(i-3)!+1/2*(i-2)!*(-4+i)!*(-1+i)!)*(i-5/2)*(-1+i)*(i^2+1/2*i+8)/(2+i)/(i+1)/(4-i)!/(i-2)!/(i-3)!/(-4+i)!/(-1+i)!*hypergeom([2, -1...
(2*(i-2)/i*(i-3)*(-4+i)*((((-4+i)!+1/6*(-1+i)!)*(i-2)!+(-4+i)!*(-1+i)!)*(i-3)!+1/2*(i-2)!*(-4+i)!*(-1+i)!)*(i-5/2)*(-1+i)*(i^2+1/2*i+8)/(2+i)/(i+1)/(4-i)!/(i-2)!/(i-3)!/(-4+i)!/(-1+i)!*hypergeom([2, -1...
(2*(i-2)/i*(i-3)*(-4+i)*((((-4+i)!+1/6*(-1+i)!)*(i-2)!+(-4+i)!*(-1+i)!)*(i-3)!+1/2*(i-2)!*(-4+i)!*(-1+i)!)*(i-5/2)*(-1+i)*(i^2+1/2*i+8)/(2+i)/(i+1)/(4-i)!/(i-2)!/(i-3)!/(-4+i)!/(-1+i)!*hypergeom([2, -1...
(2*(i-2)/i*(i-3)*(-4+i)*((((-4+i)!+1/6*(-1+i)!)*(i-2)!+(-4+i)!*(-1+i)!)*(i-3)!+1/2*(i-2)!*(-4+i)!*(-1+i)!)*(i-5/2)*(-1+i)*(i^2+1/2*i+8)/(2+i)/(i+1)/(4-i)!/(i-2)!/(i-3)!/(-4+i)!/(-1+i)!*hypergeom([2, -1...
(2*(i-2)/i*(i-3)*(-4+i)*((((-4+i)!+1/6*(-1+i)!)*(i-2)!+(-4+i)!*(-1+i)!)*(i-3)!+1/2*(i-2)!*(-4+i)!*(-1+i)!)*(i-5/2)*(-1+i)*(i^2+1/2*i+8)/(2+i)/(i+1)/(4-i)!/(i-2)!/(i-3)!/(-4+i)!/(-1+i)!*hypergeom([2, -1...
(2*(i-2)/i*(i-3)*(-4+i)*((((-4+i)!+1/6*(-1+i)!)*(i-2)!+(-4+i)!*(-1+i)!)*(i-3)!+1/2*(i-2)!*(-4+i)!*(-1+i)!)*(i-5/2)*(-1+i)*(i^2+1/2*i+8)/(2+i)/(i+1)/(4-i)!/(i-2)!/(i-3)!/(-4+i)!/(-1+i)!*hypergeom([2, -1...
(2*(i-2)/i*(i-3)*(-4+i)*((((-4+i)!+1/6*(-1+i)!)*(i-2)!+(-4+i)!*(-1+i)!)*(i-3)!+1/2*(i-2)!*(-4+i)!*(-1+i)!)*(i-5/2)*(-1+i)*(i^2+1/2*i+8)/(2+i)/(i+1)/(4-i)!/(i-2)!/(i-3)!/(-4+i)!/(-1+i)!*hypergeom([2, -1...

 (more interesting things happen for n=4)

 to actually do this sum we break it in two pieces

>    i:='i': j:='j': n:='3':

>    sum(n!/(i!*j!*(n-i-j)!),j=0..n-i);

48/i!/(3-i)!/(2^i)

>    sum(%,i=0..n);

27

 now we know how to do the general sum (1+1+1)^n = 3^n

>    i:='i': j:='j': n:='n':

>    sum(n!/(i!*j!*(n-i-j)!),j=0..n-i);

n!/i!/(n-i)!/(2^(-n+i))

>    simplify(sum(%,i=0..n));

3^n

 Conversions (convert)

>    convert([1,2,3],set);

{1, 2, 3}

>    convert({1,2,3,4},list);

[1, 2, 3, 4]

>    convert([1,1,2,3],set);

{1, 2, 3}

>    l:=seq(i^3,i=1..5);

l := 1, 8, 27, 64, 125

>    convert([l],`+`);

225

>    convert([l],`*`);

1728000

>    series(tan(x),x,9);

series(1*x+1/3*x^3+2/15*x^5+17/315*x^7+O(x^9),x,9)

>    w(%);

w(series(1*x+1/3*x^3+2/15*x^5+17/315*x^7+O(x^9),x,9))

>    convert(%%,polynom);

x+1/3*x^3+2/15*x^5+17/315*x^7

>    w(%);

w(x+1/3*x^3+2/15*x^5+17/315*x^7)

>    convert(16341,hex);

`3FD5`

>    convert(16341,binary);

11111111010101

>    convert(16341,octal);

37725

  For more general bases, use the following syntax

>    convert(13,base,6);

[1, 2]

>    convert(15,base,15);

[0, 1]

  Partial fraction decomposition of a rational function

>    r:=(3*x+1)/(x^2-4*x+4);    w(r);

r := (3*x+1)/(x^2-4*x+4)

w((3*x+1)/(x^2-4*x+4))

>    convert(r,parfrac,x);

7/(x-2)^2+3/(x-2)

 Iterative constructs (for, while)

>    for i from 1 to 4 do i^2; od;

1

4

9

16

>    l:=[];

l := []

>    for i from 1 to 4 do l:=[op(l),i^2]; od;

l := [1]

l := [1, 4]

l := [1, 4, 9]

l := [1, 4, 9, 16]

>    i:=1;

i := 1

>    while i+2 < 6 do i:=i+1: od;

i := 2

i := 3

i := 4

 Functions

>    f:=x->x+1;

f := proc (x) options operator, arrow; x+1 end proc

>    f(2);

3

>    type(f,'function');
type(cos(x),'function');

false

true

>    whattype(f);

symbol

    The operators @  and @@  are used to compose functions

>    (f@f)(2);

4

>    (f@@5)(2);

7

>    g:=y->y+2;

g := proc (y) options operator, arrow; y+2 end proc

>    (f@g)(x);

x+3

>    ((f@g)@@3)(x);

x+9

>    h:=(x,y,z)->x^2+y^2-z^2;

h := proc (x, y, z) options operator, arrow; x^2+y^2-z^2 end proc

>    h(3,4,5);

0

  Recursive functions (factorial, Fibonacci)

>    fact:=n->if n = 0 then 1; else n*f(n-1); fi;

fact := proc (n) options operator, arrow; if n = 0 then 1 else n*f(n-1) end if end proc

>    fact(3); fact(7);

9

49

>    fibo:=n->if n = 0 then 0
    elif n = 1 then 1
    else fibo(n-1)+fibo(n-2);
fi;

fibo := proc (n) options operator, arrow; if n = 0 then 0 elif n = 1 then 1 else fibo(n-1)+fibo(n-2) end if end proc

>    seq(fibo(i),i=0..12);

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144

>    phi := (1+sqrt(5))/2;  
Phi := (1-sqrt(5))/2;

phi := 1/2+1/2*5^(1/2)

Phi := 1/2-1/2*5^(1/2)

>    fibo1:=n->(phi^n-Phi^n)/(sqrt(5));

fibo1 := proc (n) options operator, arrow; (phi^n-Phi^n)/sqrt(5) end proc

>    fibo1(100);

1/5*((1/2+1/2*5^(1/2))^100-(1/2-1/2*5^(1/2))^100)*5^(1/2)

>    expand(%);

354224848179261915075

>    seq(fibo(i)-expand(fibo1(i)),i=1..10);

0, 0, 0, 0, 0, 0, 0, 0, 0, 0

>    time(fibo(25));

4.270

>    time(expand(fibo1(500)));

4.910

 Procedures

>    f:=proc(x) x^2; end;

f := proc (x) x^2 end proc

>    f(3);

9

>    ff:=proc(n,g)
   local aux,i;
      aux:=[seq(i^2,i=1..n)];
      map(g,aux);
end proc;

ff := proc (n, g) local aux, i; aux := [seq(i^2,i = 1 .. n)]; map(g,aux) end proc

>    ff(3,f);

[1, 16, 81]

  The option remember updates the remember table

>    fib := proc(n) option remember;
    if n<2 then n
           else fib(n-1)+fib(n-2) end if;
end proc;

fib := proc (n) option remember; if n < 2 then n else fib(n-1)+fib(n-2) end if end proc

>    st:=time():
fib(100);
time()-st;

354224848179261915075

0.

  Accessing the remember table

>    op(4,eval(fib));

TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...
TABLE([60 = 1548008755920, 0 = 0, 61 = 2504730781961, 1 = 1, 62 = 4052739537881, 2 = 1, 63 = 6557470319842, 3 = 2, 4 = 3, 64 = 10610209857723, 5 = 5, 65 = 17167680177565, 6 = 8, 66 = 27777890035288, 7 ...

>    op(4,eval(cos));

TABLE([0 = 1, 1/6*Pi = 1/2*3^(1/2), -x = cos(x), 1/2*Pi = 0, Pi = -1, -infinity = undefined, 1/4*Pi = 1/2*2^(1/2), 1/3*Pi = 1/2, infinity = undefined, x = cos(x), I = cosh(1)])

>    cos(Pi/12):=1/4*sqrt(6)+1/4*sqrt(2);

cos(1/12*Pi) := 1/4*6^(1/2)+1/4*2^(1/2)

>    op(4,eval(cos));

TABLE([0 = 1, 1/6*Pi = 1/2*3^(1/2), -x = cos(x), 1/2*Pi = 0, Pi = -1, -infinity = undefined, 1/12*Pi = 1/4*6^(1/2)+1/4*2^(1/2), 1/4*Pi = 1/2*2^(1/2), 1/3*Pi = 1/2, infinity = undefined, x = cos(x), I =...

>    cos(Pi/12);

1/4*6^(1/2)+1/4*2^(1/2)

 Matrices (linalg)

>    with(linalg);

Warning, the protected names norm and trace have been redefined and unprotected

[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, Wronskian, addcol, addrow, adj, adjoint, angle, augment, backsub, band, basis, bezout, blockmatrix, charmat, charpoly, cholesky, col, coldi...
[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, Wronskian, addcol, addrow, adj, adjoint, angle, augment, backsub, band, basis, bezout, blockmatrix, charmat, charpoly, cholesky, col, coldi...
[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, Wronskian, addcol, addrow, adj, adjoint, angle, augment, backsub, band, basis, bezout, blockmatrix, charmat, charpoly, cholesky, col, coldi...
[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, Wronskian, addcol, addrow, adj, adjoint, angle, augment, backsub, band, basis, bezout, blockmatrix, charmat, charpoly, cholesky, col, coldi...
[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, Wronskian, addcol, addrow, adj, adjoint, angle, augment, backsub, band, basis, bezout, blockmatrix, charmat, charpoly, cholesky, col, coldi...
[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, Wronskian, addcol, addrow, adj, adjoint, angle, augment, backsub, band, basis, bezout, blockmatrix, charmat, charpoly, cholesky, col, coldi...
[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, Wronskian, addcol, addrow, adj, adjoint, angle, augment, backsub, band, basis, bezout, blockmatrix, charmat, charpoly, cholesky, col, coldi...

>    m:=matrix([[1,2,3],[2,3,1],[3,2,1]]);

m := matrix([[1, 2, 3], [2, 3, 1], [3, 2, 1]])

>    det(m);

-12

>    im:=inverse(m);

im := matrix([[-1/12, -1/3, 7/12], [-1/12, 2/3, -5/12], [5/12, -1/3, 1/12]])

>    id3:=evalm(m&*im);

id3 := matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]])

>    eigenvalues(m);

1, -2, 6

>    map(evalf,[%]);

[1., -2., 6.]

>    charpoly(m,lambda);

lambda^3-5*lambda^2-8*lambda+12

>    eig:=[eigenvectors(m)];

eig := [[1, 1, {vector([1, -3/2, 1])}], [6, 1, {vector([1, 1, 1])}], [-2, 1, {vector([-13/11, 3/11, 1])}]]

>    eiv:=eig[2][1];

eiv := 6

>    eigm2:=op(eig[2][3]);

eigm2 := vector([1, 1, 1])

>    evalm((m-eiv*id3)&*eigm2);

vector([0, 0, 0])

>    nullspace(m);

{}

  Gaussian Elimination

>    m:=matrix([[1,2,3,-2],[2,3,1,4],[3,2,1,7]]);

m := matrix([[1, 2, 3, -2], [2, 3, 1, 4], [3, 2, 1, 7]])

>    m:=addrow(addrow(m,1,2,-2),1,3,-3);

m := matrix([[1, 2, 3, -2], [0, -1, -5, 8], [0, -4, -8, 13]])

>    m:=addrow(m,2,3,-4);

m := matrix([[1, 2, 3, -2], [0, -1, -5, 8], [0, 0, 12, -19]])

>    submatrix(m,1..2,1..3);

matrix([[1, 2, 3], [0, -1, -5]])

>    minor(m,1,2);

matrix([[0, -5, 8], [0, 12, -19]])

 Matrices (LinearAlgebra)

>    with(LinearAlgebra);

Warning, the assigned name GramSchmidt now has a global binding

[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...
[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...
[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...
[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...
[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...
[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...
[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...
[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...
[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...
[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...
[Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, ...

>    M:=Matrix(3,[[1,2,3],[-1,-3,-5],[-5,7,9]]);

M := Matrix(%id = 154901776)

>    sort(CharacteristicPolynomial(M,lambda));

lambda^3-7*lambda^2+31*lambda-10

>    Eigenvalues(M);

Vector(%id = 135892440)

>    f:= (i,j) -> x^i+y^j;

f := proc (i, j) options operator, arrow; x^i+y^j end proc

>    Matrix(3,f);

Matrix(%id = 136188444)

  Differentiation, Integration (diff, int)

>    restart;

>    diff(x^n,x);

x^n*n/x

>    diff(1/sqrt(x+1),x);

-1/(2*(x+1)^(3/2))

>    diff(x^x,x);

x^x*(ln(x)+1)

>    diff(exp(x*ln(x)),x);

(ln(x)+1)*exp(x*ln(x))

>    int(x^n,x);

x^(n+1)/(n+1)

>    int(sin(x)^5,x);

-1/5*sin(x)^4*cos(x)-4/15*sin(x)^2*cos(x)-8/15*cos(x)

>    int(exp(-x^2),x);

1/2*Pi^(1/2)*erf(x)

>    int(sqrt(1-x^2),x=-1..1);

1/2*Pi

>    int(1/(1+x^3),x=0..infinity);

2/9*Pi*3^(1/2)

>    int(1/(1+x^6),x=0..infinity);

1/3*Pi

>    int(1/(1+x^7),x=0..infinity);

1/7*Beta(1/7,6/7)

>    int(1/(1+x^n),x=0..infinity);

Definite integration: Can't determine if the integral is convergent.

Need to know the sign of --> n

Will now try indefinite integration and then take limits.

1/n*Pi*csc(Pi/n)

>    int(cos(x)/(1+x^2),x=0..infinity);

-1/2*Pi*sinh(1)+1/2*cosh(1)*Pi

>    convert(%,exp);

-1/2*Pi*(1/2*exp(1)-1/2*1/exp(1))+1/2*(1/2*exp(1)+1/2*1/exp(1))*Pi

>    simplify(%);

1/2*Pi*exp(-1)

>    int(log(x)/(1+x)^3,x=0..infinity);

-1/2

 Number Theory (numtheory)

>    with(numtheory);

Warning, the protected name order has been redefined and unprotected

[GIgcd, bigomega, cfrac, cfracpol, cyclotomic, divisors, factorEQ, factorset, fermat, imagunit, index, integral_basis, invcfrac, invphi, issqrfree, jacobi, kronecker, lambda, legendre, mcombine, mersen...
[GIgcd, bigomega, cfrac, cfracpol, cyclotomic, divisors, factorEQ, factorset, fermat, imagunit, index, integral_basis, invcfrac, invphi, issqrfree, jacobi, kronecker, lambda, legendre, mcombine, mersen...
[GIgcd, bigomega, cfrac, cfracpol, cyclotomic, divisors, factorEQ, factorset, fermat, imagunit, index, integral_basis, invcfrac, invphi, issqrfree, jacobi, kronecker, lambda, legendre, mcombine, mersen...

>    ifactor(2^(2^6)+1);

``(67280421310721)*``(274177)

>    divisors(111111);

{1, 3, 7, 11, 13, 21, 33, 37, 39, 77, 91, 111, 143, 231, 259, 273, 407, 429, 481, 777, 1001, 1221, 1443, 2849, 3003, 3367, 5291, 8547, 10101, 15873, 37037, 111111}

>    n:=2352356;

n := 2352356

>    phi(n);

1169640

>    pr:=convert(factorset(n),list);

pr := [2, 191, 3079]

>    n*product('(1-1/pr[i])','i'=1..nops(pr));

1169640

>    cyclotomic(3,x);

x^2+x+1

>    cyclotomic(6,x);

x^2-x+1

>    cyclotomic(7,x);

x^6+x^5+x^4+x^3+x^2+x+1

>    cyclotomic(8,x);

x^4+1

>    cyclotomic(100,x);

x^40-x^30+x^20-x^10+1

>    for n from 1 to 500 do
   if {coeffs(cyclotomic(n,x))} minus {1,-1} <> {}
      then lprint(n); fi;
od;

105

165

195

210

255

273

285

315

330

345

357

385

390

420

429

455

495

>    s:={};
st:=time():
for n from 1 to 2000 do
   s:=s union {coeffs(cyclotomic(n,x))}: od:
time()-st;
s;

s := {}

40.161

{-4, -3, -2, -1, 1, 2, 3, 4, 5}

  Count the number of primes less than a given limit

>    pi(100);

25

 resultants

>    p:=x^6+y^2-1;

p := x^6+y^2-1

>    q:=x^3+y^3-5;

q := x^3+y^3-5

>    factor(resultant(p,q,x));

(24+y^2-10*y^3+y^6)^3

>    solve({p,q},{x,y});

{y = RootOf(24+_Z^2-10*_Z^3+_Z^6), x = RootOf(_Z^3-5+RootOf(24+_Z^2-10*_Z^3+_Z^6)^3)}

>    with(plots):

Warning, the name changecoords has been redefined

>    p:=x^2+y^2-1;q:=x-y-1;

p := x^2+y^2-1

q := x-y-1

>    a:=implicitplot(p,x=-1..1,y=-1..1):

>    b:=implicitplot(q,x=-2..2,y=-2..1,color=green):

>    display({a,b});

[Maple Plot]

 other packages

>    with(orthopoly);

[G, H, L, P, T, U]

>    seq(T(n,x),n=1..6);

x, -1+2*x^2, 4*x^3-3*x, 1+8*x^4-8*x^2, 16*x^5-20*x^3+5*x, -1+32*x^6-48*x^4+18*x^2

>    map(factor,[%]);

[x, -1+2*x^2, x*(4*x^2-3), 1+8*x^4-8*x^2, x*(16*x^4-20*x^2+5), (-1+2*x^2)*(16*x^4-16*x^2+1)]

>    with(stats);with(stats[statplots]):

[anova, describe, fit, importdata, random, statevalf, statplots, transform]

>    mark:=rand(1..100);

mark := proc () local t; global _seed; _seed := irem(a*_seed,p); t := _seed;  to concats do _seed := irem(a*_seed,p); t := s*t+_seed end do; irem(t,divisor)+offset end proc
mark := proc () local t; global _seed; _seed := irem(a*_seed,p); t := _seed;  to concats do _seed := irem(a*_seed,p); t := s*t+_seed end do; irem(t,divisor)+offset end proc
mark := proc () local t; global _seed; _seed := irem(a*_seed,p); t := _seed;  to concats do _seed := irem(a*_seed,p); t := s*t+_seed end do; irem(t,divisor)+offset end proc
mark := proc () local t; global _seed; _seed := irem(a*_seed,p); t := _seed;  to concats do _seed := irem(a*_seed,p); t := s*t+_seed end do; irem(t,divisor)+offset end proc
mark := proc () local t; global _seed; _seed := irem(a*_seed,p); t := _seed;  to concats do _seed := irem(a*_seed,p); t := s*t+_seed end do; irem(t,divisor)+offset end proc

>    marks:=[seq(mark(),i=1..30)];

marks := [77, 39, 86, 69, 22, 10, 56, 64, 58, 61, 75, 86, 17, 62, 8, 50, 87, 99, 67, 10, 74, 82, 75, 67, 74, 43, 92, 94, 1, 12]

>    histogram(marks);

[Maple Plot]

>