"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/mapmrsin.c" (21 Feb 2010, 4809 Bytes) of package /linux/misc/old/mapm-4.9.5a.tar.gz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
For more information about "mapmrsin.c" see the
Fossies "Dox" file reference documentation.
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 */
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 */
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 /****************************************************************************/