How to calculate the sine function in Haskell?

Question

I need a Haskell function that computes an approximation of the sine of some number, using the associated Taylor series.

In C++ I wrote this:

double msin(double number, int counter = 0, double sum = 0)
{
    // sin(x) = x - (x'3 / 3!) + (x'5 / 5!) - (x'7 / 7!) + (x'9 / 9!)
    if (counter <= 20)
    {
        if (counter % 2 == 0)
            sum += mpow(number, counter * 2 + 1) / mfak(counter * 2 + 1);
        else
            sum -= mpow(number, counter * 2 + 1) / mfak(counter * 2 + 1);

        counter++;
        sum =  msin(number, counter, sum);

        return sum;
    }

    return (sum* 180.0 / _PI);
}

Now I am trying to do it in Haskell and I have no idea how. For now I was trying something like this:

This works:

mfak number = if number < 2
                    then 1
                    else number *( mfak (number -1 )) 


mpow number potenca = if potenca == 0
                        then 0
                        else if potenca == 1
                        then 1
                        else (number * (mpow number (potenca-1)))

This doesn't work:

msin :: Double -> Int -> Double -> Double                           
msin number counter sum = if counter <= 20
                                    then if counter `mod` 2==0
                                            then let sum = sum + (msin 1 (let counter = counter+1 in counter) sum) in sum 
                                            else let sum = sum + (msin 1 (let counter = counter+1 in counter) sum) in sum                                   
                                    else sum* 180.0 / 3.14

Updated, doesn't compile: "Couldn't match expected type 'Double' with actual type 'Int'"

msin :: Double -> Int -> Double -> Double                           
msin number counter sum = if counter <= 20
                                    then if counter `mod` 2==0
                                            then let sum' = sum + ((mpow number (counter*2+1))/(mfak counter*2+1)) in msin number (counter+1) sum'
                                            else let sum' = sum - ((mpow number (counter*2+1))/(mfak counter*2+1)) in msin number (counter+1) sum'                              
                                    else sum* 180.0 / 3.14

As you can see, the biggest problem is how to add something to sum, increase counter and go in recursion again with these new values...

P.S. I am new to Haskell so try to explain your solution as much as you can please. I was reading some tutorials and that, but I can't find how to save the result of some expression into a value and then continue with other code after it. It just returns my value each time I try to do that, and I don't want that.

Answer 1

The expression is of the form x*P(x²).

For maximal efficiency, the polynomial in x² must be evaluated using the Horner rule rather than computing the powers of x² separately.

The coefficient serie with the factorial values can be expressed recursively in Haskell, just like is commonly done for the Fibonacci series. Using the ghci interpreter as our testbed, we have:

$ ghci
GHCi, version 8.8.4: https://www.haskell.org/ghc/  :? for help
 λ> 
 λ> 
 λ> nextCoeffs d c = c : (nextCoeffs (d+1) ((-c)/(fromIntegral $ (2*d+2)*(2*d+3))))
 λ> 
 λ> allCoeffs = nextCoeffs 0 1.0
 λ> 

where d is the depth inside the serie and c the current coefficient.

Sanity check: the coefficient at depth 3 must be the inverse of 7!

 λ> 
 λ> 1.0 /(allCoeffs !! 3)
 -5040.0
 λ> 

The Horner rule can be rendered in Haskell thru the foldr1 :: (a -> a -> a) -> [a] -> a library function.

As is customary in Haskell, I take the liberty to put the term count as the leftmost argument because it is the one most likely to be held constant. This is for currying (partial evaluation) purposes.

So we have:

 λ> :{
|λ>  msin count x = let { s = x*x ;  cs = take count allCoeffs ;
|λ>                       stepFn c acc = acc*s + c ; }
|λ>                 in  x * (foldr1 stepFn cs)
|λ> :}

Sanity checks, taking 20 terms:

 λ> 
 λ> pi
 3.141592653589793
 λ>
 λ> msin 20 (pi/6)
 0.49999999999999994
 λ> 
 λ> msin 20 (pi/2)
 1.0
 λ> 

Side note 1: final multiplication by 180 / π is only of interest for inverse trigonometric functions.

Side note 2: in practice, to get a reasonably fast convergence, one should reduce the input variable x into the [-π,+π] interval using the periodicity of the sine function.

Answer 2

I’ve tried to follow your conventions as much as I could. For mfak and mpow, you should avoid using if as it is clearer to write them using pattern matching :

mfak :: Int -> Int
mfak 0 = 1
mfak 1 = 1
mfak n = n * mfak (n - 1)

