Skip to content

Commit 5da17f1

Browse files
Prok-AlProk-Al
authored
Merge pull request #26 from PeeachPie/feature/continued-fraction
Цепные дроби 2
2 parents 1d9a716 + c4ff178 commit 5da17f1

File tree

12 files changed

+913
-91
lines changed

12 files changed

+913
-91
lines changed

‎src/modules/continued_fraction/continued_fraction.cpp‎

Lines changed: 227 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include <stdexcept>
44
#include <cmath>
55
#include <set>
6+
#include <tuple>
67

78

89
ContinuedFraction::ContinuedFraction()
@@ -52,7 +53,7 @@ ContinuedFraction ContinuedFraction::FROM_Q_CF(const Rational& rational) {
5253

5354
// p = q, q = остаток
5455
p = Integer(q.as_string());
55-
q = Natural(remainder.ABS_Z_N().as_string());
56+
q = remainder.ABS_Z_N();
5657
}
5758

5859
return ContinuedFraction(coeffs);
@@ -221,29 +222,50 @@ Rational ContinuedFraction::APPROX_CF_Q(const Natural& max_denominator) const {
221222

222223
// CF-6 | Сравнение двух цепных дробей
223224
int ContinuedFraction::COM_CF_D(const ContinuedFraction& other) const {
224-
size_t min_len = std::min(coefficients_.size(), other.coefficients_.size());
225+
// нормализация: убираем завершающие 1 ([..., a, 1] -> [..., a+1])
226+
auto normalize = [](std::vector<Integer> coeffs) {
227+
while (coeffs.size() > 1 && coeffs.back().COM_ZZ_D(Integer("1")) == 0) {
228+
coeffs.pop_back();
229+
coeffs.back() = coeffs.back().ADD_ZZ_Z(Integer("1"));
230+
}
231+
return coeffs;
232+
};
233+
234+
std::vector<Integer> a = (period_start_ < 0) ? normalize(coefficients_) : coefficients_;
235+
std::vector<Integer> b = (other.period_start_ < 0) ? normalize(other.coefficients_) : other.coefficients_;
236+
237+
size_t min_len = std::min(a.size(), b.size());
225238

226239
for (size_t i = 0; i < min_len; ++i) {
227-
int cmp = coefficients_[i].COM_ZZ_D(other.coefficients_[i]);
240+
int cmp = a[i].COM_ZZ_D(b[i]);
228241
if (cmp != 0) {
229242
// Для чётных позиций: больший коэффициент = большая дробь
230243
// Для нечётных позиций: больший коэффициент = меньшая дробь
244+
// COM_ZZ_D возвращает: 1 если this > other, -1 если this < other
231245
if (i % 2 == 0) {
232-
return (cmp == 2) ? 1 : -1;
246+
return cmp; // на чётной позиции порядок прямой
233247
} else {
234-
return (cmp == 2) ? -1 : 1;
248+
return -cmp; // на нечётной позиции порядок обратный
235249
}
236250
}
237251
}
238252

239-
if (coefficients_.size() == other.coefficients_.size()) {
253+
if (a.size() == b.size()) {
254+
// Для периодических дробей также сравниваем начало периода
255+
if (period_start_ >= 0 && other.period_start_ >= 0) {
256+
if (period_start_ != other.period_start_) {
257+
return (period_start_ < other.period_start_) ? -1 : 1;
258+
}
259+
}
240260
return 0;
241261
}
242262

243-
if (coefficients_.size() > other.coefficients_.size()) {
244-
return (min_len % 2 == 0) ? 1 : -1;
245-
} else {
263+
// Более короткая дробь эквивалентна дроби с "бесконечностью" на позиции min_len
264+
if (a.size() > b.size()) {
246265
return (min_len % 2 == 0) ? -1 : 1;
266+
} else {
267+
// other длиннее, this короче (this имеет "∞" на позиции min_len)
268+
return (min_len % 2 == 0) ? 1 : -1;
247269
}
248270
}
249271

@@ -464,6 +486,202 @@ std::vector<Integer> ContinuedFraction::GET_APERIODIC_CF() const {
464486
return aperiodic;
465487
}
466488

489+
// CF-13 | Преобразование периодической цепной дроби в квадратичную иррациональность
490+
// Возвращает (a, D, c) такие что дробь = (a + sqrt(D)) / c
491+
std::tuple<Integer, Natural, Integer> ContinuedFraction::TO_CF_QUAD() const {
492+
// Если дробь не периодическая, возвращаем рациональное число
493+
if (!IS_PERIODIC_CF()) {
494+
Rational r = TO_CF_Q();
495+
return {r.numerator(), Natural("0"), Integer(r.denominator().as_string())};
496+
}
497+
498+
// Получаем период и непериодическую часть
499+
std::vector<Integer> period = GET_PERIOD_CF();
500+
std::vector<Integer> aperiodic = GET_APERIODIC_CF();
501+
502+
if (period.empty()) {
503+
Rational r = TO_CF_Q();
504+
return {r.numerator(), Natural("0"), Integer(r.denominator().as_string())};
505+
}
506+
507+
// Шаг 1: Вычисляем конвергенты периода [b1; b2, ..., bk]
508+
// Для уравнения x = (p_k + p_{k-1} * x) / (q_k + q_{k-1} * x)
509+
// где x = [b1; b2, ..., bk, x] - чисто периодическая часть
510+
511+
Integer p_prev = Integer("1"); // p_{-1} = 1
512+
Integer q_prev = Integer("0"); // q_{-1} = 0
513+
Integer p_curr = period[0]; // p_0 = b_1
514+
Integer q_curr = Integer("1"); // q_0 = 1
515+
516+
for (size_t i = 1; i < period.size(); ++i) {
517+
Integer p_new = period[i].MUL_ZZ_Z(p_curr).ADD_ZZ_Z(p_prev);
518+
Integer q_new = period[i].MUL_ZZ_Z(q_curr).ADD_ZZ_Z(q_prev);
519+
520+
p_prev = p_curr;
521+
q_prev = q_curr;
522+
p_curr = p_new;
523+
q_curr = q_new;
524+
}
525+
526+
// p_k = p_curr, p_{k-1} = p_prev
527+
// q_k = q_curr, q_{k-1} = q_prev
528+
529+
// Шаг 2: Решаем квадратное уравнение для чисто периодической части
530+
// x = [b1; b2, ..., bk, x]
531+
// x = b1 + 1/(b2 + 1/(...+ 1/(bk + 1/x)))
532+
// Используя конвергенты: x = (p_k * x + p_{k-1}) / (q_k * x + q_{k-1})
533+
// => x * (q_k * x + q_{k-1}) = p_k * x + p_{k-1}
534+
// => q_k * x^2 + q_{k-1} * x = p_k * x + p_{k-1}
535+
// => q_k * x^2 + (q_{k-1} - p_k) * x - p_{k-1} = 0
536+
537+
Integer coef_a = q_curr; // q_k
538+
Integer coef_b = q_prev.SUB_ZZ_Z(p_curr); // q_{k-1} - p_k
539+
Integer coef_c = p_prev.MUL_ZM_Z(); // -p_{k-1}
540+
541+
// D = b^2 - 4ac = (q_{k-1} - p_k)^2 + 4 * q_k * p_{k-1}
542+
Integer b_sq = coef_b.MUL_ZZ_Z(coef_b);
543+
Integer four = Integer("4");
544+
Integer four_ac = four.MUL_ZZ_Z(coef_a).MUL_ZZ_Z(p_prev);
545+
Integer D_int = b_sq.ADD_ZZ_Z(four_ac);
546+
547+
// D должен быть неотрицательным
548+
if (D_int.SGN_Z_D() < 0) {
549+
Rational r = TO_CF_Q();
550+
return {r.numerator(), Natural("0"), Integer(r.denominator().as_string())};
551+
}
552+
553+
Natural D_nat = D_int.ABS_Z_N();
554+
555+
// x = (-b + sqrt(D)) / (2a) = (p_k - q_{k-1} + sqrt(D)) / (2 * q_k)
556+
Integer x_num = coef_b.MUL_ZM_Z(); // -coef_b = p_k - q_{k-1}
557+
Integer x_den = Integer("2").MUL_ZZ_Z(coef_a); // 2 * q_k
558+
559+
// Если нет непериодической части, результат = x
560+
if (aperiodic.empty()) {
561+
// Нормализуем знак знаменателя
562+
if (x_den.SGN_Z_D() < 0) {
563+
x_num = x_num.MUL_ZM_Z();
564+
x_den = x_den.MUL_ZM_Z();
565+
}
566+
return {x_num, D_nat, x_den};
567+
}
568+
569+
// Шаг 3: Учитываем непериодическую часть
570+
// result = [a0; a1, ..., a_{m-1}, x]
571+
// result = (P_{m-1} * x + P_{m-2}) / (Q_{m-1} * x + Q_{m-2})
572+
// где P_i/Q_i - конвергенты непериодической части
573+
574+
Integer P_prev = Integer("1"); // P_{-1} = 1
575+
Integer Q_prev = Integer("0"); // Q_{-1} = 0
576+
Integer P_curr = aperiodic[0]; // P_0 = a_0
577+
Integer Q_curr = Integer("1"); // Q_0 = 1
578+
579+
for (size_t i = 1; i < aperiodic.size(); ++i) {
580+
Integer P_new = aperiodic[i].MUL_ZZ_Z(P_curr).ADD_ZZ_Z(P_prev);
581+
Integer Q_new = aperiodic[i].MUL_ZZ_Z(Q_curr).ADD_ZZ_Z(Q_prev);
582+
583+
P_prev = P_curr;
584+
Q_prev = Q_curr;
585+
P_curr = P_new;
586+
Q_curr = Q_new;
587+
}
588+
589+
// P_{m-1} = P_curr, P_{m-2} = P_prev
590+
// Q_{m-1} = Q_curr, Q_{m-2} = Q_prev
591+
592+
// result = (P_{m-1} * x + P_{m-2}) / (Q_{m-1} * x + Q_{m-2})
593+
// где x = (x_num + sqrt(D)) / x_den
594+
//
595+
// Подставляем:
596+
// числитель = P_{m-1} * (x_num + sqrt(D)) / x_den + P_{m-2}
597+
// = (P_{m-1} * x_num + P_{m-2} * x_den + P_{m-1} * sqrt(D)) / x_den
598+
// = (A + B * sqrt(D)) / x_den
599+
// где A = P_{m-1} * x_num + P_{m-2} * x_den, B = P_{m-1}
600+
//
601+
// знаменатель = Q_{m-1} * (x_num + sqrt(D)) / x_den + Q_{m-2}
602+
// = (Q_{m-1} * x_num + Q_{m-2} * x_den + Q_{m-1} * sqrt(D)) / x_den
603+
// = (C + E * sqrt(D)) / x_den
604+
// где C = Q_{m-1} * x_num + Q_{m-2} * x_den, E = Q_{m-1}
605+
//
606+
// result = (A + B * sqrt(D)) / (C + E * sqrt(D))
607+
//
608+
// Рационализируем знаменатель:
609+
// result = (A + B * sqrt(D)) * (C - E * sqrt(D)) / ((C + E * sqrt(D)) * (C - E * sqrt(D)))
610+
// = (A*C - B*E*D + (B*C - A*E) * sqrt(D)) / (C^2 - E^2 * D)
611+
612+
Integer A = P_curr.MUL_ZZ_Z(x_num).ADD_ZZ_Z(P_prev.MUL_ZZ_Z(x_den));
613+
Integer B = P_curr;
614+
Integer C = Q_curr.MUL_ZZ_Z(x_num).ADD_ZZ_Z(Q_prev.MUL_ZZ_Z(x_den));
615+
Integer E = Q_curr;
616+
617+
// числитель: (A*C - B*E*D) + (B*C - A*E) * sqrt(D)
618+
Integer D_as_int = Integer(D_nat.as_string());
619+
Integer AC = A.MUL_ZZ_Z(C);
620+
Integer BED = B.MUL_ZZ_Z(E).MUL_ZZ_Z(D_as_int);
621+
Integer result_const = AC.SUB_ZZ_Z(BED); // A*C - B*E*D
622+
623+
Integer BC = B.MUL_ZZ_Z(C);
624+
Integer AE = A.MUL_ZZ_Z(E);
625+
Integer result_sqrt_coef = BC.SUB_ZZ_Z(AE); // B*C - A*E
626+
627+
// знаменатель: C^2 - E^2 * D
628+
Integer C_sq = C.MUL_ZZ_Z(C);
629+
Integer E_sq_D = E.MUL_ZZ_Z(E).MUL_ZZ_Z(D_as_int);
630+
Integer result_den = C_sq.SUB_ZZ_Z(E_sq_D);
631+
632+
// Результат: (result_const + result_sqrt_coef * sqrt(D)) / result_den
633+
634+
// Упрощаем: делим всё на result_sqrt_coef (если он не 0)
635+
if (result_sqrt_coef.SGN_Z_D() == 0) {
636+
// Нет sqrt(D) в результате - это рациональное число
637+
if (result_den.SGN_Z_D() < 0) {
638+
result_const = result_const.MUL_ZM_Z();
639+
result_den = result_den.MUL_ZM_Z();
640+
}
641+
// Сокращаем
642+
Natural gcd = result_const.ABS_Z_N().GCF_NN_N(result_den.ABS_Z_N());
643+
if (gcd.NZER_N_B() && gcd.COM_NN_D(Natural("1")) != 0) {
644+
result_const = result_const.DIV_ZZ_Z(Integer(gcd.as_string()));
645+
result_den = result_den.DIV_ZZ_Z(Integer(gcd.as_string()));
646+
}
647+
return {result_const, Natural("0"), result_den};
648+
}
649+
650+
// Находим НОД всех трёх коэффициентов
651+
Natural abs_const = result_const.ABS_Z_N();
652+
Natural abs_coef = result_sqrt_coef.ABS_Z_N();
653+
Natural abs_den = result_den.ABS_Z_N();
654+
655+
Natural gcd = abs_const.GCF_NN_N(abs_coef);
656+
gcd = gcd.GCF_NN_N(abs_den);
657+
658+
if (gcd.NZER_N_B() && gcd.COM_NN_D(Natural("1")) != 0) {
659+
result_const = result_const.DIV_ZZ_Z(Integer(gcd.as_string()));
660+
result_sqrt_coef = result_sqrt_coef.DIV_ZZ_Z(Integer(gcd.as_string()));
661+
result_den = result_den.DIV_ZZ_Z(Integer(gcd.as_string()));
662+
}
663+
664+
// Нормализуем знак знаменателя
665+
if (result_den.SGN_Z_D() < 0) {
666+
result_const = result_const.MUL_ZM_Z();
667+
result_sqrt_coef = result_sqrt_coef.MUL_ZM_Z();
668+
result_den = result_den.MUL_ZM_Z();
669+
}
670+
671+
// Если коэффициент при sqrt(D) отрицательный, меняем знаки
672+
if (result_sqrt_coef.SGN_Z_D() < 0) {
673+
result_const = result_const.MUL_ZM_Z();
674+
result_sqrt_coef = result_sqrt_coef.MUL_ZM_Z();
675+
result_den = result_den.MUL_ZM_Z();
676+
}
677+
678+
// D_new = k^2 * D, где k = result_sqrt_coef
679+
Natural k_nat = result_sqrt_coef.ABS_Z_N();
680+
Natural D_new = k_nat.MUL_NN_N(k_nat).MUL_NN_N(D_nat);
681+
682+
return {result_const, D_new, result_den};
683+
}
684+
467685
// Строковое представление
468686
std::string ContinuedFraction::as_string() const {
469687
if (coefficients_.empty()) {

‎src/modules/continued_fraction/continued_fraction.hpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
#include <vector>
66
#include <string>
77
#include <optional>
8+
#include <tuple>
89

910
/*
1011
* Класс цепных дробей для рациональных чисел
@@ -72,6 +73,11 @@ class ContinuedFraction {
7273

7374
// CF-12 | Получение непериодической части
7475
[[nodiscard]] std::vector<Integer> GET_APERIODIC_CF() const;
76+
77+
// CF-13 | Преобразование периодической цепной дроби в квадратичную иррациональность
78+
// Возвращает (a, D, c) такие что дробь = (a + sqrt(D)) / c
79+
// Если дробь не периодическая, возвращает рациональное число как (p, 0, q) где p/q - значение
80+
[[nodiscard]] std::tuple<Integer, Natural, Integer> TO_CF_QUAD() const;
7581

7682
// Геттеры
7783
[[nodiscard]] const std::vector<Integer>& coefficients() const { return coefficients_; }

0 commit comments

Comments
 (0)