diff --git a/src/pkg/math/acosh.go b/src/pkg/math/acosh.go
index 13afc5a..d8067c0 100644
--- a/src/pkg/math/acosh.go
+++ b/src/pkg/math/acosh.go
@@ -42,7 +42,7 @@
 func Acosh(x float64) float64 {
 	const (
 		Ln2   = 6.93147180559945286227e-01 // 0x3FE62E42FEFA39EF
-		Large = 1 << 28                    // 2^28
+		Large = 1 << 28                    // 2**28
 	)
 	// TODO(rsc): Remove manual inlining of IsNaN
 	// when compiler does it for us
@@ -53,9 +53,9 @@
 	case x == 1:
 		return 0
 	case x >= Large:
-		return Log(x) + Ln2 // x > 2^28
+		return Log(x) + Ln2 // x > 2**28
 	case x > 2:
-		return Log(2*x - 1/(x+Sqrt(x*x-1))) // 2^28 > x > 2
+		return Log(2*x - 1/(x+Sqrt(x*x-1))) // 2**28 > x > 2
 	}
 	t := x - 1
 	return Log1p(t + Sqrt(2*t+t*t)) // 2 >= x > 1
diff --git a/src/pkg/math/asinh.go b/src/pkg/math/asinh.go
index b38bbd7..90dcd27 100644
--- a/src/pkg/math/asinh.go
+++ b/src/pkg/math/asinh.go
@@ -28,7 +28,7 @@
 //	asinh(x) := x  if  1+x*x=1,
 //	         := sign(x)*(log(x)+ln2)) for large |x|, else
 //	         := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
-//	         := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
+//	         := sign(x)*log1p(|x| + x**2/(1 + sqrt(1+x**2)))
 //
 
 // Asinh(x) calculates the inverse hyperbolic sine of x.
