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 */
61void	m_apm_gcd(M_APM r, M_APM u, M_APM v)
62{
63M_APM   tmpM, tmpN, tmpT, tmpU, tmpV;
64int	kk, kr, mm;
65long    pow_2;
66
67/* 'is_integer' will return 0 || 1 */
68
69if ((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
77if (u->m_apm_sign == 0)
78  {
79   m_apm_absolute_value(r, v);
80   return;
81  }
82
83if (v->m_apm_sign == 0)
84  {
85   m_apm_absolute_value(r, u);
86   return;
87  }
88
89tmpM = M_get_stack_var();
90tmpN = M_get_stack_var();
91tmpT = M_get_stack_var();
92tmpU = M_get_stack_var();
93tmpV = M_get_stack_var();
94
95m_apm_absolute_value(tmpU, u);
96m_apm_absolute_value(tmpV, v);
97
98/* Step B1 */
99
100kk = 0;
101
102while (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
123if (mm)
124  {
125   m_apm_negate(tmpT, tmpV);
126   goto B4;
127  }
128
129m_apm_copy(tmpT, tmpU);
130
131/* Step: */
132
133B3:
134
135m_apm_multiply(tmpN, MM_0_5, tmpT);
136m_apm_copy(tmpT, tmpN);
137
138/* Step: */
139
140B4:
141
142if (m_apm_is_even(tmpT))
143  goto B3;
144
145/* Step B5 */
146
147if (tmpT->m_apm_sign == 1)
148  m_apm_copy(tmpU, tmpT);
149else
150  m_apm_negate(tmpV, tmpT);
151
152/* Step B6 */
153
154m_apm_subtract(tmpT, tmpU, tmpV);
155
156if (tmpT->m_apm_sign != 0)
157  goto B3;
158
159/*
160 *  result = U * 2 ^ kk
161 */
162
163if (kk == 0)
164   m_apm_copy(r, tmpU);
165else
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
212M_restore_stack(5);
213}
214/****************************************************************************/
215/*
216 *                      u * v
217 *      LCM(u,v)  =  ------------
218 *                     GCD(u,v)
219 */
220
221void	m_apm_lcm(M_APM r, M_APM u, M_APM v)
222{
223M_APM   tmpN, tmpG;
224
225tmpN = M_get_stack_var();
226tmpG = M_get_stack_var();
227
228m_apm_multiply(tmpN, u, v);
229m_apm_gcd(tmpG, u, v);
230m_apm_integer_divide(r, tmpN, tmpG);
231
232M_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 */
256void	m_apm_gcd_traditional(M_APM r, M_APM u, M_APM v)
257{
258M_APM   tmpD, tmpN, tmpU, tmpV;
259
260tmpD = M_get_stack_var();
261tmpN = M_get_stack_var();
262tmpU = M_get_stack_var();
263tmpV = M_get_stack_var();
264
265m_apm_absolute_value(tmpU, u);
266m_apm_absolute_value(tmpV, v);
267
268while (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
278m_apm_copy(r, tmpU);
279M_restore_stack(4);
280}
281/****************************************************************************/
282
283#endif
284
285