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 */
M_raw_sin(M_APM rr,int places,M_APM xx)63 void M_raw_sin(M_APM rr, int places, M_APM xx)
64 {
65 M_APM sum, term, tmp2, tmp7, tmp8;
66 int tolerance, flag, local_precision, dplaces;
67 long m1, m2;
68
69 sum = M_get_stack_var();
70 term = M_get_stack_var();
71 tmp2 = M_get_stack_var();
72 tmp7 = M_get_stack_var();
73 tmp8 = M_get_stack_var();
74
75 m_apm_copy(sum, xx);
76 m_apm_copy(term, xx);
77 m_apm_multiply(tmp8, xx, xx);
78 m_apm_round(tmp2, (places + 6), tmp8);
79
80 dplaces = (places + 8) - xx->m_apm_exponent;
81 tolerance = xx->m_apm_exponent - (places + 4);
82
83 m1 = 2L;
84 flag = 0;
85
86 while (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
118 m_apm_round(rr, places, sum);
119 M_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 */
M_raw_cos(M_APM rr,int places,M_APM xx)127 void M_raw_cos(M_APM rr, int places, M_APM xx)
128 {
129 M_APM sum, term, tmp7, tmp8, tmp9;
130 int tolerance, flag, local_precision, prev_exp;
131 long m1, m2;
132
133 sum = M_get_stack_var();
134 term = M_get_stack_var();
135 tmp7 = M_get_stack_var();
136 tmp8 = M_get_stack_var();
137 tmp9 = M_get_stack_var();
138
139 m_apm_copy(sum, MM_One);
140 m_apm_copy(term, MM_One);
141
142 m_apm_multiply(tmp8, xx, xx);
143 m_apm_round(tmp9, (places + 6), tmp8);
144
145 local_precision = places + 8;
146 tolerance = -(places + 4);
147 prev_exp = 0;
148
149 m1 = 1L;
150 flag = 0;
151
152 while (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
188 m_apm_round(rr, places, sum);
189 M_restore_stack(5);
190 }
191 /****************************************************************************/
192