diff --git a/changelog b/changelog
index c7165ea..eea576d 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,4 @@
+20080221 wxh src/interp/sfsfun.boot increase precision of PI (236)
 20080221 tpd src/input/gamma.input investigate complex gamma issues
 20080220 tpd src/hyper/bookvol11 add additional hyperdoc page translations
 20080219 tpd src/hyper/Makefile handle hyperdoc bitmaps properly
diff --git a/src/interp/sfsfun.boot.pamphlet b/src/interp/sfsfun.boot.pamphlet
index 50bd7b5..6536e2e 100644
--- a/src/interp/sfsfun.boot.pamphlet
+++ b/src/interp/sfsfun.boot.pamphlet
@@ -105,6 +105,11 @@ negintp(x) ==
         else
                 false
 
+-- Lisp PI is a long float and causes type errors, here we give
+-- enough digits to have double accuracy even after conversion
+-- to binary
+DEFCONSTANT(dfPi,3.14159265358979323846264338328)
+
 --- Small float implementation of Gamma function.  Valid for
 --- real arguments.  Maximum error (relative) seems to be
 --- 2-4 ulps for real x 2<x<9, and up to ten ulps for larger x
@@ -132,7 +137,7 @@ gammaStirling(x) ==
        EXP(lnrgamma(x))
 
 lnrgammaRatapprox(x) ==
-       (x-.5)*LOG(x) - x + LOG(SQRT(2.0*PI)) + phiRatapprox(x)
+       (x-.5)*LOG(x) - x + LOG(SQRT(2.0*dfPi)) + phiRatapprox(x)
 
 phiRatapprox(x) ==
         arg := 1/(x**2)
@@ -167,7 +172,7 @@ gammaRatapprox (x) ==
                      prod := */[x+i for i in 0..n-1]
                      result := gammaRatapprox(reducedarg)/prod
                  else
-                        Pi := PI
+                        Pi := dfPi
                         lx := MULTIPLE_-VALUE_-LIST(FLOOR(x))
                         intpartx := CAR(lx)+1
                         restx := CADR(lx)
@@ -247,12 +252,12 @@ PiMinusLogSinPi(z1,z2,z) ==
         cgammaG(z1,z2)  - logH(z1,z2,z)
 
 cgammaG(z1,z2) ==
-        LOG(2*PI) + PI*z2 - COMPLEX(0.0,1.0)*PI*(z1-.5)
+        LOG(2*dfPi) + dfPi*z2 - COMPLEX(0.0,1.0)*dfPi*(z1-.5)
 
 logH(z1,z2,z) ==
         z1bar := CADR(MULTIPLE_-VALUE_-LIST(FLOOR(z1))) ---frac part of z1
-        piz1bar := PI*z1bar
-        piz2 := PI*z2
+        piz1bar := dfPi*z1bar
+        piz2 := dfPi*z2
         twopiz2 := 2.0*piz2
         i := COMPLEX(0.0,1.0)
         part2 := EXP(twopiz2)*(2.0*SIN(piz1bar)**2 + SIN(2.0*piz1bar)*i)
@@ -284,7 +289,7 @@ logS(z1,z2,z,n,zpn) ==
 --- adjusted by 2 Pi if it is negative.
 cgammaAdjust(z) ==
         if IMAGPART(z)<0.0
-        then result := z + COMPLEX(0.0, 2.0*PI)
+        then result := z + COMPLEX(0.0, 2.0*dfPi)
         else result := z
         result
 
@@ -292,7 +297,7 @@ clngammacase3(z) ==
         (z- .5)*LOG(z) - z + cgammaBernsum(z)
 
 cgammaBernsum (z) ==
