Ticket #3156: poly.hs

File poly.hs, 2.5 KB (added by zachk, 5 years ago)
Line 
1{-AUTHOR ZachK-}
2{-LICENSE BSD3-}
3
4{-This is a polynomial tester-}
5
6module Main where
7
8import IO
9import System
10import Data.List
11import Control.Monad
12
13fact 0=1 --factorial
14fact m@(n+1)=m*fact n
15
16t 0 _ _ = 0
17t a p x = a*x^p --polynomial term
18
19summa _ 0 = 0 --sums k^p over 1..n
20summa p n = sum.map (^p) $ enumFromTo 1 n
21
22summo p=(map (summa p) [0..] !!!) --memoizes the sum function (: HOPEFULLY :) 
23
24m = enumFromTo 0 --sets up polynomial coeffecients
25m2 n = [2,4..n] --gives us even numbers
26
27ri string = read string::Integer
28
29(!!!)::[a]->Integer->a
30(!!!) list n=list!!(fromInteger n) 
31
32makePoly :: Integer -> [([(Integer,Integer)],Integer)]
33makePoly maxP = do
34 let nMax = 10 * maxP
35 let dMax = 10 * maxP --fact (maxP+1) --it is too large of an upper bound
36 a  <- genList nMax (maxP + 2)  [] --generates list of coeffecients
37 let xs = [1..5] --different numbers to test the poly on
38 base <- m2 dMax --the denominator of the poly
39 let terms=map (\u->(a!!!(u-1),u)) $ reverse $ enumFromTo 1 (maxP+1) --generate tuples of the monomials in the polynomial
40 let poly y=sum $ map (\u->t (fst u) (snd u) y) $ terms --turn the tuples into a function
41 guard ((sum a==base)  &&& (foldl' (&&&) True $ map (\x->poly x==base*summo maxP x) xs)) --check
42 return (terms,base)
43 
44
45genList _ 0 as = return (as) --generates list of coeffecients for the polynomial
46genList max n as = do
47 a <- enumFromTo 0 max
48 genList max (n-1) (a:as)
49
50(&&&) True True = True --shortcircuit AND operator
51(&&&) False  _  = False
52(&&&) _  False  = False
53
54main=do
55 maxP <- fmap ri $ fmap head $ getArgs
56 putStrLn "Testing polynomials..."
57 run maxP
58 putStrLn "Done."
59
60run :: Integer -> IO () 
61run maxP = do 
62 let poly = head $ makePoly maxP --take first poly that solves and output it
63 print poly
64 putStrLn $ polyShow poly
65
66var="n" --what we use to represent the variable of the polynomial in the output
67
68termShow (a,p) = --shows single term of polynomial
69 if a==0 
70  then "" 
71  else if a==0 
72   then "" 
73   else (if a==1 
74    then "+" 
75    else "+" ++ (show a)) ++ var ++ (if p==1 
76     then "" 
77     else ("^" ++ (show p)))
78
79--shows the polynomial in all its glory
80polyShow (terms,base) = "(" ++ (pshow terms []) ++ ")/" ++ (show base)
81
82pshow [] a = a --this should never be called unless you pass in a null list
83pshow (t:[]) a = (drop 1 a) ++ (termShow t) --drop the leading +
84pshow (t:t2:ts) a = pshow (t2:ts) (a ++ (termShow t)) 
85--recurses and builds up a string representing the polynomial
86