Skip to content

Commit

Permalink
Update math.hpp
Browse files Browse the repository at this point in the history
  • Loading branch information
steve02081504 committed Oct 28, 2023
1 parent 682768f commit b3f673c
Showing 1 changed file with 56 additions and 44 deletions.
100 changes: 56 additions & 44 deletions parts/header_file/files/elc/_files/base_defs/math.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -910,7 +910,8 @@ namespace math{

//求出最大公约数
template<integer_type T>
[[nodiscard]]inline auto base_gcd(T x, T y)noexcept{
[[nodiscard]]inline auto binary_gcd(T x, T y)noexcept{
static_assert(BIT_POSSIBILITY==2);
size_t shift = 0;
while(y){
//都是偶数时,shift++并位移
Expand All @@ -927,64 +928,75 @@ namespace math{
// 返回 x 左移 shift 位的结果
return move(x) << shift;
}
//求出最大公约数
//求出最大公约数:欧几里得算法
template<integer_type T>
[[nodiscard]]inline auto euclidean_gcd(T x, T y)noexcept{
if constexpr(is_signed<T>)
if(is_negative(x) && is_negative(y))
return euclidean_gcd(abs(x), abs(y));
while(y){
if(x < y)swap(x, y);
x%=y;
}
return x;
}
template<unsigned_big_integer_type T>
[[nodiscard]]inline auto gcd(T x, T y)noexcept{
[[nodiscard]]inline auto lehmer_gcd(T x, T y)noexcept{
//假设x≥y
if(x < y)swap(x, y);
//如果y只包含一位数,则使用其他方法求出结果
if (y.size_in_base_type() == 1)
return base_gcd(x, y);

//如果x和y的数位长度不同,使x和y的长度相等于m
size_t m=exλ{
size_t m1=x.size_in_base_type(),m2=y.size_in_base_type();
if(m1>m2)
x = x.right_shift_in_base_type(m1-m2);
//else if(m1<m2):不可能,因为x≥y
return m2;
}();

size_t m=y.size_in_base_type();
if(m <= 1)return euclidean_gcd(x, y);//如果y只包含一位数(或根本没有),则使用其他方法求出结果
T A, B, C, D;
do{
while(x && y) {//外循环
--m;//提高精度
T a = x.right_shift_in_base_type(m), b = y.right_shift_in_base_type(m);//取最显著数位
/*初始化矩阵[A B a]为扩展的同一矩阵[1 0 a]
[C D b] [0 1 b]
*/
A = 1u, B = T{}, C = T{}, D = 1u;
//对(a+A,b+C)和(a+B,b+D)同时执行欧几里得算法,直到商不同为止
do{
const auto w1 = (a+A) / (b+C);
outer_loop:
while(x && y) {//外循环
--m;//提高精度
T a = x.right_shift_in_base_type(m), b = y.right_shift_in_base_type(m);//取最显著数位
/*初始化矩阵[A B a]为扩展的同一矩阵[1 0 a]
[C D b] [0 1 b]
*/
A = D = 1u; B = C = T{};
//对(a+A,b+C)和(a+B,b+D)同时执行欧几里得算法,直到商不同为止
do{//内循环
const auto w = (a+A) / (b+C);//w为欧几里得算法长除法链中当前长除法的商
{
const auto w2 = (a+B) / (b+D);
if(w1 != w2)break;
const auto& w = w1;//w为欧几里得算法长除法链中当前长除法的商
/*根据扩展欧几里得算法的矩阵表述替换当前矩阵[A B a]到矩阵积[0 1] [A B a] [C D a ]
[C D b] ⋅ =
[1 −w] [C D b] [A−wC B−wD a−wb]
*/
A-=w*C;swap(A,C);
B-=w*D;swap(B,D);
a-=w*b;swap(a,b);
}while(B);//B≠0则继续
//B=0意味着陷入了死锁;执行欧几里得算法的正常步骤,然后重新开始外循环
x = swap(y, x%y);
}
if(w != w2)
if(!B){//B=0意味着陷入了死锁;执行欧几里得算法的正常步骤,然后重新开始外循环
x = swap(y, x%y);
goto outer_loop;
}else break;
}
/*根据扩展欧几里得算法的矩阵表述替换当前矩阵[A B a]到矩阵积[0 1] [A B a] [C D a ]
[C D b] ⋅ =
[1 −w] [C D b] [A−wC B−wD a−wb]
*/
A-=w*C;swap(A,C);
B-=w*D;swap(B,D);
a-=w*b;swap(a,b);
}while(B);//B≠0则继续
//将对压缩形式的前数位执行的欧几里得算法步骤应用到x和y上
auto new_x = A*x+B*y, new_y = C*x+D*y;
x = move(new_x),y = move(new_y);
}while(y);//如果y≠0,则转到外循环的起点。
}
return x;
}
template<basic_integer_type T>
[[nodiscard]]force_inline auto gcd(T x, T y)noexcept{
return base_gcd(x, y);
[[nodiscard]]force_inline auto lehmer_gcd(T x, T y)noexcept{
return euclidean_gcd(move(x), move(y));
}
template<signed_big_integer_type T>
[[nodiscard]]force_inline auto gcd(T x, T y)noexcept{
[[nodiscard]]force_inline auto lehmer_gcd(T x, T y)noexcept{
const bool sign = is_negative(x) && is_negative(y);
return copy_as_negative(gcd(abs(x), abs(y)), sign);
return copy_as_negative(gcd(abs(move(x)), abs(move(y))), sign);
}

template<integer_type T>
[[nodiscard]]inline auto gcd(T x, T y)noexcept{
if constexpr(BIT_POSSIBILITY==2)
return binary_gcd(move(x), move(y));
else
return lehmer_gcd(move(x), move(y));
}
}
using math::gcd;
Expand Down

0 comments on commit b3f673c

Please sign in to comment.