Square form factorization: Difference between revisions
New post.
(New post.) |
(New post.) |
||
Line 360:
f = 26675843 N/f = 57626245319
...
</pre>
=={{header|C++}}==
<syntaxhighlight lang="c++">
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <random>
uint64_t test_value = 0;
uint64_t sqrt_test_value = 0;
class BQF { // Binary quadratic form
public:
BQF(const uint64_t& a, const uint64_t& b, const uint64_t& c) : a(a), b(b), c(c) {
q = ( sqrt_test_value + b ) / c;
bb = q * c - b;
}
BQF rho() {
return BQF(c, bb, a + q * ( b - bb ));
}
BQF rho_inverse() {
return BQF(c, bb, ( test_value - bb * bb ) / c);
}
uint64_t a, b, c;
private:
uint64_t q, bb;
};
uint64_t squfof(const uint64_t& number) {
const uint32_t sqrt = std::sqrt(number);
if ( sqrt * sqrt == number ) {
return sqrt;
}
test_value = number;
sqrt_test_value = std::sqrt(test_value);
// Principal form
BQF form(0, sqrt_test_value, 1);
form = form.rho_inverse();
// Search principal cycle
for ( uint32_t i = 0; i < 4 * std::sqrt(2 * sqrt_test_value); i += 2 ) {
// Even step
form = form.rho();
uint64_t sqrt_c = std::sqrt(form.c);
if ( sqrt_c * sqrt_c == form.c ) { // Square form found
// Inverse square root
BQF form_inverse(0, -form.b, sqrt_c);
form_inverse = form_inverse.rho_inverse();
// Search ambiguous cycle
uint64_t previous_b = 0;
do {
previous_b = form_inverse.b;
form_inverse = form_inverse.rho();
} while ( form_inverse.b != previous_b );
// Symmetry point
const uint64_t g = std::gcd(number, form_inverse.a);
if ( g != 1 ) {
return g;
}
}
// Odd step
form = form.rho();
}
if ( number % 2 == 0 ) {
return 2;
}
return 0; // Failed to factorise
}
int main() {
std::random_device random;
std::mt19937 generator(random());
const uint64_t lower_limit = 100'000'000'000'000'000;
std::uniform_int_distribution<uint64_t> distribution(lower_limit, 10 * lower_limit);
for ( uint32_t i = 0; i < 20; ++i ) {
uint64_t test = distribution(random);
std::cout << "N = " << test;
uint64_t factor = squfof(test);
if ( factor == 0 ) {
std::cout << " Failed to factorise" << std::endl;
} else {
std::cout << " = " << factor << " * " << test / factor << std::endl;
}
std::cout << std::endl;
}
}
</syntaxhighlight>
{{ out }}
<pre>
N = 822140815871714649 = 141 * 5830785928168189
N = 473377979025428817 = 3 * 157792659675142939
N = 482452941918160803 = 4410431 * 109389069213
N = 165380937127655630 = 65438 * 2527292049385
N = 191677853606692475 = 7589219 * 25256598025
N = 480551815975206727 = 2843 * 169029833265989
N = 178710207362206205 = 5 * 35742041472441241
N = 484660189375949842 = 1094 * 443016626486243
N = 758704390319635770 = 1605 * 472713015775474
N = 820453356193182720 = 97280 * 8433936638499
N = 706982627912630220 = 121273 * 5829678724140
N = 614913973550671312 = 437204432 * 1406467841
N = 601482456081568543 = 131 * 4591469130393653
N = 610533314488947626 = 14 * 43609522463496259
N = 336343281182924332 = 70108 * 4797502156429
N = 308127213282933401 = 7 * 44018173326133343
N = 582455924775519843 = 3 * 194151974925173281
N = 694215100094443276 = 32070628 * 21646445467
N = 398821795604697523 = 181 * 2203435334832583
N = 477964959783291032 = 517029608 * 924444079
</pre>
| |||