OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
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__
20Int PB_Cgcd( Int M, Int N )
21#else
22Int 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}
Int
#define Int
Definition Bconfig.h:22
PB_Cgcd
Int PB_Cgcd()
N
#define N
Definition example_user.c:10