プログラミング初心者でもペル方程式の最小解を計算したい - 音楽室と化学室と美術室とPC室の融合部屋　所謂自由室

平方数でない $D$ について
$x^2 - Dy^2 = \pm 1$
についての整数解 $(x,y)$ は必ず存在し，この解は連分数展開を用いて機械的に求めることができる．
自明な解は $(1,0)$ で，非自明な解はその他の整数解を指す．

#include <iostream>
#include <string>
#include <cmath>
#include <vector>
#include <boost/multiprecision/cpp_int.hpp>

#define MYINT boost::multiprecision::cpp_int

/************************************************************************************
数値入力部
************************************************************************************/
int input(void) {
int num;
//ちゃんと数値が入るまで
while (true) {
std::cout << "\nx^2 - by^2 = 1　のbの値を入力 > ";
//数値入ったかどうか
if (std::cin >> num) {
//正の数値のときにループから出る
if (num > 0) { break; }
else { std::cout << "正の整数を入力してください\n"; }
}
else {
//数字じゃなかった場合
std::cout << "数字じゃないじゃん！\n";
}
//正しく入力されなかった場合入力をリセット
std::cin.clear();
std::cin.ignore(100, '\n');
}
return num;
}

/************************************************************************************
x^2 - Ay^2 = 1の解を計算
アルゴリズムについては
ペル方程式
http://izumi-math.jp/sanae/MyText/number/number_05.pdf
を参照
************************************************************************************/
bool calc(std::vector ans, int num, MYINT *a, MYINT *b) {
MYINT pn, qn, pn1 = 1, qn1 = 0;
//周期計算
int loop, cycle = ans.size() - 1;
//周期が奇数であるとき = 1の解は2m-1計算しないといけない．偶数であるときはm-1でよい
loop = ((cycle % 2) == 1) ? 2 * cycle - 1 : cycle - 1;
for (int i = loop - 1; i >= 0; i--) {
pn = ans[i%cycle + 1] * pn1 + qn1;
qn = pn1;
pn1 = pn;
qn1 = qn;
}
//最後に整数部分
*a = ans[0] * pn1 + qn1;
*b = pn1;
//x^2 - Ay^2 = -1の解が存在するかどうか
return (cycle % 2) == 1;
}

/************************************************************************************
x^2 - Ay^2 = -1の解を計算
アルゴリズムについては
ペル方程式
http://izumi-math.jp/sanae/MyText/number/number_05.pdf
を参照
************************************************************************************/
void calc_other(std::vector ans, int num, MYINT *a, MYINT *b) {
MYINT pn, qn, pn1, qn1;
int loop, cycle = ans.size() - 1;
loop = cycle - 1;
pn1 = 1;
qn1 = 0;
for (int i = loop - 1; i >= 0; i--) {
pn = ans[i%cycle + 1] * pn1 + qn1;
qn = pn1;
pn1 = pn;
qn1 = qn;
}
*a = ans[0] * pn1 + qn1;
*b = pn1;
}

/************************************************************************************
/連分数展開
アルゴリズムについては
Lagrange による連分数展開のアルゴリズムの一般化の試み
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1199-22.pdf
を参照
************************************************************************************/
std::vector continued_fraction(int num) {
int Pk = 0, Pk1, Qk = 1, Qk1, ak, TP, TQ, rootnum;
std::vector ans;
bool f = true;
//整数部分を計算
rootnum = (int)sqrt(num);
ak = rootnum;
//数値保存
ans.push_back(ak);
//計算部
while (true) {
Pk1 = ak * Qk - Pk;
Qk1 = (num - Pk1 * Pk1) / Qk;
if (f) {
//循環する場合最初に戻るので最初の値を保存
TP = Pk1, TQ = Qk1;
f = false;
}
else {
//循環した場合ループから抜ける
if (TP == Pk1 && TQ == Qk1) break;
}
ak = (int)((Pk1 + rootnum) / Qk1);
Pk = Pk1;
Qk = Qk1;
//数値保存
ans.push_back(ak);
}
return ans;
}

//計算本体
void pell(int num) {
std::vector ans;
//連分数を計算
ans = continued_fraction(num);
MYINT pn, qn;
//ペル方程式の解を計算
bool other = calc(ans, num, &pn, &qn);
std::cout << "x^2 - " + std::to_string(num) + "y^2 = 1の非自明最小解(x,y)は\n(" + pn.str() + " , " + qn.str() + ")\nです" << std::endl;
//x^2 - Ay^2 = -1の解が存在する場合計算
if (other) {
calc_other(ans, num, &pn, &qn);
std::cout << "また，x^2 - " + std::to_string(num) + "y^2 = -1の非自明最小解(x,y)は\n(" + pn.str() + " , " + qn.str() + ")\nです" << std::endl;
}
}

//計算
void solution(int num) {
//平方根の整数部分
int c = (int)sqrt(num);
if (c*c == num) {
//入力値が平方数の場合
std::cout << "非自明解はありません" << std::endl;
} else {
//入力値が平方数でない場合
pell(num);
}
std::cout << "\n";
}

int main(void) {
while (true) {
int num = input();
if (num < 0)break;
solution(num);
}
return 0;
}

以下実行結果

