1 2 /* 3 * M_APM - mapm_div.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: mapm_div.c,v 1.12 2007/12/03 01:35:18 mike Exp $ 23 * 24 * This file contains the basic division functions 25 * 26 * $Log: mapm_div.c,v $ 27 * Revision 1.12 2007/12/03 01:35:18 mike 28 * Update license 29 * 30 * Revision 1.11 2003/07/21 20:07:13 mike 31 * Modify error messages to be in a consistent format. 32 * 33 * Revision 1.10 2003/03/31 22:09:18 mike 34 * call generic error handling function 35 * 36 * Revision 1.9 2002/11/03 22:06:50 mike 37 * Updated function parameters to use the modern style 38 * 39 * Revision 1.8 2001/07/16 19:03:22 mike 40 * add function M_free_all_div 41 * 42 * Revision 1.7 2001/02/11 22:30:42 mike 43 * modify parameters to REALLOC 44 * 45 * Revision 1.6 2000/09/23 19:07:17 mike 46 * change _divide to M_apm_sdivide function name 47 * 48 * Revision 1.5 2000/04/11 18:38:55 mike 49 * use new algorithm to determine q-hat. uses more digits of 50 * the numerator and denominator. 51 * 52 * Revision 1.4 2000/02/03 22:45:08 mike 53 * use MAPM_* generic memory function 54 * 55 * Revision 1.3 1999/06/23 01:10:49 mike 56 * use predefined constant for '15' 57 * 58 * Revision 1.2 1999/06/23 00:55:09 mike 59 * change mult factor to 15 60 * 61 * Revision 1.1 1999/05/10 20:56:31 mike 62 * Initial revision 63 */ 64 65 #include "m_apm_lc.h" 66 67 static M_APM M_div_worka; 68 static M_APM M_div_workb; 69 static M_APM M_div_tmp7; 70 static M_APM M_div_tmp8; 71 static M_APM M_div_tmp9; 72 73 static int M_div_firsttime = TRUE; 74 75 /****************************************************************************/ 76 void M_free_all_div() 77 { 78 if (M_div_firsttime == FALSE) 79 { 80 m_apm_free(M_div_worka); 81 m_apm_free(M_div_workb); 82 m_apm_free(M_div_tmp7); 83 m_apm_free(M_div_tmp8); 84 m_apm_free(M_div_tmp9); 85 86 M_div_firsttime = TRUE; 87 } 88 } 89 /****************************************************************************/ 90 void m_apm_integer_div_rem(M_APM qq, M_APM rr, M_APM aa, M_APM bb) 91 { 92 m_apm_integer_divide(qq, aa, bb); 93 m_apm_multiply(M_div_tmp7, qq, bb); 94 m_apm_subtract(rr, aa, M_div_tmp7); 95 } 96 /****************************************************************************/ 97 void m_apm_integer_divide(M_APM rr, M_APM aa, M_APM bb) 98 { 99 /* 100 * we must use this divide function since the 101 * faster divide function using the reciprocal 102 * will round the result (possibly changing 103 * nnm.999999... --> nn(m+1).0000 which would 104 * invalidate the 'integer_divide' goal). 105 */ 106 107 M_apm_sdivide(rr, 4, aa, bb); 108 109 if (rr->m_apm_exponent <= 0) /* result is 0 */ 110 { 111 M_set_to_zero(rr); 112 } 113 else 114 { 115 if (rr->m_apm_datalength > rr->m_apm_exponent) 116 { 117 rr->m_apm_datalength = rr->m_apm_exponent; 118 M_apm_normalize(rr); 119 } 120 } 121 } 122 /****************************************************************************/ 123 void M_apm_sdivide(M_APM r, int places, M_APM a, M_APM b) 124 { 125 int j, k, m, b0, sign, nexp, indexr, icompare, iterations; 126 long trial_numer; 127 void *vp; 128 129 if (M_div_firsttime) 130 { 131 M_div_firsttime = FALSE; 132 133 M_div_worka = m_apm_init(); 134 M_div_workb = m_apm_init(); 135 M_div_tmp7 = m_apm_init(); 136 M_div_tmp8 = m_apm_init(); 137 M_div_tmp9 = m_apm_init(); 138 } 139 140 sign = a->m_apm_sign * b->m_apm_sign; 141 142 if (sign == 0) /* one number is zero, result is zero */ 143 { 144 if (b->m_apm_sign == 0) 145 { 146 M_apm_log_error_msg(M_APM_RETURN, "\'M_apm_sdivide\', Divide by 0"); 147 } 148 149 M_set_to_zero(r); 150 return; 151 } 152 153 /* 154 * Knuth step D1. Since base = 100, base / 2 = 50. 155 * (also make the working copies positive) 156 */ 157 158 if (b->m_apm_data[0] >= 50) 159 { 160 m_apm_absolute_value(M_div_worka, a); 161 m_apm_absolute_value(M_div_workb, b); 162 } 163 else /* 'normal' step D1 */ 164 { 165 k = 100 / (b->m_apm_data[0] + 1); 166 m_apm_set_long(M_div_tmp9, (long)k); 167 168 m_apm_multiply(M_div_worka, M_div_tmp9, a); 169 m_apm_multiply(M_div_workb, M_div_tmp9, b); 170 171 M_div_worka->m_apm_sign = 1; 172 M_div_workb->m_apm_sign = 1; 173 } 174 175 /* setup trial denominator for step D3 */ 176 177 b0 = 100 * (int)M_div_workb->m_apm_data[0]; 178 179 if (M_div_workb->m_apm_datalength >= 3) 180 b0 += M_div_workb->m_apm_data[1]; 181 182 nexp = M_div_worka->m_apm_exponent - M_div_workb->m_apm_exponent; 183 184 if (nexp > 0) 185 iterations = nexp + places + 1; 186 else 187 iterations = places + 1; 188 189 k = (iterations + 1) >> 1; /* required size of result, in bytes */ 190 191 if (k > r->m_apm_malloclength) 192 { 193 if ((vp = MAPM_REALLOC(r->m_apm_data, (k + 32))) == NULL) 194 { 195 /* fatal, this does not return */ 196 197 M_apm_log_error_msg(M_APM_FATAL, "\'M_apm_sdivide\', Out of memory"); 198 } 199 200 r->m_apm_malloclength = k + 28; 201 r->m_apm_data = (UCHAR *)vp; 202 } 203 204 /* clear the exponent in the working copies */ 205 206 M_div_worka->m_apm_exponent = 0; 207 M_div_workb->m_apm_exponent = 0; 208 209 /* if numbers are equal, ratio == 1.00000... */ 210 211 if ((icompare = m_apm_compare(M_div_worka, M_div_workb)) == 0) 212 { 213 iterations = 1; 214 r->m_apm_data[0] = 10; 215 nexp++; 216 } 217 else /* ratio not 1, do the real division */ 218 { 219 if (icompare == 1) /* numerator > denominator */ 220 { 221 nexp++; /* to adjust the final exponent */ 222 M_div_worka->m_apm_exponent += 1; /* multiply numerator by 10 */ 223 } 224 else /* numerator < denominator */ 225 { 226 M_div_worka->m_apm_exponent += 2; /* multiply numerator by 100 */ 227 } 228 229 indexr = 0; 230 m = 0; 231 232 while (TRUE) 233 { 234 /* 235 * Knuth step D3. Only use the 3rd -> 6th digits if the number 236 * actually has that many digits. 237 */ 238 239 trial_numer = 10000L * (long)M_div_worka->m_apm_data[0]; 240 241 if (M_div_worka->m_apm_datalength >= 5) 242 { 243 trial_numer += 100 * M_div_worka->m_apm_data[1] 244 + M_div_worka->m_apm_data[2]; 245 } 246 else 247 { 248 if (M_div_worka->m_apm_datalength >= 3) 249 trial_numer += 100 * M_div_worka->m_apm_data[1]; 250 } 251 252 j = (int)(trial_numer / b0); 253 254 /* 255 * Since the library 'normalizes' all the results, we need 256 * to look at the exponent of the number to decide if we 257 * have a lead in 0n or 00. 258 */ 259 260 if ((k = 2 - M_div_worka->m_apm_exponent) > 0) 261 { 262 while (TRUE) 263 { 264 j /= 10; 265 if (--k == 0) 266 break; 267 } 268 } 269 270 if (j == 100) /* qhat == base ?? */ 271 j = 99; /* if so, decrease by 1 */ 272 273 m_apm_set_long(M_div_tmp8, (long)j); 274 m_apm_multiply(M_div_tmp7, M_div_tmp8, M_div_workb); 275 276 /* 277 * Compare our q-hat (j) against the desired number. 278 * j is either correct, 1 too large, or 2 too large 279 * per Theorem B on pg 272 of Art of Compter Programming, 280 * Volume 2, 3rd Edition. 281 * 282 * The above statement is only true if using the 2 leading 283 * digits of the numerator and the leading digit of the 284 * denominator. Since we are using the (3) leading digits 285 * of the numerator and the (2) leading digits of the 286 * denominator, we eliminate the case where our q-hat is 287 * 2 too large, (and q-hat being 1 too large is quite remote). 288 */ 289 290 if (m_apm_compare(M_div_tmp7, M_div_worka) == 1) 291 { 292 j--; 293 m_apm_subtract(M_div_tmp8, M_div_tmp7, M_div_workb); 294 m_apm_copy(M_div_tmp7, M_div_tmp8); 295 } 296 297 /* 298 * Since we know q-hat is correct, step D6 is unnecessary. 299 * 300 * Store q-hat, step D5. Since D6 is unnecessary, we can 301 * do D5 before D4 and decide if we are done. 302 */ 303 304 r->m_apm_data[indexr++] = (UCHAR)j; /* j == 'qhat' */ 305 m += 2; 306 307 if (m >= iterations) 308 break; 309 310 /* step D4 */ 311 312 m_apm_subtract(M_div_tmp9, M_div_worka, M_div_tmp7); 313 314 /* 315 * if the subtraction yields zero, the division is exact 316 * and we are done early. 317 */ 318 319 if (M_div_tmp9->m_apm_sign == 0) 320 { 321 iterations = m; 322 break; 323 } 324 325 /* multiply by 100 and re-save */ 326 M_div_tmp9->m_apm_exponent += 2; 327 m_apm_copy(M_div_worka, M_div_tmp9); 328 } 329 } 330 331 r->m_apm_sign = sign; 332 r->m_apm_exponent = nexp; 333 r->m_apm_datalength = iterations; 334 335 M_apm_normalize(r); 336 } 337 /****************************************************************************/ 338