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