x^2 - by^2 = 1　のbの値を入力 > 5
x^2 - 5y^2 = 1の非自明最小解(x,y)は
(9 , 4)
です
また，x^2 - 5y^2 = -1の非自明最小解(x,y)は
(2 , 1)
です

x^2 - by^2 = 1　のbの値を入力 > 114514
x^2 - 114514y^2 = 1の非自明最小解(x,y)は
(3058389164815894335086675882217709431950420307140756009821362546111334285928768064662409120517323199 , 9037815138660369922198555785216162916412331641365948545459353586895717702576049626533527779108680)
です

x^2 - by^2 = 1　のbの値を入力 > 334
x^2 - 334y^2 = 1の非自明最小解(x,y)は
(63804373719695 , 3491219999244)
です

x^2 - by^2 = 1　のbの値を入力 > 226
x^2 - 226y^2 = 1の非自明最小解(x,y)は
(451 , 30)
です
また，x^2 - 226y^2 = -1の非自明最小解(x,y)は
(15 , 1)
です

x^2 - by^2 = 1　のbの値を入力 > 485
x^2 - 485y^2 = 1の非自明最小解(x,y)は
(969 , 44)
です
また，x^2 - 485y^2 = -1の非自明最小解(x,y)は
(22 , 1)
です

x^2 - by^2 = 1　のbの値を入力 > 4649
x^2 - 4649y^2 = 1の非自明最小解(x,y)は
(9420195137863692801 , 138159300628153120)
です
また，x^2 - 4649y^2 = -1の非自明最小解(x,y)は
(2170275920 , 31829893)
です

x^2 - by^2 = 1　のbの値を入力 > 31415
x^2 - 31415y^2 = 1の非自明最小解(x,y)は
(3460685671 , 19525116)
です

x^2 - by^2 = 1　のbの値を入力 > 3141592
x^2 - 3141592y^2 = 1の非自明最小解(x,y)は
(1756577403460059873785963109897142608287956736268848789675455684752093429444073349144761124407145075120858664137427130562438631838919558791344584853540997297670271941098249246336600609712623050975568565075702908922524161480792156247335314319726249011086988377484344170180137699355300671038415360631964473689826373416427212991353645732880462725443628592488305697508079797981278132447858933249950819070987738570235904884557114824417277249478291790330076515459063777262004807039765026590467370947039575371713555438759329583205704602640485797 , 991042776817833278738581646937092191647848193349824120779873611896871122319520430520077057844928034595214923753466819580676669874741905302802504079975844026481601886524137653816046717133314294156801127677737933086175774330421477761764291656881581817410335592200457372287959215612365706928260924100158215328958698484052055807151645954316415388569900162977578210755097130815965504322338604299720082609596463775009954764162947334320725374116620651931874810653207281293318815223995642121493057852768934184209873625056855789617869150545593)
です

x^2 - by^2 = 1　のbの値を入力 > 1024
非自明解はありません

x^2 - by^2 = 1　のbの値を入力 > 2048
x^2 - 2048y^2 = 1の非自明最小解(x,y)は
(886731088897 , 19594173939)
です

x^2 - by^2 = 1　のbの値を入力 > 893
x^2 - 893y^2 = 1の非自明最小解(x,y)は
(6091434999 , 203842100)
です

x^2 - by^2 = 1　のbの値を入力 > 141421356
x^2 - 141421356y^2 = 1の非自明最小解(x,y)は
(4578929726620897898623228667405046142991532861062394252908117068128473249652338180999568586493350792881956799995596457063408973390356507946700539299479807186516914141263310105646285045628698802057346685151171364668675187454463238036251592004478870187200786413954913902491773456399178783768068639976888009439812338935518927835030073134124637121796871119325785965201499843243697929219657053963128787552664924867736924364717256506938246047773265483952892867549237886852160984188773100343923911076416229567699943726807919493670981396285213132706602279346175852428974632053163823308285562470831370437610839462540244771864986580936194097052614323899161136957898085970489552696796829142941767873398469659323629162778846292353592314711482790641419939134587376117856242838882447857005932406327159087905651878683729549525155438714265186838888338542003497392583424805908035867677620612657408740665417187557782298921404248986397940433562598892792901895658782456159403260094469918950260597243143622042253692118943242417876593241536240565697845116242607464859307477168648376199 , 385040559604473530575023943446816858012800613980368960980295128072110959893390479012368186598828957639354389680129370950594550393689543594955832200697022500457423065980019278092108933444059135941624679752023471470951269300136074695589984776607883342494227592728757287316615506971284102230958216030906168369735111390121812979739354800647236905142065600923667936685502288606685730143880684378137983216590775589367936104520504280974171057664524187221488147186060568096871332283895497662975533299319571983317021845185422144474057304206631677983552123091532537834137643666171056952179530988729058253795875839402473090358300341175895306312850476337660814301572478630181006705416390498348147797270761818904241580521096617352248929817534367482150442926192953026930699210256851163216640411967900002837171622112708738765835136610569258115306770783550510189121088760972988391603122170801838050745380584534545563282887922370633264034314830136685893034054133133111746289684800950866705041235404460861992008971167085663856753695469706713142080497998282718365788233548770890)
です