Randomness

gistfile1.m

// generator returns values in range [0, 1)
typedef double (^CWRandomGenerator)(void);
CWRandomGenerator CWMakeRandomGenerator(void);

// generate a normally-distributed number with the given mean and standard deviation
// see http://en.wikipedia.org/wiki/Normal_distribution
double CWNormalDistribution(CWRandomGenerator rnd, double mean, double stddev);




CWRandomGenerator CWMakeRandomGenerator(void)
{
    static const double randmax = 0x80000000;
    
    NSMutableData *seedData = [NSMutableData dataWithLength: 3 * sizeof(unsigned short)];
    int ret = SecRandomCopyBytes(kSecRandomDefault, [seedData length], [seedData mutableBytes]);
    if(ret == -1)
        NSLog(@"SecRandomCopyBytes failed with error %d (%s)", errno, strerror(errno));
    NSLog(@"Using random seed %@", seedData);
    
    CWRandomGenerator generator = ^{
        return nrand48([seedData mutableBytes]) / randmax;
    };
    
    return [[generator copy] autorelease];
}

double CWStandardNormalDistribution(CWRandomGenerator rnd)
{
    // formula from http://en.wikipedia.org/wiki/Box–Muller_transform
    // this can generate two random variables from two inputs
    // because we are lazy, we just toss one
    
    // u_1 and u_2 are supposed to be in the range (0, 1]
    // rnd gives [0, 1), so convert
    // (this is extremely unlikely to matter)
    double u_1 = 1.0 - rnd();
    double u_2 = 1.0 - rnd();
    
    double z_0 = sqrt(-2 * log(u_1)) * cos(2 * M_PI * u_2);
    
    return z_0;
}

double CWNormalDistribution(CWRandomGenerator rnd, double mean, double stddev)
{
    // per Wikipedia: "The parameter σ^2 is called the variance...."
    // "The square root of σ^2 is called the standard deviation...."
    // "...a N(µ, σ^2) can be generated as X = µ + σZ, where Z is standard normal."
    
    double standardNormal = CWStandardNormalDistribution(rnd);
    return mean + stddev * standardNormal;
}