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