Matlab/Octave codes for the fast construction of good lattice rules
Warning: this code needs some documentation...
For the moment these snippets are literally what can be found in my
PhD thesis, as such they are actually heavily commented.
In my PhD thesis you can find example code as well.
All codes fit on one page and have very fast running times.
This makes the code easy to study and experiment with.
However, the code is not necessarily easy to read: it is high level
code which implements the ideas of my thesis as close as possible.
All this is following the principles set forth by Nick Trefethen in his
Ten Digit Algorithms — “Ten digits, five seconds, and just one page”
A must read!
Please note however that the code will not give you 10 digits...
It will give you useful numerical results, but the numerics could be
In fact this is a long standing issue in the calculation of all kind of
discrepancies where most authors tell you to use as high precision as
Since the fast component-by-component algorithm replaces
an O(n2) calculation by
an O(n log(n)) one the numerical qualities of the
algorithm are slightly improved, but the intrinsic problem still remains.
Fast construction codes
The first code actually has some comments in the header for the
Matlab/Octave help command. The other codes are meant to be used in a
- fastrank1pt.m: component-by-component construction of rank-1 lattice rules with prime number of points and with product weights;
- fastrank1od.m: component-by-component construction of rank-1 lattice rules with prime number of points and with order dependent weights;
- fastrank1expt.m: component-by-component construction of rank-1 lattice sequences with prime power of points and with product weights, this can also be used for fast construction of a rank-1 lattice rule with a prime power of points;
- fastrank1exod.m: component-by-component construction of rank-1 lattice sequences with prime power of points and with order dependent weights, same remark as above.
- fastpoly2rank1pt.m: component-by-component construction of polynomial rank-1 lattice sequences in base 2 and with product weights. (The cludge with getfield is to please Matlab.)
- generatorp.m: find a primitive root modulo n for n prime;
- generator.m: find a primitive root modulo n for any suitable n;
- powmod.m: modular exponentiation
(for Octave you need to make a function usejava which returns false, or just do usejava = @(x) false; or download OCTAVE ONLY: usejava.m);
- vm.m: mapping function for polynomial lattice rules modulo an irreducible polynomial. (The double(w) at the end is to please Matlab.)