/usr/share/axiom-20170501/src/algebra/IROOT.spad is in axiom-source 20170501-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 | )abbrev package IROOT IntegerRoots
++ Author: Michael Monagan
++ Date Created: November 1987
++ Description:
++ The \spadtype{IntegerRoots} package computes square roots and
++ nth roots of integers efficiently.
IntegerRoots(I) : SIG == CODE where
I : IntegerNumberSystem
NNI ==> NonNegativeInteger
SIG ==> with
perfectNthPower? : (I, NNI) -> Boolean
++ \spad{perfectNthPower?(n,r)} returns true if n is an \spad{r}th
++ power and false otherwise
perfectNthRoot : (I,NNI) -> Union(I,"failed")
++ \spad{perfectNthRoot(n,r)} returns the \spad{r}th root of n if n
++ is an \spad{r}th power and returns "failed" otherwise
perfectNthRoot : I -> Record(base:I, exponent:NNI)
++ \spad{perfectNthRoot(n)} returns \spad{[x,r]}, where \spad{n = x\^r}
++ and r is the largest integer such that n is a perfect \spad{r}th power
approxNthRoot : (I,NNI) -> I
++ \spad{approxRoot(n,r)} returns an approximation x
++ to \spad{n**(1/r)} such that \spad{-1 < x - n**(1/r) < 1}
perfectSquare? : I -> Boolean
++ \spad{perfectSquare?(n)} returns true if n is a perfect square
++ and false otherwise
perfectSqrt : I -> Union(I,"failed")
++ \spad{perfectSqrt(n)} returns the square root of n if n is a
++ perfect square and returns "failed" otherwise
approxSqrt : I -> I
++ \spad{approxSqrt(n)} returns an approximation x
++ to \spad{sqrt(n)} such that \spad{-1 < x - sqrt(n) < 1}.
++ Compute an approximation s to \spad{sqrt(n)} such that
++ \spad{-1 < s - sqrt(n) < 1}
++ A variable precision Newton iteration is used.
++ The running time is \spad{O( log(n)**2 )}.
CODE ==> add
import IntegerPrimesPackage(I)
resMod144: List I := [0::I,1::I,4::I,9::I,16::I,25::I,36::I,49::I,_
52::I,64::I,73::I,81::I,97::I,100::I,112::I,121::I]
two := 2::I
perfectSquare? a == (perfectSqrt a) case I
perfectNthPower?(b, n) == perfectNthRoot(b, n) case I
perfectNthRoot n == -- complexity (log log n)**2 (log n)**2
m:NNI
(n = 1) or zero? n or n = -1 => [n, 1]
e:NNI := 1
p:NNI := 2
while p::I <= length(n) + 1 repeat
for m in 0.. while (r := perfectNthRoot(n, p)) case I repeat
n := r::I
e := e * p ** m
p := convert(nextPrime(p::I))@Integer :: NNI
[n, e]
approxNthRoot(a, n) == -- complexity (log log n) (log n)**2
zero? n => error "invalid arguments"
(n = 1) => a
n=2 => approxSqrt a
negative? a =>
odd? n => - approxNthRoot(-a, n)
0
zero? a => 0
(a = 1) => 1
-- quick check for case of large n
((3*n) quo 2)::I >= (l := length a) => two
-- the initial approximation must be >= the root
y := max(two, shift(1, (n::I+l-1) quo (n::I)))
z:I := 1
n1:= (n-1)::NNI
while z > 0 repeat
x := y
xn:= x**n1
y := (n1*x*xn+a) quo (n*xn)
z := x-y
x
perfectNthRoot(b, n) ==
(r := approxNthRoot(b, n)) ** n = b => r
"failed"
perfectSqrt a ==
a < 0 or not member?(a rem (144::I), resMod144) => "failed"
(s := approxSqrt a) * s = a => s
"failed"
approxSqrt a ==
a < 1 => 0
if (n := length a) > (100::I) then
-- variable precision newton iteration
n := n quo (4::I)
s := approxSqrt shift(a, -2 * n)
s := shift(s, n)
return ((1 + s + a quo s) quo two)
-- initial approximation for the root is within a factor of 2
(new, old) := (shift(1, n quo two), 1)
while new ^= old repeat
(new, old) := ((1 + new + a quo new) quo two, new)
new
|