"Fossies" - the Fresh Open Source Software Archive 
Member "mapm_4.9.5a/mapm_gcd.c" (21 Feb 2010, 5530 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 "mapm_gcd.c" see the
Fossies "Dox" file reference documentation.
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