[project @ 1999-11-02 12:03:11 by simonpj]
authorsimonpj <unknown>
Tue, 2 Nov 1999 12:03:11 +0000 (12:03 +0000)
committersimonpj <unknown>
Tue, 2 Nov 1999 12:03:11 +0000 (12:03 +0000)
Add Doug McIlroys power-series program

spectral/power/Main.hs [new file with mode: 0644]
spectral/power/Makefile [new file with mode: 0644]
spectral/power/power.stdout [new file with mode: 0644]

diff --git a/spectral/power/Main.hs b/spectral/power/Main.hs
new file mode 100644 (file)
index 0000000..6a522be
--- /dev/null
@@ -0,0 +1,138 @@
+{-     Taken from Doug McIlroy's paper
+       "Power series, power serious"
+       JFP 9(3) May 1999
+-}
+
+module Main where
+
+import IO
+import Ratio
+
+infixl 7 .*
+infixr 5 :+: 
+
+default (Integer, Rational, Double)
+
+main = do { putStrLn (show (extract 50 (sinx - sqrt (1-cosx^2)))) ;
+           putStrLn (show (extract 50 (sinx/cosx - revert (integral (1/(1+x^2)))))) ;
+           putStrLn (show (extract 50 ts)) ;
+           putStrLn (show (extract 50 tree)) 
+         }
+       
+
+-- From Section 6
+tree   = 0 :+: forest
+forest = compose list tree
+list   = 1 :+: list
+
+ts = 1 :+: ts^2
+       
+
+-- The main implementation follows
+data Ps a = Pz | a :+: Ps a
+
+extract :: Int -> Ps a -> [a]
+extract 0 ps        = []
+extract n Pz        = []
+extract n (x :+: ps) = x : extract (n-1) ps
+
+deriv:: Num a => Ps a -> Ps a
+integral:: Fractional a => Ps a -> Ps a
+compose:: Num a => Ps a -> Ps a -> Ps a
+revert:: Fractional a => Ps a -> Ps a
+toList:: Num a => Ps a -> [a]
+takePs:: Num a => Int -> Ps a -> [a]
+(.*):: Num a => a -> Ps a -> Ps a
+x:: Num a => Ps a
+expx, sinx, cosx:: Fractional a => Ps a
+
+c .* Pz = Pz
+c .* (f :+: fs) = c*f :+: c.*fs
+
+x = 0 :+: 1 :+: Pz
+
+toList Pz = []                                         --(0)
+toList (f :+: fs) = f : (toList fs)
+
+takePs n fs = take n (toList fs)
+
+instance Num a => Eq (Ps a) where                      --(1)
+       Pz == Pz = True
+       Pz == (f :+: fs) = f==0 && Pz==fs
+       fs == Pz = Pz==fs
+       (f :+: fs) == (g :+: gs) = f==g && fs==gs
+
+instance Num a => Show (Ps a) where                    --(2)
+       showsPrec p Pz = showsPrec p [0]
+       showsPrec p fs = showsPrec p (toList fs)
+
+instance Num a => Num (Ps a) where
+       negate Pz = Pz
+       negate (f :+: fs) = -f :+: -fs
+
+       Pz + gs = gs
+       fs + Pz = fs
+       (f :+: fs) + (g :+: gs) = f+g :+: fs+gs
+
+       Pz * _ = Pz
+       _ * Pz = Pz
+       (f :+: fs) * (g :+: gs) =
+               f*g :+: f.*gs + g.*fs + x*fs*gs         --(3)
+
+       fromInteger 0 = Pz
+       fromInteger c = fromInteger c :+: Pz
+
+instance Fractional a => Fractional (Ps a) where
+       recip fs = 1/fs
+
+       Pz/Pz = error "power series 0/0"
+       Pz / (0 :+: gs) = Pz / gs
+       Pz / _ = Pz
+       (0 :+: fs) / (0 :+: gs) = fs / gs
+       (f :+: fs) / (g :+: gs) = let q = f/g in
+               q :+: (fs - q.*gs)/(g :+: gs)
+
+compose Pz _ = Pz
+compose (f :+: _) Pz = f :+: Pz
+compose (f :+: fs) (0 :+: gs) = f :+: gs*(compose fs (0 :+: gs))
+compose (f :+: fs) gs = (f :+: Pz) + gs*(compose fs gs)        --(4)
+
+revert (0 :+: fs) = rs where
+       rs = 0 :+: 1/(compose fs rs)
+revert (f0 :+: f1 :+: Pz) = -1/f1 :+: 1/f1 :+: Pz      --(5)
+
+deriv Pz = Pz
+deriv (_ :+: fs) = deriv1 fs 1 where
+       deriv1 Pz _ = Pz
+       deriv1 (f :+: fs) n = n*f :+: (deriv1 fs (n+1))
+
+integral fs = 0 :+: (int1 fs 1) where                  --(6)
+       int1 Pz _ = Pz
+       int1 (f :+: fs) n = f/n :+: (int1 fs (n+1))
+
+instance Fractional a => Floating (Ps a) where
+       sqrt Pz = Pz
+       sqrt (0 :+: 0 :+: fs) = 0 :+: (sqrt fs)
+       sqrt (1 :+: fs) = qs where
+               qs = 1 + integral((deriv (1:+:fs))/(2.*qs))
+
+expx = 1 + (integral expx)
+sinx = integral cosx
+cosx  = 1 - (integral sinx)
+
+--(0) Convert power series to a list; used in printing.
+--(1) Equality works on polynomials; diverges otherwise.
+--(2) Specifies how to print the new data type.
+--(3) x*fs*gs replaces 0:fs*gs to avoid extra zero
+--    at end of product of polynomials; it works
+--    because x is a finite series.
+--(4) This extra production works for the composition
+--    of polynomials with non-zero constant term,
+--    but not for infinite series.
+--(5) Special case for reverting a linear function.
+--(6) There is no special case for (integral Pz)
+--    because this would defeat the property that
+--    (integral) emits one term before evaluating
+--    its operand--a property used in solving
+--    differential equations.
+
diff --git a/spectral/power/Makefile b/spectral/power/Makefile
new file mode 100644 (file)
index 0000000..cf44ea3
--- /dev/null
@@ -0,0 +1,5 @@
+TOP = ../..
+include $(TOP)/mk/boilerplate.mk
+-include opts.mk
+include $(TOP)/mk/target.mk
+
diff --git a/spectral/power/power.stdout b/spectral/power/power.stdout
new file mode 100644 (file)
index 0000000..548408e
--- /dev/null
@@ -0,0 +1,4 @@
+[0.0,0.0,0.0,0.0,0.0,1.734723475976807e-18,-0.0,5.963111948670274e-19,0.0,8.300922883092143e-20,-0.0,9.072516958481608e-21,0.0,9.344297488493562e-22,-0.0,9.51040609914865e-23,0.0,9.649022765261823e-24,-0.0,9.780102876153644e-25,0.0,9.910257086498755e-26,-0.0,1.0041431087911443e-26,0.0,1.0174158482263249e-27,-0.0,1.0308593527566662e-28,0.0,1.0444792983506079e-29,-0.0,1.0582788897503616e-30,0.0,1.0722607231550364e-31,-0.0,1.0864272629796521e-32,0.0,1.1007809639231849e-33,-0.0,1.115324302376871e-34,0.0,1.1300597847253922e-35,-0.0,1.1449899497768106e-36,0.0,1.1601173697065218e-37,-0.0,1.175444650628751e-38,0.0,1.190974433078049e-39]
+[0.0,0.0,0.0,5.551115123125783e-17,0.0,5.551115123125783e-17,0.0,6.938893903907228e-17,0.0,6.938893903907228e-17,0.0,5.898059818321144e-17,0.0,2.6020852139652106e-17,0.0,-2.710505431213761e-17,0.0,-8.174884380540703e-17,0.0,-1.0551997643715172e-16,0.0,-6.343937961755808e-17,0.0,7.108300493358088e-17,0.0,3.119791751327039e-16,0.0,6.743576576599859e-16,0.0,1.2085991251851655e-15,0.0,2.0415929248551994e-15,0.0,3.405919057393068e-15,0.0,5.651485695109e-15,0.0,9.260817345739914e-15,0.0,1.4906042891650998e-14,0.0,2.3575628773554824e-14,0.0,3.6776045155667545e-14,0.0,5.68217272527916e-14,0.0,8.730113111004884e-14,0.0,1.3392789172916284e-13]
+[1.0,1.0,2.0,5.0,14.0,42.0,132.0,429.0,1430.0,4862.0,16796.0,58786.0,208012.0,742900.0,2674440.0,9694845.0,3.535767e7,1.2964479e8,4.776387e8,1.76726319e9,6.56412042e9,2.446626702e10,9.148256364e10,3.4305961365e11,1.289904147324e12,4.861946401452e12,1.8367353072152e13,6.9533550916004e13,2.6374795175036e14,1.002242216651368e15,3.814986502092304e15,1.4544636039226908e16,5.553406487704819e16,2.123361304122431e17,8.129440421497308e17,3.1162854949073014e18,1.1959798385860452e19,4.595080432462175e19,1.767338627870067e20,6.804253717299758e20,2.6221270422764923e21,1.0113918591637899e22,3.904442991190445e22,1.5085347920508534e23,5.833001195929967e23,2.257117854077248e24,8.740328711533173e24,3.3868773757191048e25,1.3132789824216935e26,5.0955224517961704e26]
+[0.0,1.0,1.0,2.0,5.0,14.0,42.0,132.0,429.0,1430.0,4862.0,16796.0,58786.0,208012.0,742900.0,2674440.0,9694845.0,3.535767e7,1.2964479e8,4.776387e8,1.76726319e9,6.56412042e9,2.446626702e10,9.148256364e10,3.4305961365e11,1.289904147324e12,4.861946401452e12,1.8367353072152e13,6.9533550916004e13,2.6374795175036e14,1.002242216651368e15,3.814986502092304e15,1.4544636039226908e16,5.553406487704819e16,2.123361304122431e17,8.129440421497308e17,3.1162854949073014e18,1.1959798385860452e19,4.595080432462175e19,1.767338627870067e20,6.804253717299758e20,2.6221270422764923e21,1.0113918591637899e22,3.904442991190445e22,1.5085347920508534e23,5.833001195929967e23,2.257117854077248e24,8.740328711533173e24,3.3868773757191048e25,1.3132789824216935e26]