What’s the point of spending hundreds of billions of pounds to renew the Trident nuclear deterrent to protect against hostile nuclear states, if you’re going to accept billions of pounds from China to build an untested nuclear reactor in Essex?

# Month: September 2016

## Counting semiprimes

The prime-counting function returns the number of prime numbers less than or equal to . The prime number theorem gives the asymptotic behaviour .

We shall call a number *semiprime* if it is the product of two primes. Let denote the semiprime-counting function. What is the asymptotic behaviour of ? I claim that

**Edit (17 September 2016) **Based on experiments in Mathematica, the above seems to be a lower bound, while an upper bound is given by

I obtain these by the following heuristic, probabilistic argument. The probability that a random integer chosen between 1 and is ; loosely, this means that the probability of a randomly chosen large number is approximately . Similarly, is the probability that a randomly-chosen is semiprime.

An exact formula for is

The sum over is trivially . Hence

We split the sum at into two parts, and assume that the first part is dominated by the second. In the second part, is large and we invoke and approximate the sum as an integral, to get

Evaluating the integral and keeping only the highest-order terms, we get the desired formula.

The semiprime-counting function can be implemented in Mathematica by

SemiprimeQ[n_] := SemiprimeQ[n] = Module[{f}, f = FactorInteger[n]; ( ((Length[f] == 2) && (Total[Last /@ f] == 2)) || ((Length[f] == 1) && f[[1]][[2]] == 2) ) ]; w[n_] := w[n] = Length[Select[Range[n], SemiprimeQ]]

(where we memoise the results to speed things up). The following code produces a graph showing the validity of these approximation:

W1[n_] := n/Log[n] Log[Log[n]]; W2[n_] := n/Log[n] Log[Log[n]/Log[2]]; nn = 150000; nj = 500; Show[ ListPlot[Table[{n, w[n]/W1[n]}, {n, 1, nn, nj}], PlotStyle -> Red], ListPlot[Table[{n, w[n]/W2[n]}, {n, 1, nn, nj}], PlotStyle -> Green] ]

## ‘Buoyancy frequency’ is such a boring name

Revising my notes on internal gravity waves. Obsessively trying to pronounce ‘Väisälä’. Feeling ve’y silly.

## Task scheduling in MPI

I wrote this simple parallel task-scheduling system in C++/MPI, following discussions with Juha Jäyykä. It seems to work *most* of the time, but fails if jobs finish so quickly or so close to each other that a node does not have the change to update `nexttask_ind`

.

I can’t work out the exact failure criterion, and would welcome any other advice with the code.

/* A simple task scheduler in C++ and MPI. * Usage: * $ mpic++ taskscheduler.cpp -o taskscheduler * $ mpirun -np 4 ./taskscheduler Ntasks */ // see http://stackoverflow.com/questions/11180624/mpi-task-scheduling // and // http://stackoverflow.com/questions/12810391/mpi-asynchronous-broadcast-gather#12810617 #include<iostream> #include<mpi.h> #include<unistd.h> #include<assert.h> #include<time.h> #include<stdlib.h> int main (int argc, char** argv) { MPI_Init (&argc, &argv); /* starts MPI */ MPI_Request mpireq; MPI_Status* mpistatus; assert(argc > 1); int Ntasks = atoi(argv[1]); /* get number of nodes */ int nnodes; MPI_Comm_size (MPI_COMM_WORLD, &nnodes); /* get my node rank */ int rank; MPI_Comm_rank (MPI_COMM_WORLD, &rank); /* get my node hostname */ char hostname[MPI_MAX_PROCESSOR_NAME]; int name_len; MPI_Get_processor_name(hostname, &name_len); /* This is the task that we are currently working on. */ int mytask_ind; /* This buffer stores the next task that needs to be done. */ int nexttask_ind = nnodes; /* We start off by giving task i to node i. */ mytask_ind = rank; srandom(rank); printf("Node %d signing on.\n", rank); while (true) { /* Do our task, if there is any more work to do. */ if (mytask_ind >= Ntasks) break; /* Do our task. */ printf("Node rank %d on %s is now doing task %d.\n", rank, hostname, mytask_ind); /* Whatever we're doing takes a while... */ int usec = 0; while(true) { if (random() % 100000 == 0) break; usec++; usleep(1); /* Every now and then, check for updates to nexttask_ind from other * nodes. */ if (usec % 100 == 0) for (int node = 0; node < nnodes; node++) if (node != rank) { int flag = 0; MPI_Iprobe(node, MPI_ANY_TAG, MPI_COMM_WORLD, &flag, mpistatus); if (flag) { MPI_Irecv( &nexttask_ind, 1, MPI_INT, node, MPI_ANY_TAG, MPI_COMM_WORLD, &mpireq ); } } } /* Advertise the fact that we have started work on this task, by * updating nexttask_ind. */ mytask_ind = nexttask_ind; nexttask_ind++; for (int node = 0; node < nnodes; node++) if (node != rank) MPI_Isend( &nexttask_ind, 1, MPI_INT, node, 0, MPI_COMM_WORLD, &mpireq); } MPI_Finalize(); return 0; }

## Not really qualified to judge…

LinkedIn sometimes doesn’t realise the relationship between you and your head of department: