"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/mapmasin.c" (21 Feb 2010, 12646 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 "mapmasin.c" see the
Fossies "Dox" file reference documentation.
1
2 /*
3 * M_APM - mapmasin.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: mapmasin.c,v 1.28 2007/12/03 01:49:10 mike Exp $
23 *
24 * This file contains the 'ARC' family of functions; ARC-SIN, ARC-COS,
25 * ARC-TAN, and ARC-TAN2.
26 *
27 * $Log: mapmasin.c,v $
28 * Revision 1.28 2007/12/03 01:49:10 mike
29 * Update license
30 *
31 * Revision 1.27 2003/07/24 16:34:02 mike
32 * update arctan_large_input
33 *
34 * Revision 1.26 2003/07/21 20:27:48 mike
35 * Modify error messages to be in a consistent format.
36 *
37 * Revision 1.25 2003/07/21 19:19:26 mike
38 * add new arctan with large input value
39 *
40 * Revision 1.24 2003/05/01 21:58:49 mike
41 * remove math.h
42 *
43 * Revision 1.23 2003/04/09 21:43:00 mike
44 * optimize iterative asin & acos with lessons learned
45 * from the new log function
46 *
47 * Revision 1.22 2003/03/31 21:58:11 mike
48 * call generic error handling function
49 *
50 * Revision 1.21 2002/11/03 21:41:54 mike
51 * Updated function parameters to use the modern style
52 *
53 * Revision 1.20 2001/02/07 19:07:07 mike
54 * eliminate MM_skip_limit_PI_check
55 *
56 * Revision 1.19 2001/02/06 21:50:56 mike
57 * don't display accuracy when iteration count maxes out
58 *
59 * Revision 1.18 2000/12/02 20:10:09 mike
60 * add calls to more efficient functions if
61 * the input args are close to zero
62 *
63 * Revision 1.17 2000/09/05 22:18:02 mike
64 * re-arrange code to eliminate goto from atan2
65 *
66 * Revision 1.16 2000/05/28 23:58:41 mike
67 * minor optimization to arc-tan2
68 *
69 * Revision 1.15 2000/05/19 17:13:29 mike
70 * use local copies of PI variables & recompute
71 * on the fly as needed
72 *
73 * Revision 1.14 2000/03/27 21:43:23 mike
74 * dtermine how many iterations should be required at
75 * run time for arc-sin and arc-cos
76 *
77 * Revision 1.13 1999/09/21 21:00:33 mike
78 * make sure the sign of 'sin' from M_cos_to_sin is non-zero
79 * before assigning it from the original angle.
80 *
81 * Revision 1.12 1999/07/21 03:05:06 mike
82 * added some comments
83 *
84 * Revision 1.11 1999/07/19 02:33:39 mike
85 * reset local precision again
86 *
87 * Revision 1.10 1999/07/19 02:18:05 mike
88 * more fine tuning of local precision
89 *
90 * Revision 1.9 1999/07/19 00:08:34 mike
91 * adjust local precision during iterative loops
92 *
93 * Revision 1.8 1999/07/18 22:35:56 mike
94 * make arc-sin and arc-cos use dynamically changing
95 * precision to speed up iterative routines for large N
96 *
97 * Revision 1.7 1999/07/09 22:52:00 mike
98 * skip limit PI check when not needed
99 *
100 * Revision 1.6 1999/07/09 00:10:39 mike
101 * use better method for arc sin and arc cos
102 *
103 * Revision 1.5 1999/07/08 22:56:20 mike
104 * replace local MAPM constant with a global
105 *
106 * Revision 1.4 1999/06/20 16:55:01 mike
107 * changed local static variables to MAPM stack variables
108 *
109 * Revision 1.3 1999/05/15 02:10:27 mike
110 * add check for number of decimal places
111 *
112 * Revision 1.2 1999/05/10 21:10:21 mike
113 * added some comments
114 *
115 * Revision 1.1 1999/05/10 20:56:31 mike
116 * Initial revision
117 */
118
119 #include "m_apm_lc.h"
120
121 /****************************************************************************/
122 void m_apm_arctan2(M_APM rr, int places, M_APM yy, M_APM xx)
123 {
124 M_APM tmp5, tmp6, tmp7;
125 int ix, iy;
126
127 iy = yy->m_apm_sign;
128 ix = xx->m_apm_sign;
129
130 if (ix == 0) /* x == 0 */
131 {
132 if (iy == 0) /* y == 0 */
133 {
134 M_apm_log_error_msg(M_APM_RETURN, "\'m_apm_arctan2\', Both Inputs = 0");
135 M_set_to_zero(rr);
136 return;
137 }
138
139 M_check_PI_places(places);
140 m_apm_round(rr, places, MM_lc_HALF_PI);
141 rr->m_apm_sign = iy;
142 return;
143 }
144
145 if (iy == 0)
146 {
147 if (ix == 1)
148 {
149 M_set_to_zero(rr);
150 }
151 else
152 {
153 M_check_PI_places(places);
154 m_apm_round(rr, places, MM_lc_PI);
155 }
156
157 return;
158 }
159
160 /*
161 * the special cases have been handled, now do the real work
162 */
163
164 tmp5 = M_get_stack_var();
165 tmp6 = M_get_stack_var();
166 tmp7 = M_get_stack_var();
167
168 m_apm_divide(tmp6, (places + 6), yy, xx);
169 m_apm_arctan(tmp5, (places + 6), tmp6);
170
171 if (ix == 1) /* 'x' is positive */
172 {
173 m_apm_round(rr, places, tmp5);
174 }
175 else /* 'x' is negative */
176 {
177 M_check_PI_places(places);
178
179 if (iy == 1) /* 'y' is positive */
180 {
181 m_apm_add(tmp7, tmp5, MM_lc_PI);
182 m_apm_round(rr, places, tmp7);
183 }
184 else /* 'y' is negative */
185 {
186 m_apm_subtract(tmp7, tmp5, MM_lc_PI);
187 m_apm_round(rr, places, tmp7);
188 }
189 }
190
191 M_restore_stack(3);
192 }
193 /****************************************************************************/
194 /*
195 Calculate arctan using the identity :
196
197 x
198 arctan (x) == arcsin [ --------------- ]
199 sqrt(1 + x^2)
200
201 */
202 void m_apm_arctan(M_APM rr, int places, M_APM xx)
203 {
204 M_APM tmp8, tmp9;
205
206 if (xx->m_apm_sign == 0) /* input == 0 ?? */
207 {
208 M_set_to_zero(rr);
209 return;
210 }
211
212 if (xx->m_apm_exponent <= -4) /* input close to 0 ?? */
213 {
214 M_arctan_near_0(rr, places, xx);
215 return;
216 }
217
218 if (xx->m_apm_exponent >= 4) /* large input */
219 {
220 M_arctan_large_input(rr, places, xx);
221 return;
222 }
223
224 tmp8 = M_get_stack_var();
225 tmp9 = M_get_stack_var();
226
227 m_apm_multiply(tmp9, xx, xx);
228 m_apm_add(tmp8, tmp9, MM_One);
229 m_apm_sqrt(tmp9, (places + 6), tmp8);
230 m_apm_divide(tmp8, (places + 6), xx, tmp9);
231 m_apm_arcsin(rr, places, tmp8);
232 M_restore_stack(2);
233 }
234 /****************************************************************************/
235 /*
236
237 for large input values use :
238
239 arctan(x) = (PI / 2) - arctan(1 / |x|)
240
241 and sign of result = sign of original input
242
243 */
244 void M_arctan_large_input(M_APM rr, int places, M_APM xx)
245 {
246 M_APM tmp1, tmp2;
247
248 tmp1 = M_get_stack_var();
249 tmp2 = M_get_stack_var();
250
251 M_check_PI_places(places);
252
253 m_apm_divide(tmp1, (places + 6), MM_One, xx); /* 1 / xx */
254 tmp1->m_apm_sign = 1; /* make positive */
255 m_apm_arctan(tmp2, (places + 6), tmp1);
256 m_apm_subtract(tmp1, MM_lc_HALF_PI, tmp2);
257 m_apm_round(rr, places, tmp1);
258
259 rr->m_apm_sign = xx->m_apm_sign; /* fix final sign */
260
261 M_restore_stack(2);
262 }
263 /****************************************************************************/
264 void m_apm_arcsin(M_APM r, int places, M_APM x)
265 {
266 M_APM tmp0, tmp1, tmp2, tmp3, current_x;
267 int ii, maxiter, maxp, tolerance, local_precision;
268
269 current_x = M_get_stack_var();
270 tmp0 = M_get_stack_var();
271 tmp1 = M_get_stack_var();
272 tmp2 = M_get_stack_var();
273 tmp3 = M_get_stack_var();
274
275 m_apm_absolute_value(tmp0, x);
276
277 ii = m_apm_compare(tmp0, MM_One);
278
279 if (ii == 1) /* |x| > 1 */
280 {
281 M_apm_log_error_msg(M_APM_RETURN, "\'m_apm_arcsin\', |Argument| > 1");
282 M_set_to_zero(r);
283 M_restore_stack(5);
284 return;
285 }
286
287 if (ii == 0) /* |x| == 1, arcsin = +/- PI / 2 */
288 {
289 M_check_PI_places(places);
290 m_apm_round(r, places, MM_lc_HALF_PI);
291 r->m_apm_sign = x->m_apm_sign;
292
293 M_restore_stack(5);
294 return;
295 }
296
297 if (m_apm_compare(tmp0, MM_0_85) == 1) /* check if > 0.85 */
298 {
299 M_cos_to_sin(tmp2, (places + 4), x);
300 m_apm_arccos(r, places, tmp2);
301 r->m_apm_sign = x->m_apm_sign;
302
303 M_restore_stack(5);
304 return;
305 }
306
307 if (x->m_apm_sign == 0) /* input == 0 ?? */
308 {
309 M_set_to_zero(r);
310 M_restore_stack(5);
311 return;
312 }
313
314 if (x->m_apm_exponent <= -4) /* input close to 0 ?? */
315 {
316 M_arcsin_near_0(r, places, x);
317 M_restore_stack(5);
318 return;
319 }
320
321 tolerance = -(places + 4);
322 maxp = places + 8 - x->m_apm_exponent;
323 local_precision = 20 - x->m_apm_exponent;
324
325 /*
326 * compute the maximum number of iterations
327 * that should be needed to calculate to
328 * the desired accuracy. [ constant below ~= 1 / log(2) ]
329 */
330
331 maxiter = (int)(log((double)(places + 2)) * 1.442695) + 3;
332
333 if (maxiter < 5)
334 maxiter = 5;
335
336 M_get_asin_guess(current_x, x);
337
338 /* Use the following iteration to solve for arc-sin :
339
340 sin(X) - N
341 X = X - ------------
342 n+1 cos(X)
343 */
344
345 ii = 0;
346
347 while (TRUE)
348 {
349 M_4x_cos(tmp1, local_precision, current_x);
350
351 M_cos_to_sin(tmp2, local_precision, tmp1);
352 if (tmp2->m_apm_sign != 0)
353 tmp2->m_apm_sign = current_x->m_apm_sign;
354
355 m_apm_subtract(tmp3, tmp2, x);
356 m_apm_divide(tmp0, local_precision, tmp3, tmp1);
357
358 m_apm_subtract(tmp2, current_x, tmp0);
359 m_apm_copy(current_x, tmp2);
360
361 if (ii != 0)
362 {
363 if (((2 * tmp0->m_apm_exponent) < tolerance) || (tmp0->m_apm_sign == 0))
364 break;
365 }
366
367 if (++ii == maxiter)
368 {
369 M_apm_log_error_msg(M_APM_RETURN,
370 "\'m_apm_arcsin\', max iteration count reached");
371 break;
372 }
373
374 local_precision *= 2;
375
376 if (local_precision > maxp)
377 local_precision = maxp;
378 }
379
380 m_apm_round(r, places, current_x);
381 M_restore_stack(5);
382 }
383 /****************************************************************************/
384 void m_apm_arccos(M_APM r, int places, M_APM x)
385 {
386 M_APM tmp0, tmp1, tmp2, tmp3, current_x;
387 int ii, maxiter, maxp, tolerance, local_precision;
388
389 current_x = M_get_stack_var();
390 tmp0 = M_get_stack_var();
391 tmp1 = M_get_stack_var();
392 tmp2 = M_get_stack_var();
393 tmp3 = M_get_stack_var();
394
395 m_apm_absolute_value(tmp0, x);
396
397 ii = m_apm_compare(tmp0, MM_One);
398
399 if (ii == 1) /* |x| > 1 */
400 {
401 M_apm_log_error_msg(M_APM_RETURN, "\'m_apm_arccos\', |Argument| > 1");
402 M_set_to_zero(r);
403 M_restore_stack(5);
404 return;
405 }
406
407 if (ii == 0) /* |x| == 1, arccos = 0, PI */
408 {
409 if (x->m_apm_sign == 1)
410 {
411 M_set_to_zero(r);
412 }
413 else
414 {
415 M_check_PI_places(places);
416 m_apm_round(r, places, MM_lc_PI);
417 }
418
419 M_restore_stack(5);
420 return;
421 }
422
423 if (m_apm_compare(tmp0, MM_0_85) == 1) /* check if > 0.85 */
424 {
425 M_cos_to_sin(tmp2, (places + 4), x);
426
427 if (x->m_apm_sign == 1)
428 {
429 m_apm_arcsin(r, places, tmp2);
430 }
431 else
432 {
433 M_check_PI_places(places);
434 m_apm_arcsin(tmp3, (places + 4), tmp2);
435 m_apm_subtract(tmp1, MM_lc_PI, tmp3);
436 m_apm_round(r, places, tmp1);
437 }
438
439 M_restore_stack(5);
440 return;
441 }
442
443 if (x->m_apm_sign == 0) /* input == 0 ?? */
444 {
445 M_check_PI_places(places);
446 m_apm_round(r, places, MM_lc_HALF_PI);
447 M_restore_stack(5);
448 return;
449 }
450
451 if (x->m_apm_exponent <= -4) /* input close to 0 ?? */
452 {
453 M_arccos_near_0(r, places, x);
454 M_restore_stack(5);
455 return;
456 }
457
458 tolerance = -(places + 4);
459 maxp = places + 8;
460 local_precision = 18;
461
462 /*
463 * compute the maximum number of iterations
464 * that should be needed to calculate to
465 * the desired accuracy. [ constant below ~= 1 / log(2) ]
466 */
467
468 maxiter = (int)(log((double)(places + 2)) * 1.442695) + 3;
469
470 if (maxiter < 5)
471 maxiter = 5;
472
473 M_get_acos_guess(current_x, x);
474
475 /* Use the following iteration to solve for arc-cos :
476
477 cos(X) - N
478 X = X + ------------
479 n+1 sin(X)
480 */
481
482 ii = 0;
483
484 while (TRUE)
485 {
486 M_4x_cos(tmp1, local_precision, current_x);
487
488 M_cos_to_sin(tmp2, local_precision, tmp1);
489 if (tmp2->m_apm_sign != 0)
490 tmp2->m_apm_sign = current_x->m_apm_sign;
491
492 m_apm_subtract(tmp3, tmp1, x);
493 m_apm_divide(tmp0, local_precision, tmp3, tmp2);
494
495 m_apm_add(tmp2, current_x, tmp0);
496 m_apm_copy(current_x, tmp2);
497
498 if (ii != 0)
499 {
500 if (((2 * tmp0->m_apm_exponent) < tolerance) || (tmp0->m_apm_sign == 0))
501 break;
502 }
503
504 if (++ii == maxiter)
505 {
506 M_apm_log_error_msg(M_APM_RETURN,
507 "\'m_apm_arccos\', max iteration count reached");
508 break;
509 }
510
511 local_precision *= 2;
512
513 if (local_precision > maxp)
514 local_precision = maxp;
515 }
516
517 m_apm_round(r, places, current_x);
518 M_restore_stack(5);
519 }
520 /****************************************************************************/