mpow :: Double -> Int -> Double
mpow _ 0 = 1
mpow x 1 = x
mpow x p = x * mpow x (p - 1)

Before calculating the sinus, we create a list of coefficients [(sign, power, factorial)] :

x - (x^3 / 3!) + (x^5 / 5!) - (x^7 / 7!) + (x^9 / 9!)
→ [(1,1,1), (-1,3,6), (1,5,120), (-1,7,5040), (1,9,362880)]

The list is created infinite by a list comprehension. First we zip the lists [1,-1,1,-1,1,-1...] and [1,3,5,7,9,11...]. This gives us the list [(1,1), (-1,3), (1,5), (-1,7)...]. From this list, we create the final list [(1,1,1), (-1,3,6), (1,5,120), (-1,7,5040)...]:

sinCoeff :: [(Double, Int, Double)]
sinCoeff = [ (fromIntegral s, i, fromIntegral $ mfak i)
           | (s, i) <- zip (cycle [1, -1]) [1,3..]]

(cycle repeats a list indefinitely, [1,3..] creates an infinite list which starts at 1 with a step of 2)

Finally, the msin function is near the definition. It also uses a list comprehension to achieve its goeal (note that I kept the * 180 / pi though I’m not sure it should be there. Haskell knows pi).

msin :: Int -> Double -> Double
msin n x = 180 * sum [ s * mpow x p / f | (s, p, f) <- take n sinCoeff] / pi

(take n sinCoeff returns the first n elements of a list)

You may try the previous code with the following :

main = do
    print $ take 10 sinCoeff
    print $ msin 5 0.5
    print $ msin 10 0.5

Answer 3

You may also use recursive lists definitions to get [x, x^3, x^5 ...] and [1, 1/3!, 1/5! ...] infinite sequences. When they are done, the rest is to multiply their items each by other and take the sum.

sinus count x = sum (take count $ zipWith (*) ifactorials xpowers)
    where xpowers     = x : map ((x*x)*) xpowers 
          ifactorials = 1 : zipWith (/) ifactorials [i*(i+1) | i <- [2, 4 .. ]]

Also, it would be better to define xpowers = iterate ((x*x)*) x, as it seems to be much more readable.

Answer 4

I would rework the algorithm a bit. First we can define the list of factorial inverses:

factorialInv :: [Double]
factorialInv = scanl (/) 1 [1..]  -- 1/0! , 1/1! , 1/2! , 1/3! , ...

Then, we follow with the sine coefficients:

sineCoefficients :: [Double]
sineCoefficients = 0 : 1 : 0 : -1 : sineCoefficients

Then, given x, we multiply both the above lists with the powers of x, pointwise:

powerSeries :: [Double]   -- ^ Coefficients
            -> Double     -- ^ Point x on which to compute the series
            -> [Double]   -- ^ Series terms
powerSeries cs x = zipWith3 (\a b c -> a * b * c) cs powers factorialInv
   where powers = iterate (*x) 1   -- 1 , x , x^2 , x^3 , ...

Finally, we take the first 20 terms and sum them up.

sine :: Double -> Double                   
sine = sum . take 20 . powerSeries sineCoefficients
 -- i.e., sine x = sum (take 20 (powerSeries sineCoefficients x))

Answer 5

As you can see the biggest problem is how to add something to "vsota",

In a functional language you would use recursion here - the variable vstota is implemented as a function parameter which is passed from call to call as a list is processed.

For example, to sum a list of numbers, we would write something like:

sum xs = go 0 xs
  where go total [] = total
        go total (x:xs) = go (total+x) xs

In an imperative language total would be a variable which gets updated. Here is is a function parameter which gets passed to the next recursive call to go.

In your case, I would first write a function which generates the terms of the power series:

sinusTerms n x = ... -- the first n terms of x - (x'3 / 3!) + (x'5 / 5!) - (x'7 / 7!) ...

and then use the sum function above:

sinus n x = sum (sinusTerms n x)

Answer 6

The problem is expressions like let stevec = stevec+1 in stevec. Haskell is not an imperative language. This does not add one to stevec. Instead it defines stevec to be a number that is one more than itself. No such number exists, thus you will get an infinite loop or, if you are lucky, a crash.

Instead of

stevec++;
vsota =  msin(stevilo, stevec, vsota);

You should use something like

let stevec' = stevec + 1
in  msin stevilo stevec' vsota

or just

msin stevilo (stevec + 1) vsota

(There's also something here that I don't understand. You are going to need mpow and mfak. Where are they?)