lnpi.hs

lnpi.hs

import qualified Data.Map.Strict as M                                                        
import qualified Data.Set        as S                                                        
import Control.Arrow ((&&&))                                                                 
import Control.Monad (guard)                                                                 
import Data.List (findIndices, tails)                                                        
import Data.Maybe (mapMaybe)                                                                 
import Data.Ratio (Ratio, (%), numerator, denominator)                                       
                                                                                             
-------------------------------------------------------------------------------              
                                                                                             
type Val  = Ratio Integer                                                                    
data Fn   = Id | Fact                     deriving (Bounded, Enum)                           
data Op   = Plus | Sub | Mult | Div | Exp deriving (Bounded, Enum)                           
data Expr = V Val | E1 Fn Expr | E2 Op (Expr,Expr)                                           
                                                                                             
instance Show Op where                                                                       
  show Plus = "+"                                                                            
  show Sub  = "-"                                                                            
  show Mult = "*"                                                                            
  show Div  = "/"                                                                            
  show Exp  = "^"                                                                            
                                                                                             
instance Show Expr where                                                                     
  show (V a)          = show $ round a                                                       
  show (E1 Id   e)    = show e                                                               
  show (E1 Fact e)    = show e ++ "!"                                                        
  show (E2 o (e1,e2)) = "(" ++ show e1 ++ show o ++ show e2 ++ ")"                           
                                                                                             
-------------------------------------------------------------------------------              
                                                                                             
functions  = [minBound..maxBound] :: [Fn]                                                    
operations = [minBound..maxBound] :: [Op]                                                    
                                                                                             
facts :: [Val]                                                                               
facts = map (\n -> product [1..n]) [1..100]                                                  
                                                                                             
isInt :: (Num a, Eq a) => Ratio a -> Bool                                                    
isInt x = 1 == denominator x                                                                 
                                                                                             
mkPartitions :: [a] -> [([a],[a])]                                                           
mkPartitions xs = map (\n -> (take n &&& drop n) xs) [1..(length xs-1)]                                                                                                                   
-------------------------------------------------------------------------------              
                                                                                             
eval :: Expr -> Maybe Val                                                                    
eval (V a)           = Just a                                                                
eval (E1 f e)        = calc1 f =<< eval e                                                    
eval (E2 o (e1, e2)) = calc2 o =<< (,) <$> eval e1 <*> eval e2                               
                                                                                             
calc1 :: Fn -> Val -> Maybe Val                                                              
calc1 Id   a = Just a                                                                        
calc1 Fact a = guard (isInt a && a<100 && a>0) >> Just f                                     
  where f = facts !! (pred . fromIntegral . numerator) a                                     

calc2 :: Op -> (Val, Val) -> Maybe Val                                                       
calc2 Plus (a,b) = Just $ a + b                                                              
calc2 Sub  (a,b) = Just $ a - b                                                              
calc2 Mult (a,b) = Just $ a * b                                                              
calc2 Div  (a,b) = guard (b/=0) >> Just (a/b)                                                
calc2 Exp  (a,b) = makeExp                                                                   
  where makeExp                                                                              
          | a == 0        = Just 0                                                           
          | a == 1        = Just 1                                                           
          | not (isInt b) = Nothing                                                          
          | abs b > 1023  = Nothing                                                          
          | otherwise     = Just $ a ^^ numerator b                                          
                                                                                             
valuesFrom :: [Val] -> S.Set Val                                                             
valuesFrom range = M.findWithDefault S.empty range                                           
                 $ mkMap (length range)                                                      
  where                                                                                      
    mkMap :: Int -> M.Map [Val] (S.Set Val)                                                  
    mkMap 0 = M.empty                                                                        
    mkMap n = M.union prevMap thisMap                                                        
      where                                                                                  
        prevMap = mkMap $ pred n                                                             
        thisMap = foldr (\sq -> M.insert sq $ mkSet sq) M.empty $ subsequences n             
        mkSet   = S.fromList . mapMaybe eval . exprsFrom prevMap                             
                                                                                             
    exprsFrom :: M.Map [Val] (S.Set Val) -> [Val] -> [Expr]                                  
    exprsFrom prevMap [x] = [ E1 f (V x) | f <- functions ]                                  
    exprsFrom prevMap xs  = [ E1 f $ E2 o (V va, V vb)                                       
                            |        f <- functions                                          
                            ,        o <- operations                                         
                            , (as, bs) <- mkPartitions xs                                    
                            ,       va <- valsFrom as                                        
                            ,       vb <- valsFrom bs                                        
                            ]                                                                
      where valsFrom xs = S.toList $ M.findWithDefault S.empty xs prevMap                    
                                                                                             
    subsequences n = filter ((== n) . length) $ map (take n) $ tails range                   
                                                                                             
--------------------------------------------------------------------------------             
                                                                                             
lnpi :: [Val] -> Integer                                                                     
lnpi = numerator . firstGap . S.toList . S.filter (>0) . S.filter isInt                      
     . valuesFrom                                                                            
  where firstGap :: (Num a, Eq a) => [a] -> a                                                
        firstGap ls  = 1 + ls !! (head . findIndices (/=1) . gaps) ls                        
                                                                                             
        gaps :: Num t => [t] -> [t]                                                          
        gaps [x] = []                                                                        
        gaps (x:xs) = (head xs - x) : gaps xs                                                
                                                                                             
main = print $ lnpi [1..5]