MuPAD

There are a few basic number theory functions in the kernel; others are contained in the numlib library.

>> isprime(997);
                                   TRUE
>> Factor(2^67-1);
                          193707721 761838257287
>> nextprime(1000000);
                                  1000003
>> powermod(9382471,322973,1298377);
                                  880825
>> phi(nextprime(2^20)-1);
                                  498400

Here phi is Euler's totient function returning the number of integers less than and relatively prime to its argument. These functions allow us to perform simple RSA encryption and decryption. Suppose we choose two primes and compute their product:

>> p:=nextprime(5678);
                                   5683
>> q:=nextprime(6789);
                                   6791
>> N:=p*q;
                                 38593253

>> e:=17:

>> d:=1/e mod (p-1)*(q-1);
                                  6808373

>> M:=24367139;

>> M1:=powermod(M,e,N);
                                 18476508

>> powermod(M1,d,N);
                                 24367139

We have seen a glimpse of MuPAD's symbolic abilities in the equation solving above. But MuPAD can do much more than this: all manner of algebraic simplification; rewriting in a different form; partial fractions; and so on.

>> expand((x+2*y-3*z)^4);
 4       4       4         3      3            3       3            3
x  + 16 y  + 81 z  + 32 x y  + 8 x  y - 108 x z  - 12 x  z - 216 y z  -
       3              2          2         2           2  2       2  2
   96 y  z + 216 x y z  - 144 x y  z - 72 x  y z + 24 x  y  + 54 x  z  +
        2  2
   216 y  z
>> Factor(%);
                                            4
                             (x + 2 y - 3 z)
>> sum(1/(k*(k+2)*(k+4)),k=1..n);
                                  2        3       4
                     310 n + 337 n  + 110 n  + 11 n
                 ----------------------------------------
                                2        3       4
                 4800 n + 3360 n  + 960 n  + 96 n  + 2304
>> partfrac(%);
               1           1           1           1
           --------- - --------- - --------- + --------- + 11/96
           8 (n + 3)   8 (n + 1)   8 (n + 2)   8 (n + 4)
>> normal(%);
                                  2        3       4
                     310 n + 337 n  + 110 n  + 11 n
                 ----------------------------------------
                                2        3       4
                 4800 n + 3360 n  + 960 n  + 96 n  + 2304
>> Factor(%);
                                            2
                      n (n + 5) (55 n + 11 n  + 62)
                    ----------------------------------
                    96 (n + 1) (n + 2) (n + 3) (n + 4)
>> sum(sin(k*x),k=1..n);
(I exp(-I x) - I exp(I x) + I exp(-I n x) - I exp(I n x) -
   I exp(- I x - I n x) + I exp(I x + I n x)) /
   4 - 2 exp(-I x) - 2 exp(I x)
>> rewrite(%,sincos);
(2 sin(x) + 2 sin(n x) + I cos(x + n x) - 2 sin(x + n x) - I cos(- x - n x)
   ) / 4 - 4 cos(x)

MuPAD's calculus skills are excellent. It implements very powerful algorithms for differentiation, integration, and limits.

>> diff(x^3,x);
                                      2
                                   3 x
>> diff(exp(exp(x)),x$4);
                             2                       3
exp(x) exp(exp(x)) + 7 exp(x)  exp(exp(x)) + 6 exp(x)  exp(exp(x)) +
         4
   exp(x)  exp(exp(x))

The dollar operator, $, is MuPAD's sequencing operator. As with most operators, it is overloaded; in the context of a derivative it is interpreted as a multiple derivative. We can of course also perform partial differentiation.

>> int(sec(x),x);
                    ln(2 sin(x) + 2)   ln(2 - 2 sin(x))
                    ---------------- - ----------------
                           2                  2
>> int(cos(x)^3, x=-PI/4..PI/3);
                                 1/2      1/2
                              5 2      3 3
                              ------ + ------
                                12       8
>> int(E^(-x^2),x=0..0.5);
                           /    1                 \
                        int| --------, x = 0..0.5 |
                           |        2             |
                           |       x              |
                           \ exp(1)               /
>> float(%);
0.461281006412792448755702936740453103083759088964291146680472565934983884\
2952938567126622486999424745

>> hold(int(exp(-x^2),x=0..0.5));