1
2/*
3 *  M_APM  -  mapmrsin.c
4 *
5 *  Copyright (C) 1999 - 2007   Michael C. Ring
6 *
7 *  Permission to use, copy, and distribute this software and its
8 *  documentation for any purpose with or without fee is hereby granted,
9 *  provided that the above copyright notice appear in all copies and
10 *  that both that copyright notice and this permission notice appear
11 *  in supporting documentation.
12 *
13 *  Permission to modify the software is granted. Permission to distribute
14 *  the modified code is granted. Modifications are to be distributed by
15 *  using the file 'license.txt' as a template to modify the file header.
16 *  'license.txt' is available in the official MAPM distribution.
17 *
18 *  This software is provided "as is" without express or implied warranty.
19 */
20
21/*
22 *      $Id: mapmrsin.c,v 1.8 2007/12/03 01:57:00 mike Exp $
23 *
24 *      This file contains the basic series expansion functions for
25 *	the SIN / COS functions.
26 *
27 *      $Log: mapmrsin.c,v $
28 *      Revision 1.8  2007/12/03 01:57:00  mike
29 *      Update license
30 *
31 *      Revision 1.7  2003/06/08 18:22:09  mike
32 *      optimize the raw sin algorithm
33 *
34 *      Revision 1.6  2002/11/03 21:58:27  mike
35 *      Updated function parameters to use the modern style
36 *
37 *      Revision 1.5  2001/07/10 22:14:43  mike
38 *      optimize raw_sin & cos by using fewer digits
39 *      as subsequent terms get smaller
40 *
41 *      Revision 1.4  2000/03/30 21:53:48  mike
42 *      change compare to terminate series expansion using ints instead
43 *      of MAPM numbers
44 *
45 *      Revision 1.3  1999/06/20 16:23:10  mike
46 *      changed local static variables to MAPM stack variables
47 *
48 *      Revision 1.2  1999/05/12 21:06:36  mike
49 *      changed global var names
50 *
51 *      Revision 1.1  1999/05/10 20:56:31  mike
52 *      Initial revision
53 */
54
55#include "m_apm_lc.h"
56
57/****************************************************************************/
58/*
59                                  x^3     x^5     x^7     x^9
60		 sin(x)  =  x  -  ---  +  ---  -  ---  +  ---  ...
61                                   3!      5!      7!      9!
62*/
63void	M_raw_sin(M_APM rr, int places, M_APM xx)
64{
65M_APM	sum, term, tmp2, tmp7, tmp8;
66int     tolerance, flag, local_precision, dplaces;
67long	m1, m2;
68
69sum  = M_get_stack_var();
70term = M_get_stack_var();
71tmp2 = M_get_stack_var();
72tmp7 = M_get_stack_var();
73tmp8 = M_get_stack_var();
74
75m_apm_copy(sum, xx);
76m_apm_copy(term, xx);
77m_apm_multiply(tmp8, xx, xx);
78m_apm_round(tmp2, (places + 6), tmp8);
79
80dplaces   = (places + 8) - xx->m_apm_exponent;
81tolerance = xx->m_apm_exponent - (places + 4);
82
83m1   = 2L;
84flag = 0;
85
86while (TRUE)
87  {
88   m_apm_multiply(tmp8, term, tmp2);
89
90   if ((tmp8->m_apm_exponent < tolerance) || (tmp8->m_apm_sign == 0))
91     break;
92
93   local_precision = dplaces + term->m_apm_exponent;
94
95   if (local_precision < 20)
96     local_precision = 20;
97
98   m2 = m1 * (m1 + 1);
99   m_apm_set_long(tmp7, m2);
100
101   m_apm_divide(term, local_precision, tmp8, tmp7);
102
103   if (flag == 0)
104     {
105      m_apm_subtract(tmp7, sum, term);
106      m_apm_copy(sum, tmp7);
107     }
108   else
109     {
110      m_apm_add(tmp7, sum, term);
111      m_apm_copy(sum, tmp7);
112     }
113
114   m1 += 2;
115   flag = 1 - flag;
116  }
117
118m_apm_round(rr, places, sum);
119M_restore_stack(5);
120}
121/****************************************************************************/
122/*
123                                  x^2     x^4     x^6     x^8
124		 cos(x)  =  1  -  ---  +  ---  -  ---  +  ---  ...
125                                   2!      4!      6!      8!
126*/
127void	M_raw_cos(M_APM rr, int places, M_APM xx)
128{
129M_APM	sum, term, tmp7, tmp8, tmp9;
130int     tolerance, flag, local_precision, prev_exp;
131long	m1, m2;
132
133sum  = M_get_stack_var();
134term = M_get_stack_var();
135tmp7 = M_get_stack_var();
136tmp8 = M_get_stack_var();
137tmp9 = M_get_stack_var();
138
139m_apm_copy(sum, MM_One);
140m_apm_copy(term, MM_One);
141
142m_apm_multiply(tmp8, xx, xx);
143m_apm_round(tmp9, (places + 6), tmp8);
144
145local_precision = places + 8;
146tolerance       = -(places + 4);
147prev_exp        = 0;
148
149m1   = 1L;
150flag = 0;
151
152while (TRUE)
153  {
154   m2 = m1 * (m1 + 1);
155   m_apm_set_long(tmp7, m2);
156
157   m_apm_multiply(tmp8, term, tmp9);
158   m_apm_divide(term, local_precision, tmp8, tmp7);
159
160   if (flag == 0)
161     {
162      m_apm_subtract(tmp7, sum, term);
163      m_apm_copy(sum, tmp7);
164     }
165   else
166     {
167      m_apm_add(tmp7, sum, term);
168      m_apm_copy(sum, tmp7);
169     }
170
171   if ((term->m_apm_exponent < tolerance) || (term->m_apm_sign == 0))
172     break;
173
174   if (m1 != 1L)
175     {
176      local_precision = local_precision + term->m_apm_exponent - prev_exp;
177
178      if (local_precision < 20)
179        local_precision = 20;
180     }
181
182   prev_exp = term->m_apm_exponent;
183
184   m1 += 2;
185   flag = 1 - flag;
186  }
187
188m_apm_round(rr, places, sum);
189M_restore_stack(5);
190}
191/****************************************************************************/
192