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 /****************************************************************************/ 521