Skip to content

Instantly share code, notes, and snippets.

@marionette-of-u

marionette-of-u/integer.hpp

Last active August 29, 2015 14:26
Show Gist options
  • Save marionette-of-u/5825c4a7c2b542e21f67 to your computer and use it in GitHub Desktop.
Save marionette-of-u/5825c4a7c2b542e21f67 to your computer and use it in GitHub Desktop.
Download ZIP
多倍精度整数を使った有理数
Raw
integer.hpp
#ifndef MULTI_PRECISION_INCLUDE_INTEGER_HPP
#define MULTI_PRECISION_INCLUDE_INTEGER_HPP
#include <algorithm>
#include <vector>
#include <iostream>
#include <string>
#include <iterator>
#include <utility>
#include <random>
#include <type_traits>
#include <cstdint>
#include <climits>
#include <cmath>
namespace multi_precision{
template<class UInt = std::uint16_t, class DoubleUInt = std::uint32_t, class DoubleInt = std::int32_t, DoubleUInt BitNum = sizeof(UInt) * CHAR_BIT>
class integer{
static_assert(sizeof(UInt) * 2 <= sizeof(DoubleUInt), "sizeof(UInt) * 2 <= sizeof(DoubleUInt)");
static_assert(std::is_signed<DoubleInt>::value, "std::is_signed<DoubleInt>::value");
public:
static const DoubleUInt base_mask = (static_cast<DoubleUInt>(1) << BitNum) - 1;
static const UInt half = static_cast<UInt>(static_cast<DoubleUInt>(1) << (BitNum - 1));
using data_type = std::vector<UInt>;
data_type data;
DoubleInt sign = +1;
integer() = default;
integer(const integer&) = default;
integer(integer &&other) : data(std::move(other.data)), sign(other.sign){}
integer(int x){
if(x == 0){ return; }
data.resize(1);
data[0] = std::abs(x);
sign = x >= 0 ? +1 : - 1;
}
integer(UInt x){
if(x == 0){ return; }
data.resize(1);
data[0] = x;
sign = +1;
}
integer(const char *str){
build_from_str(str);
}
integer(const std::string &str){
build_from_str(str.begin());
}
template<class StrIter>
void build_from_str(StrIter iter){
if(*iter == '-'){
sign = - 1;
++iter;
}
char buff[2] = { 0 };
while(*iter){
buff[0] = *iter;
int a = std::atoi(buff);
*this *= 10;
unsigned_single_add(a);
++iter;
}
}
~integer() = default;
integer &operator =(const integer &other){
data = other.data;
sign = other.sign;
return *this;
}
integer &operator =(integer &&other){
data = std::move(other.data);
sign = std::move(other.sign);
return *this;
}
operator bool() const{
return !data.empty();
}
bool operator <(const integer &other) const{
if(data.empty() && other.data.empty()){ return false; }
if(data.size() > 0 && other.data.empty()){ return sign == -1; }
if(data.empty() && other.data.size() > 0){ return other.sign == +1; }
if(sign < other.sign){ return true; }
if(sign > other.sign){ return false; }
if(data.size() < other.data.size()){ return sign > 0; }
if(data.size() > other.data.size()){ return sign < 0; }
std::size_t i = data.size() - 1;
do{
if(data[i] < other.data[i]){ return true; }
} while(i-- > 0);
return false;
}
bool operator >(const integer &other) const{
return other < *this;
}
bool operator <=(const integer &other) const{
if(sign < other.sign){ return true; }
if(sign > other.sign){ return false; }
if(data.size() < other.data.size()){ return sign > 0; }
if(data.size() > other.data.size()){ return sign < 0; }
std::size_t i = data.size() - 1;
do{
if(data[i] > other.data[i]){ return false; }
} while(i-- > 0);
return true;
}
bool operator >=(const integer &other) const{
return other <= *this;
}
bool operator ==(const integer &other) const{
return sign == other.sign && data == other.data;
}
bool operator !=(const integer &other) const{
return !(*this == other);
}
static bool unsigned_less(const data_type &lhs, const data_type &rhs){
if(lhs.size() < rhs.size()){ return true; }
if(lhs.size() > rhs.size()){ return false; }
for(std::size_t i = lhs.size() - 1; i + 1 > 0; --i){
if(lhs[i] < rhs[i]){ return true; }
if(lhs[i] > rhs[i]){ return false; }
}
return false;
}
static bool unsigned_less_eq(const data_type &lhs, const data_type &rhs){
if(lhs.size() < rhs.size()){ return true; }
if(lhs.size() > rhs.size()){ return false; }
for(std::size_t i = lhs.size() - 1; i + 1 > 0; --i){
if(lhs[i] < rhs[i]){ return true; }
if(lhs[i] > rhs[i]){ return false; }
}
return true;
}
void add(const integer &other, std::size_t n = 0){
if(sign == other.sign){
unsigned_add(other, n);
}else{
if(unsigned_less(data, other.data)){
integer temp(*this);
data = other.data;
unsigned_sub(temp);
sign = other.sign;
}else{
unsigned_sub(other);
}
}
normalize_data_size();
}
void sub(const integer &other, std::size_t n = 0){
if(sign != other.sign){
unsigned_add(other, n);
}else{
if(unsigned_less(data, other.data)){
integer temp(*this);
data = other.data;
unsigned_sub(temp);
sign *= -1;
}else{
unsigned_sub(other);
}
}
normalize_data_size();
}
void unsigned_add(const integer &other, std::size_t n = 0){
if(data.size() < other.data.size() + n){
data.resize(other.data.size() + n);
}
DoubleUInt c = 0;
std::size_t i;
for(i = 0; i < other.data.size(); ++i){
DoubleUInt v = static_cast<DoubleUInt>(data[i + n]) + static_cast<DoubleUInt>(other.data[i]) + c;
data[i + n] = static_cast<UInt>(v & base_mask);
c = v >> BitNum;
}
for(; c > 0; ++i){
if(i >= data.size()){ data.resize(data.size() + 1); }
DoubleUInt v = static_cast<DoubleUInt>(data[i + n]) + c;
data[i + n] = static_cast<UInt>(v & base_mask);
c = v >> BitNum;
}
}
integer &operator ++(){
if(sign == +1){
unsigned_single_add(1);
}else{
unsigned_single_sub(1);
}
return *this;
}
integer operator ++(int){
integer r;
++*this;
return r;
}
integer &operator --(){
if(sign == -1){
unsigned_single_add(1);
}else{
unsigned_single_sub(1);
}
return *this;
}
integer operator --(int){
integer r;
--*this;
return r;
}
void unsigned_single_add(UInt c, std::size_t i = 0){
if(i >= data.size()){ data.resize(i + 1); }
for(; c > 0; ++i){
DoubleUInt v = static_cast<DoubleUInt>(data[i]) + c;
data[i] = static_cast<UInt>(v & base_mask);
c = v >> BitNum;
}
}
void unsigned_sub(const integer &other, std::size_t i = 0){
unsigned_sub(other.data.begin(), other.data.end(), i);
}
void unsigned_sub(const typename data_type::const_iterator &other_first, const typename data_type::const_iterator &other_last, std::size_t i = 0){
typename data_type::const_iterator iter = other_first;
UInt c = 0;
for(std::size_t n = std::distance(other_first, other_last); i < n; ++i, ++iter){
const UInt &other_value(*iter);
UInt t = data[i] - (other_value + c);
if(data[i] < other_value + c){
c = 1;
}else{
c = 0;
}
data[i] = t;
}
for(; c > 0; ++i){
UInt t = data[i] - c;
if(data[i] < c){
c = 1;
}else{
c = 0;
}
data[i] = t;
}
normalize_data_size();
}
void unsigned_single_sub(UInt c, std::size_t i = 0){
for(; c > 0; ++i){
UInt t = data[i] - c;
if(data[i] < c){
c = 1;
}else{
c = 0;
}
data[i] = t;
}
}
static void unsigned_mul(integer &r, const integer &lhs, const integer &rhs){
std::size_t s = lhs.data.size() + rhs.data.size() - 1;
r.data.resize(s + 1);
for(std::size_t i = 0; i < lhs.data.size(); ++i){
for(std::size_t j = 0; j < rhs.data.size(); ++j){
DoubleUInt c = static_cast<DoubleUInt>(lhs.data[i]) * static_cast<DoubleUInt>(rhs.data[j]);
for(std::size_t k = 0; i + j + k < s + 1; ++k){
std::size_t u = i + j + k;
DoubleUInt v = static_cast<DoubleUInt>(r.data[u]) + c;
UInt a = static_cast<UInt>(v & base_mask);
r.data[u] = a;
c = v >> BitNum;
}
}
}
r.normalize_data_size();
}
static integer mul(const integer &lhs, const integer &rhs){
integer r;
r.kar_mul(lhs, rhs);
return r;
}
static integer classic_mul(const integer &lhs, const integer &rhs){
integer r;
unsigned_mul(r, lhs, rhs);
return r;
}
static void unsigned_single_mul(integer &r, const integer &lhs, UInt rhs){
std::size_t s = lhs.data.size() + 1;
r.data.resize(s + 1);
for(std::size_t i = 0; i < lhs.data.size(); ++i){
DoubleUInt c = static_cast<DoubleUInt>(lhs.data[i]) * static_cast<DoubleUInt>(rhs);
for(std::size_t k = 0; i + k < s + 1; ++k){
std::size_t u = i + k;
DoubleUInt v = static_cast<DoubleUInt>(r.data[u]) + c;
UInt a = static_cast<UInt>(v & base_mask);
r.data[u] = a;
c = v >> BitNum;
}
}
r.normalize_data_size();
}
struct kar_range{
kar_range() = default;
kar_range(typename data_type::const_iterator first, typename data_type::const_iterator last) : first(first), last(last), size(last - first){}
kar_range(typename data_type::const_iterator first, typename data_type::const_iterator last, std::size_t size) : first(first), last(last), size(size){}
typename data_type::const_iterator first, last;
std::size_t size;
};
void kar_mul(const integer &lhs, const integer &rhs){
kar_rec_mul(kar_range(lhs.data.begin(), lhs.data.end()), kar_range(rhs.data.begin(), rhs.data.end()));
sign = lhs.sign * rhs.sign;
}
void kar_mul_range(const kar_range &lhs, const kar_range &rhs){
if(lhs.size == 0 || rhs.size == 0){
data.clear();
sign = +1;
return;
}
data.resize(lhs.size + rhs.size);
if(data.empty()){
normalize_data_size();
return;
}
std::size_t n = 0;
for(typename data_type::const_iterator rhs_iter = rhs.first; rhs_iter != rhs.last; ++rhs_iter, ++n){
UInt rhs_value = *rhs_iter, lhs_value;
std::size_t m = 0;
for(typename data_type::const_iterator lhs_iter = lhs.first; lhs_iter != lhs.last; ++lhs_iter, ++m){
lhs_value = *lhs_iter;
DoubleUInt a = static_cast<DoubleUInt>(rhs_value) * static_cast<DoubleUInt>(lhs_value);
unsigned_single_add(static_cast<UInt>(a & base_mask), n + m);
UInt temp = static_cast<UInt>((a >> BitNum) & base_mask);
if(temp > 0){
unsigned_single_add(temp, n + m + 1);
}
}
}
}
bool kar_less(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last) const{
return kar_less_impl<false>(rhs_first, rhs_last);
}
bool kar_less_eq(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last) const{
return kar_less_impl<true>(rhs_first, rhs_last);
}
template<bool Final>
bool kar_less_impl(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last) const{
return kar_less_impl<Final>(data.begin(), data.end(), rhs_first, rhs_last);
}
template<bool Final>
static bool kar_less_impl(
typename data_type::const_iterator lhs_first, typename data_type::const_iterator lhs_last,
typename data_type::const_iterator rhs_first, typename data_type::const_iterator rhs_last
){
std::size_t lhs_n = std::distance(lhs_first, lhs_last), rhs_n = std::distance(rhs_first, rhs_last);
if(lhs_n == 0 && rhs_n == 0){
return false;
}else if(lhs_n == 0 && rhs_n > 0){
return true;
}else if(rhs_n == 0 && lhs_n > 0){
return false;
}
typename data_type::const_iterator lhs_iter = lhs_first, rhs_iter = rhs_first;
if(lhs_n < rhs_n){
std::size_t n = rhs_n - lhs_n;
for(std::size_t i = 0; i < n; ++i){
if(*(rhs_iter + (rhs_n - i - 1)) > 0){
return true;
}
}
rhs_iter += lhs_n - 1;
lhs_iter = lhs_first + (lhs_n - 1);
}else if(lhs_n > rhs_n){
std::size_t n = lhs_n - rhs_n;
for(std::size_t i = 0; i < n; ++i){
if(*(lhs_iter + (lhs_n - i - 1)) > 0){
return false;
}
lhs_iter += rhs_n - 1;
rhs_iter = rhs_iter + (rhs_n - 1);
}
}else{
lhs_iter += lhs_n - 1;
rhs_iter += rhs_n - 1;
}
for(; ; ){
UInt l = *lhs_iter, r = *rhs_iter;
if(l < r){ return true; }
if(l > r){ return false; }
if(lhs_iter == lhs_first){
break;
}
--lhs_iter, --rhs_iter;
}
return Final;
}
void kar_add(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last, std::size_t n = 0){
std::size_t rhs_size = std::distance(rhs_first, rhs_last);
if(rhs_size == 0){
return;
}
rhs_size += n;
if(data.size() < rhs_size){
data.resize(rhs_size);
}
typename data_type::iterator operand_iter = data.begin() + n;
DoubleUInt c = 0;
typename data_type::const_iterator iter = rhs_first, end = rhs_last;
if(iter != end){
--end;
}
for(; iter != end; ++iter, ++operand_iter){
UInt &operand(*operand_iter);
DoubleUInt s = operand + *iter + c;
operand = static_cast<UInt>(s & base_mask);
c = s >> BitNum;
}
if(operand_iter != data.end()){
UInt &operand(*operand_iter);
DoubleUInt s = static_cast<DoubleUInt>(operand) + static_cast<DoubleUInt>(*iter) + c;
operand = static_cast<UInt>(s & base_mask);
c = s >> BitNum;
if(c > 0){
data.push_back(static_cast<UInt>(c & base_mask));
}
}
}
void kar_signed_add(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last){
if(sign < 0){
if(kar_less_eq(rhs_first, rhs_last)){
integer s(std::move(*this));
data.assign(rhs_first, rhs_last);
unsigned_sub(s);
sign = +1;
}else{
unsigned_sub(rhs_first, rhs_last);
}
}else{
kar_add(rhs_first, rhs_last);
}
}
void kar_sub(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last){
if(kar_less_eq(rhs_first, rhs_last)){
integer s(std::move(*this));
data.assign(rhs_first, rhs_last);
unsigned_sub(s);
sign *= -1;
}else{
std::size_t n = std::distance(rhs_first, rhs_last);
if(n > data.size()){
unsigned_sub(rhs_first, rhs_first + data.size());
}else{
unsigned_sub(rhs_first, rhs_last);
}
}
normalize_data_size();
}
void kar_rec_mul(const kar_range &x, const kar_range &y){
std::size_t n = ceil_pow2_single((std::max)(x.size, y.size));
if(n < 2){
if(x.size > 0 && y.size > 0){
kar_mul_range(x, y);
normalize_data_size();
}
return;
}
n /= 2;
std::size_t xn = x.size < n ? x.size : n, yn = y.size < n ? y.size : n;
kar_range x0, y0, x1, y1;
x0.first = x.first, x0.last = x.first + xn, x0.size = xn;
x1.first = x.first + xn, x1.last = x.last, x1.size = x.size - xn;
y0.first = y.first, y0.last = y.first + yn, y0.size = yn;
y1.first = y.first + yn, y1.last = y.last, y1.size = y.size - yn;
integer z2, z0;
{
integer tx, ty;
if(x1.first != x1.last){
tx.data.assign(x1.first, x1.last);
tx.kar_sub(x0.first, x0.last);
}else{
tx.data.assign(x0.first, x0.last);
tx.sign = - 1;
}
if(y1.first != y1.last){
ty.data.assign(y1.first, y1.last);
ty.kar_sub(y0.first, y0.last);
}else{
ty.data.assign(y0.first, y0.last);
ty.sign = - 1;
}
kar_rec_mul(kar_range(tx.data.begin(), tx.data.end()), kar_range(ty.data.begin(), ty.data.end()));
sign = tx.sign != ty.sign ? +1 : - 1;
}
z0.kar_rec_mul(x0, y0);
add(z0);
if(x.size >= n && y.size >= n){
z2.kar_rec_mul(x1, y1);
add(z2);
}
elemental_shift(n);
add(z0);
if(x.size >= n && y.size >= n){
add(z2, n * 2);
}
normalize_data_size();
}
struct quo_rem{
integer quo, rem;
quo_rem() = default;
quo_rem(const quo_rem&) = default;
quo_rem(quo_rem &&other) : quo(std::move(other.quo)), rem(std::move(other.rem)){}
~quo_rem() = default;
};
static quo_rem unsigned_div(const integer &lhs, const integer &rhs){
return unsigned_div(lhs.data.begin(), lhs.data.end(), rhs.data.begin(), rhs.data.end());
}
static quo_rem unsigned_div(
typename data_type::const_iterator l_first,
typename data_type::const_iterator l_last,
typename data_type::const_iterator r_first,
typename data_type::const_iterator r_last
){
quo_rem qr;
std::size_t ln = std::distance(l_first, l_last), rn = std::distance(r_first, r_last);
qr.quo.data.reserve(ln);
qr.rem.data.reserve(rn);
for(std::size_t i = ln * static_cast<std::size_t>(BitNum) - 1; i + 1 > 0; --i){
if(!qr.rem.data.empty()){
qr.rem <<= 1;
}else{
qr.rem.data.push_back(0);
}
qr.rem.data[0] |= (l_first[i / static_cast<std::size_t>(BitNum)] >> (i % static_cast<std::size_t>(BitNum))) & 1;
if(kar_less_impl<true>(r_first, r_last, qr.rem.data.begin(), qr.rem.data.end())){
qr.rem.unsigned_sub(r_first, r_last);
std::size_t t = i / static_cast<std::size_t>(BitNum) + 1;
if(qr.quo.data.size() < t){
qr.quo.data.resize(t);
}
qr.quo.data[t - 1] |= 1 << (i % static_cast<std::size_t>(BitNum));
}
}
return qr;
}
static quo_rem kar_rec_div(kar_range lhs, kar_range rhs){
quo_rem qr;
std::size_t n = (rhs.size + (rhs.size % 2)) / 2;
if(n <= 2 || lhs.size > n * 4){
qr = unsigned_div(lhs.first, lhs.last, rhs.first, rhs.last);
qr.quo.normalize_data_size();
qr.rem.normalize_data_size();
return qr;
}
{
kar_range u[4], v[2];
{
std::size_t a = static_cast<std::size_t>(lhs.last - lhs.first);
for(std::size_t i = 0; i < 4; ++i){
std::size_t p = i * n;
u[i].first = lhs.first + std::min(p, a);
u[i].last = u[i].first + std::min(n, static_cast<std::size_t>(lhs.last - u[i].first));
u[i].size = u[i].last - u[i].first;
}
a = static_cast<std::size_t>(rhs.last - rhs.first);
for(std::size_t i = 0; i < 2; ++i){
std::size_t p = i * n;
v[i].first = rhs.first + std::min(p, a);
v[i].last = v[i].first + std::min(n, static_cast<std::size_t>(rhs.last - v[i].first));
v[i].size = v[i].last - v[i].first;
}
}
quo_rem qr1;
{
kar_range u3_u2(u[2].first, u[3].last);
qr1 = kar_rec_div(u3_u2, v[1]);
qr1.rem.elemental_shift(n);
qr1.rem.kar_add(u[1].first, u[1].last);
{
integer vw;
vw.kar_mul_range(v[0], kar_range(qr1.quo.data.begin(), qr1.quo.data.end()));
qr1.rem.kar_sub(vw.data.begin(), vw.data.end());
}
if(qr1.rem.sign < 0){
integer &r1(qr1.rem), prod;
quo_rem mod = kar_rec_div(kar_range(r1.data.begin(), r1.data.end()), rhs);
if(!mod.rem.data.empty()){
mod.quo.unsigned_single_add(1);
}
prod.kar_mul_range(kar_range(mod.quo.data.begin(), mod.quo.data.end()), rhs);
r1.kar_signed_add(prod.data.begin(), prod.data.end());
qr1.quo.unsigned_sub(mod.quo);
}
}
quo_rem qr2;
{
kar_range r1(qr1.rem.data.begin(), qr1.rem.data.end());
qr2 = kar_rec_div(r1, v[1]);
integer &r0(qr2.rem);
r0.elemental_shift(n);
r0.kar_add(u[0].first, u[0].last);
{
integer vw;
vw.kar_mul_range(v[0], kar_range(qr2.quo.data.begin(), qr2.quo.data.end()));
qr2.rem.kar_sub(vw.data.begin(), vw.data.end());
}
if(r0.sign < 0){
integer prod;
quo_rem mod = kar_rec_div(kar_range(r0.data.begin(), r0.data.end()), rhs);
if(!mod.rem.data.empty()){
mod.quo.unsigned_single_add(1);
}
prod.kar_mul_range(kar_range(mod.quo.data.begin(), mod.quo.data.end()), rhs);
r0.kar_signed_add(prod.data.begin(), prod.data.end());
qr2.quo.unsigned_sub(mod.quo);
}
}
qr.quo = std::move(qr2.quo);
qr.quo.kar_add(qr1.quo.data.begin(), qr1.quo.data.end(), n);
qr.rem = std::move(qr2.rem);
}
return qr;
}
#if 0
static quo_rem kar_large_div(kar_range lhs, kar_range rhs){
quo_rem qr;
integer &quo(qr.quo), &rem(qr.rem);
quo.data.reserve(lhs.size);
rem.data.assign(lhs.first, lhs.last);
using iter_type = typename data_type::const_iterator;
auto compare = [](iter_type rem_first, iter_type rem_last, iter_type rhs_iter, iter_type rhs_last){
for(auto iter = rem_last - 1, jter = rhs_last - 1; ; --iter, --jter){
if(*iter > *jter){
return true;
}else if(*iter == *jter){
if(iter == rem_last){
break;
}
}else{
return false;
}
}
return true;
};
std::size_t i = rem.data.size() - 1;
for(; rem.data.size() > rhs.size || (rem.data.size() == rhs.size && compare(rem.data.begin(), rem.data.end(), rhs.first, rhs.last)); --i){
kar_range rem_range(rem.data.begin() + i, rem.data.end());
quo_rem st = kar_rec_div(rem_range, rhs);
quo.kar_add(st.quo.data.begin(), st.quo.data.end(), i);
rem.data.resize(i + st.rem.data.size());
std::copy(st.rem.data.begin(), st.rem.data.end(), rem.data.begin() + i);
}
return qr;
}
#endif
static quo_rem div(const integer &lhs, const integer &rhs){
quo_rem qr = unsigned_div(lhs, rhs);
qr.quo.sign = qr.rem.sign = lhs.sign == rhs.sign ? +1 : - 1;
return qr;
}
static quo_rem modern_div(const integer &lhs, const integer &rhs){
quo_rem qr = kar_rec_div(kar_range(lhs.data.begin(), lhs.data.end()), kar_range(rhs.data.begin(), rhs.data.end()));
qr.quo.sign = qr.rem.sign = lhs.sign * rhs.sign;
return qr;
}
// Extended Euclidean Algorithm.
struct eea_result{
DoubleInt result, s, t;
};
static eea_result eea(DoubleInt a, DoubleInt b){
eea_result p;
if(b == 0){
p.result = a;
p.s = 1;
p.t = 0;
return p;
}
p = eea(b, std::abs(a % b));
eea_result q;
q.result = p.result;
q.s = p.t;
q.t = p.s - (a / b) * p.t;
return q;
}
// FFT.
static data_type fft(const data_type &arg_data, int coe){
// arg_data の内部表現データのサイズを 2 の冪乗に繰り上げて得る.
std::size_t n = ceil_pow2_single(static_cast<UInt>(arg_data.size()));
// log_2(n) を得る.
std::size_t lg = bit_num(n) - 1;
// arg_data 内部表現を並べ替えて a にコピーする.
std::vector<DoubleInt> a(n), b(n, 1);
for(std::size_t i = 0; i < arg_data.size(); ++i){
a[bit_reverse(i, lg)] = static_cast<DoubleInt>(arg_data[i]);
}
// n ごとに素数 p = k * 2^n + 1 の最小の k を示すテーブル.
// 実用上は 2^64 以下までで十分.
// enumeration_prime_k 関数によって計算された.
static const UInt k_table[] = {
1, 1, 2, 1, 3, 3, 2, 1,
15, 12, 6, 3, 5, 4, 2, 1,
6, 3, 11, 7, 33, 25, 20, 10,
5, 7, 15, 12, 6, 3, 35, 18,
9, 12, 6, 3, 20, 10, 5, 6,
3, 9, 9, 15, 35, 19, 27, 15,
14, 7, 14, 7, 20, 10, 5, 27,
29, 54, 27, 31
};
for(std::size_t i = 1; i <= lg; ++i){
std::size_t m = 1 << i;
UInt q = k_table[i - 1], p = q * static_cast<UInt>(m) + 1;
DoubleInt zeta = eea(coe * static_cast<DoubleInt>(q), p).s, omega = 1;
for(std::size_t j = 0; j < m / 2; ++j){
for(std::size_t k = j; k < n; k += m){
std::size_t l = k + m / 2;
DoubleInt t = mod(omega * b[l], p), u = b[k];
b[k] = mod(t + u, p);
if(b[k] == 0){
a[k] = a[l] + a[k];
}
b[l] = mod(t - u, p);
if(b[l] == 0){
a[l] = a[l] - a[k];
}
}
omega = mod(omega * zeta, p);
}
}
data_type result;
result.resize(a.size());
for(std::size_t i = 0; i < a.size(); ++i){
result[i] = static_cast<UInt>(a[i]);
}
return result;
}
// Bit Reverse.
static std::size_t bit_reverse(UInt k, UInt n){
std::size_t l = 0u;
for(UInt i = 0u; i < n; ++i){
l += ((k & (1 << i)) ? 1 : 0) << (n - i - 1);
}
return l;
}
// prime = k * n + 1 になる k を探す.
static int search_prime_k(const integer &n){
int k = 1;
while(!miller_rabin_test(k * n + 1, 100)){
k += 1;
}
return k;
}
// prime = k * n + 1 になる k を列挙する.
template<class OStream>
static void enumeration_prime_k(OStream &os){
integer bound = 1;
for(int i = 0; i < 64; ++i){
bound <<= 1;
}
integer p;
for(int a = 1; ; ++a){
p = 1;
p <<= a;
if(p > bound){
break;
}
int k = search_prime_k(p);
os << k << (a % 8 == 0 ? ",\n" : ",\t");
if(k * p + 1 > bound){
break;
}
}
}
// Primality Test (2).
static bool miller_rabin_test(const integer &n, std::size_t s){
static std::mt19937 g(0xABABABABu);
integer a, m = n - 2;
for(std::size_t i = 0; i < s; ++i){
if(witness_test(m.random(g) + 1, n)){
return false;
}
}
return true;
}
// Primality Test (1).
static bool witness_test(const integer &a, const integer &n){
integer b = n - 1, d = 1, x, t;
for(int i = static_cast<int>(b.bit_num()) - 1; i >= 0; --i){
x = d;
d = (d * d) % n;
if(d == 1 && x != 1 && x != b){
return true;
}
t = b >> static_cast<UInt>(i);
if(t > 0 && t.data[0] % 2 == 1){
d = (d * a) % n;
}
}
return d != 1;
}
// Random Number
template<class Generator>
integer random(Generator &g) const{
std::size_t s = g() % data.size() + 1;
integer x;
x.data.resize(s);
for(std::size_t i = 0; i < s - 1; ++i){
x.data[i] = g() & base_mask;
}
if(s == data.size()){
x.data.back() = g() % data.back();
}
x.normalize_data_size();
return x;
}
// 一般化逆数.
static integer generalized_inverse(const integer &n){
std::size_t b = (n - 1).bit_num();
integer r = 1;
{
integer s;
r <<= b;
s = r;
for(; ; ){
r = 2 * r - ((((r * r) >> b) * n) >> b);
if(r <= s){
break;
}
s = r;
}
}
integer y = (integer(1) << (2 * b)) - n * r;
while(y < 0){
r -= 1;
y += n;
}
return r;
}
static quo_rem barrett_div(integer x, const integer &n){
integer r = generalized_inverse(n);
std::size_t s = 2 * (r.bit_num() - 1);
integer quo;
for(; ; ){
integer d = (x * r) >> s;
x -= n * d;
if(x >= n){
x -= n;
d += 1;
}
quo += d;
if(x < n){ break; }
}
quo_rem result;
result.quo = std::move(quo);
result.quo.sign = x.sign * n.sign;
result.rem = std::move(x);
return result;
}
static integer pow(const integer &x, const integer &y){
integer z = x;
for(std::size_t j = y.bit_num() - 2; ; ){
z = z * z;
if(y.nth_bit(j) == 1){
z = z * x;
}
if(j == 0){
break;
}else{
--j;
}
}
return z;
}
static DoubleInt mod(DoubleInt x, DoubleInt y){
x = x % y;
return x >= 0 ? x : (y + x);
}
static UInt ceil_pow2_single(UInt x){
std::size_t n = 0, m = 0;
while(x > 0){
if((x & 1) > 0){
++n;
}
++m;
x >>= 1;
}
if(n == 1){
return static_cast<UInt>(1 << (m - 1));
}
return static_cast<UInt>(1 << m);
}
void elemental_shift(std::size_t n){
data.insert(data.begin(), n, 0);
}
std::size_t bit_num() const{
return bit_num(data.back()) + BitNum * (data.size() - 1);
}
static std::size_t bit_num(UInt x){
std::size_t n = 0;
while(x > 0){
++n;
x >>= 1;
}
return n;
}
UInt nth_bit(std::size_t i) const{
return (data[i / BitNum] >> (i % BitNum)) & 1;
}
std::size_t finite_bit_shift_lsr(std::size_t n){
std::size_t
digit = n / static_cast<std::size_t>(BitNum),
shift = n % static_cast<std::size_t>(BitNum),
rev_shift = static_cast<std::size_t>(BitNum) - shift;
UInt c = 0;
for(std::size_t i = 0; i < data.size(); ++i){
UInt x = data[i];
data[i] = (x << shift) | c;
if(rev_shift < BitNum){
c = x >> rev_shift;
}else{
c = 0;
}
}
if(c > 0){ data.push_back(c); }
return digit;
}
void bit_shift_lsr(std::size_t n){
elemental_shift(finite_bit_shift_lsr(n));
}
integer operator <<(const std::size_t n) const{
integer r = *this;
r.bit_shift_lsr(n);
return r;
}
integer &operator <<=(const std::size_t n){
bit_shift_lsr(n);
return *this;
}
void bit_shift_rsl(std::size_t n){
std::size_t
digit = n / static_cast<std::size_t>(BitNum),
shift = n % static_cast<std::size_t>(BitNum),
rev_shift = BitNum - shift,
size = data.size();
if(digit > 0){
for(std::size_t i = 0, length = size - digit; i < length; ++i){
data[i] = data[i + digit];
}
for(std::size_t i = 0; i < digit; ++i){
data[size - i - 1] = 0;
}
}
UInt c = 0;
for(std::size_t i = 0; i < data.size(); ++i){
std::size_t j = size - i - 1;
UInt x = data[j];
data[j] = (x >> shift) | c;
if(rev_shift < BitNum){
c = x << rev_shift;
}else{
c = 0;
}
}
normalize_data_size();
}
integer operator >>(const std::size_t n) const{
integer r = *this;
r.bit_shift_rsl(n);
return r;
}
integer &operator >>=(const std::size_t n){
bit_shift_rsl(n);
return *this;
}
void normalize_data_size(){
if(data.empty()){
sign = +1;
return;
}
std::size_t i = data.size() - 1;
for(; ; ){
if(i == 0 && data.front() == 0){
data.clear();
sign = +1;
return;
}
if(data[i] == 0){
--i;
}else{
break;
}
}
data.resize(i + 1);
}
static typename data_type::const_iterator normalize_range(typename data_type::const_iterator first, typename data_type::const_iterator last){
do{
--last;
if(last == first){
return last;
}
}while(*last == 0);
return last + 1;
}
#define MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(assign_op, op, type) \
integer &operator assign_op(const type &rhs){ *this = *this op integer(rhs); return *this; }
MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(+=, +, int);
MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(-=, -, int);
MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(*=, *, int);
MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(/=, /, int);
MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(%=, %, int);
MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(+=, +, unsigned int);
MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(-=, -, unsigned int);
MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(*=, *, unsigned int);
MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(/=, /, unsigned int);
MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(%=, %, unsigned int);
};
#define MULTI_PRECISION_COMPARE_OPERATOR(op, type) \
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum> \
bool operator op(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const type &rhs){ \
return lhs op integer<UInt, DoubleUInt, DoubleInt, BitNum>(rhs); \
} \
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum> \
bool operator op(const type &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){ \
return integer<UInt, DoubleUInt, DoubleInt, BitNum>(lhs) op rhs; \
}
MULTI_PRECISION_COMPARE_OPERATOR(==, int);
//MULTI_PRECISION_COMPARE_OPERATOR(!=, int);
MULTI_PRECISION_COMPARE_OPERATOR(<=, int);
MULTI_PRECISION_COMPARE_OPERATOR(>=, int);
MULTI_PRECISION_COMPARE_OPERATOR(<, int);
MULTI_PRECISION_COMPARE_OPERATOR(>, int);
MULTI_PRECISION_COMPARE_OPERATOR(==, unsigned int);
//MULTI_PRECISION_COMPARE_OPERATOR(!=, unsigned int);
MULTI_PRECISION_COMPARE_OPERATOR(<=, unsigned int);
MULTI_PRECISION_COMPARE_OPERATOR(>=, unsigned int);
MULTI_PRECISION_COMPARE_OPERATOR(<, unsigned int);
MULTI_PRECISION_COMPARE_OPERATOR(>, unsigned int);
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> operator +(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
integer<UInt, DoubleUInt, DoubleInt, BitNum> t(lhs);
t.add(rhs);
return t;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> operator -(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
integer<UInt, DoubleUInt, DoubleInt, BitNum> t(lhs);
t.sub(rhs);
return t;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> operator *(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
return integer<UInt, DoubleUInt, DoubleInt, BitNum>::mul(lhs, rhs);
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> operator /(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
return integer<UInt, DoubleUInt, DoubleInt, BitNum>::div(lhs, rhs).quo;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> operator %(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
return integer<UInt, DoubleUInt, DoubleInt, BitNum>::div(lhs, rhs).rem;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> operator <<(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, std::size_t n){
integer<UInt, DoubleUInt, DoubleInt, BitNum> x(lhs);
x.bit_shift_lsr(n);
return x;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> operator >>(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, std::size_t n){
integer<UInt, DoubleUInt, DoubleInt, BitNum> x(lhs);
x.bit_shift_rsl(n);
return x;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> &operator +=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
lhs.add(rhs);
return lhs;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> &operator -=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
lhs.sub(rhs);
return lhs;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> &operator *=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
lhs = integer<UInt, DoubleUInt, DoubleInt, BitNum>::mul(lhs, rhs);
return lhs;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> &operator /=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
lhs = integer<UInt, DoubleUInt, DoubleInt, BitNum>::div(lhs, rhs).quo;
return lhs;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> &operator %=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
lhs = integer<UInt, DoubleUInt, DoubleInt, BitNum>::div(lhs, rhs).rem;
return lhs;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> &operator <<=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, std::size_t n){
lhs.bit_shift_lsr(n);
return lhs;
}
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
integer<UInt, DoubleUInt, DoubleInt, BitNum> &operator >>=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, std::size_t n){
lhs.bit_shift_rsl(n);
return lhs;
}
#define MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(op, type) \
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum> \
integer<UInt, DoubleUInt, DoubleInt, BitNum> operator op(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const type &rhs){ return lhs op integer<UInt, DoubleUInt, DoubleInt, BitNum>(rhs); } \
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum> \
integer<UInt, DoubleUInt, DoubleInt, BitNum> operator op(const type &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){ return integer<UInt, DoubleUInt, DoubleInt, BitNum>(lhs) op rhs; }
MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(+, int);
MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(-, int);
MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(*, int);
MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(/, int);
MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(%, int);
MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(+, unsigned int);
MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(-, unsigned int);
MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(*, unsigned int);
MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(/, unsigned int);
MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(%, unsigned int);
template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
std::ostream &operator <<(
std::ostream &os,
integer<UInt, DoubleUInt, DoubleInt, BitNum> rhs
){
using integer = integer<UInt, DoubleUInt, DoubleInt, BitNum>;
std::vector<std::string> rseq;
integer lo(10);
for(; rhs.data.size() > 0; rhs /= lo){
integer temp = rhs % lo;
if(temp != 0){
rseq.push_back(std::to_string(temp.data[0]));
}else{
rseq.push_back(std::to_string(0));
}
}
os << (rhs.sign > 0 ? "" : "-");
if(!rseq.empty()){
for(auto iter = rseq.rbegin(); iter != rseq.rend(); ++iter){
os << *iter;
}
}else{
os << "0";
}
return os;
}
}
#endif // MULTI_PRECISION_INCLUDE_INTEGER_HPP
Raw
rational.hpp
#ifndef RATIONAL_HPP
#define RATIONAL_HPP
#include <iostream> // for std::istream and std::ostream
#include <ios> // for std::noskipws
#include <stdexcept> // for std::domain_error
#include <string> // for std::string implicit constructor
#include <cstdlib> // for std::abs
#include <limits> // for std::numeric_limits
#include <cstdint> // for std::int64_t
#include "integer.hpp"
//----------------
// rational (from boost 1.53)
namespace rational_impl{
using integer = multi_precision::integer<>;
class bad_rational : public std::domain_error{
public:
explicit bad_rational() : std::domain_error("bad rational: zero denominator") {}
};
template<class Signature>
inline integer gcd_binary(integer u, integer v){
if(u && v){
std::size_t shifts = 0;
integer one = 1;
while(!(u & one) && !(v & one)){
++shifts;
u >>= 1;
v >>= 1;
}
integer r[] = { u, v };
integer which = static_cast<bool>(u & one);
do{
while(!(r[which] & 1u)){
r[which] >>= 1;
}
if(r[!which] > r[which]){
which ^= 1u;
}
r[which] -= r[!which];
}while(r[which]);
return r[!which] << shifts;
}else{
return u + v;
}
}
template<bool ImplicitExponentialNumerator = true>
class rational{
private:
struct helper{ integer parts[2]; };
typedef integer (helper::* bool_type)[2];
public:
inline rational() : num(0), den(1){}
inline rational(const rational &other) : num(other.num), den(other.den){}
inline rational(rational &&other) : num(std::move(other.num)), den(std::move(other.den)){}
inline rational(const integer &n) : num(n), den(1){}
inline rational(const integer &n, const integer &d) : num(n), den(d){}
inline rational &operator =(const integer &n){ return assign(n, 1); }
rational &assign(const integer &n, const integer &d);
inline const integer &numerator () const{ return num; }
inline const integer &denominator () const{ return den; }
rational &operator +=(const rational &r);
rational &operator -=(const rational &r);
rational &operator *=(const rational &r);
rational &operator /=(const rational &r);
rational &operator +=(const integer &i);
rational &operator -=(const integer &i);
rational &operator *=(const integer &i);
rational &operator /=(const integer &i);
rational operator +(const rational &r) const;
rational operator -(const rational &r) const;
rational operator *(const rational &r) const;
rational operator /(const rational &r) const;
const rational &operator ++();
const rational &operator --();
bool operator !() const{ return !num; }
bool operator <(const rational &r) const;
bool operator >(const rational &r) const;
bool operator ==(const rational &r) const;
bool operator !=(const rational &r) const;
bool operator <=(const rational &r) const;
bool operator >=(const rational &r) const;
bool operator <(const integer &i) const;
bool operator >(const integer &i) const;
bool operator ==(const integer &i) const;
void normalize();
private:
integer num, den;
bool test_invariant() const;
};
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator> &rational<ImplicitExponentialNumerator>::assign(const integer &n, const integer &d){
num = n;
den = d;
if(ImplicitExponentialNumerator){
normalize();
}
return *this;
}
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator> operator +(const rational<ImplicitExponentialNumerator> &r){
return r;
}
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator> operator -(const rational<ImplicitExponentialNumerator> &r){
return rational<ImplicitExponentialNumerator>(-r.numerator(), r.denominator());
}
template<bool ImplicitExponentialNumerator>
rational<ImplicitExponentialNumerator> &rational<ImplicitExponentialNumerator>::operator+= (const rational<ImplicitExponentialNumerator> &r){
integer r_num = r.num;
integer r_den = r.den;
integer g = gcd_binary<void>(den, r_den);
den /= g;
num = num * (r_den / g) + r_num * den;
g = gcd_binary<void>(num, g);
num /= g;
den *= r_den/g;
return *this;
}
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator> rational<ImplicitExponentialNumerator>::operator+ (const rational<ImplicitExponentialNumerator> &r) const{
rational a = *this;
a += r;
return a;
}
template<bool ImplicitExponentialNumerator>
rational<ImplicitExponentialNumerator> &rational<ImplicitExponentialNumerator>::operator-= (const rational<ImplicitExponentialNumerator> &r){
integer r_num = r.num;
integer r_den = r.den;
integer g = gcd_binary<void>(den, r_den);
den /= g;
num = num * (r_den / g) - r_num * den;
g = gcd_binary<void>(num, g);
num /= g;
den *= r_den/g;
return *this;
}
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator> rational<ImplicitExponentialNumerator>::operator- (const rational<ImplicitExponentialNumerator> &r) const{
rational a = *this;
a -= r;
return a;
}
template<bool ImplicitExponentialNumerator>
rational<ImplicitExponentialNumerator> &rational<ImplicitExponentialNumerator>::operator*= (const rational<ImplicitExponentialNumerator> &r){
integer r_num = r.num;
integer r_den = r.den;
integer gcd1 = gcd_binary<void>(num, r_den);
integer gcd2 = gcd_binary<void>(r_num, den);
num = (num/gcd1) * (r_num/gcd2);
den = (den/gcd2) * (r_den/gcd1);
return *this;
}
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator> rational<ImplicitExponentialNumerator>::operator* (const rational<ImplicitExponentialNumerator> &r) const{
rational a = *this;
a *= r;
return a;
}
template<bool ImplicitExponentialNumerator>
rational<ImplicitExponentialNumerator> &rational<ImplicitExponentialNumerator>::operator/= (const rational<ImplicitExponentialNumerator> &r){
integer r_num = r.num;
integer r_den = r.den;
integer zero(0);
if(r_num == zero)
throw bad_rational();
if(num == zero)
return *this;
integer gcd1 = gcd_binary<void>(num, r_num);
integer gcd2 = gcd_binary<void>(r_den, den);
num = (num/gcd1) * (r_den/gcd2);
den = (den/gcd2) * (r_num/gcd1);
if(den < zero){
num = -num;
den = -den;
}
return *this;
}
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator> rational<ImplicitExponentialNumerator>::operator/ (const rational<ImplicitExponentialNumerator> &r) const{
rational a = *this;
a += r;
return a;
}
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator>& rational<ImplicitExponentialNumerator>::operator+= (const integer &i){
return operator+= (rational<ImplicitExponentialNumerator>(i));
}
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator>& rational<ImplicitExponentialNumerator>::operator -=(const integer &i){
return operator-= (rational<ImplicitExponentialNumerator>(i));
}
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator>& rational<ImplicitExponentialNumerator>::operator*= (const integer &i){
return operator*= (rational<ImplicitExponentialNumerator>(i));
}
template<bool ImplicitExponentialNumerator>
inline rational<ImplicitExponentialNumerator>& rational<ImplicitExponentialNumerator>::operator/= (const integer &i){
return operator/= (rational<ImplicitExponentialNumerator>(i));
}
template<bool ImplicitExponentialNumerator>
inline const rational<ImplicitExponentialNumerator> &rational<ImplicitExponentialNumerator>::operator ++(){
num += den;
return *this;
}
template<bool ImplicitExponentialNumerator>
inline const rational<ImplicitExponentialNumerator> &rational<ImplicitExponentialNumerator>::operator --(){
num -= den;
return *this;
}
template<bool ImplicitExponentialNumerator>
bool rational<ImplicitExponentialNumerator>::operator <(const rational<ImplicitExponentialNumerator> &r) const{
integer const zero(0);
struct { integer n, d, q, r; } ts = { this->num, this->den, this->num / this->den, this->num % this->den }, rs = { r.num, r.den, r.num / r.den, r.num % r.den };
unsigned reverse = 0u;
while(ts.r < zero){ ts.r += ts.d; --ts.q; }
while(rs.r < zero){ rs.r += rs.d; --rs.q; }
while(true){
if(ts.q != rs.q){
return reverse ? ts.q > rs.q : ts.q < rs.q;
}
reverse ^= 1u;
if((ts.r == zero) || (rs.r == zero)){
break;
}
ts.n = ts.d; ts.d = ts.r;
ts.q = ts.n / ts.d; ts.r = ts.n % ts.d;
rs.n = rs.d; rs.d = rs.r;
rs.q = rs.n / rs.d; rs.r = rs.n % rs.d;
}
if(ts.r == rs.r){
return false;
}else{
return (ts.r != zero) != (reverse != 0);
}
}
template<bool ImplicitExponentialNumerator>
bool rational<ImplicitExponentialNumerator>::operator >(const rational<ImplicitExponentialNumerator> &r) const{
return r.operator <(*this);
}
template<bool ImplicitExponentialNumerator>
bool rational<ImplicitExponentialNumerator>::operator <=(const rational<ImplicitExponentialNumerator> &r) const{
integer const zero(0);
struct { integer n, d, q, r; } ts = { this->num, this->den, this->num / this->den, this->num % this->den }, rs = { r.num, r.den, r.num / r.den, r.num % r.den };
unsigned reverse = 0u;
while(ts.r < zero){ ts.r += ts.d; --ts.q; }
while(rs.r < zero){ rs.r += rs.d; --rs.q; }
while(true){
if(ts.q != rs.q){
return reverse ? ts.q > rs.q : ts.q < rs.q;
}
reverse ^= 1u;
if((ts.r == zero) || (rs.r == zero)){
break;
}
ts.n = ts.d; ts.d = ts.r;
ts.q = ts.n / ts.d; ts.r = ts.n % ts.d;
rs.n = rs.d; rs.d = rs.r;
rs.q = rs.n / rs.d; rs.r = rs.n % rs.d;
}
if(ts.r == rs.r){
return true;
}else{
return (ts.r != zero) != (reverse != 0);
}
}
template<bool ImplicitExponentialNumerator>
bool rational<ImplicitExponentialNumerator>::operator >=(const rational<ImplicitExponentialNumerator> &r) const{
return r.operator <=(*this);
}
template<bool ImplicitExponentialNumerator>
bool rational<ImplicitExponentialNumerator>::operator <(const integer &i) const{
integer const zero(0);
integer q = this->num / this->den, r = this->num % this->den;
while(r < zero){ r += this->den; --q; }
return q < i;
}
template<bool ImplicitExponentialNumerator>
bool rational<ImplicitExponentialNumerator>::operator> (const integer &i) const{
if(num == i && den == integer(1))
return false;
return !operator <(i);
}
template<bool ImplicitExponentialNumerator>
inline bool rational<ImplicitExponentialNumerator>::operator ==(const rational<ImplicitExponentialNumerator> &r) const{
return ((num == r.num) && (den == r.den));
}
template<bool ImplicitExponentialNumerator>
inline bool rational<ImplicitExponentialNumerator>::operator !=(const rational<ImplicitExponentialNumerator> &r) const{
return !(*this == r);
}
template<bool ImplicitExponentialNumerator>
inline bool rational<ImplicitExponentialNumerator>::operator ==(const integer &i) const{
return ((den == integer(1)) && (num == i));
}
template<bool ImplicitExponentialNumerator>
inline bool rational<ImplicitExponentialNumerator>::test_invariant() const{
return (this->den > integer(0)) && (gcd_binary<void>(this->num, this->den) ==
integer(1));
}
template<bool ImplicitExponentialNumerator>
void rational<ImplicitExponentialNumerator>::normalize(){
integer zero(0);
if(den == zero)
throw bad_rational();
if(num == zero){
den = integer(1);
return;
}
integer g = gcd_binary<void>(num, den);
num /= g;
den /= g;
if(den < zero){
num = -num;
den = -den;
}
}
template<bool ImplicitExponentialNumerator>
std::ostream& operator<< (std::ostream& os, const rational<ImplicitExponentialNumerator> &r){
os << r.numerator() << '/' << r.denominator();
return os;
}
} // namespace rational_impl
typedef rational_impl::rational<true> rational;
typedef rational_impl::rational<false> explicit_exponential_rational;
template<class T, bool Cond>
inline T rational_cast(const rational_impl::rational<Cond> &src){
return static_cast<T>(src.numerator()) / static_cast<T>(src.denominator());
}
template<bool Cond>
inline rational_impl::rational<Cond> abs(const rational_impl::rational<Cond> &r){
if(r.numerator() >= multi_precision::integer<>(0)){
return r;
}
return rational_impl::rational<Cond>(-r.numerator(), r.denominator());
}
#endif // RATIONAL_HPP
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment