5
\$\begingroup\$

I'm very new to Haskell as was hoping to get some feedback on my code AND I have some specific questions. I've posted code below or you can see it here.

I'd welcome ideas on how better to calculate the Maclaurin series. Ideas on how better to do the math itself are welcome but I had in mind what Haskell ideas am I missing out on.

My specific questions are:

  1. Are (almost) all the calculations for Factorial going to be done every time for each value of my Maclaurin series? Or is there some sort of caching going on? One can appreciate that calculating 100! from scratch seems crazy if you've just done 99!

  2. Notice the signature at end: maclaurin :: ( Floating b) => Int-> b -> b I tried using maclaurin :: ( Integral a, Floating b) => a-> b -> b , thinking that would be best. Didn't work. I stumbled on this scheme by removing signature and typing :t maclaurin in WinGHci. I'm surprised it worked! I thought we had to declare each parameter?

--import Data.List(genericTake)

 -- MacClaurin series for Sin(x) is f(0) + f'(0)(x)/1! + f''(0)(x)^2 / 2!  + 
  ... + f^n (0) x^n / n!
  -- or in the case of Sin(x):  0 + 1x + 0 + -1x^3 / 3! + 0 + 1 x^5 / 5! + 0 
  + -1x^7/7! ...
   -- See Wikipedia https://en.wikipedia.org/wiki/Taylor_series




-- return list of alterating 0 and 1
    -- not needed for maclaurin or macResult
ymod2 :: Integral a =>  [a]
ymod2  = map (\y->(y `mod` 2)) [0..]    

factorial :: (Integral a) => a -> a  
factorial 0 = 1  
factorial n = n * factorial (n - 1)



powersMac :: (Integral a) => a -> a
powersMac n
    | n == 0 = 0
    | n == 1 = 1
    | ( (n `mod` 2) == 0 ) = 0
    | otherwise =  ( (-1) * powersMac (n-2) )

-- return list of 0, 1, 0, -1, 0, 1, 0, -1 ...
powersMacList :: (Integral a) => [a]
powersMacList  = map (\y->(powersMac y)) [0..]


-- return list of x^0, x^1, x^2, ... x^n
powersOfx :: (Enum b, Floating b) =>  b -> [b]
powersOfx  x = map (\y->(x^^y)) [0..]

-- return list of 0, 1, 0, -1, 0, ... to n 
firstLst :: Integral a =>  [a]
firstLst   = zipWith (*) (ymod2 ) (powersMacList)

-- return list of 0x, 1x, 0x^2, -1x^3, ... to n
secondLst :: (Enum b, Floating b) =>  b -> [b]
secondLst  x = zipWith (*) (map fromIntegral(powersMacList )) (powersOfx  x)

-- return list of 0, 1x, 0 , -1x^3 / 3! , 0 , 1 x^5 / 5! , 0 ... to n
thirdLst :: ( Enum b, Floating b) =>  b -> [b]
thirdLst  x = zipWith (/) (secondLst  x) (map fromIntegral ( (map(\y->(factorial y) ) [0..])) )

-- sum the list to get result
macResult :: (  Enum b, Floating b) =>  Int -> b -> b
macResult  n x = sum (take n (thirdLst  x))

maclaurin :: (  Floating b) => Int-> b -> b
maclaurin  n x = sum ( take n (zipWith (*) (map(\y->(x^^y))[0..])(zipWith (/) (map fromIntegral (powersMacList )) (map fromIntegral (map(\y->(factorial y) )[0..]) )) ))
\$\endgroup\$
3
\$\begingroup\$

There's a few builtins you may have missed:

ymod2 = cycle [0,1]

If I read powersMac correctly, it basically goes along [0, 1, 0, -1] in cycles as well. So you can reformulate it as:

powersMac n = [0, 1, 0, -1] !! (n `mod` 4)

Considering that you only directly use it as powersMacList, we can again reformulate with cycle:

powersMacList = cycle [0, 1, 0, -1]

BTW: While we're on these two methods:

map (\y -> f y)

is equivalent to:

map f

Noticing that we can rewrite powersOfx as:

powersOfx = map (x^^) [0..]

firstLst does nothing. When there's a 0 in the one list, there's always a 0 in the other list as well.

The formatting in secondLst is somewhat confusing and difficult to read. Consider:

secondLst x = zipWith (*) (map fromIntegral powersMacList) (powersOfx x)

