1 2 /* 3 * M_APM - mapm_gcd.c 4 * 5 * Copyright (C) 2001 - 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_gcd.c,v 1.8 2007/12/03 01:41:05 mike Exp $ 23 * 24 * This file contains the GCD and LCM functions 25 * 26 * $Log: mapm_gcd.c,v $ 27 * Revision 1.8 2007/12/03 01:41:05 mike 28 * Update license 29 * 30 * Revision 1.7 2003/07/21 20:16:43 mike 31 * Modify error messages to be in a consistent format. 32 * 33 * Revision 1.6 2003/03/31 22:12:33 mike 34 * call generic error handling function 35 * 36 * Revision 1.5 2002/11/03 22:44:21 mike 37 * Updated function parameters to use the modern style 38 * 39 * Revision 1.4 2002/05/17 22:17:47 mike 40 * fix comment 41 * 42 * Revision 1.3 2002/05/17 22:16:36 mike 43 * move even/odd util functions to mapmutl2 44 * 45 * Revision 1.2 2001/07/15 20:55:20 mike 46 * add comments 47 * 48 * Revision 1.1 2001/07/15 20:48:27 mike 49 * Initial revision 50 */ 51 52 #include "m_apm_lc.h" 53 54 /****************************************************************************/ 55 /* 56 * From Knuth, The Art of Computer Programming: 57 * 58 * This is the binary GCD algorithm as described 59 * in the book (Algorithm B) 60 */ 61 void m_apm_gcd(M_APM r, M_APM u, M_APM v) 62 { 63 M_APM tmpM, tmpN, tmpT, tmpU, tmpV; 64 int kk, kr, mm; 65 long pow_2; 66 67 /* 'is_integer' will return 0 || 1 */ 68 69 if ((m_apm_is_integer(u) + m_apm_is_integer(v)) != 2) 70 { 71 M_apm_log_error_msg(M_APM_RETURN, "\'m_apm_gcd\', Non-integer input"); 72 73 M_set_to_zero(r); 74 return; 75 } 76 77 if (u->m_apm_sign == 0) 78 { 79 m_apm_absolute_value(r, v); 80 return; 81 } 82 83 if (v->m_apm_sign == 0) 84 { 85 m_apm_absolute_value(r, u); 86 return; 87 } 88 89 tmpM = M_get_stack_var(); 90 tmpN = M_get_stack_var(); 91 tmpT = M_get_stack_var(); 92 tmpU = M_get_stack_var(); 93 tmpV = M_get_stack_var(); 94 95 m_apm_absolute_value(tmpU, u); 96 m_apm_absolute_value(tmpV, v); 97 98 /* Step B1 */ 99 100 kk = 0; 101 102 while (TRUE) 103 { 104 mm = 1; 105 if (m_apm_is_odd(tmpU)) 106 break; 107 108 mm = 0; 109 if (m_apm_is_odd(tmpV)) 110 break; 111 112 m_apm_multiply(tmpN, MM_0_5, tmpU); 113 m_apm_copy(tmpU, tmpN); 114 115 m_apm_multiply(tmpN, MM_0_5, tmpV); 116 m_apm_copy(tmpV, tmpN); 117 118 kk++; 119 } 120 121 /* Step B2 */ 122 123 if (mm) 124 { 125 m_apm_negate(tmpT, tmpV); 126 goto B4; 127 } 128 129 m_apm_copy(tmpT, tmpU); 130 131 /* Step: */ 132 133 B3: 134 135 m_apm_multiply(tmpN, MM_0_5, tmpT); 136 m_apm_copy(tmpT, tmpN); 137 138 /* Step: */ 139 140 B4: 141 142 if (m_apm_is_even(tmpT)) 143 goto B3; 144 145 /* Step B5 */ 146 147 if (tmpT->m_apm_sign == 1) 148 m_apm_copy(tmpU, tmpT); 149 else 150 m_apm_negate(tmpV, tmpT); 151 152 /* Step B6 */ 153 154 m_apm_subtract(tmpT, tmpU, tmpV); 155 156 if (tmpT->m_apm_sign != 0) 157 goto B3; 158 159 /* 160 * result = U * 2 ^ kk 161 */ 162 163 if (kk == 0) 164 m_apm_copy(r, tmpU); 165 else 166 { 167 if (kk == 1) 168 m_apm_multiply(r, tmpU, MM_Two); 169 170 if (kk == 2) 171 m_apm_multiply(r, tmpU, MM_Four); 172 173 if (kk >= 3) 174 { 175 mm = kk / 28; 176 kr = kk % 28; 177 pow_2 = 1L << kr; 178 179 if (mm == 0) 180 { 181 m_apm_set_long(tmpN, pow_2); 182 m_apm_multiply(r, tmpU, tmpN); 183 } 184 else 185 { 186 m_apm_copy(tmpN, MM_One); 187 m_apm_set_long(tmpM, 0x10000000L); /* 2 ^ 28 */ 188 189 while (TRUE) 190 { 191 m_apm_multiply(tmpT, tmpN, tmpM); 192 m_apm_copy(tmpN, tmpT); 193 194 if (--mm == 0) 195 break; 196 } 197 198 if (kr == 0) 199 { 200 m_apm_multiply(r, tmpU, tmpN); 201 } 202 else 203 { 204 m_apm_set_long(tmpM, pow_2); 205 m_apm_multiply(tmpT, tmpN, tmpM); 206 m_apm_multiply(r, tmpU, tmpT); 207 } 208 } 209 } 210 } 211 212 M_restore_stack(5); 213 } 214 /****************************************************************************/ 215 /* 216 * u * v 217 * LCM(u,v) = ------------ 218 * GCD(u,v) 219 */ 220 221 void m_apm_lcm(M_APM r, M_APM u, M_APM v) 222 { 223 M_APM tmpN, tmpG; 224 225 tmpN = M_get_stack_var(); 226 tmpG = M_get_stack_var(); 227 228 m_apm_multiply(tmpN, u, v); 229 m_apm_gcd(tmpG, u, v); 230 m_apm_integer_divide(r, tmpN, tmpG); 231 232 M_restore_stack(2); 233 } 234 /****************************************************************************/ 235 236 #ifdef BIG_COMMENT_BLOCK 237 238 /* 239 * traditional GCD included for reference 240 * (also useful for testing ...) 241 */ 242 243 /* 244 * From Knuth, The Art of Computer Programming: 245 * 246 * To compute GCD(u,v) 247 * 248 * A1: 249 * if (v == 0) return (u) 250 * A2: 251 * t = u mod v 252 * u = v 253 * v = t 254 * goto A1 255 */ 256 void m_apm_gcd_traditional(M_APM r, M_APM u, M_APM v) 257 { 258 M_APM tmpD, tmpN, tmpU, tmpV; 259 260 tmpD = M_get_stack_var(); 261 tmpN = M_get_stack_var(); 262 tmpU = M_get_stack_var(); 263 tmpV = M_get_stack_var(); 264 265 m_apm_absolute_value(tmpU, u); 266 m_apm_absolute_value(tmpV, v); 267 268 while (TRUE) 269 { 270 if (tmpV->m_apm_sign == 0) 271 break; 272 273 m_apm_integer_div_rem(tmpD, tmpN, tmpU, tmpV); 274 m_apm_copy(tmpU, tmpV); 275 m_apm_copy(tmpV, tmpN); 276 } 277 278 m_apm_copy(r, tmpU); 279 M_restore_stack(4); 280 } 281 /****************************************************************************/ 282 283 #endif 284 285