1
2 /*
3 * M_APM - mapmasn0.c
4 *
5 * Copyright (C) 2000 - 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: mapmasn0.c,v 1.8 2007/12/03 01:49:49 mike Exp $
23 *
24 * This file contains the 'ARC' family of functions; ARC-SIN,
25 * ARC-COS, ARC-TAN when the input arg is very close to 0 (zero).
26 *
27 * $Log: mapmasn0.c,v $
28 * Revision 1.8 2007/12/03 01:49:49 mike
29 * Update license
30 *
31 * Revision 1.7 2003/06/02 16:51:13 mike
32 * *** empty log message ***
33 *
34 * Revision 1.6 2003/06/02 16:49:48 mike
35 * tweak the decimal places
36 *
37 * Revision 1.5 2003/06/02 16:47:39 mike
38 * tweak arctan algorithm some more
39 *
40 * Revision 1.4 2003/05/31 22:38:07 mike
41 * optimize arctan by using fewer digits as subsequent
42 * terms get smaller
43 *
44 * Revision 1.3 2002/11/03 21:36:43 mike
45 * Updated function parameters to use the modern style
46 *
47 * Revision 1.2 2000/12/02 20:11:37 mike
48 * add comments
49 *
50 * Revision 1.1 2000/12/02 20:08:27 mike
51 * Initial revision
52 */
53
54 #include "m_apm_lc.h"
55
56 /****************************************************************************/
57 /*
58 Calculate arcsin using the identity :
59
60 x
61 arcsin (x) == arctan [ --------------- ]
62 sqrt(1 - x^2)
63
64 */
M_arcsin_near_0(M_APM rr,int places,M_APM aa)65 void M_arcsin_near_0(M_APM rr, int places, M_APM aa)
66 {
67 M_APM tmp5, tmp6;
68
69 tmp5 = M_get_stack_var();
70 tmp6 = M_get_stack_var();
71
72 M_cos_to_sin(tmp5, (places + 8), aa);
73 m_apm_divide(tmp6, (places + 8), aa, tmp5);
74 M_arctan_near_0(rr, places, tmp6);
75
76 M_restore_stack(2);
77 }
78 /****************************************************************************/
79 /*
80 Calculate arccos using the identity :
81
82 arccos (x) == PI / 2 - arcsin (x)
83
84 */
M_arccos_near_0(M_APM rr,int places,M_APM aa)85 void M_arccos_near_0(M_APM rr, int places, M_APM aa)
86 {
87 M_APM tmp1, tmp2;
88
89 tmp1 = M_get_stack_var();
90 tmp2 = M_get_stack_var();
91
92 M_check_PI_places(places);
93 M_arcsin_near_0(tmp1, (places + 4), aa);
94 m_apm_subtract(tmp2, MM_lc_HALF_PI, tmp1);
95 m_apm_round(rr, places, tmp2);
96
97 M_restore_stack(2);
98 }
99 /****************************************************************************/
100 /*
101 calculate arctan (x) with the following series:
102
103 x^3 x^5 x^7 x^9
104 arctan (x) = x - --- + --- - --- + --- ...
105 3 5 7 9
106
107 */
M_arctan_near_0(M_APM rr,int places,M_APM aa)108 void M_arctan_near_0(M_APM rr, int places, M_APM aa)
109 {
110 M_APM tmp0, tmp2, tmpR, tmpS, digit, term;
111 int tolerance, dplaces, local_precision;
112 long m1;
113
114 tmp0 = M_get_stack_var();
115 tmp2 = M_get_stack_var();
116 tmpR = M_get_stack_var();
117 tmpS = M_get_stack_var();
118 term = M_get_stack_var();
119 digit = M_get_stack_var();
120
121 tolerance = aa->m_apm_exponent - (places + 4);
122 dplaces = (places + 8) - aa->m_apm_exponent;
123
124 m_apm_copy(term, aa);
125 m_apm_copy(tmpS, aa);
126 m_apm_multiply(tmp0, aa, aa);
127 m_apm_round(tmp2, (dplaces + 8), tmp0);
128
129 m1 = 1L;
130
131 while (TRUE)
132 {
133 /*
134 * do the subtraction term
135 */
136
137 m_apm_multiply(tmp0, term, tmp2);
138
139 if ((tmp0->m_apm_exponent < tolerance) || (tmp0->m_apm_sign == 0))
140 {
141 m_apm_round(rr, places, tmpS);
142 break;
143 }
144
145 local_precision = dplaces + tmp0->m_apm_exponent;
146
147 if (local_precision < 20)
148 local_precision = 20;
149
150 m1 += 2;
151 m_apm_set_long(digit, m1);
152 m_apm_round(term, local_precision, tmp0);
153 m_apm_divide(tmp0, local_precision, term, digit);
154 m_apm_subtract(tmpR, tmpS, tmp0);
155
156 /*
157 * do the addition term
158 */
159
160 m_apm_multiply(tmp0, term, tmp2);
161
162 if ((tmp0->m_apm_exponent < tolerance) || (tmp0->m_apm_sign == 0))
163 {
164 m_apm_round(rr, places, tmpR);
165 break;
166 }
167
168 local_precision = dplaces + tmp0->m_apm_exponent;
169
170 if (local_precision < 20)
171 local_precision = 20;
172
173 m1 += 2;
174 m_apm_set_long(digit, m1);
175 m_apm_round(term, local_precision, tmp0);
176 m_apm_divide(tmp0, local_precision, term, digit);
177 m_apm_add(tmpS, tmpR, tmp0);
178 }
179
180 M_restore_stack(6); /* restore the 6 locals we used here */
181 }
182 /****************************************************************************/
183