P-Adic square roots: Difference between revisions
Content added Content deleted
m (Improved code.) |
(New post.) |
||
Line 36: | Line 36: | ||
__TOC__ |
__TOC__ |
||
=={{header|C++}}== |
|||
This example displays p-adic numbers in standard mathematical format, consisting of a possibly infinite list of digits extending leftwards from the p-adic point. |
|||
<syntaxhighlight lang="c++"> |
|||
#include <algorithm> |
|||
#include <cmath> |
|||
#include <cstdint> |
|||
#include <iostream> |
|||
#include <stdexcept> |
|||
#include <vector> |
|||
class P_adic_square_root { |
|||
public: |
|||
// Create a P_adic_square_root number, with p = 'prime', from the given rational 'numerator' / 'denominator'. |
|||
P_adic_square_root(const uint32_t& prime, const uint32_t& precision, int32_t numerator, int32_t denominator) |
|||
: prime(prime), precision(precision) { |
|||
if ( denominator == 0 ) { |
|||
throw std::invalid_argument("Denominator cannot be zero"); |
|||
} |
|||
order = 0; |
|||
// Process rational zero |
|||
if ( numerator == 0 ) { |
|||
digits.assign(digits_size, 0); |
|||
order = ORDER_MAX; |
|||
return; |
|||
} |
|||
// Remove multiples of 'prime' and adjust the order of the P_adic_square_root number accordingly |
|||
while ( modulo(numerator, prime) == 0 ) { |
|||
numerator /= static_cast<int32_t>(prime); |
|||
order += 1; |
|||
} |
|||
while ( modulo(denominator, prime) == 0 ) { |
|||
denominator /= static_cast<int32_t>(prime); |
|||
order -= 1; |
|||
} |
|||
if ( ( order & 1 ) != 0 ) { |
|||
throw std::invalid_argument("Number does not have a square root in " + std::to_string(prime) + "-adic"); |
|||
} |
|||
order >>= 1; |
|||
if ( prime == 2 ) { |
|||
square_root_even_prime(numerator, denominator); |
|||
} else { |
|||
square_root_odd_prime(numerator, denominator); |
|||
} |
|||
} |
|||
// Return the additive inverse of this P_adic_square_root number. |
|||
P_adic_square_root negate() { |
|||
if ( digits.empty() ) { |
|||
return *this; |
|||
} |
|||
std::vector<uint32_t> negated = digits; |
|||
negate_digits(negated); |
|||
return P_adic_square_root(prime, precision, negated, order); |
|||
} |
|||
// Return the product of this P_adic_square_root number and the given P_adic_square_root number. |
|||
P_adic_square_root multiply(P_adic_square_root other) { |
|||
if ( prime != other.prime ) { |
|||
throw std::invalid_argument("Cannot multiply p-adic's with different primes"); |
|||
} |
|||
if ( digits.empty() || other.digits.empty() ) { |
|||
return P_adic_square_root(prime, precision, 0 , 1); |
|||
} |
|||
return P_adic_square_root(prime, precision, multiply(digits, other.digits), order + other.order); |
|||
} |
|||
// Return a string representation of this P_adic_square_root as a rational number. |
|||
std::string convertToRational() { |
|||
std::vector<uint32_t> numbers = digits; |
|||
if ( numbers.empty() ) { |
|||
return "0 / 1"; |
|||
} |
|||
// Lagrange lattice basis reduction in two dimensions |
|||
int64_t series_sum = numbers.front(); |
|||
int64_t maximum_prime = 1; |
|||
for ( uint32_t i = 1; i < precision; ++i ) { |
|||
maximum_prime *= prime; |
|||
series_sum += numbers[i] * maximum_prime; |
|||
} |
|||
std::vector<int64_t> one = { maximum_prime, series_sum }; |
|||
std::vector<int64_t> two = { 0, 1 }; |
|||
int64_t previous_norm = series_sum * series_sum + 1; |
|||
int64_t current_norm = previous_norm + 1; |
|||
uint32_t i = 0; |
|||
uint32_t j = 1; |
|||
while ( previous_norm < current_norm ) { |
|||
int64_t numerator = one[i] * one[j] + two[i] * two[j]; |
|||
int64_t denominator = previous_norm; |
|||
current_norm = std::floor(static_cast<double>(numerator) / denominator + 0.5); |
|||
one[i] -= current_norm * one[j]; |
|||
two[i] -= current_norm * two[j]; |
|||
current_norm = previous_norm; |
|||
previous_norm = one[i] * one[i] + two[i] * two[i]; |
|||
if ( previous_norm < current_norm ) { |
|||
std::swap(i, j); |
|||
} |
|||
} |
|||
int64_t x = one[j]; |
|||
int64_t y = two[j]; |
|||
if ( y < 0 ) { |
|||
y = -y; |
|||
x = -x; |
|||
} |
|||
if ( std::abs(one[i] * y - x * two[i]) != maximum_prime ) { |
|||
throw std::invalid_argument("Rational reconstruction failed."); |
|||
} |
|||
for ( int32_t k = order; k < 0; ++k ) { |
|||
y *= prime; |
|||
} |
|||
for ( int32_t k = order; k > 0; --k ) { |
|||
x *= prime; |
|||
} |
|||
return std::to_string(x) + " / " + std::to_string(y); |
|||
} |
|||
// Return a string representation of this P_adic_square_root. |
|||
std::string to_string() { |
|||
std::vector<uint32_t> numbers = digits; |
|||
pad_with_zeros(numbers); |
|||
std::string result = ""; |
|||
for ( int64_t i = numbers.size() - 1; i >= 0; --i ) { |
|||
result += std::to_string(digits[i]); |
|||
} |
|||
if ( order >= 0 ) { |
|||
for ( int32_t i = 0; i < order; ++i ) { |
|||
result += "0"; |
|||
} |
|||
result += ".0"; |
|||
} else { |
|||
result.insert(result.length() + order, "."); |
|||
while ( result[result.length() - 1] == '0' ) { |
|||
result = result.substr(0, result.length() - 1); |
|||
} |
|||
} |
|||
return " ..." + result.substr(result.length() - precision - 1); |
|||
} |
|||
private: |
|||
/** |
|||
* Create a P_adic_square_root, with p = 'prime', directly from a vector of digits. |
|||
* |
|||
* For example: with 'order' = 0, the vector [1, 2, 3, 4, 5] creates the P_adic_square_root ...54321.0, |
|||
* 'order' > 0 shifts the vector 'order' places to the left and |
|||
* 'order' < 0 shifts the vector 'order' places to the right. |
|||
*/ |
|||
P_adic_square_root( |
|||
const uint32_t& prime, const uint32_t& precision, const std::vector<uint32_t>& digits, const int32_t& order) |
|||
: prime(prime), precision(precision), digits(digits), order(order) { |
|||
} |
|||
// Create a 2-adic number which is the square root of the rational 'numerator' / 'denominator'. |
|||
void square_root_even_prime(const int32_t& numerator, const int32_t& denominator) { |
|||
if ( modulo(numerator * denominator, 8) != 1 ) { |
|||
throw std::invalid_argument("Number does not have a square root in 2-adic"); |
|||
} |
|||
// First digit |
|||
uint64_t sum = 1; |
|||
digits.emplace_back(sum); |
|||
// Further digits |
|||
while ( digits.size() < digits_size ) { |
|||
int64_t factor = denominator * sum * sum - numerator; |
|||
uint32_t valuation = 0; |
|||
while ( modulo(factor, 2) == 0 ) { |
|||
factor /= 2; |
|||
valuation += 1; |
|||
} |
|||
sum += std::pow(2, valuation - 1); |
|||
for ( uint32_t i = digits.size(); i < valuation - 1; ++i ) { |
|||
digits.emplace_back(0); |
|||
} |
|||
digits.emplace_back(1); |
|||
} |
|||
} |
|||
// Create a p-adic number, with an odd prime number, p = 'prime', |
|||
// which is the p-adic square root of the given rational 'numerator' / 'denominator'. |
|||
void square_root_odd_prime(const int32_t& numerator, const int32_t& denominator) { |
|||
// First digit |
|||
int32_t first_digit = 0; |
|||
for ( int32_t i = 1; i < prime && first_digit == 0; ++i ) { |
|||
if ( modulo(denominator * i * i - numerator, prime) == 0 ) { |
|||
first_digit = i; |
|||
} |
|||
} |
|||
if ( first_digit == 0 ) { |
|||
throw std::invalid_argument("Number does not have a square root in " + std::to_string(prime) + "-adic"); |
|||
} |
|||
digits.emplace_back(first_digit); |
|||
// Further digits |
|||
const uint64_t coefficient = modulo_inverse(modulo(2 * denominator * first_digit, prime)); |
|||
uint64_t sum = first_digit; |
|||
for ( uint32_t i = 2; i < digits_size; ++i ) { |
|||
int64_t next_sum = sum - ( coefficient * ( denominator * sum * sum - numerator ) ); |
|||
next_sum = modulo(next_sum, static_cast<uint64_t>(std::pow(prime, i))); |
|||
next_sum -= sum; |
|||
sum += next_sum; |
|||
const uint32_t digit = next_sum / std::pow(prime, i - 1); |
|||
digits.emplace_back(digit); |
|||
} |
|||
} |
|||
// Return the list obtained by multiplying the digits of the given two lists, |
|||
// where the digits in each list are regarded as forming a single number in reverse. |
|||
// For example 12 * 13 = 156 is computed as [2, 1] * [3, 1] = [6, 5, 1]. |
|||
std::vector<uint32_t> multiply(const std::vector<uint32_t>& one, const std::vector<uint32_t>& two) { |
|||
std::vector<uint32_t> product(one.size() + two.size(), 0); |
|||
for ( uint32_t b = 0; b < two.size(); ++b ) { |
|||
uint32_t carry = 0; |
|||
for ( uint32_t a = 0; a < one.size(); ++a ) { |
|||
product[a + b] += one[a] * two[b] + carry; |
|||
carry = product[a + b] / prime; |
|||
product[a + b] %= prime; |
|||
} |
|||
product[b + one.size()] = carry; |
|||
} |
|||
return std::vector(product.begin(), product.begin() + digits_size); |
|||
} |
|||
// Return the multiplicative inverse of the given number modulo 'prime'. |
|||
uint32_t modulo_inverse(const uint32_t& number) const { |
|||
uint32_t inverse = 1; |
|||
while ( modulo(inverse * number, prime) != 1 ) { |
|||
inverse += 1; |
|||
} |
|||
return inverse; |
|||
} |
|||
// Return the given number modulo 'prime' in the range 0..'prime' - 1. |
|||
int32_t modulo(const int64_t& number, const int64_t& modulus) const { |
|||
const int32_t div = static_cast<int32_t>(number % modulus); |
|||
return ( div >= 0 ) ? div : div + modulus; |
|||
} |
|||
// Transform the given vector of digits representing a p-adic number |
|||
// into a vector which represents the negation of the p-adic number. |
|||
void negate_digits(std::vector<uint32_t>& numbers) { |
|||
numbers[0] = modulo(prime - numbers[0], prime); |
|||
for ( uint64_t i = 1; i < numbers.size(); ++i ) { |
|||
numbers[i] = prime - 1 - numbers[i]; |
|||
} |
|||
} |
|||
// The given vector is padded on the right by zeros up to a maximum length of 'DIGITS_SIZE'. |
|||
void pad_with_zeros(std::vector<uint32_t>& vector) { |
|||
while ( vector.size() < digits_size ) { |
|||
vector.emplace_back(0); |
|||
} |
|||
} |
|||
const int32_t prime; |
|||
const uint32_t precision; |
|||
const uint32_t digits_size = precision + 5; |
|||
std::vector<uint32_t> digits; |
|||
int32_t order; |
|||
static const uint32_t ORDER_MAX = 1'000; |
|||
}; |
|||
int main() { |
|||
std::vector<std::vector<int32_t>> tests = { { 2, 20, 497, 10496 }, { 5, 14, 86, 25 }, { 7, 10, -19, 1 } }; |
|||
for ( const std::vector<int32_t>& test : tests ) { |
|||
std::cout << "Number: " << test[2] << " / " << test[3] << " in " << test[0] << "-adic" << std::endl; |
|||
P_adic_square_root square_root(test[0], test[1], test[2], test[3]); |
|||
std::cout << "The two square roots are:" << std::endl; |
|||
std::cout << " " << square_root.to_string() << std::endl; |
|||
std::cout << " " << square_root.negate().to_string() << std::endl; |
|||
P_adic_square_root square = square_root.multiply(square_root); |
|||
std::cout << "The p-adic value is " << square.to_string() << std::endl; |
|||
std::cout << "The rational value is " << square.convertToRational() << std::endl; |
|||
std::cout << std::endl; |
|||
} |
|||
} |
|||
</syntaxhighlight> |
|||
{{ out }} |
|||
<pre> |
|||
Number: 497 / 10496 in 2-adic |
|||
The two square roots are: |
|||
...0101011010110011.1101 |
|||
...1010100101001100.0011 |
|||
The p-adic value is ...111000001100.10001001 |
|||
The rational value is 497 / 10496 |
|||
Number: 86 / 25 in 5-adic |
|||
The two square roots are: |
|||
...0104011102411.1 |
|||
...4340433342033.4 |
|||
The p-adic value is ...000000000003.21 |
|||
The rational value is 86 / 25 |
|||
Number: -19 / 1 in 7-adic |
|||
The two square roots are: |
|||
...045320643.0 |
|||
...621346024.0 |
|||
The p-adic value is ...666666642.0 |
|||
The rational value is -19 / 1 |
|||
</pre> |
|||
=={{header|FreeBASIC}}== |
=={{header|FreeBASIC}}== |