Conference rusure::math

    Re .0:
    
    That program does not calculate the sine of multiples of 2*pi.  It
    computes multiples of numbers that are _approximately_ 2*pi, since 2*pi
    cannot be completely represented in the machine representation.  As the
    number gets larger, you lose more bits off the end, and your
    approximation gets further away from multiples of 2*pi.
    
    To use an example, suppose the machine had only three decimal digits,
    so your approximation were 6.28 -- you are .0032 away from 2*pi, so you
    get the sine of -.0032.  Then when you try 62.8, you are .032 off.
    
    One way to normalize this is by taking 62.8, dividing by 2*pi, getting
    something around 10, computing 10*2*pi, and subtracting that from 62.8.
    Since the machine only has three digits, when it computes 20*pi, it
    gets 62.8, and when it subtracts 62.8 from 62.8, it gets zero. 
    
    This is wrong.  The sine of 62.8 is not the sine of zero.
    
    The normalization routines in the library take 62.8 and do some work to
    subtract 20*pi as exactly as they can, getting the correct answer of
    -.032.
    
    The math routines correctly calculate the sine of the number you give
    them with good accuracy given the limitations of the machine.  The
    numbers you are getting back are very close to the sines of the numbers
    you are really giving the routine, but the numbers you are giving the
    routine are getting further and further from multiples of 2*pi as the
    numbers get bigger. 
                                                             
    
    				-- edp

    Roundoff error is indeed the problem as stated in .1

    However, there exist library routines designed to minimize this
    problem for situations where you want angles to be geometrically
    accurate.

    See the Runtime Library Routines manual, part II.  The functions
    you want are called MTH$SIND, MTH$COSD, and so on, and take angles
    in degrees, instead of radians.  These are probably what the customer
    will want to use...

    - Jim

    {originally posted in PASCAL, topic 1373 but is better suited for here}
    
    I have an assignment to calculate sine using the approximation formula.
    
                        x - x^3 + x^5 - x^7 ...
              sin(x) = --   ---   ---   ---
                        1!   3!     5!   7!
    
    As X increases (> 5 in my program), the output becomes *MUCH* larger
    than using the Pascal function SIN(X).  Does anyone know the formula or
    algorithm Pascal uses to calculate sine?  Do I have to convert the input 
    to radians (* pi/180) to correct this?
    
    thanks,
    Mike

        That's the correct series when x is given in radians, so if
        your input is in degrees then you would need to convert it
        first. Even for large x that series will converge to six x. 
        But that is assuming "real" arithmetic, i.e., all infinitely
        many decimal places. :-)  With limited precision floating
        point arithmetic on a computer, you have the problem of
        subtracting large individual terms of the sum and being left
        with a remainder which doesn't have very many bits of
        precision left.  So for larger x the result computed the most
        straightforward way starts to drift from the correct answer.
        
	Dan
        

    Well I convert the input to radians and now they start "drifting" from
    the actual SIN(X) results after 121 degrees.  Up to then it's
    relatively close but never exact.

        If you have an input in degrees, then there is a "reduced"
        value between -90 and 90 degrees inclusive such that it and
        the original angle have the same sine.  If you convert that
        redcued angle to radians and plug it into the series, you
        should get a result which is closer to the correct answer.
        
        Dan

    It seems the instructor opened up a hornets nest on this approximation
    formula assignment.  Once you get to the 5th term/level, the factorial of
    9 exceeds MaxInt on the PC's.  I think the VAX handled up to around the
    20th level.
    
    I'd still like to see an algorithm on how the SIN(X) function does it
    in Pascal.
    
    Mike

            <<< Note 762.7 by FRETZ::HEISER "light without heat" >>>
                          -< not a good idea to use >-

<    It seems the instructor opened up a hornets nest on this approximation
<    formula assignment.  Once you get to the 5th term/level, the factorial of
<    9 exceeds MaxInt on the PC's.  I think the VAX handled up to around the
<    20th level.

There's no reason to calculate the factorial function explicitly.  Note that

	sin(x) = sum(i=0..., x**(2*i+1)/(2*i+1)!)

and thus the quotient between term i and term i+1 is

	x**2/((2*i+2)*(2*i+3))

This makes it irrelevant how big the factorial function grows.  Incidentally,
this makes the calculation quicker, too.

 <<< Note 762.8 by AMCCXN::BERGH "Peter Bergh, (719) 592-5036, DTN 592-5036" >>>
                      -< No need to evaluate factorials >-

< There's no reason to calculate the factorial function explicitly.  Note that

<	sin(x) = sum(i=0..., x**(2*i+1)/(2*i+1)!)

Sorry  --  forgot the alternating signs.  The argument against explicit
evaluation of factorials will be the same when using the correct formula.

on systems such as vaxen and ibm mainframes there is a math library.
all of the languages call the sin/cos function.  on a pc or mac, where
each compiler is a complete is a complete kit, the sin/cos function
is included with each compiler.  but for a given vendor, it is likely
that the function is the same in all their language kits.

typically, the first step is to reduce the argument to the range -pi to +pi
for sin, and 0 to pi for cos.  this is harder than it sounds without losing 
some accuracy.  this is just an application of sin(2pi + x) = sin(x).
then, since sin(-x) = -sin(x), they usually remember the sign and take the
absolute value of x.  this insures that the series expansion is an
alternating series, which is a good thing. the error is less than the
last term used.  it also provides a way to tell when to stop.

then they usually look at x.  if x is now small enough, they use the
series you gave, a Maclaurin (sp?) series. For larger values they may
use a Taylor series which is similar.

When evaluating the terms of the series, they never compute x^n or 
n!, instead they multiply the previos term by (x*x)/((n-1)n).
the numerical analyst who wrote the function probably spent a few hours
deciding what order to do those operations: div, mul, div, mul or
maybe mul, div, mul, div or maybe something else.

on vms, the goal of the math library is to be good to the last bit.
on some other systems the goal is to be fast.

the whole area of numerical analysis is much more difficult than
it appears on the surface.  the vms rtl group has a document with
the detailed analysis of all the math functions.