ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
PB_Cgcd.c
Go to the documentation of this file.
1 /* ---------------------------------------------------------------------
2 *
3 * -- PBLAS auxiliary routine (version 2.0) --
4 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5 * and University of California, Berkeley.
6 * April 1, 1998
7 *
8 * ---------------------------------------------------------------------
9 */
10 /*
11 * Include files
12 */
13 #include "../pblas.h"
14 #include "../PBpblas.h"
15 #include "../PBtools.h"
16 #include "../PBblacs.h"
17 #include "../PBblas.h"
18
19 #ifdef __STDC__
20 int PB_Cgcd( int M, int N )
21 #else
22 int PB_Cgcd( M, N )
23 /*
24 * .. Scalar Arguments ..
25 */
26  int M, N;
27 #endif
28 {
29 /*
30 * Purpose
31 * =======
32 *
33 * PB_Cgcd computes and returns the Greatest Common Divisor (GCD) of two
34 * positive integers M and N using a binary gcd algorithm.
35 *
36 * Arguments
37 * =========
38 *
39 * M (input) INTEGER
40 * On entry, M must be at least zero.
41 *
42 * N (input) INTEGER
43 * On entry, N must be at least zero.
44 *
45 * -- Written on April 1, 1998 by
46 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
47 *
48 * ---------------------------------------------------------------------
49 */
50 /*
51 * .. Local Scalars ..
52 */
53  int gcd=1, m_val, n_val, t;
54 /* ..
55 * .. Executable Statements ..
56 *
57 */
58  if( M > N ) { m_val = N; n_val = M; }
59  else { m_val = M; n_val = N; }
60
61  while( m_val > 0 )
62  {
63  while( !( m_val & 1 ) )
64  {
65 /*
66 * m is even
67 */
68  m_val >>= 1; /* if n is odd, gcd( m, n ) = gcd( m / 2, n ) */
69 /*
70 * if n is odd, gcd( m, n ) = gcd( m / 2, n )
71 */
72  if( !( n_val & 1 ) )
73  {
74 /*
75 * otherwise gcd( m, n ) = 2 * gcd( m / 2, n / 2 )
76 */
77  n_val >>= 1;
78  gcd <<= 1;
79  }
80  }
81 /*
82 * m is odd now. If n is odd, gcd( m, n ) = gcd( m, ( m - n ) / 2 ).
83 * Otherwise, gcd( m, n ) = gcd ( m, n / 2 ).
84 */
85  n_val -= ( n_val & 1 ) ? m_val : 0;
86  n_val >>= 1;
87  while( n_val >= m_val )
88  {
89 /*
90 * If n is odd, gcd( m, n ) = gcd( m, ( m - n ) / 2 ).
91 * Otherwise, gcd( m, n ) = gcd ( m, n / 2 )
92 */
93  n_val -= ( n_val & 1 ) ? m_val : 0;
94  n_val >>= 1;
95  }
96 /*
97 * n < m, gcd( m, n ) = gcd( n, m )
98 */
99  t = n_val;
100  n_val = m_val;
101  m_val = t;
102  }
103  return ( n_val * gcd );
104 /*
105 * End of PB_Cgcd
106 */
107 }
PB_Cgcd
int PB_Cgcd(int M, int N)
Definition: PB_Cgcd.c:22