P-Adic square roots: Difference between revisions
m
Minor formatting of code.
m (Minor formatting of code.) |
|||
| (12 intermediate revisions by 7 users not shown) | |||
Line 1:
{{
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}}==
<
' ***********************************************
'subject: p-adic square roots, Hensel lifting.
Line 432 ⟶ 770:
end
</syntaxhighlight>
{{out|Examples}}
Line 606 ⟶ 944:
-255/256
</pre>
=={{header|Haskell}}==
Uses solution for [[P-Adic_numbers,_basic#Haskell]]
<syntaxhighlight lang="haskell">{-# LANGUAGE KindSignatures, DataKinds #-}
import Data.Ratio
import Data.List (find)
import GHC.TypeLits
import Padic
pSqrt :: KnownNat p => Rational -> Padic p
pSqrt r = res
where
res = maybe Null mkUnit series
(a, b) = (numerator r, denominator r)
series = case modulo res of
2 | eqMod 4 a 3 -> Nothing
| not (eqMod 8 a 1) -> Nothing
| otherwise -> Just $ 1 : 0 : go 8 1
where
go pk x =
let q = ((b*x*x - a) `div` pk) `mod` 2
in q : go (2*pk) (x + q * (pk `div` 2))
p -> do
y <- find (\x -> eqMod p (b*x*x) a) [1..p-1]
df <- recipMod p (2*b*y)
let go pk x =
let f = (b*x*x - a) `div` pk
d = (f * (p - df)) `mod` p
in x `div` (pk `div` p) : go (p*pk) (x + d*pk)
Just $ go p y
eqMod :: Integral a => a -> a -> a -> Bool
eqMod p a b = a `mod` p == b `mod` p</syntaxhighlight>
<pre>λ> :set -XDataKinds
λ> pSqrt 2 :: Padic 7
7-adic: 12011266421216213.0
λ> pSqrt 39483088 :: Padic 7
7-adic: 22666526525245311.0
λ> pSqrt 5 :: Padic 7
Null
λ> pSqrt 9 :: Padic 2
2-adic: 11111111111111101.0
λ> it ^ 2
2-adic: 1001.0
λ> toRational it
9 % 1
λ> pSqrt 17 ^ 2 == (17 :: Padic 2)
True
λ> pSqrt 50 ^ 2 == (50 :: Padic 7)
True
λ> pSqrt 52 ^ 2 == (52 :: Padic 7)
False
λ> pSqrt 52 :: Padic 7
Null</pre>
=={{header|Java}}==
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. p-adic numbers and rational reconstructions are given correct to O(prime^20).
<syntaxhighlight lang="java">
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.stream.Collectors;
public final class PAdicSquareRoots {
public static void main(String[] args) {
List<List<Integer>> tests = List.of( List.of( 2, 497, 10496 ),
List.of( 3, 15403, 26685 ),
List.of( 7, -19, 1 ) );
for ( List<Integer> test : tests ) {
System.out.println("Number: " + test.get(1) + " / " + test.get(2) + " in " + test.get(0) + "-adic");
PadicSquareRoot squareRoot = new PadicSquareRoot(test.get(0), test.get(1), test.get(2));
System.out.println("The two square roots are:");
System.out.println(" " + squareRoot);
System.out.println(" " + squareRoot.negate());
PadicSquareRoot square = squareRoot.multiply(squareRoot);
System.out.println("The p-adic value is " + square);
System.out.println("The rational value is " + square.convertToRational());
System.out.println();
}
}
}
final class PadicSquareRoot {
/**
* Create a PadicSquareRoot number, with p = 'aPrime',
* which is the p-adic square root of the given rational 'aNumerator' / 'aDenominator'.
*/
public PadicSquareRoot(int aPrime, int aNumerator, int aDenominator) {
if ( aDenominator == 0 ) {
throw new IllegalArgumentException("Denominator cannot be zero");
}
prime = aPrime;
digits = new ArrayList<Integer>(DIGITS_SIZE);
order = 0;
// Process rational zero
if ( aNumerator == 0 ) {
order = MAX_ORDER;
return;
}
// Remove multiples of 'prime' and adjust the order of the p-adic number accordingly
while ( Math.floorMod(aNumerator, prime) == 0 ) {
aNumerator /= prime;
order += 1;
}
while ( Math.floorMod(aDenominator, prime) == 0 ) {
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 ) {
squareRootEvenPrime(aNumerator, aDenominator);
} else {
squareRootOddPrime(aNumerator, aDenominator);
}
}
/**
* Return the additive inverse of this PadicSquareRoot number.
*/
public PadicSquareRoot negate() {
if ( digits.isEmpty() ) {
return this;
}
List<Integer> negated = new ArrayList<Integer>(digits);
negateDigits(negated);
return new PadicSquareRoot(prime, negated, order);
}
/**
* Return the product of this PadicSquareRoot number and the given PadicSquareRoot number.
*/
public PadicSquareRoot multiply(PadicSquareRoot aOther) {
if ( prime != aOther.prime ) {
throw new IllegalArgumentException("Cannot multiply p-adic's with different primes");
}
if ( digits.isEmpty() || aOther.digits.isEmpty() ) {
return new PadicSquareRoot(prime, 0 , 1);
}
return new PadicSquareRoot(prime, multiply(digits, aOther.digits), order + aOther.order);
}
/**
* Return a string representation of this PadicSquareRoot as a rational number.
*/
public String convertToRational() {
List<Integer> numbers = new ArrayList<Integer>(digits);
if ( numbers.isEmpty() ) {
return "0 / 1";
}
// Lagrange lattice basis reduction in two dimensions
long seriesSum = numbers.getFirst();
long maximumPrime = 1;
for ( int i = 1; i < PRECISION; i++ ) {
maximumPrime *= prime;
seriesSum += numbers.get(i) * maximumPrime;
}
final MathContext mathContext = new MathContext(PRECISION, RoundingMode.HALF_UP);
final BigDecimal primeBig = BigDecimal.valueOf(prime);
final BigDecimal maximumPrimeBig = BigDecimal.valueOf(maximumPrime);
final BigDecimal seriesSumBig = BigDecimal.valueOf(seriesSum);
BigDecimal[] one = new BigDecimal[] { maximumPrimeBig, seriesSumBig };
BigDecimal[] two = new BigDecimal[] { BigDecimal.ZERO, BigDecimal.ONE };
BigDecimal previousNorm = BigDecimal.valueOf(seriesSum).pow(2).add(BigDecimal.ONE);
BigDecimal currentNorm = previousNorm.add(BigDecimal.ONE);
int i = 0;
int j = 1;
while ( previousNorm.compareTo(currentNorm) < 0 ) {
currentNorm = one[i].multiply(one[j]).add(two[i].multiply(two[j])).divide(previousNorm, mathContext);
currentNorm = currentNorm.setScale(0, RoundingMode.HALF_UP);
one[i] = one[i].subtract(currentNorm.multiply(one[j]));
two[i] = two[i].subtract(currentNorm.multiply(two[j]));
currentNorm = previousNorm;
previousNorm = one[i].multiply(one[i]).add(two[i].multiply(two[i]));
if ( previousNorm.compareTo(currentNorm) < 0 ) {
final int temp = i; i = j; j = temp;
}
}
BigDecimal x = one[j];
BigDecimal y = two[j];
if ( y.signum() == -1 ) {
y = y.negate();
x = x.negate();
}
if ( ! one[i].multiply(y).subtract(x.multiply(two[i])).abs().equals(maximumPrimeBig) ) {
throw new AssertionError("Rational reconstruction failed.");
}
for ( int k = order; k < 0; k++ ) {
y = y.multiply(primeBig);
}
for ( int k = order; k > 0; k-- ) {
x = x.multiply(primeBig);
}
return x + " / " + y;
}
/**
* Return a string representation of this PadicSquareRoot number.
*/
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");
}
builder.append(".0");
} else {
builder.insert(builder.length() + order, ".");
while ( builder.toString().endsWith("0") ) {
builder.deleteCharAt(builder.length() - 1);
}
}
return " ..." + builder.toString().substring(builder.length() - PRECISION - 1);
}
// PRIVATE //
/**
* Create a PadicSquareRoot, with p = 'aPrime', directly from a list of digits.
*
* With 'aOrder' = 0, the list [1, 2, 3, 4, 5] creates the p-adic ...54321.0
* 'aOrder' > 0 shifts the list 'aOrder' places to the left and
* 'aOrder' < 0 shifts the list 'aOrder' places to the right.
*/
private PadicSquareRoot(int aPrime, List<Integer> aDigits, int aOrder) {
prime = aPrime;
digits = new ArrayList<Integer>(aDigits);
order = aOrder;
}
/**
* Create a 2-adic number which is the square root of the rational 'aNumerator' / 'aDenominator'.
*/
private void squareRootEvenPrime(int aNumerator, int aDenominator) {
if ( Math.floorMod(aNumerator * aDenominator, 8) != 1 ) {
throw new AssertionError("Number does not have a square root in 2-adic");
}
// First digit
BigInteger sum = BigInteger.ONE;
digits.addLast(sum.intValue());
// 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);
int valuation = 0;
while ( factor.mod(BigInteger.TWO).signum() == 0 ) {
factor = factor.shiftRight(1);
valuation += 1;
}
sum = sum.add(BigInteger.TWO.pow(valuation - 1));
for ( int i = digits.size(); i < valuation - 1; i++ ) {
digits.addLast(0);
}
digits.addLast(1);
}
}
/**
* Create a p-adic number, with an odd prime number, p = 'prime',
* which is the p-adic square root of the given rational 'aNumerator' / 'aDenominator'.
*/
private void squareRootOddPrime(int aNumerator, int aDenominator) {
// First digit
int firstDigit = 0;
for ( int i = 1; i < prime && firstDigit == 0; i++ ) {
if ( ( aDenominator * i * i - aNumerator ) % prime == 0 ) {
firstDigit = i;
}
}
if ( firstDigit == 0 ) {
throw new IllegalArgumentException("Number does not have a square root in " + prime + "-adic");
}
digits.addLast(firstDigit);
// 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);
final BigInteger coefficient =
denominator.multiply(firstDigitBig).shiftLeft(1).mod(primeBig).modInverse(primeBig);
BigInteger sum = firstDigitBig;
for ( int i = 2; i < DIGITS_SIZE; i++ ) {
BigInteger nextSum =
sum.subtract(coefficient.multiply(denominator.multiply(sum).multiply(sum).subtract(numerator)));
nextSum = nextSum.mod(primeBig.pow(i));
nextSum = nextSum.subtract(sum);
sum = sum.add(nextSum);
final int digit = nextSum.divideAndRemainder(primeBig.pow(i - 1))[0].intValueExact();
digits.addLast(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].
*/
private List<Integer> multiply(List<Integer> aOne, List<Integer> aTwo) {
List<Integer> product = new ArrayList<Integer>(Collections.nCopies(aOne.size() + aTwo.size(), 0));
for ( int b = 0; b < aTwo.size(); b++ ) {
int carry = 0;
for ( int a = 0; a < aOne.size(); a++ ) {
product.set(a + b, product.get(a + b) + aOne.get(a) * aTwo.get(b) + carry);
carry = product.get(a + b) / prime;
product.set(a + b, product.get(a + b) % prime);
}
product.set(b + aOne.size(), carry);
}
return product.subList(0, DIGITS_SIZE);
}
/**
* Transform the given list of digits representing a p-adic number
* into a list which represents the negation of the p-adic number.
*/
private void negateDigits(List<Integer> aDigits) {
aDigits.set(0, Math.floorMod(prime - aDigits.get(0), prime));
for ( int i = 1; i < aDigits.size(); i++ ) {
aDigits.set(i, prime - 1 - aDigits.get(i));
}
}
/**
* The given list is padded on the right by zeros up to a maximum length of 'DIGITS_SIZE'.
*/
private static void padWithZeros(List<Integer> aList) {
while ( aList.size() < DIGITS_SIZE ) {
aList.addLast(0);
}
}
private List<Integer> digits;
private int order;
private final int prime;
private static final int MAX_ORDER = 1_000;
private static final int PRECISION = 20;
private static final int DIGITS_SIZE = PRECISION + 5;
}
</syntaxhighlight>
{{ out }}
<pre>
umber: 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: 15403 / 26685 in 3-adic
The two square roots are:
...2112112212211122110.1
...0110110010011100112.2
The p-adic value is ...111212012201101222.01
The rational value is 15403 / 26685
Number: -19 / 1 in 7-adic
The two square roots are:
...0126260226045320643.0
...6540406440621346024.0
The p-adic value is ...6666666666666666642.0
The rational value is -19 / 1
</pre>
=={{header|Julia}}==
<
set_printing_mode(FlintPadicField, :terse)
Line 622 ⟶ 1,384:
while dot(a1, a1) > dot(a2, a2)
q = dot(a1, a2) // dot(a2, a2)
a1, a2 = a2, a1 -
end
if dot(a1, a1) < N
Line 654 ⟶ 1,416:
for (num1, den1, P, K) in DATA
Qp = PadicField(P, K)
a = Qp(QQ
c = sqrt(a)
r = toRational(c * c)
println(a, "\nsqrt +/-\n", dstring(c), "\n", dstring(-c), "\nCheck sqrt^2:\n", dstring(c * c), "\n", r, "\n")
end
</
<pre>
121 + O(2^7)
Line 757 ⟶ 1,519:
-255//256
</pre>
=={{header|Nim}}==
{{trans|FreeBASIC}}
{{trans|Wren}}
<syntaxhighlight lang="nim">import strformat
const
Emx = 64 # Exponent maximum.
Amx = 6000 # Argument maximum.
PMax = 32749 # Prime maximum.
type
Ratio = tuple[a, b: int]
Padic = object
p: int # Prime.
k: int # Precision.
v: int
d: array[-Emx..(Emx-1), int]
PadicError = object of ValueError
proc sqrt(pa: var Padic; q: Ratio; sw: bool) =
## Return the p-adic square root of q = a/b. Set sw to print.
var (a, b) = q
var i, x: int
if b == 0:
raise newException(PadicError, &"Wrong rational: {a}/{b}" )
if b < 0:
b = -b
a = -a
if pa.p < 2:
raise newException(PadicError, &"Wrong value for p: {pa.p}")
if pa.k < 1:
raise newException(PadicError, &"Wrong value for k: {pa.k}")
pa.p = min(pa.p, PMax) # Maximum short prime.
if sw: echo &"{a}/{b} + 0({pa.p}^{pa.k})"
# Initialize.
pa.v = 0
pa.d.reset()
if a == 0: return
# Valuation.
while b mod pa.p == 0:
b = b div pa.p
dec pa.v
while a mod pa.p == 0:
a = a div pa.p
inc pa.v
if (pa.v and 1) != 0:
# Odd valuation.
raise newException(PadicError, &"Non-residue mod {pa.p}.")
# Maximum array length.
pa.k = min(pa.k + pa.v, Emx - 1) - pa.v
pa.v = pa.v shr 1
if abs(a) > Amx or b > Amx:
raise newException(PadicError, &"Rational exceeding limits: {a}/{b}.")
if pa.p == 2:
# 1 / b = b (mod 8); a / b = 1 (mod 8).
if (a * b and 7) - 1 != 0:
raise newException(PadicError, "Non-residue mod 8.")
# Initialize.
x = 1
pa.d[pa.v] = 1
pa.d[pa.v + 1] = 0
var pk = 4
i = pa.v + 2
while i < pa.k + pa.v:
pk *= 2
let f = b * x * x - a
let q = f div pk
if f != q * pk: break # Overflow.
# Next digit.
pa.d[i] = if (q and 1) != 0: 1 else: 0
# Lift "x".
x += pa.d[i] * (pk shr 1)
inc i
else:
# Find root for small "p".
var r = 1
while r < pa.p:
if (b * r * r - a) mod pa.p == 0: break
inc r
if r == pa.p:
raise newException(PadicError, &"Non-residue mod {pa.p}.")
let t = (b * r shl 1) mod pa.p
var s = 0
# Modular inverse for small "p".
var f1 = 1
while f1 < pa.p:
inc s, t
if s >= pa.p: dec s, pa.p
if s == 1: break
inc f1
if f1 == pa.p:
raise newException(PadicError, "Impossible to compute inverse modulo")
f1 = pa.p - f1
x = r
pa.d[pa.v] = x
var pk = 1
i = pa.v + 1
while i < pa.k + pa.v:
pk *= pa.p
let f = b * x * x - a
let q = f div pk
if f != q * pk: break # Overflow.
pa.d[i] = q * f1 mod pa.p
if pa.d[i] < 0: pa.d[i] += pa.p
x += pa.d[i] * pk
inc i
pa.k = i - pa.v
if sw: echo &"lift: {x} mod {pa.p}^{pa.k}"
proc crat(pa: Padic; sw: bool): Ratio =
## Rational reconstruction.
# Weighted digit sum.
var
s = 0
pk = 1
for i in min(pa.v, 0)..<(pa.k + pa.v):
let pm = pk
pk *= pa.p
if pk div pm - pa.p != 0:
# Overflow.
pk = pm
break
s += pa.d[i] * pm
# Lattice basis reduction.
var
m = [pk, s]
n = [0, 1]
i = 0
j = 1
s = s * s + 1
# Lagrange's algorithm.
while true:
# Euclidean step.
var q = ((m[i] * m[j] + n[i] * n[j]) / s).toInt
m[i] -= q * m[j]
n[i] -= q * n[j]
q = s
s = m[i] * m[i] + n[i] * n[i]
# Compare norms.
if s < q: swap i, j # Interchange vectors.
else: break
var x = m[j]
var y = n[j]
if y < 0:
y = -y
x = -x
# Check determinant.
if abs(m[i] * y - x * n[i]) != pk:
raise newException(PadicError, "Rational reconstruction failed.")
# Negative powers.
for i in pa.v..(-1): y *= pa.p
if sw: echo x, if y > 1: '/' & $y else: ""
result = (x, y)
func cmpt(pa: Padic): Padic =
## Return the complement.
result = Padic(p: pa.p, k: pa.k, v: pa.v)
var c = 1
for i in pa.v..(pa.k + pa.v):
inc c, pa.p - 1 - pa.d[i]
if c >= pa.p:
result.d[i] = c - pa.p
c = 1
else:
result.d[i] = c
c = 0
func sqr(pa: Padic): Padic =
## Return the square of a P-adic number.
result = Padic(p: pa.p, k: pa.k, v: pa.v * 2)
var c = 0
for i in 0..pa.k:
for j in 0..i:
c += pa.d[pa.v + j] * pa.d[pa.v + i - j]
# Euclidean step.
let q = c div pa.p
result.d[result.v + i] = c - q * pa.p
c = q
func `$`(pa: Padic): string =
## String representation.
let t = min(pa.v, 0)
for i in countdown(pa.k - 1 + t, t):
result.add $pa.d[i]
if i == 0 and pa.v < 0: result.add "."
result.add " "
when isMainModule:
const Data = [[-7, 1, 2, 7],
[9, 1, 2, 8],
[17, 1, 2, 9],
[497, 10496, 2, 18],
[10496, 497, 2, 19],
[3141, 5926, 3, 15],
[2718, 281, 3, 13],
[-1, 1, 5, 8],
[86, 25, 5, 8],
[2150, 1, 5, 8],
[2,1, 7, 8],
[-2645, 28518, 7, 9],
[3029, 4821, 7, 9],
[379, 449, 7, 8],
[717, 8, 11, 7],
[1414, 213, 41, 5],
[-255, 256, 257, 3]]
for d in Data:
try:
let q: Ratio = (d[0], d[1])
var a = Padic(p: d[2], k: d[3])
a.sqrt(q, true)
echo "sqrt +/-"
echo "...", a
a = a.cmpt()
echo "...", a
let c = sqr(a)
echo "sqrt^2"
echo " ", c
let r = c.crat(true)
if q.a * r.b - r.a * q.b != 0:
echo "fail: sqrt^2"
echo ""
except PadicError:
echo getCurrentExceptionMsg()</syntaxhighlight>
{{out}}
<pre>-7/1 + 0(2^7)
lift: 53 mod 2^7
sqrt +/-
...0 1 1 0 1 0 1
...1 0 0 1 0 1 1
sqrt^2
1 1 1 1 0 0 1
-7
9/1 + 0(2^8)
lift: 253 mod 2^8
sqrt +/-
...1 1 1 1 1 1 0 1
...0 0 0 0 0 0 1 1
sqrt^2
0 0 0 0 1 0 0 1
9
17/1 + 0(2^9)
lift: 233 mod 2^9
sqrt +/-
...0 1 1 1 0 1 0 0 1
...1 0 0 0 1 0 1 1 1
sqrt^2
0 0 0 0 1 0 0 0 1
17
497/10496 + 0(2^18)
lift: 92989 mod 2^18
sqrt +/-
...0 1 0 1 1 0 1 0 1 1 0 0 1 1. 1 1 0 1
...1 0 1 0 0 1 0 1 0 0 1 1 0 0. 0 0 1 1
sqrt^2
1 0 0 0 0 0 1 1 0 0. 1 0 0 0 1 0 0 1
497/10496
10496/497 + 0(2^19)
lift: 148501 mod 2^19
sqrt +/-
...1 0 0 0 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0
...0 1 1 1 0 1 1 1 1 1 0 1 0 1 1 0 0 0 0
sqrt^2
0 0 1 1 0 1 1 1 0 0 1 0 0 0 0 0 0 0 0
10496/497
3141/5926 + 0(3^15)
lift: 3406966 mod 3^15
sqrt +/-
...2 0 1 0 2 0 0 2 1 1 0 2 2 1 0
...0 2 1 2 0 2 2 0 1 1 2 0 0 2 0
sqrt^2
2 1 1 1 2 2 0 0 0 2 0 1 1 0 0
3141/5926
2718/281 + 0(3^13)
lift: 693355 mod 3^13
sqrt +/-
...0 2 2 0 2 0 0 0 2 2 1 1 0
...2 0 0 2 0 2 2 2 0 0 1 2 0
sqrt^2
2 2 0 0 1 2 2 0 2 2 1 0 0
2718/281
-1/1 + 0(5^8)
lift: 280182 mod 5^8
sqrt +/-
...3 2 4 3 1 2 1 2
...1 2 0 1 3 2 3 3
sqrt^2
4 4 4 4 4 4 4 4
-1
86/25 + 0(5^8)
lift: 95531 mod 5^8
sqrt +/-
...1 1 0 2 4 1 1. 1
...3 3 4 2 0 3 3. 4
sqrt^2
0 0 0 0 0 3. 2 1
86/25
2150/1 + 0(5^8)
lift: 95531 mod 5^8
sqrt +/-
...1 0 2 4 1 1 1 0
...3 4 2 0 3 3 4 0
sqrt^2
0 0 0 3 2 1 0 0
2150
2/1 + 0(7^8)
lift: 1802916 mod 7^8
sqrt +/-
...2 1 2 1 6 2 1 3
...4 5 4 5 0 4 5 4
sqrt^2
0 0 0 0 0 0 0 2
2
-2645/28518 + 0(7^9)
lift: 39082527 mod 7^9
sqrt +/-
...6 5 3 1 2 4 1 4. 1
...0 1 3 5 4 2 5 2. 6
sqrt^2
1 1 3 3 4 4 5. 1 1
-2645/28518
3029/4821 + 0(7^9)
lift: 31104424 mod 7^9
sqrt +/-
...5 2 5 2 4 5 3 1 1
...1 4 1 4 2 1 3 5 6
sqrt^2
6 5 6 4 1 3 0 2 1
3029/4821
379/449 + 0(7^8)
lift: 555087 mod 7^8
sqrt +/-
...0 4 5 0 1 2 2 1
...6 2 1 6 5 4 4 6
sqrt^2
5 3 1 2 4 1 4 1
379/449
717/8 + 0(11^7)
lift: 5280093 mod 11^7
sqrt +/-
...2 10 8 7 0 1 5
...8 0 2 3 10 9 6
sqrt^2
1 4 1 4 2 1 3
717/8
1414/213 + 0(41^5)
lift: 25564902 mod 41^5
sqrt +/-
...9 1 38 6 8
...31 39 2 34 33
sqrt^2
12 36 31 15 23
1414/213
-255/256 + 0(257^3)
lift: 8867596 mod 257^3
sqrt +/-
...134 66 68
...122 190 189
sqrt^2
255 255 255
-255/256</pre>
=={{header|Phix}}==
{{trans|FreeBASIC}}
<
constant DMX = 1e5 // approximation loop maximum
constant AMX = 700000 // argument maximum
Line 1,095 ⟶ 2,266:
end if
printf(1,"\n")
end for</
{{out}}
<pre>
Line 1,120 ⟶ 2,291:
{{trans|FreeBASIC}}
{{libheader|Wren-dynamic}}
{{libheader|Wren-big}}
<
import "./
// constants
Line 1,160 ⟶ 2,329:
}
if (P < 2 || K < 1) return 1
P =
if (sw != 0) {
System.write("%(a)/%(b) + ") // numerator, denominator
Line 1,188 ⟶ 2,357:
return -1
}
K =
_v = (_v/2).truncate
Line 1,282 ⟶ 2,451:
// rational reconstruction
crat(sw) {
var t =
// weighted digit sum
var s = 0
Line 1,354 ⟶ 2,523:
// print expansion
printf(sw) {
var t =
for (i in K - 1 + t..t) {
System.write(_d[i + EMX])
Line 1,446 ⟶ 2,615:
System.print()
}
}</
{{out}}
| |||