P-Adic square roots: Difference between revisions
m
Minor formatting of code.
(New post.) |
m (Minor formatting of code.) |
||
| (3 intermediate revisions by the same user not shown) | |||
Line 36:
__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}}==
Line 685 ⟶ 1,023:
List<List<Integer>> tests = List.of( List.of( 2, 497, 10496 ),
List.of( 3, 15403, 26685 ),
List.of(
for ( List<Integer> test : tests ) {
Line 697 ⟶ 1,035:
System.out.println("The rational value is " + square.convertToRational());
System.out.println();
}
}
Line 704 ⟶ 1,042:
final class PadicSquareRoot {
/**
* Create a
* which is the p-adic square root of the given rational 'aNumerator' / 'aDenominator'.
*/
Line 714 ⟶ 1,052:
prime = aPrime;
digits = new ArrayList<Integer>(DIGITS_SIZE);
order = 0;
// Process rational zero
Line 731 ⟶ 1,069:
aDenominator /= prime;
order -= 1;
}
if ( ( order & 1 ) != 0 ) {
throw new AssertionError("Number does not have a square root in " + prime + "-adic");
}
order >>= 1;
if ( prime == 2 ) {
Line 749 ⟶ 1,083:
}
* Return the additive inverse of this
*/
public PadicSquareRoot negate() {
Line 757 ⟶ 1,091:
}
List<Integer>
negateDigits(
return new PadicSquareRoot(prime,
}
/**
* Return the product of this
*/
public PadicSquareRoot multiply(PadicSquareRoot aOther) {
Line 779 ⟶ 1,113:
/**
* Return a string representation of this
*/
public String convertToRational() {
Line 785 ⟶ 1,119:
if ( numbers.isEmpty() ) {
return "
}
Line 797 ⟶ 1,131:
}
final MathContext mathContext = new MathContext(PRECISION, RoundingMode.HALF_UP);
final BigDecimal primeBig = BigDecimal.valueOf(prime);
final BigDecimal maximumPrimeBig = BigDecimal.valueOf(maximumPrime);
Line 847 ⟶ 1,181:
/**
* Return a string representation of this
*/
public String toString() {
List<Integer> numbers = new ArrayList<Integer>(digits);
padWithZeros(numbers);
Collections.reverse(numbers);
String numberString = numbers.stream().map(String::valueOf).collect(Collectors.joining());
StringBuilder builder = new StringBuilder(numberString);
if ( order >= 0 ) {
for ( int i = 0; i < order; i++ ) {
builder.append("0");
}
Line 865 ⟶ 1,198:
} else {
builder.insert(builder.length() + order, ".");
while ( builder.toString().endsWith("0") ) {
builder.deleteCharAt(builder.length() - 1);
}
}
Line 873 ⟶ 1,210:
/**
* Create a
*
* With 'aOrder' = 0, the list [1, 2, 3, 4, 5] creates the p-adic ...54321.0
Line 890 ⟶ 1,227:
private void squareRootEvenPrime(int aNumerator, int aDenominator) {
if ( Math.floorMod(aNumerator * aDenominator, 8) != 1 ) {
throw new AssertionError(
}
Line 898 ⟶ 1,235:
// Further digits
final BigInteger numerator = BigInteger.valueOf(aNumerator);
final BigInteger denominator = BigInteger.valueOf(aDenominator);
while ( digits.size() < DIGITS_SIZE ) {
BigInteger factor = denominator.multiply(sum.multiply(sum)).subtract(numerator);
Line 929 ⟶ 1,269:
if ( firstDigit == 0 ) {
throw new IllegalArgumentException("Number does not have a square root in " + prime + "-adic");
}
Line 936 ⟶ 1,275:
// Further digits
final BigInteger numerator = BigInteger.valueOf(aNumerator);
final BigInteger denominator = BigInteger.valueOf(aDenominator);
final BigInteger firstDigitBig = BigInteger.valueOf(firstDigit);
final BigInteger primeBig = BigInteger.valueOf(prime);
Line 942 ⟶ 1,283:
BigInteger sum = firstDigitBig;
for ( int i = 2; i <
BigInteger nextSum =
sum.subtract(coefficient.multiply(denominator.multiply(sum).multiply(sum).subtract(numerator)));
Line 986 ⟶ 1,327:
/**
* The given list is padded on the right by zeros up to a maximum length of '
*/
private static void padWithZeros(List<Integer> aList) {
Line 996 ⟶ 1,337:
private List<Integer> digits;
private int order;
private final int prime;
Line 1,009 ⟶ 1,348:
{{ out }}
<pre>
The two square roots are:
...0101011010110011.1101
Line 1,023 ⟶ 1,362:
The rational value is 15403 / 26685
Number: -19 / 1 in
The two square roots are:
...
...
The p-adic value is ...
The rational value is -19 / 1
</pre>
| |||