P-Adic square roots
You are encouraged to solve this task according to the task description, using any language you may know.
- Task.
Convert rational a/b to its approximate p-adic square root. To check the result, square the root and construct rational m/n to compare with radicand a/b.
For rational reconstruction Lagrange's lattice basis reduction algorithm is used.
Recipe: find root x1 modulo p and build a sequence of solutions
f(xk) ≡ 0 (mod pk),
using the lifting equation
xk+1 = xk + dk * pk
with dk = –(f(xk) / pk) /
f ′(x1) (mod p).
The multipliers dk are the successive p-adic digits to find.
If evaluation of f(x) = bx2 – a overflows, the expansion is cut off and might be too short to retrieve the radicand. Setting a higher precision won't help, using a programming language with built-in large integer support will.
- Related task.
- Reference.
[1] Solving x2 ≡ a (mod n)
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.
#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;
}
}
- Output:
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
FreeBASIC
' ***********************************************
'subject: p-adic square roots, Hensel lifting.
'tested : FreeBasic 1.07.0
'The root is squared, approximated by a
'rational, and compared with radicand a/b.
const emx = 48
'exponent maximum
const amx = 700000
'tentative argument maximum
'------------------------------------------------
const Mxd = cdbl(2)^53 - 1
'max. float64 integer
const Pmax = 32749
'max. prime < 2^15
type ratio
as longint a, b
end type
type padic
declare function sqrt (byref q as ratio, byval sw as integer) as integer
'p-adic square root of q = a/b, set sw to print
declare sub printf (byval sw as integer)
'print expansion, set sw to print rational
declare function crat (byval sw as integer) as ratio
'rational reconstruction
declare sub cmpt (byref a as padic)
'let self:= complement_a
declare sub sqr (byref a as padic)
'let self:= a ^ 2
as long d(-emx to emx - 1)
as integer v
end type
'global variables
dim shared as long p1, p = 7
'default prime
dim shared as integer k = 11
'precision
#define min(a, b) iif((a) > (b), b, a)
'------------------------------------------------
'p-adic square root of g = a/b
function padic.sqrt (byref g as ratio, byval sw as integer) as integer
dim as longint a = g.a, b = g.b
dim as longint q, x, pk, pm
dim as long f1, r, s, t
dim i as integer, f as double
sqrt = 0
if b = 0 then return 1
if b < 0 then b = -b: a = -a
if p < 2 or k < 1 then return 1
'max. short prime
p = min(p, Pmax)
if sw then
'echo numerator, denominator,
print a;"/";str(b);" + ";
'prime and precision
print "O(";str(p);"^";str(k);")"
end if
'initialize
v = 0
p1 = p - 1
for i = -emx to emx - 1
d(i) = 0: next
if a = 0 then return 0
'valuation
do until b mod p
b \= p: v -= 1
loop
do until a mod p
a \= p: v += 1
loop
if (v and 1) = 1 then
'odd valuation
print "non-residue mod"; p
return -1
end if
'max. array length
k = min(k + v, emx - 1)
k -= v: v shr= 1
if abs(a) > amx or b > amx then return -1
if p = 2 then
'1 / b = b (mod 8)
'a / b = 1 (mod 8)
t = a * b
if (t and 7) - 1 then
print "non-residue mod 8"
return -1
end if
else
'find root for small p
for r = 1 to p1
q = b * r * r - a
if q mod p = 0 then exit for
next r
if r = p then
print "non-residue mod"; p
return -1
end if
'f'(r) = 2br
t = b * r shl 1
s = 0
t mod= p
'modular inverse for small p
for f1 = 1 to p1
s += t
if s > p1 then s -= p
if s = 1 then exit for
next f1
if f1 = p then
print "impossible inverse mod"
return -1
end if
end if
'evaluate f(x)
#macro evalf(x)
f = b * x * cdbl(x / pk)
f -= cdbl(a / pk)
'overflow
if f > Mxd then exit for
q = clngint(f)
#endmacro
if p = 2 then
'initialize
x = 1
d(v) = 1
d(v + 1) = 0
pk = 4
for i = v + 2 to k - 1 + v
pk shl= 1
'2-power overflow
if pk < 1 then exit for
evalf(x)
'next digit
d(i) = iif(q and 1, 1, 0)
'lift x
x += d(i) * (pk shr 1)
next i
else
'-1 / f'(x) mod p
f1 = p - f1
x = r
d(v) = x
pk = 1
for i = v + 1 to k - 1 + v
pm = pk: pk *= p
if pk \ pm - p then exit for
evalf(x)
d(i) = q * f1 mod p
if d(i) < 0 then d(i) += p
x += d(i) * pk
next i
end if
k = i - v
if sw then print "lift:";x;" mod";p;"^";str(k)
end function
'------------------------------------------------
'rational reconstruction
function padic.crat (byval sw as integer) as ratio
dim as integer i, j, t = min(v, 0)
dim as longint s, pk, pm
dim as long q, x, y
dim as double f, h
dim r as ratio
'weighted digit sum
s = 0: pk = 1
for i = t to k - 1 + v
pm = pk: pk *= p
if pk \ pm - p then
'overflow
pk = pm: exit for
end if
s += d(i) * pm '(mod pk)
next i
'lattice basis reduction
dim as longint m(1) = {pk, s}
dim as longint n(1) = {0, 1}
'norm(v)^2
h = cdbl(s) * s + 1
i = 0: j = 1
'Lagrange's algorithm
do
f = m(i) * (m(j) / h)
f += n(i) * (n(j) / h)
'Euclidean step
q = int(f +.5)
m(i) -= q * m(j)
n(i) -= q * n(j)
f = h
h = cdbl(m(i)) * m(i)
h += cdbl(n(i)) * n(i)
'compare norms
if h < f then
'interchange vectors
swap i, j
else
exit do
end if
loop
x = m(j): y = n(j)
if y < 0 then y = -y: x = -x
'check determinant
t = abs(m(i) * y - x * n(i)) = pk
if t = 0 then
print "crat: fail"
x = 0: y = 1
else
'negative powers
for i = v to -1
y *= p: next
if sw then
print x;
if y > 1 then print "/";str(y);
print
end if
end if
r.a = x: r.b = y
return r
end function
'print expansion
sub padic.printf (byval sw as integer)
dim as integer i, t = min(v, 0)
for i = k - 1 + t to t step -1
print d(i);
if i = 0 andalso v < 0 then print ".";
next i
print
'rational approximation
if sw then crat(sw)
end sub
'------------------------------------------------
'let self:= complement_a
sub padic.cmpt (byref a as padic)
dim i as integer, r as padic
dim as long c = 1
with r
.v = a.v
for i = .v to k +.v
c += p1 - a.d(i)
'carry
if c > p1 then
.d(i) = c - p: c = 1
else
.d(i) = c: c = 0
end if
next i
end with
this = r
end sub
'let self:= a ^ 2
sub padic.sqr (byref a as padic)
dim as long ptr rp, ap = @a.d(a.v)
dim as longint q, c = 0
dim as integer i, j
dim r as padic
with r
.v = a.v shl 1
rp = @.d(.v)
for i = 0 to k
for j = 0 to i
c += ap[j] * ap[i - j]
next j
'Euclidean step
q = c \ p
rp[i] = c - q * p
c = q
next i
end with
this = r
end sub
'main
'------------------------------------------------
dim as integer sw
dim as padic a, c
dim as ratio q, r
width 64, 30
cls
' -7 + O(2^7)
data -7,1, 2,7
data 9,1, 2,8
data 17,1, 2,9
data 497,10496, 2,18
data 10496,497, 2,19
data -577215,664901, 3,23
data 15403,26685, 3,18
data -1,1, 5,8
data 86,25, 5,8
data 2150,1, 5,8
data 2,1, 7,8
data 11696,621467, 7,11
data -27764,11521, 7,11
data -27584,12953, 7,11
data -166420,135131, 11,11
data 14142,135623, 5,15
data -255,256, 257,3
data 0,0, 0,0
print
do
read q.a,q.b, p,k
sw = a.sqrt(q, 1)
if sw = 1 then exit do
if sw then ? : continue do
print "sqrt +/-"
print "...";
a.printf(0)
a.cmpt(a)
print "...";
a.printf(0)
c.sqr(a)
print "sqrt^2"
print " ";
c.printf(0)
r = c.crat(1)
'{r = q}
if q.a * r.b - r.a * q.b then
print "fail: sqrt^2"
end if
print : ?
loop
end- Examples:
-7/1 + O(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 + O(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 + O(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 + O(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 + O(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
-577215/664901 + O(3^23)
lift: 44765673211 mod 3^23
sqrt +/-
... 1 0 2 1 1 1 2 2 1 0 1 1 1 2 1 1 0 1 0 2 0 1 0
... 1 2 0 1 1 1 0 0 1 2 1 1 1 0 1 1 2 1 2 0 2 2 0
sqrt^2
0 1 2 1 1 0 0 2 1 1 0 0 1 1 0 2 0 1 1 0 1 0 0
-577215/664901
15403/26685 + O(3^18)
lift: 238961782 mod 3^18
sqrt +/-
... 1 2 1 1 2 2 1 2 2 1 1 1 2 2 1 1 0. 1
... 1 0 1 1 0 0 1 0 0 1 1 1 0 0 1 1 2. 2
sqrt^2
1 2 1 2 0 1 2 2 0 1 1 0 1 2 2 2. 0 1
15403/26685
-1/1 + O(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 + O(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 + O(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 + O(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
11696/621467 + O(7^11)
lift: 967215488 mod 7^11
sqrt +/-
... 3 2 6 5 3 1 2 4 1 4. 1
... 3 4 0 1 3 5 4 2 5 2. 6
sqrt^2
3 3 1 1 3 3 4 4 5. 1 1
11696/621467
-27764/11521 + O(7^11)
lift: 453209984 mod 7^11
sqrt +/-
... 1 4 1 4 2 1 3 5 6 2 3
... 5 2 5 2 4 5 3 1 0 4 4
sqrt^2
5 1 0 2 5 5 5 3 6 6 2
-27764/11521
-27584/12953 + O(7^11)
lift: 1239987302 mod 7^11
sqrt +/-
... 4 2 5 0 4 5 0 1 2 2 1
... 2 4 1 6 2 1 6 5 4 4 6
sqrt^2
3 2 6 5 3 1 2 4 1 4 1
-27584/12953
-166420/135131 + O(11^11)
lift: 34219617965 mod 11^11
sqrt +/-
... 1 3 5 7 0 0 9 10 5 0 5
... 9 7 5 3 10 10 1 0 5 10 6
sqrt^2
1 4 1 4 2 1 3 5 6 2 3
-166420/135131
14142/135623 + O(5^15)
lift: 236397452 mod 5^15
sqrt +/-
... 0 0 0 4 4 1 0 0 4 2 0 4 3 0 2
... 4 4 4 0 0 3 4 4 0 2 4 0 1 4 3
sqrt^2
4 2 1 3 3 4 3 4 3 4 2 3 2 0 4
14142/135623
-255/256 + O(257^3)
lift: 8867596 mod 257^3
sqrt +/-
... 134 66 68
... 122 190 189
sqrt^2
255 255 255
-255/256
Haskell
Uses solution for P-Adic_numbers,_basic#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
λ> :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
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).
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;
}
- Output:
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
Julia
using Nemo, LinearAlgebra
set_printing_mode(FlintPadicField, :terse)
""" convert to Rational (rational reconstruction) """
function toRational(pa::padic)
rat = lift(QQ, pa)
r, den = BigInt(numerator(rat)), Int(denominator(rat))
p, k = Int(prime(parent(pa))), Int(precision(pa))
N = BigInt(p^k)
a1, a2 = [N, 0], [r, 1]
while dot(a1, a1) > dot(a2, a2)
q = dot(a1, a2) // dot(a2, a2)
a1, a2 = a2, a1 - round(q) * a2
end
if dot(a1, a1) < N
return (Rational{Int}(a1[1]) // Rational{Int}(a1[2])) // Int(den)
else
return Int(r) // den
end
end
function dstring(pa::padic)
u, v, n, p, k = pa.u, pa.v, pa.N, pa.parent.p, pa.parent.prec_max
d = digits(v > 0 ? u * p^v : u, base=p, pad=k)
return prod([i == k + v && v != 0 ? "$x . " : "$x " for (i, x) in enumerate(reverse(d))])
end
const DATA = [
[-7, 1, 2, 7],
[9, 1, 2, 8],
[17, 1, 2, 9],
[-1, 1, 5, 8],
[86, 25, 5, 8],
[2150, 1, 5, 8],
[2, 1, 7, 8],
[3029, 4821, 7, 9],
[379, 449, 7, 8],
[717, 8, 11, 7],
[1414, 213, 41, 5],
[-255, 256, 257, 3]
]
for (num1, den1, P, K) in DATA
Qp = PadicField(P, K)
a = Qp(QQ(num1 // den1))
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
- Output:
121 + O(2^7) sqrt +/- 1 1 1 0 1 0 1 0 0 0 1 0 1 1 Check sqrt^2: 1 1 1 1 0 0 1 -7//1 9 + O(2^8) sqrt +/- 1 1 1 1 1 1 0 1 0 0 0 0 0 0 1 1 Check sqrt^2: 0 0 0 0 1 0 0 1 9//1 17 + O(2^9) sqrt +/- 0 1 1 1 0 1 0 0 1 1 0 0 0 1 0 1 1 1 Check sqrt^2: 0 0 0 0 1 0 0 0 1 17//1 390624 + O(5^8) sqrt +/- 3 2 4 3 1 2 1 2 1 2 0 1 3 2 3 3 Check sqrt^2: 4 4 4 4 4 4 4 4 -1//1 86/25 + O(5^6) sqrt +/- 1 1 0 2 4 1 1 . 1 3 3 4 2 0 3 3 . 4 Check sqrt^2: 0 0 0 0 0 3 . 2 1 86//25 2150 + O(5^10) sqrt +/- 1 1 0 2 4 1 1 1 0 . 3 3 4 2 0 3 3 4 0 . Check sqrt^2: 0 0 0 3 2 1 0 0 2150//1 2 + O(7^8) sqrt +/- 2 1 2 1 6 2 1 3 4 5 4 5 0 4 5 4 Check sqrt^2: 0 0 0 0 0 0 0 2 2//1 39483088 + O(7^9) sqrt +/- 5 2 5 2 4 5 3 1 1 1 4 1 4 2 1 3 5 6 Check sqrt^2: 6 5 6 4 1 3 0 2 1 3029//4821 4493721 + O(7^8) sqrt +/- 0 4 5 0 1 2 2 1 6 2 1 6 5 4 4 6 Check sqrt^2: 5 3 1 2 4 1 4 1 379//449 2435986 + O(11^7) sqrt +/- 2 10 8 7 0 1 5 8 0 2 3 10 9 6 Check sqrt^2: 1 4 1 4 2 1 3 717//8 36443037 + O(41^5) sqrt +/- 9 1 38 6 8 31 39 2 34 33 Check sqrt^2: 12 36 31 15 23 1414//213 16908285 + O(257^3) sqrt +/- 134 66 68 122 190 189 Check sqrt^2: 255 255 255 -255//256
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()
- Output:
-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
Phix
constant EMX = 48 // exponent maximum (if indexing starts at -EMX)
constant DMX = 1e5 // approximation loop maximum
constant AMX = 700000 // argument maximum
constant PMAX = 32749 // prime maximum
// global variables
integer p1 = 0
integer p = 7 // default prime
integer k = 11 // precision
type Ratio(sequence r)
return length(r)=2 and integer(r[1]) and integer(r[2])
end type
class Padic
integer v = 0
sequence d = repeat(0,EMX*2)
function square_root(Ratio g, integer sw)
-- p-adic square root of g = a/b
integer {a,b} = g
atom f, q, pk, x
integer f1, r, s, t, i, res = 0
if b = 0 then return 1 end if
if b < 0 then
b = -b
a = -a
end if
if p < 2 or k < 1 then return 1 end if
-- max. short prime
p = min(p, PMAX)
if sw then
-- numerator, denominator, prime, precision
printf(1,"%d/%d + O(%d^%d)\n",{a,b,p,k})
end if
-- initialize
v = 0
p1 = p - 1
sequence ntd = repeat(0,2*EMX) -- (new this.d)
if a = 0 then return 0 end if
-- valuation
while remainder(b,p)=0 do
b /= p
v -= 1
end while
while remainder(a,p)=0 do
a /= p
v += 1
end while
if remainder(v,2) then
-- odd valuation
printf(1,"(1)non-residue mod %d\n",p)
return -1
end if
-- max. array length
k = min(k + v, EMX - 1) - v
v = floor(v/2)
if abs(a) > AMX or b > AMX then return -1 end if
if p = 2 then
--1 / b = b (mod 8)
--a / b = 1 (mod 8)
t = a * b
if mod(t,8)-1 then
printf(1,"(2)non-residue mod 8\n")
return -1
end if
else
-- find root for small p
for r = 1 to p1 do
f = b * r * r - a
if mod(f,p) = 0 then exit end if
end for
if r = p then
printf(1,"(3)non-residue mod %d\n", p)
return -1
end if
-- f'(r) = 2br
t = b * r * 2
s = 0
t = mod(t,p)
-- modular inverse for small p
for f1 = 1 to p1 do
s += t
if s > p1 then s -= p end if
if s = 1 then exit end if
end for
if f1 = p then
printf(1,"impossible inverse mod\n")
return -1
end if
end if
if p = 2 then
-- initialize
x = 1
ntd[v+EMX+1] = 1
ntd[v+EMX+2] = 0
pk = 4
for i = v+2 to k-1+v do
pk *= 2
f = b * x * x - a
q = floor(f/pk)
-- overflow
if f != q * pk then exit end if
-- next digit
ntd[i+EMX+1] = and_bits(q,1)
-- lift x
x += ntd[i+EMX+1] * floor(pk/2)
end for
else
-- -1 / f'(x) mod p
f1 = p - f1
x = r
ntd[v+EMX+1] = x
pk = 1
for i = v+1 to k-1 do
pk *= p
f = b * x * x - a
q = floor(f/pk)
-- overflow
if f - q * pk then exit end if
r = mod(q*f1,p)
if r < 0 then r += p end if
ntd[i+EMX+1] = r
x += r * pk
end for
end if
this.d = ntd
k = i-v
if sw then
printf(1,"lift: %d mod %d^%d\n",{x,p,k})
end if
return 0
end function
function square()
integer c = 0
Padic r = new()
r.v = this.v * 2
sequence td = this.d,
rd = r.d
for i=0 to k do
for j=0 to i do
c += td[v+j+EMX+1] * td[v+i-j+EMX+1]
end for
// Euclidean step
integer q = floor(c/p)
rd[r.v+i+EMX+1] = c - q*p
c = q
end for
r.d = rd
return r
end function
function complement()
integer c = 1
Padic r = new({v})
sequence rd = r.d
for i=v to k+v do
integer dx = i+EMX+1
c += p1 - this.d[dx]
if c>p1 then
rd[dx] = c - p
c = 1
else
rd[dx] = c
c = 0
end if
end for
r.d = rd
return r
end function
function crat(integer sw)
-- rational reconstruction
integer i, j, t = min(v, 0)
Ratio r
atom f
integer x, y
atom p1,pk, q, s
-- weighted digit sum
s = 0
pk = 1
for i = t to k-1+v do
p1 = pk
pk *= p
if floor(pk/p1) - p then
-- overflow
pk = p1
exit
end if
s += d[i+EMX+1] * p1 --(mod pk)
end for
-- lattice basis reduction
sequence m = {pk, s},
n = {0, 1}
i = 1
j = 2
s = s * s + 1 -- norm(v)^2
-- Lagrange's algorithm
while true do
f = (m[i] * m[j] + n[i] * n[j]) / s
-- Euclidean step
q = floor(f +.5)
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 then
-- interchange vectors
{i,j} = {j,i}
else
exit
end if
end while
x = m[j]
y = n[j]
if y < 0 then
y = -y
x = -x
end if
-- check determinant
t = abs(m[i] * y - x * n[i]) == pk
if not t then
printf(1,"crat: fail\n")
x = 0
y = 1
else
-- negative powers
for i = v to -1 do
y *= p
end for
if sw then
-- printf(1,iff(y=1?"%d":"%d/%d"),{x*sgn,y})
printf(1,iff(y=1?"%d\n":"%d/%d\n"),{x,y})
end if
end if
r = {x,y}
return r
end function
procedure prntf(bool sw)
-- print expansion
integer t = min(v, 0)
for i=k-1+t to t by -1 do
printf(1,"%d",d[i+EMX+1])
printf(1,iff(i=0 and v<0?". ":" "))
end for
printf(1,"\n")
// rational approximation
if sw then crat(sw) end if
end procedure
end class
constant tests = {
{{-7,1},2,7},
--/*
{{9,1},2,8},
{{17,1},2,9},
{{497,10496},2,18},
{{10496,497},2,19},
{{3141,5926},3,17},
{{2718,281},3,15},
{{-1,1},5,8},
{{86,25},5,8},
{{2150,1},5,10},
{{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}
}
Padic a = new(), c
Ratio q, r
for i=1 to length(tests) do
{q,p,k} = tests[i]
integer sw = a.square_root(q, 1)
if sw=1 then exit end if
if sw=0 then
printf(1,"square_root +/-\n")
printf(1,"... ")
a.prntf(0)
a = a.complement()
printf(1,"... ")
a.prntf(0)
c = a.square()
printf(1,"square_root^2\n")
printf(1," ")
c.prntf(0)
r = c.crat(1)
if q[1] * r[2] - r[1] * q[2] then
printf(1,"fail: square_root^2\n")
end if
end if
printf(1,"\n")
end for
- Output:
-7/1 + O(2^7)
lift: 53 mod 2^7
square_root +/-
... 0 1 1 0 1 0 1
... 1 0 0 1 0 1 1
square_root^2
1 1 1 1 0 0 1
-7
-255/256 + O(257^3)
lift: 8867596 mod 257^3
square_root +/-
... 134 66 68
... 122 190 189
square_root^2
255 255 255
-255/256
Wren
import "./dynamic" for Struct
import "./big" for BigInt
// constants
var EMX = 64 // exponent maximum (if indexing starts at -EMX)
var AMX = 6000 // argument maximum
var PMAX = 32749 // prime maximum
// global variables
var P1 = 0
var P = 7 // default prime
var K = 11 // precision
var Ratio = Struct.create("Ratio", ["a", "b"])
class Padic {
// uninitialized
construct new() {
_v = 0
_d = List.filled(2 * EMX, 0) // add EMX to index to be consistent wih FB
}
// properties
v { _v }
v=(o) { _v = o }
d { _d }
// (re)initialize 'this' to the square root of a Ratio, set 'sw' to print
sqrt(g, sw) {
var a = g.a
var b = g.b
if (b == 0) return 1
if (b < 0) {
b = -b
a = -a
}
if (P < 2 || K < 1) return 1
P = P.min(PMAX) // maximum short prime
if (sw != 0) {
System.write("%(a)/%(b) + ") // numerator, denominator
System.print("0(%(P)^%(K))") // prime, precision
}
// (re)initialize
_v = 0
P1 = P - 1
_d = List.filled(2 * EMX, 0)
if (a == 0) return 0
//valuation
while (b%P== 0) {
b = (b/P).truncate
_v = _v - 1
}
while (a%P == 0) {
a = (a/P).truncate
_v = _v + 1
}
if ((_v & 1) == 1) {
// odd valuation
System.print("non-residue mod %(P)")
return -1
}
K = (K + _v).min(EMX - 1) - _v // maximum array length
_v = (_v/2).truncate
if (a.abs > AMX || b > AMX) return -1
var bb = BigInt.new(b) // to avoid overflowing 'f(x) = b * x * x – a'
var r
var s
var t
var f
var f1
if (P == 2) {
t = a * b
if ((t & 7) - 1 != 0) {
System.print("non-residue mod 8")
return -1
}
} else {
// find root for small P
r = 1
while (r <= P1) {
f = bb * r * r - a
if ((f % P) == 0) break
r = r + 1
}
if (r == P) {
System.print("non-residue mod %(P)")
return -1
}
t = 2 * b * r
s = 0
t = t % P
// modular inverse for small P
f1 = 1
while (f1 <= P1) {
s = s + t
if (s > P1) s = s - P
if (s == 1) break
f1 = f1 + 1
}
if (f1 == P) {
System.print("impossible inverse mod")
return -1
}
}
var x
var pk
var q
var i
if (P == 2) {
// initialize
x = 1
_d[_v+EMX] = 1
_d[_v+1+EMX] = 0
pk = 4
i = _v + 2
while (i <= K - 1 + _v) {
pk = pk * 2
f = bb * x * x - a
q = f / pk
// overflow
if (f != q * pk) break
// next digit
_d[i+EMX] = ((q & 1) != 0) ? 1 : 0
// lift x
x = x + _d[i+EMX]*(pk >> 1)
i = i + 1
}
} else {
f1 = P - f1
x = r
_d[_v+EMX] = x
pk = 1
i = _v + 1
while (i <= K - 1 + _v) {
pk = pk * P
f = bb * x * x - a
q = f / pk
// overflow
if (f != q * pk) break
_d[i+EMX] = q.toSmall * f1 % P
if (_d[i+EMX] < 0) _d[i+EMX] = _d[i+EMX] + P
x = x + _d[i+EMX]*pk
i = i + 1
}
}
K = i - _v
if (sw != 0) System.print("lift: %(x) mod %(P)^%(K)")
return 0
}
// rational reconstruction
crat(sw) {
var t = _v.min(0)
// weighted digit sum
var s = 0
var pk = 1
for (i in t..K-1+_v) {
P1 = pk
pk = pk * P
if (((pk/P1).truncate - P) != 0) {
// overflow
pk = p1
break
}
s = s + _d[i+EMX]*P1
}
// lattice basis reduction
var m = [pk, s]
var n = [0, 1]
var i = 0
var j = 1
s = s * s + 1
// Lagrange's algorithm
while (true) {
var f = (m[i] * m[j] + n[i] * n[j]) / s
// Euclidean step
var q = (f + 0.5).floor
m[i] = m[i] - q*m[j]
n[i] = n[i] - q*n[j]
q = s
s = m[i] * m[i] + n[i] * n[i]
// compare norms
if (s < q) {
// interchange vectors
var z = i
i = j
j = z
} else {
break
}
}
var x = m[j]
var y = n[j]
if (y < 0) {
y = -y
x = -x
}
// check determinant
t = (m[i]*y - x*n[i]).abs == pk
if (!t) {
System.print("crat: fail")
x = 0
y = 1
} else {
// negative powers
var i = _v
while (i <= -1) {
y = y * P
i = i + 1
}
if (sw != 0) {
System.write(x)
if (y > 1) System.write("/%(y)")
System.print()
}
}
return Ratio.new(x, y)
}
// print expansion
printf(sw) {
var t = _v.min(0)
for (i in K - 1 + t..t) {
System.write(_d[i + EMX])
if (i == 0 && _v < 0) System.write(".")
System.write(" ")
}
System.print()
// rational approximation
if (sw != 0) crat(sw)
}
// complement
cmpt {
var c = 1
var r = Padic.new()
r.v = _v
for (i in r.v..K + r.v) {
c = c + P1 - _d[i+EMX]
if (c > P1) {
r.d[i+EMX] = c - P
c = 1
} else {
r.d[i+EMX] = c
c = 0
}
}
return r
}
// square
sqr {
var c = 0
var r = Padic.new()
r.v = _v * 2
for (i in 0..K) {
for (j in 0..i) c = c + _d[_v+j+EMX] * _d[_v+i-j+EMX]
// Euclidean step
var q = (c/P).truncate
r.d[r.v+i+EMX] = c - q*P
c = q
}
return r
}
}
var 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]
]
var sw = 0
var a = Padic.new()
var c = Padic.new()
for (d in data) {
var q = Ratio.new(d[0], d[1])
P = d[2]
K = d[3]
sw = a.sqrt(q, 1)
if (sw == 1) break
if (sw == 0) {
System.print("sqrt +/-")
System.write("...")
a.printf(0)
a = a.cmpt
System.write("...")
a.printf(0)
c = a.sqr
System.print("sqrt^2")
System.write(" ")
c.printf(0)
var r = c.crat(1)
if (q.a * r.b - r.a * q.b != 0) {
System.print("fail: sqrt^2")
}
System.print()
}
}
- Output:
-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