
Sattolo's Algorithm (2017) - laybak
https://danluu.com/sattolo/
======
Jtsummers
The motivation for writing this amuses me. The Wikipedia prose "proof" is
frustratingly familiar to many math textbooks I had when getting my BS many
years ago. Fortunately, I had a couple professors who insisted on properly
structuring proofs so that they were, you know, actually readable and
understandable.

That Wikipedia prose is like the lazy author's approach, it's a brain dump and
not a real attempt to inform or educate. It is nice that there is Python code
in the relevant section on the algorithm, that's informative (and why I like
Python, it's executable pseudocode). But it would also be nice to have
properly structured proofs to go along with the awful prose section (and
either remove it or clean it up).

Also, here's the discussion from when this was first posted:
[https://news.ycombinator.com/item?id=14967697](https://news.ycombinator.com/item?id=14967697)

~~~
GuB-42
I don't know if it is math but formal scientific writing can make even simple
concepts hard to grasp for me, a developer.

During my studies, I had to assist to a conference by a mathematician
explaining some image processing algorithm. It took me half an hour to realize
that all that was a good old "make each pixel the weighted average of its
surroundings" blur, also known as a convolution filter. And I spent the other
half hour completely lost trying to understand what was special about his
blur. I think no one actually understood anything, in fact, I probably did
better than most because I was familiar with convolution filters from personal
experience.

Had he explained it in terms of pixels and loops instead of what was probably
(for us) a Greek letter soup, I am sure we would have understood.

~~~
qppo
There is a chronic problem in audio/image processing where domain experts have
developed their own notation and terminology and do a really poor job
connecting it to the simple concepts they represent.

> Had he explained it in terms of pixels and loops instead of what was
> probably (for us) a Greek letter soup, I am sure we would have understood.

Unfortunately concretizing the math to realizations leaves you with a cookbook
for implementing the least interesting bits of signal processing. Looking at
some matrix multiplication as for loops doesn't teach you how a filter works.

Like just to give an example, I don't think it's possible to come up with an
explanation of the Remez Exchange Algorithm in any physically meaningful
context since it's fundamentally about recognizing the mathematical properties
of a higher order abstraction of the system designed with it. But we can be a
bit smarter about how we write, since "an iterative algorithm used to find
simple approximations to functions, specifically, approximations by functions
in a Chebyshev space that are the best in the uniform norm L∞ sense" is
Wikipedia word soup that belies the concepts and connections that make it
useful. It's just hard to do it in a sentence and not a semester.

------
kuratkull
Weird, I was actually searching for a function that would do a "permutation of
a list with exactly one cycle" an hour ago. I found Sattolo's, but found its
use of random() disagreeable - it would depend on external/stdlib code and
generate a new sequence every time if the seed is not fixed. Instead I looked
into key schedules for cryptographic ciphers and found an easy solution - the
key schedule function for RC4 for example. (My use does not care about the
flaws/biases). It's a couple of lines of code, generates the same permutation
for the same input+key - just use a couple of arbitrary hardcoded bytes as the
key and you are done.

~~~
GuB-42
If you want a quick, easy and good quality random number generator, you might
want to look at xorshifts.

There are a lot of variations, but I particularly like the 128 bit version
from the original paper [1].

Here is the entire code

    
    
      unsigned long xor128() {
        static unsigned long x=123456789,y=362436069,z=521288629,w=88675123;
        unsigned long t;
        t=(xˆ(x<<11));x=y;y=z;z=w; 
        return( w=(wˆ(w>>19))ˆ(tˆ(t>>8)) );
      }
    

Note: Do not use for security applications! It is not cryptographically
secure. But if you need that, or if you are running Monte Carlo simulations
and want nothing but the best, you probably should use specialized libraries
anyways.

[1]
[https://www.jstatsoft.org/index.php/jss/article/view/v008i14...](https://www.jstatsoft.org/index.php/jss/article/view/v008i14/xorshift.pdf)

Edit: And since we are in the subject of picking random numbers for Sattolo's
algorithm. Using a simple modulo operator for get a range from a single random
number going from 0 to 2^n-1 (like most RNGs produce) can introduce a bias. If
the list is short and the numbers are large, it doesn't make that big a
difference, but it is something to consider.

~~~
KMag
In case anyone is wondering how to do an unbiased choice of the numbers 0 to N
given an RNG that generates numbers 0 to 2^Y-1, where 2^Y-1 >= N, there are
two ways, both rejecting some samples:

    
    
       1. Find the smallest Z s.t. 2^Z-1 <= N, generate a random number, bitwise-and with 2^Z-1, accept the result if the result <= N, otherwise repeat.  The expected number of times through the loop depends on Y and N, but is strictly less than 2.
    
       2. Let W = floor((2^Z-1) / N), generate a random number, if the number < N*W, return the number mod N, otherwise repeat.  The expected number of times through the loop depends on N and Y, but is no more than the first algorithm.  The downside is that division and modular division are much slower than bitwise-and on most processors.
    

If division or modular division (on many processors, the same instruction) is
much slower than bitwise-and, and your PRNG is pretty fast (such as xorshift /
xorishiro), you're better off with the first algorithm.

Note that finding the mask from the first algorithm on a 64-bit machine is as
simple as:

    
    
      mask = N; for( i = 1; i < 64; i *= 2) { mask |= mask >> i; }  // Check if your compiler unrolls this loop

------
yissp
Does anyone have an example for when you'd want to use this algorithm?

~~~
BeeOnRope
Sometimes you'd like to visit all elements of some logical array in a random
order, with each element pointing to the next, e.g., to benchmark the speed of
pointer chasing. Using Sattolo, you can set up such a full-sized linked list
fairly simply, e.g.:

[https://github.com/travisdowns/MemoryLanes/blob/b6cd7d2ceaa0...](https://github.com/travisdowns/MemoryLanes/blob/b6cd7d2ceaa0fc36bd02ddb15544792c4c616973/consoleversion/testingmlp.cpp#L110)

I have used it a half dozen times in more or less that specific context.

