### One simple way to do simple floating point arithmetics in shell

#!/bin/tcsh set n = 128 set beta = echo "(-0.6)*($n)*($n)" | bc -l set ppn = 16 set np = 31 set nodes = echo "if ($np%$ppn==0) {$np/$ppn;} else {$np/$ppn+1}" | bc

### Useful Linux command—paste

I have for long ignored blogging, but today I came across the useful linux command “paste” and I thought I should write it down. In the past I have used quite a few smart tricks to handle data or output files, but I failed to document these tricks. I really should do it this time.

So here is the story. I have an output file test.o1551739 which contains lines such as “iter = 125” scattering everywhere. I need to get the sum of all the numbers in these lines. How can I achieve this without writing a program? The following is a little thought process.

First, I am so used to grep that I can grab all such lines by doing:

grep iter test.o1551739

Then, use sed to get rid of the suffix “iter = “:

sed s/"iter = "//

Now that I have a bunch of numbers, one in a line. How can I sum them up? Yah, here comes the usage of paste. The command paste allows me to merge the lines and put a delimiter in between. I chose to use the plus sign, because then it gives me an arithmetic expression. Try this:

paste -sd+

Once I have the arithmetic expression, I’ll just use bc to calculate the result. So piping all these steps, here is the one line command:

grep iter test.o1551739 | sed s/"iter = "// | paste -sd+ | bc

Nice. I realize that there must be tons of ways to do the same job. But I really love my solution.

### The cases environment in latex

The cases environment defined by the amsmath package sets each case to be inline math style. This will look ugly when the formulas in the cases are complicated, such as fractions or integrals. Consider the following sample code:

\usepackage{amsmath} $f(x)= \begin{cases} \frac{1}{2(x-1)}, & x>1 \\ \frac{\Gamma(x)}{2^{x-1}}, & 0<x<1 \end{cases}$ 

The result looks like this:

The solution is to use the dcases environment provided by the mathtools package instead. The prefix d’ means display’. It will set the cases in displayed math style (exactly the same as \frac versus \dfrac). So if we write the following code

\usepackage{amsmath,mathtools} $f(x)= \begin{dcases} \frac{1}{2(x-1)}, & x>1 \\ \frac{\Gamma(x)}{2^{x-1}}, & 0<x<1 \end{dcases}$ 

then we will see the result which looks like this:

### Syracuse PSTricks Page

http://melusine.eu.org/syracuse/pstricks/

Many fantastic PSTricks examples.

### Bocher’s formula

A colleague of mine mentioned to me today the Bocher’s formula for computing the coefficients of the characteristic polynomial of a matrix. It seems that this formula does not appear too often in textbooks or literature. I’ll just write down the formula and the idea of a simple proof here.

Let the characteristic polynomial of a matrix $A$ be

$\displaystyle{p(\lambda)=\lambda^n+a_1\lambda^{n-1}+\cdots+a_n.}$

Then the coefficients can be computed by
$a_1=-tr(A),$
$a_2=-\frac{1}{2}\left(a_1tr(A)+tr(A^2)\right),$
$a_3=-\frac{1}{3}\left(a_2tr(A)+a_1tr(A^2)+tr(A^3)\right),$
$\vdots$
$a_n=-\frac{1}{n}\left(a_{n-1}tr(A)+\cdots+a_1tr(A^{n-1})+tr(A^n)\right).$

To prove the formula, note that the coefficient $a_j$ is the summation of all possible products of j eigenvalues, i.e.,

$\displaystyle{a_j=(-1)^j\sum_{\{t_1\cdots t_j\}\in C_n^j}\lambda_{t_1}\cdots\lambda_{t_j},}$

where $C_n^j$ denotes the j-combination of numbers from 1 to n, and the trace of $A^i$ is the sum of the $j$th power of the eigenvalues, i.e.,

$\displaystyle{tr(A^i)=\sum_{t=1}^n\lambda_t^i.}$

$\displaystyle{\left(\sum_{\{t_1\cdots t_j\}\in C_n^j}\lambda_{t_1}\cdots\lambda_{t_j}\right)\left(\sum_{t_{j+1}=1}^n\lambda_{t_{j+1}}^i\right)=\sum_{\{t_1\cdots t_{j+1}\}\in C_n^{j+1}}\lambda_{t_1}\cdots\lambda_{t_j}\lambda_{t_{j+1}}^i+\sum_{\{t_1\cdots t_j\}\in C_n^j}\lambda_{t_1}\cdots\lambda_{t_{j-1}}\lambda_{t_j}^{i+1}.}$