When you use factorial and run a map factorial [0..] you actually want the list of factorials. I'm stealing this excellent SO answer:

factorials = 1 : zipWith (*) factorials [1..]

which turns thirdList into:

thirdList x = zipWith (/) (secondLst x) factorials

applying everything we learned to maclaurin after discarding the fromIntegral calls (that are most likely a result of overspecific type signatures):

maclaurin n x = sum (take n (zipWith (*) (map (x^^) [0..]) (zipWith (/) powersMacList factorial)

Now putting it all together again, after giving out a bit more generous names than firstLst, secondLst and thirdLst:

factorials = 1 : zipWith (*) factorials [1..]

sinZeroDerivations = cycle [0,1,0,-1]
macRawElements x = zipWith (*) sinZeroDerivations (map (x^^) [0..])
macDiscounted x = zipWith (*) (macRawElements x) factorials

macResult n x = sum (take n (macDiscounted x))

Now if you wanted to be less maintainable you could inline the definitions again, but that's not really necessary, is it?

There's something about macDiscounted that annoys me a bit, namely that it only passes x through to macRawElements. We can fix that by swapping those two around like:

sinMacFactors = zipWith (/) sinZeroDerivations factorials

That allows us to specify a bit cleaner what we want:

macResult n x = sum (zipWith (*) (take n (map (x^^) [0..])) sinMacFactors)

notice the second argument to zipWith. That wants to be written as:

macResult n x = sum (zipWith (*) (map (x^^) [0..n]) sinMacFactors)
\$\endgroup\$
  • \$\begingroup\$ Thanks Vogel612! I'm still digesting your detailed answer, but I don't quite understand the factorial improvement (I did read the link and had trouble with that too!) zipWith zips two lists. But isn't the list factorials and [1..] going to be different lengths? It works, so it's a problem with my understanding! Can you explain a little more how this works or give a reference? Also, any pointers on my second question about the type signature for maclaurin? I suspect I don't understand type signatures that well. Fun stuff though! \$\endgroup\$ – Dave May 17 '18 at 16:58
  • \$\begingroup\$ note that factorials references itself in it's definition. Since haskell is lazily evaluated, this means that factorials computes itself on-the-fly (which is the best equivalent to caching a function, btw.). Both factorials and [1..] are infinitely long. the result of zipWith-ing them together also is infinitely long. It's not calculated though. \$\endgroup\$ – Vogel612 May 17 '18 at 17:02
  • 1
    \$\begingroup\$ When using the "final" code in my answer, ghci gives me the typeSignature macResult :: (Enum a, Fractional a, Integral b) => b -> a -> a. Note the difference between Fractional and Floating. Fractionals are all rational numbers \$\mathbb{R}\$, and Floating only encompasses binary floating point numbers (Float, Double, ...) At the start it's much easier to have haskell infer the types for you and just run with that. It's great that you're thinking about types, but the compiler is really good at inferring most cases. \$\endgroup\$ – Vogel612 May 17 '18 at 17:07
  • \$\begingroup\$ Thanks. I think my best bet on the factorials is to go think about it for a while. My brain is in knots right now! I'm trying to think about "take 3 factorials" and writing down the steps, but it doesn't seem like the lists are ever the same length. Point taken about the "types". \$\endgroup\$ – Dave May 17 '18 at 18:29
2
\$\begingroup\$

I am not a Haskell expert. Few words from a numerical analyst point of view.

Your concern about recalculating factorial is very well founded. I recommend to take one more step and realize that factorials are not needed here at all. Rewrite

\$ x - \dfrac{x^3}{3!} + \dfrac{x^5}{5!} - \dfrac{x^7}{7!}\$

as

\$x(1 - \dfrac{x^2}{2*3}(1 - \dfrac{x^2}{4*5}(1 - \dfrac{x^2}{6*7})))\$

(if you want more terms, the pattern shall be clear), and compute the result inside out.


If you are not required to go Maclauren path, considers a natural recursion for a sine from the half angle formula:

\$ \sin x = 2 \sin\frac{x}{2} \cos\frac{x}{2}\$

which combined with a cosine half angle formula:

\$\cos x = \cos^2\frac{x}{2} - \sin^2\frac{x}{2}\$

leads to (x > 0 obviously)

sincos :: Double -> Double -> (Double, Double)

sincos x eps =
    if x < eps then (x, 1 - x*x/2)
    else (2 * s * c, c * c - s * s)
        where (s, c) = sincos (x/2) eps
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.