small algorithm to approximate an integral

combinators.hs

@@ -13,7 +13,7 @@ instance Show Value where
 
 main :: IO ()
 main = do
-  putStr "δ> "
+  putStr "λ> "
   hFlush stdout
   input <- getLine
   stepEval $ parse input
@@ -32,17 +32,24 @@ stepEval (hed:rest) = do
         Atom combinator ->
           case combinator of
             'I' -> stepEval rest
-            'S' -> stepEval (s rest)
-            'K' -> stepEval (k rest)
-            'B' -> stepEval (b rest)
-            'C' -> stepEval (c rest)
-            'W' -> stepEval (w rest)
-            'M' -> stepEval (m rest)
-            'T' -> stepEval (t rest)
+            'S' -> apply s rest
+            'K' -> apply k rest
+            'B' -> apply b rest
+            'C' -> apply c rest
+            'W' -> apply w rest
+            'M' -> apply m rest
+            'T' -> apply t rest
             _   -> return (hed:rest)
         Compound expr -> do
           stepEval $ expr ++ rest
 
+apply :: Combinator -> [Value] -> IO [Value]
+apply f vals = case f vals of
+  Just xs -> stepEval xs
+  Nothing -> do
+    hPutStrLn stderr "Invalid application of combinator."
+    stepEval []
+
 -- Parsing
 
 parse s = parseAux s []
@@ -64,24 +71,32 @@ parseAux (x:rest) v
 
 -- Combinators
 
+type Combinator = [Value] -> Maybe [Value]
+
 -- Sxyz ~> xz(yz)
-s (x:y:z:rest) = x:z:Compound [y,z]:rest
+s (x:y:z:rest) = Just $ x:z:Compound [y,z]:rest
+s _ = Nothing
 
 -- Kxy  ~> y
-k (x:y:rest) = x:rest
-k _ = error "bad call to K"
+k (x:y:rest) = Just $ x:rest
+k _ = Nothing
 
 -- Bxyz ~> x(yz)
-b (x:y:z:rest) = x:Compound [y,z]:rest
+b (x:y:z:rest) = Just $ x:Compound [y,z]:rest
+b _ = Nothing
 
 -- Cxyz ~> xzy
-c (x:y:z:rest) = x:z:y:rest
+c (x:y:z:rest) = Just $ x:z:y:rest
+c _ = Nothing
 
 -- Wxy  ~> xyy
-w (x:y:rest) = x:y:y:rest
+w (x:y:rest) = Just $ x:y:y:rest
+w _ = Nothing
 
 -- Mx   ~> xx
-m (x:rest) = x:x:rest
+m (x:rest) = Just $ x:x:rest
+m _ = Nothing
 
 -- Txy  ~> yx
-t (x:y:rest) = y:x:rest
+t (x:y:rest) = Just $ y:x:rest
+t _ = Nothing

integral-approx.hs

@@ -0,0 +1,7 @@
+f x = x * x
+
+n = 10000
+
+integral f startX endX =
+  sum $ map (\x -> squareLen * f (x * squareLen)) [1..n]
+  where squareLen = (endX - startX) / n