The above indicates that the first part of $(-1)^ja_jtr(A^i)$ cancels the second part of $(-1)^{j-1}a_{j-1}tr(A^{i+1})$, whereas the second part of $(-1)^ja_jtr(A^i)$ cancels the first part of $(-1)^{j+1}a_{j+1}tr(A^{i-1}).$ The rest of proof becomes obvious now.

### Mess up Matlab codes and outputs in LaTeX

You do not want to mess up, right? When writing a LaTeX document, you may once in a while want to include some Matlab codes and/or outputs (preferably typeset using typewriter font if you have the same taste as me) during the course of your writing. What I used to do was to copy and paste the Matlab codes into my LaTeX file, execute the codes in Matlab, then do another copy and paste to place the results in my LaTeX file, and finally decorate them in a verbatim block or something like that. Guess what, Matlab provides a command, called publish, that helps you do all these in a simpler way.

In a nutshell, the way to use publish is to first type in the texts (as your usually LaTeX editing), including Matlab codes, in a single .m file. Let’s say the file name is example.m. Then, in Matlab, you issue the command

publish('example.m', struct('format','latex','outputDir','ltx-src'));

It means that you want Matlab to process example.m and output a LaTeX file example.tex (that you can compile to get pdf) in the sub-directory ltx-src. This is it. Instead of writing example.tex, you write a file example.m.

So, how should I write example.m? It is best to give an example. See the following:

%% % <latex> % The eigenvalues of a circulant matrix can be % obtained by performing FFT on the first column % of the matrix. First, let us construct a % $5\times5$ circulant matrix \verb|C| whose first % column \verb|c| is generated with random input: % </latex> 

c = rand(5,1); % sad that Matlab does not provide a circulant() % command... C = toeplitz(c, c([1 end:-1:2]))

%% % <latex> % The eigenvalues of \verb|C| are nothing but % </latex>

lambda = fft(c)

%% % <latex> % Check it out! The output is the same as using % the \verb|eig| command: % </latex>

eig(C)

%% % Fun, isn't it? 

It is nothing but a script file that Matlab can execute, right? The trick part is that all the texts and LaTeX markups are buried in comment blocks. How the Matlab command publish makes a LaTeX output is that whenever it meets a whole block of comments starting with ‘%%‘, it strips the comment signs and decorates the whole block using the pair \begin{par} and \end{par}. On the other hand, whenever it meets a block of codes that does not start with ‘%%‘, Matlab knows that they are executable commands. Matlab uses \begin{verbatim} and \end{verbatim} to typeset these command texts, and automatically add the Matlab outputs of the commands, which are also decorated by the \begin{verbatim} and \end{verbatim} pair, in the LaTeX file. Something I am not satisfied is that Matlab does not recognize LaTeX commands such as \verb||. I have to put <latex></latex> so that Matlab can do a verbatim copy of \verb||, instead of expand the text \verb|| in some weird way, in the output LaTeX file.

It is time to try the above example yourself. Have fun.

### Call by value or call by reference?

In C, if you pass in a pointer to a function, the content of the involved variable may be changed. In Matlab, however, if you pass in an array as an argument, the content of the array will never be changed after the function call, even if the variable with the same name is modified within the function, unless the variable is also output. Ok, now comes the question: What if the function called by Matlab is written in C as a mex function? Haha, have fun with the following three tests, t1, t2, and t3.

/* File: t1.c */ #include 

void foo(double *a) { a[0] = 1; }

 void main(void) { double a[2] = {0}; foo(a); printf("first element of array a = %g\n", a[0]); } 

% File: t2.m function t2()

a = zeros(2,1); foo(a); fprintf(1, 'first element of array a = %g\n', a(1));

function foo(a)

 a(1) = 1; 

% File: t3.m function t3()

 a = zeros(2,1); foo_mex(a); fprintf(1, 'first element of array a = %g\n', a(1)); 

/* File: foo_mex.c */ #include "mex.h"

 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *a = mxGetPr(prhs[0]); a[0] = 1; } 

### A sample MATLAB code for defining functions with cases

The code should be self-explaning..

function y = f(x) y1 = x.* (x<0); y2 = x.^2 .* (x>=0) .* (x<2); y3 = 4 .* (x>=2); y = y1 + y2 + y3; 

• 270,266 hits