@@ -40,8 +40,8 @@
 func Asinh(x float64) float64 {
 	const (
 		Ln2      = 6.93147180559945286227e-01 // 0x3FE62E42FEFA39EF
-		NearZero = 1.0 / (1 << 28)            // 2^-28
-		Large    = 1 << 28                    // 2^28
+		NearZero = 1.0 / (1 << 28)            // 2**-28
+		Large    = 1 << 28                    // 2**28
 	)
 	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
 	// when compiler does it for us
@@ -57,13 +57,13 @@
 	var temp float64
 	switch {
 	case x > Large:
-		temp = Log(x) + Ln2 // |x| > 2^28
+		temp = Log(x) + Ln2 // |x| > 2**28
 	case x > 2:
-		temp = Log(2*x + 1/(Sqrt(x*x+1)+x)) // 2^28 > |x| > 2.0
+		temp = Log(2*x + 1/(Sqrt(x*x+1)+x)) // 2**28 > |x| > 2.0
 	case x < NearZero:
-		temp = x // |x| < 2^-28
+		temp = x // |x| < 2**-28
 	default:
-		temp = Log1p(x + x*x/(1+Sqrt(1+x*x))) // 2.0 > |x| > 2^-28
+		temp = Log1p(x + x*x/(1+Sqrt(1+x*x))) // 2.0 > |x| > 2**-28
 	}
 	if sign {
 		temp = -temp
diff --git a/src/pkg/math/atanh.go b/src/pkg/math/atanh.go
index 72ae2a6..6aecb7b 100644
--- a/src/pkg/math/atanh.go
+++ b/src/pkg/math/atanh.go
@@ -45,7 +45,7 @@
 //	Atanh(-1) = -Inf
 //	Atanh(NaN) = NaN
 func Atanh(x float64) float64 {
-	const NearZero = 1.0 / (1 << 28) // 2^-28
+	const NearZero = 1.0 / (1 << 28) // 2**-28
 	// TODO(rsc): Remove manual inlining of IsNaN
 	// when compiler does it for us
 	// special cases
diff --git a/src/pkg/math/const.go b/src/pkg/math/const.go
index 965a440..6a78d00 100644
--- a/src/pkg/math/const.go
+++ b/src/pkg/math/const.go
@@ -27,11 +27,11 @@
 // Max is the largest finite value representable by the type.
 // Min is the smallest nonzero value representable by the type.
 const (
-	MaxFloat32 = 3.40282346638528859811704183484516925440e+38  /* 2^127 * (2^24 - 1) / 2^23 */
-	MinFloat32 = 1.401298464324817070923729583289916131280e-45 /* 1 / 2^(127 - 1 + 23) */
+	MaxFloat32 = 3.40282346638528859811704183484516925440e+38  /* 2**127 * (2**24 - 1) / 2**23 */
+	MinFloat32 = 1.401298464324817070923729583289916131280e-45 /* 1 / 2**(127 - 1 + 23) */
 
-	MaxFloat64 = 1.797693134862315708145274237317043567981e+308 /* 2^1023 * (2^53 - 1) / 2^52 */
-	MinFloat64 = 4.940656458412465441765687928682213723651e-324 /* 1 / 2^(1023 - 1 + 52) */
+	MaxFloat64 = 1.797693134862315708145274237317043567981e+308 /* 2**1023 * (2**53 - 1) / 2**52 */
+	MinFloat64 = 4.940656458412465441765687928682213723651e-324 /* 1 / 2**(1023 - 1 + 52) */
 )
 
 // Integer limit values.
diff --git a/src/pkg/math/erf.go b/src/pkg/math/erf.go
index f01f806..b608999 100644
--- a/src/pkg/math/erf.go
+++ b/src/pkg/math/erf.go
@@ -39,7 +39,7 @@
 //
 // Method:
 //      1. For |x| in [0, 0.84375]
-//          erf(x)  = x + x*R(x^2)
+//          erf(x)  = x + x*R(x**2)
 //          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
 //                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
 //         where R = P/Q where P is an odd poly of degree 8 and
@@ -49,7 +49,7 @@
 //
 //
 //         Remark. The formula is derived by noting
-//          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
+//          erf(x) = (2/sqrt(pi))*(x - x**3/3 + x**5/10 - x**7/42 + ....)
 //         and that
 //          2/sqrt(pi) = 1.128379167095512573896158903121545171688
 //         is close to one. The interval is chosen because the fix
@@ -77,7 +77,7 @@
 //              erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
 //              erf(x)  = 1 - erfc(x)
 //         where
-//              R1(z) = degree 7 poly in z, (z=1/x^2)
+//              R1(z) = degree 7 poly in z, (z=1/x**2)
 //              S1(z) = degree 8 poly in z
 //
 //      4. For x in [1/0.35,28]
@@ -87,7 +87,7 @@
 //              erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6, else
 //              erf(x)  = sign(x)*(1.0 - tiny)
 //         where
-//              R2(z) = degree 6 poly in z, (z=1/x^2)
+//              R2(z) = degree 6 poly in z, (z=1/x**2)
 //              S2(z) = degree 7 poly in z
 //
 //      Note1:
@@ -99,10 +99,10 @@
 //      Note2:
 //         Here 4 and 5 make use of the asymptotic series
 //                        exp(-x*x)
-//              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
+//              erfc(x) ~ ---------- * ( 1 + Poly(1/x**2) )
 //                        x*sqrt(pi)
 //         We use rational approximation to approximate
-//              g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
+//              g(s)=f(1/x**2) = log(erfc(x)*x) - x*x + 0.5625
 //         Here is the error bound for R1/S1 and R2/S2
 //              |R1/S1 - f(x)|  < 2**(-62.57)
 //              |R2/S2 - f(x)|  < 2**(-61.52)
@@ -189,7 +189,7 @@
 func Erf(x float64) float64 {
 	const (
 		VeryTiny = 2.848094538889218e-306 // 0x0080000000000000
-		Small    = 1.0 / (1 << 28)        // 2^-28
+		Small    = 1.0 / (1 << 28)        // 2**-28
 	)
 	// special cases
 	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
@@ -209,7 +209,7 @@
 	}
 	if x < 0.84375 { // |x| < 0.84375
 		var temp float64
-		if x < Small { // |x| < 2^-28
+		if x < Small { // |x| < 2**-28
 			if x < VeryTiny {
 				temp = 0.125 * (8.0*x + efx8*x) // avoid underflow
 			} else {
@@ -266,7 +266,7 @@
 //	Erfc(-Inf) = 2
 //	Erfc(NaN) = NaN
 func Erfc(x float64) float64 {
-	const Tiny = 1.0 / (1 << 56) // 2^-56
+	const Tiny = 1.0 / (1 << 56) // 2**-56
 	// special cases
 	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
 	// when compiler does it for us
@@ -285,7 +285,7 @@
 	}
 	if x < 0.84375 { // |x| < 0.84375
 		var temp float64
-		if x < Tiny { // |x| < 2^-56
+		if x < Tiny { // |x| < 2**-56
 			temp = x
 		} else {
 			z := x * x
diff --git a/src/pkg/math/exp.go b/src/pkg/math/exp.go
index 18ec684..b801d58 100644
--- a/src/pkg/math/exp.go
+++ b/src/pkg/math/exp.go
@@ -59,7 +59,7 @@
 //
 //   3. Scale back to obtain exp(x):
 //      From step 1, we have
-//         exp(x) = 2^k * exp(r)
+//         exp(x) = 2**k * exp(r)
 //
 // Special cases:
 //      exp(INF) is INF, exp(NaN) is NaN;
@@ -81,7 +81,7 @@
 // compiler will convert from decimal to binary accurately enough
 // to produce the hexadecimal values shown.
 
-// Exp returns e^x, the base-e exponential of x.
+// Exp returns e**x, the base-e exponential of x.
 //
 // Special cases are:
 //	Exp(+Inf) = +Inf
@@ -101,7 +101,7 @@
 
 		Overflow  = 7.09782712893383973096e+02
 		Underflow = -7.45133219101941108420e+02
-		NearZero  = 1.0 / (1 << 28) // 2^-28
+		NearZero  = 1.0 / (1 << 28) // 2**-28
 	)
 
 	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
@@ -140,7 +140,7 @@
 	return Ldexp(y, k)
 }
 
-// Exp2 returns 2^x, the base-2 exponential of x.
+// Exp2 returns 2**x, the base-2 exponential of x.
 //
 // Special cases are the same as Exp.
 func Exp2(x float64) float64 { return Exp(x * Ln2) }
diff --git a/src/pkg/math/expm1.go b/src/pkg/math/expm1.go
index 9e8ae3f..35100ca 100644
--- a/src/pkg/math/expm1.go
+++ b/src/pkg/math/expm1.go
@@ -34,13 +34,13 @@
 //   2. Approximating expm1(r) by a special rational function on
 //      the interval [0,0.34658]:
 //      Since
-//          r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
+//          r*(exp(r)+1)/(exp(r)-1) = 2+ r**2/6 - r**4/360 + ...
 //      we define R1(r*r) by
-//          r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
+//          r*(exp(r)+1)/(exp(r)-1) = 2+ r**2/6 * R1(r*r)
 //      That is,
 //          R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
 //                   = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
-//                   = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
+//                   = 1 - r**2/60 + r**4/2520 - r**6/100800 + ...
 //      We use a special Reme algorithm on [0,0.347] to generate
 //      a polynomial of degree 5 in r*r to approximate R1. The
 //      maximum error of this polynomial approximation is bounded
@@ -79,20 +79,20 @@
 //                 = r - E
 //   3. Scale back to obtain expm1(x):
 //      From step 1, we have
-//         expm1(x) = either 2^k*[expm1(r)+1] - 1
-//                  = or     2^k*[expm1(r) + (1-2^-k)]
+//         expm1(x) = either 2**k*[expm1(r)+1] - 1
+//                  = or     2**k*[expm1(r) + (1-2**-k)]
 //   4. Implementation notes:
 //      (A). To save one multiplication, we scale the coefficient Qi
-//           to Qi*2^i, and replace z by (x^2)/2.
+//           to Qi*2**i, and replace z by (x**2)/2.
 //      (B). To achieve maximum accuracy, we compute expm1(x) by
 //        (i)   if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
 //        (ii)  if k=0, return r-E
 //        (iii) if k=-1, return 0.5*(r-E)-0.5
 //        (iv)  if k=1 if r < -0.25, return 2*((r+0.5)- E)
 //                     else          return  1.0+2.0*(r-E);
-//        (v)   if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
-//        (vi)  if k <= 20, return 2^k((1-2^-k)-(E-r)), else
-//        (vii) return 2^k(1-((E+2^-k)-r))
+//        (v)   if (k<-2||k>56) return 2**k(1-(E-r)) - 1 (or exp(x)-1)
+//        (vi)  if k <= 20, return 2**k((1-2**-k)-(E-r)), else
+//        (vii) return 2**k(1-((E+2**-k)-r))
 //
 // Special cases:
 //      expm1(INF) is INF, expm1(NaN) is NaN;
@@ -114,7 +114,7 @@
 // to produce the hexadecimal values shown.
 //
 
-// Expm1 returns e^x - 1, the base-e exponential of x minus 1.
+// Expm1 returns e**x - 1, the base-e exponential of x minus 1.
 // It is more accurate than Exp(x) - 1 when x is near zero.
 //
 // Special cases are:
@@ -131,7 +131,7 @@
 		Ln2Hi      = 6.93147180369123816490e-01 // 0x3fe62e42fee00000
 		Ln2Lo      = 1.90821492927058770002e-10 // 0x3dea39ef35793c76
 		InvLn2     = 1.44269504088896338700e+00 // 0x3ff71547652b82fe
-		Tiny       = 1.0 / (1 << 54)            // 2^-54 = 0x3c90000000000000
+		Tiny       = 1.0 / (1 << 54)            // 2**-54 = 0x3c90000000000000
 		// scaled coefficients related to expm1
 		Q1 = -3.33333333333331316428e-02 // 0xBFA11111111110F4
 		Q2 = 1.58730158725481460165e-03  // 0x3F5A01A019FE5585
@@ -194,7 +194,7 @@
 		}
 		x = hi - lo
 		c = (hi - x) - lo
-	} else if absx < Tiny { // when |x| < 2^-54, return x
+	} else if absx < Tiny { // when |x| < 2**-54, return x
 		return x
 	} else {
 		k = 0
@@ -223,12 +223,12 @@
 			return y - 1
 		}
 		if k < 20 {
-			t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2^-k
+			t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2**-k
 			y := t - (e - x)
 			y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
 			return y
 		}
-		t := Float64frombits(uint64((0x3ff - k) << 52)) // 2^-k
+		t := Float64frombits(uint64((0x3ff - k) << 52)) // 2**-k
 		y := x - (e + t)
 		y += 1
 		y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
diff --git a/src/pkg/math/frexp.go b/src/pkg/math/frexp.go
index ab226e7..1e497fd 100644
--- a/src/pkg/math/frexp.go
+++ b/src/pkg/math/frexp.go
@@ -6,7 +6,7 @@
 
 // Frexp breaks f into a normalized fraction
 // and an integral power of two.
-// It returns frac and exp satisfying f == frac × 2^exp,
+// It returns frac and exp satisfying f == frac × 2**exp,
 // with the absolute value of frac in the interval [½, 1).
 func Frexp(f float64) (frac float64, exp int) {
 	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
diff --git a/src/pkg/math/j0.go b/src/pkg/math/j0.go
index 8a6db3b..5aaf4ab 100644
--- a/src/pkg/math/j0.go
+++ b/src/pkg/math/j0.go
@@ -25,11 +25,11 @@
 // __ieee754_j0(x), __ieee754_y0(x)
 // Bessel function of the first and second kinds of order zero.
 // Method -- j0(x):
-//      1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
+//      1. For tiny x, we use j0(x) = 1 - x**2/4 + x**4/64 - ...
 //      2. Reduce x to |x| since j0(x)=j0(-x),  and
 //         for x in (0,2)
-//              j0(x) = 1-z/4+ z^2*R0/S0,  where z = x*x;
-//         (precision:  |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
+//              j0(x) = 1-z/4+ z**2*R0/S0,  where z = x*x;
+//         (precision:  |j0-1+z/4-z**2R0/S0 |<2**-63.67 )
 //         for x in (2,inf)
 //              j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
 //         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
@@ -50,13 +50,13 @@
 // Method -- y0(x):
 //      1. For x<2.
 //         Since
-//              y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
+//              y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x**2/4 - ...)
 //         therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
 //         We use the following function to approximate y0,
-//              y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
+//              y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x**2
 //         where
-//              U(z) = u00 + u01*z + ... + u06*z^6
-//              V(z) = 1  + v01*z + ... + v04*z^4
+//              U(z) = u00 + u01*z + ... + u06*z**6
+//              V(z) = 1  + v01*z + ... + v04*z**4
 //         with absolute approximation error bounded by 2**-72.
 //         Note: For tiny x, U/V = u0 and j0(x)~1, hence
 //              y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
@@ -232,11 +232,11 @@
 }
 
 // The asymptotic expansions of pzero is
-//      1 - 9/128 s^2 + 11025/98304 s^4 - ...,	where s = 1/x.
+//      1 - 9/128 s**2 + 11025/98304 s**4 - ..., where s = 1/x.
 // For x >= 2, We approximate pzero by
 // 	pzero(x) = 1 + (R/S)
-// where  R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
-// 	  S = 1 + pS0*s^2 + ... + pS4*s^10
+// where  R = pR0 + pR1*s**2 + pR2*s**4 + ... + pR5*s**10
+// 	  S = 1 + pS0*s**2 + ... + pS4*s**10
 // and
 //      | pzero(x)-1-R/S | <= 2  ** ( -60.26)
 
@@ -331,13 +331,13 @@
 }
 
 // For x >= 8, the asymptotic expansions of qzero is
-//      -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
+//      -1/8 s + 75/1024 s**3 - ..., where s = 1/x.
 // We approximate pzero by
-// 	qzero(x) = s*(-1.25 + (R/S))
-// where  R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
-// 	  S = 1 + qS0*s^2 + ... + qS5*s^12
+//      qzero(x) = s*(-1.25 + (R/S))
+// where R = qR0 + qR1*s**2 + qR2*s**4 + ... + qR5*s**10
+//       S = 1 + qS0*s**2 + ... + qS5*s**12
 // and
-//      | qzero(x)/s +1.25-R/S | <= 2  ** ( -61.22)
+//      | qzero(x)/s +1.25-R/S | <= 2**(-61.22)
 
 // for x in [inf, 8]=1/[0,0.125]
 var q0R8 = [6]float64{
diff --git a/src/pkg/math/j1.go b/src/pkg/math/j1.go
index 5c7b799..278162e 100644
--- a/src/pkg/math/j1.go
+++ b/src/pkg/math/j1.go
@@ -25,7 +25,7 @@
 // __ieee754_j1(x), __ieee754_y1(x)
 // Bessel function of the first and second kinds of order one.
 // Method -- j1(x):
-//      1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
+//      1. For tiny x, we use j1(x) = x/2 - x**3/16 + x**5/384 - ...
 //      2. Reduce x to |x| since j1(x)=-j1(-x),  and
 //         for x in (0,2)
 //              j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
@@ -52,13 +52,13 @@
 //      1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
 //      2. For x<2.
 //         Since
-//              y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
+//              y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x**3-...)
 //         therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
 //         We use the following function to approximate y1,
-//              y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
+//              y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x**2
 //         where for x in [0,2] (abs err less than 2**-65.89)
-//              U(z) = U0[0] + U0[1]*z + ... + U0[4]*z^4
-//              V(z) = 1  + v0[0]*z + ... + v0[4]*z^5
+//              U(z) = U0[0] + U0[1]*z + ... + U0[4]*z**4
+//              V(z) = 1  + v0[0]*z + ... + v0[4]*z**5
 //         Note: For tiny x, 1/x dominate y1 and hence
 //              y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
 //      3. For x>=2.
@@ -225,11 +225,11 @@
 }
 
 // For x >= 8, the asymptotic expansions of pone is
-//      1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x.
+//      1 + 15/128 s**2 - 4725/2**15 s**4 - ..., where s = 1/x.
 // We approximate pone by
 //      pone(x) = 1 + (R/S)
-// where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
-//       S = 1 + ps0*s^2 + ... + ps4*s^10
+// where R = pr0 + pr1*s**2 + pr2*s**4 + ... + pr5*s**10
+//       S = 1 + ps0*s**2 + ... + ps4*s**10
 // and
 //      | pone(x)-1-R/S | <= 2**(-60.06)
 
@@ -324,11 +324,11 @@
 }
 
 // For x >= 8, the asymptotic expansions of qone is
-//      3/8 s - 105/1024 s^3 - ..., where s = 1/x.
+//      3/8 s - 105/1024 s**3 - ..., where s = 1/x.
 // We approximate qone by
 //      qone(x) = s*(0.375 + (R/S))
-// where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
-//       S = 1 + qs1*s^2 + ... + qs6*s^12
+// where R = qr1*s**2 + qr2*s**4 + ... + qr5*s**10
+//       S = 1 + qs1*s**2 + ... + qs6*s**12
 // and
 //      | qone(x)/s -0.375-R/S | <= 2**(-61.13)
 
diff --git a/src/pkg/math/jn.go b/src/pkg/math/jn.go
index ecd7ab6..7d31743 100644
--- a/src/pkg/math/jn.go
+++ b/src/pkg/math/jn.go
@@ -64,7 +64,7 @@
 	case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0):
 		return 0
 	}
-	// J(-n, x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
+	// J(-n, x) = (-1)**n * J(n, x), J(n, -x) = (-1)**n * J(n, x)
 	// Thus, J(-n, x) = J(n, -x)
 
 	if n == 0 {
@@ -125,7 +125,7 @@
 	} else {
 		if x < TwoM29 { // x < 2**-29
 			// x is tiny, return the first Taylor expansion of J(n,x)
-			// J(n,x) = 1/n!*(x/2)^n  - ...
+			// J(n,x) = 1/n!*(x/2)**n  - ...
 
 			if n > 33 { // underflow
 				b = 0
@@ -135,13 +135,13 @@
 				a := float64(1)
 				for i := 2; i <= n; i++ {
 					a *= float64(i) // a = n!
-					b *= temp       // b = (x/2)^n
+					b *= temp       // b = (x/2)**n
 				}
 				b /= a
 			}
 		} else {
 			// use backward recurrence
-			//                      x      x^2      x^2
+			//                      x      x**2      x**2
 			//  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 			//                      2n  - 2(n+1) - 2(n+2)
 			//
@@ -187,7 +187,7 @@
 			}
 			a := t
 			b = 1
-			//  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
+			//  estimate log((2/x)**n*n!) = n*log(2/x)+n*ln(n)
 			//  Hence, if n*(log(2n/x)) > ...
 			//  single 8.8722839355e+01
 			//  double 7.09782712893383973096e+02
diff --git a/src/pkg/math/ldexp.go b/src/pkg/math/ldexp.go
index dc6b3a8..ab8392a 100644
--- a/src/pkg/math/ldexp.go
+++ b/src/pkg/math/ldexp.go
@@ -5,7 +5,7 @@
 package math
 
 // Ldexp is the inverse of Frexp.
-// It returns frac × 2^exp.
+// It returns frac × 2**exp.
 func Ldexp(frac float64, exp int) float64 {
 	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
 	// when compiler does it for us
diff --git a/src/pkg/math/lgamma.go b/src/pkg/math/lgamma.go
index 38558f3..dc31be9 100644
--- a/src/pkg/math/lgamma.go
+++ b/src/pkg/math/lgamma.go
@@ -39,7 +39,7 @@
 //      minimum (ymin=1.461632144968362245) to maintain monotonicity.
 //      On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
 //              Let z = x-ymin;
-//              lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
+//              lgamma(x) = -1.214862905358496078218 + z**2*poly(z)
 //              poly(z) is a 14 degree polynomial.
 //   2. Rational approximation in the primary interval [2,3]
 //      We use the following approximation:
diff --git a/src/pkg/math/log.go b/src/pkg/math/log.go
index 1727c77..02e767b 100644
--- a/src/pkg/math/log.go
+++ b/src/pkg/math/log.go
@@ -27,7 +27,7 @@
 //
 // Method :
 //   1. Argument Reduction: find k and f such that
-//			x = 2^k * (1+f),
+//			x = 2**k * (1+f),
 //	   where  sqrt(2)/2 < 1+f < sqrt(2) .
 //
 //   2. Approximation of log(1+f).
diff --git a/src/pkg/math/log1p.go b/src/pkg/math/log1p.go
index 87a8b0e..e1fc275 100644
--- a/src/pkg/math/log1p.go
+++ b/src/pkg/math/log1p.go
@@ -24,7 +24,7 @@
 //
 // Method :
 //   1. Argument Reduction: find k and f such that
-//                      1+x = 2^k * (1+f),
+//                      1+x = 2**k * (1+f),
 //         where  sqrt(2)/2 < 1+f < sqrt(2) .
 //
 //      Note. If k=0, then f=x is exact. However, if k!=0, then f
@@ -96,9 +96,9 @@
 	const (
 		Sqrt2M1     = 4.142135623730950488017e-01  // Sqrt(2)-1 = 0x3fda827999fcef34
 		Sqrt2HalfM1 = -2.928932188134524755992e-01 // Sqrt(2)/2-1 = 0xbfd2bec333018866
-		Small       = 1.0 / (1 << 29)              // 2^-29 = 0x3e20000000000000
-		Tiny        = 1.0 / (1 << 54)              // 2^-54
-		Two53       = 1 << 53                      // 2^53
+		Small       = 1.0 / (1 << 29)              // 2**-29 = 0x3e20000000000000
+		Tiny        = 1.0 / (1 << 54)              // 2**-54
+		Two53       = 1 << 53                      // 2**53
 		Ln2Hi       = 6.93147180369123816490e-01   // 3fe62e42fee00000
 		Ln2Lo       = 1.90821492927058770002e-10   // 3dea39ef35793c76
 		Lp1         = 6.666666666666735130e-01     // 3FE5555555555593
@@ -131,8 +131,8 @@
 	var iu uint64
 	k := 1
 	if absx < Sqrt2M1 { //  |x| < Sqrt(2)-1
-		if absx < Small { // |x| < 2^-29
-			if absx < Tiny { // |x| < 2^-54
+		if absx < Small { // |x| < 2**-29
+			if absx < Tiny { // |x| < 2**-54
 				return x
 			}
 			return x - x*x*0.5
@@ -175,7 +175,7 @@
 	}
 	hfsq := 0.5 * f * f
 	var s, R, z float64
-	if iu == 0 { // |f| < 2^-20
+	if iu == 0 { // |f| < 2**-20
 		if f == 0 {
 			if k == 0 {
 				return 0
diff --git a/src/pkg/math/pow.go b/src/pkg/math/pow.go
index ab8bdb6..259deb1 100644
--- a/src/pkg/math/pow.go
+++ b/src/pkg/math/pow.go
@@ -82,11 +82,11 @@
 		return Exp(y * Log(x))
 	}
 
-	// ans = a1 * 2^ae (= 1 for now).
+	// ans = a1 * 2**ae (= 1 for now).
 	a1 := float64(1)
 	ae := 0
 
-	// ans *= x^yf
+	// ans *= x**yf
 	if yf != 0 {
 		if yf > 0.5 {
 			yf--
@@ -95,7 +95,7 @@
 		a1 = Exp(yf * Log(x))
 	}
 
-	// ans *= x^yi
+	// ans *= x**yi
 	// by multiplying in successive squarings
 	// of x according to bits of yi.
 	// accumulate powers of two into exp.
@@ -113,7 +113,7 @@
 		}
 	}
 
-	// ans = a1*2^ae
+	// ans = a1*2**ae
 	// if flip { ans = 1 / ans }
 	// but in the opposite order
 	if flip {
diff --git a/src/pkg/math/sqrt_port.go b/src/pkg/math/sqrt_port.go
index 125afcd..c818834 100644
--- a/src/pkg/math/sqrt_port.go
+++ b/src/pkg/math/sqrt_port.go
@@ -31,8 +31,8 @@
 //   Bit by bit method using integer arithmetic. (Slow, but portable)
 //   1. Normalization
 //      Scale x to y in [1,4) with even powers of 2:
-//      find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
-//              sqrt(x) = 2^k * sqrt(y)
+//      find an integer k such that  1 <= (y=x*2**(2k)) < 4, then
+//              sqrt(x) = 2**k * sqrt(y)
 //   2. Bit by bit computation
 //      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
 //           i                                                   0