diff --git a/programming/random/cpp/sh.h b/programming/random/cpp/sh.h
new file mode 100644
index 0000000..9160697
--- /dev/null
+++ b/programming/random/cpp/sh.h
@@ -0,0 +1,72 @@
+/* Solving the associated legendre polynomials iteratively. Most of this code was taken from
+ * the work described by Robin Green from Sony Computer Entertainment America:
+ *
+ * 	R. Green, Spherical Harmonic Lighting: The Gritty Details, Sony Computer Entertainment America, 2003
+ * 	http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.pdf
+ */
+
+
+/// This function computes the factorial n!
+///		where n! = n * (n-1) * (n-2) * (n-3) * ...... * 2 * 1
+unsigned long factorial(unsigned int n)
+{
+    if (n == 0)						//the standard definition of a factorial specifies that 0! = 1
+       return 1;
+    return n * factorial(n - 1);	//recursively call this function
+}
+
+/// Evaluate an Associated Legendre Polynomial P(l,m,x) at x
+/// This implements evaluation of the associated Legendre polynomials using
+///		a recursive method with the following rules:
+///
+///		1) (l - m)P(l, m) = x(2l-1)P(l-1, m) - (l + m - 1)P(l-2, m)
+///		2) P(m, m)        = (-1)^m (2m - 1)!! (1 - x^2)^{m/2}
+///		3) P(m+1, m)      = x (2m + 1)P(m, m)
+double P(int l,int m,double x)
+{
+
+	double pmm = 1.0;							//first iteration, P(0, 0) = 1
+	if(m>0) {									//solve for positive values of m
+		double somx2 = sqrt((1.0-x)*(1.0+x));
+		double fact = 1.0;
+		for(int i=1; i<=m; i++) {
+			pmm *= (-fact) * somx2;
+			fact += 2.0;
+		}
+	}
+	if(l==m) return pmm;						//if the degree equals the order, the solution is found above
+	double pmmp1 = x * (2.0*m+1.0) * pmm;
+	if(l==m+1) return pmmp1;
+	double pll = 0.0;
+	for(int ll=m+2; ll<=l; ++ll) {
+		pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
+		pmm = pmmp1;
+		pmmp1 = pll;
+	}
+	return pll;
+}
+
+/// Calculate the scaling factor for a normalized Legendre polynomial
+/// l is the degree, ranging from [0..N]
+/// m is the order in the range [-l..l]
+double K(int l, int m)
+{
+	static const double K_PI = 3.14159;
+	// renormalisation constant for SH function
+	double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*K_PI*factorial(l+m));
+	return sqrt(temp);
+}
+
+/// return a point sample of a Spherical Harmonic basis function
+/// l is the degree, range [0..N]
+/// m is the order in the range [-l..l]
+/// theta in the range [0..Pi]
+/// phi in the range [0..2*Pi]
+double SH(int l, int m, double theta, double phi){
+
+	const double sqrt2 = sqrt(2.0);
+	if(m==0) return K(l,0)*P(l,m,cos(theta));
+	else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
+	else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
+}
+
--
libgit2 0.21.4