-        sum := LOG(2.0*PI)/2.0
+        sum := LOG(2.0*dfPi)/2.0
         zterm := z
         zsquaredinv := 1.0/(z*z)
         l:= [.083333333333333333333, -.0027777777777777777778,_
@@ -365,7 +370,7 @@ rPsi(n,x) ==
                                 skipit := 1
                         else
                                 skipit := 0
-                        sign*((PI**(n+1))*cotdiffeval(n,PI*(-x),skipit) + rPsi(n,1.0-x))
+                        sign*((dfPi**(n+1))*cotdiffeval(n,dfPi*(-x),skipit) + rPsi(n,1.0-x))
         else if n=0
         then
                 - rPsiW(n,x)
@@ -508,7 +513,7 @@ PsiXotic(n,result) ==
 cpsireflect(n,z) ==
         m := MOD(n,2)
         sign := (-1)**m
-        sign*PI**(n+1)*cotdiffeval(n,PI*z,0) + cPsi(n,1.0-z)
+        sign*dfPi**(n+1)*cotdiffeval(n,dfPi*z,0) + cPsi(n,1.0-z)
 
 --- c parameter to 0F1, possibly complex
 --- z argument to 0F1
@@ -600,7 +605,7 @@ f01(c,z)==
 ---             t := SQRT(-z)
 ---             c1 := c-1.0
 ---             p := PHASE(c)
----             if ABS(c)>10.0*ABS(t) and p>=0.0 and PHASE(c)<.90*PI
+---             if ABS(c)>10.0*ABS(t) and p>=0.0 and PHASE(c)<.90*dfPi
 ---             then BesselJAsymptOrder(c1,2*t)*cgamma(c/(t**(c1)))
 ---             else if ABS(t)>10.0*ABS(c) and ABS(t)>50.0
 ---             then BesselJAsympt(c1,2*t)*cgamma(c/(t**(c1)))
@@ -904,7 +909,7 @@ BesselasymptB(mu,z,zsqr,zfth) ==
 --- Asymptotic series only works when |v| < |z|.
 
 BesselJAsympt (v,z) ==
-        pi := PI
+        pi := dfPi
         mu := 4.0*v*v
         zsqr := z*z
         zfth := zsqr*zsqr
@@ -931,9 +936,9 @@ BesselIAsympt(v,z,n) ==
                         (float(r)*eight*z)
                 sum1 := sum1 + term1
                 sum2 := sum2 + ABS(term1)
-        sqrttwopiz := SQRT(two*PI*z)
+        sqrttwopiz := SQRT(two*dfPi*z)
         EXP(z)/sqrttwopiz*(1.0 + sum1 ) +_
-                EXP(-(float(n)+.5)*PI*i)*EXP(-z)/sqrttwopiz*(1.0+ sum2)
+                EXP(-(float(n)+.5)*dfPi*i)*EXP(-z)/sqrttwopiz*(1.0+ sum2)
 
 
 ---Asymptotic formula for BesselJ when order is large comes from
@@ -951,7 +956,7 @@ BesselJAsymptOrder(v,z) ==
     --  cothalpha := 1.0/tanhalpha
         ca := 1.0/tanhalpha
 
-        Pi := PI
+        Pi := dfPi
 	ca2:=ca*ca
 	ca4:=ca2*ca2
 	ca8:=ca4*ca4
@@ -974,7 +979,7 @@ BesselJAsymptOrder(v,z) ==
 ---  See Olver, p. 376-382.
 BesselIAsymptOrder(v,vz) ==
         z := vz/v
-        Pi := PI
+        Pi := dfPi
 ---     Use reflection formula (Atlas, p. 492)  if v not in right half plane;  Is this always accurate?
         if REALPART(v)<0.0
         then return BesselIAsymptOrder(-v,vz) + 2.0/Pi*SIN(-v*Pi)*BesselKAsymptOrder(-v,vz)
@@ -1015,7 +1020,7 @@ BesselKAsymptOrder (v,vz) ==
         u4p := (3675.0/32768.0+(-96833.0/40960.0+(144001.0/16384.0+(-7436429.0/663552.0+37182145.0/7962624.0*p2)*p2)*p2)*p2)*p4
         u5p := ((59535.0/262144.0+(-67608983.0/9175040.0+(250881631.0/5898240.0+(-108313205.0/1179648.0+(5391411025.0/63700992.0-5391411025.0/191102976.0*p2)*p2)*p2)*p2)*p2)*p4*p)*(-1.0)
         hornerresult := horner([u5p,u4p,u3p,u2p,u1p,u0p],vinv)
-        SQRT(PI/(2.0*v))*EXP(-v*eta)/(SQRT(opzsqroh))*hornerresult
+        SQRT(dfPi/(2.0*v))*EXP(-v*eta)/(SQRT(opzsqroh))*hornerresult