/usr/share/axiom-20170501/src/algebra/INTFACT.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 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 | )abbrev package INTFACT IntegerFactorizationPackage
++ Description:
++ This Package contains basic methods for integer factorization.
++ The factor operation employs trial division up to 10,000. It
++ then tests to see if n is a perfect power before using Pollards
++ rho method. Because Pollards method may fail, the result
++ of factor may contain composite factors. We should also employ
++ Lenstra's eliptic curve method.
IntegerFactorizationPackage(I) : SIG == CODE where
I: IntegerNumberSystem
B ==> Boolean
FF ==> Factored I
NNI ==> NonNegativeInteger
LMI ==> ListMultiDictionary I
FFE ==> Record(flg:Union("nil","sqfr","irred","prime"),
fctr:I, xpnt:Integer)
SIG ==> with
factor : I -> FF
++ factor(n) returns the full factorization of integer n
squareFree : I -> FF
++ squareFree(n) returns the square free factorization of integer n
BasicMethod : I -> FF
++ BasicMethod(n) returns the factorization
++ of integer n by trial division
PollardSmallFactor : I -> Union(I,"failed")
++ PollardSmallFactor(n) returns a factor
++ of n or "failed" if no one is found
CODE ==> add
import IntegerRoots(I)
BasicSieve: (I, I) -> FF
squareFree(n:I):FF ==
u:I
if n<0 then (m := -n; u := -1)
else (m := n; u := 1)
(m > 1) and ((v := perfectSqrt m) case I) =>
for rec in (l := factorList(sv := squareFree(v::I))) repeat
rec.xpnt := 2 * rec.xpnt
makeFR(u * unit sv, l)
-- avoid using basic sieve when the lim is too big
-- we know the sieve constants up to sqrt(100000000)
lim := 1 + approxSqrt(m)
lim > (100000000::I) => makeFR(u, factorList factor m)
x := BasicSieve(m, lim)
y :=
((m:= unit x) = 1) => factorList x
(v := perfectSqrt m) case I =>
concat_!(factorList x, ["sqfr",v,2]$FFE)
concat_!(factorList x, ["sqfr",m,1]$FFE)
makeFR(u, y)
PollardSmallFactor(n:I):Union(I,"failed") ==
-- Use the Brent variation
x0 := random()$I
m := 100::I
y := x0 rem n
r:I := 1
q:I := 1
G:I := 1
until G > 1 repeat
x := y
for i in 1..convert(r)@Integer repeat
y := (y*y+5::I) rem n
k:I := 0
until (k>=r) or (G>1) repeat
ys := y
for i in 1..convert(min(m,r-k))@Integer repeat
y := (y*y+5::I) rem n
q := q*abs(x-y) rem n
G := gcd(q,n)
k := k+m
r := 2*r
if G=n then
until G>1 repeat
ys := (ys*ys+5::I) rem n
G := gcd(abs(x-ys),n)
G=n => "failed"
G
BasicSieve(n, lim) ==
p:=primes(1::I,lim::I)$IntegerPrimesPackage(I)
l:List(I) := append([first p],reverse rest p)
ls := empty()$List(FFE)
for d in l repeat
if n<d*d then
if n>1 then ls := concat_!(ls, ["prime",n,1]$FFE)
return makeFR(1, ls)
for m in 0.. while zero?(n rem d) repeat n := n quo d
if m>0 then ls := concat_!(ls, ["prime",d,convert m]$FFE)
makeFR(n,ls)
BasicMethod n ==
u:I
if n<0 then (m := -n; u := -1)
else (m := n; u := 1)
x := BasicSieve(m, 1 + approxSqrt m)
makeFR(u, factorList x)
factor m ==
u:I
zero? m => 0
if negative? m then (n := -m; u := -1)
else (n := m; u := 1)
b := BasicSieve(n, 10000::I)
flb := factorList b
((n := unit b) = 1) => makeFR(u, flb)
a:LMI := dictionary() -- numbers yet to be factored
b:LMI := dictionary() -- prime factors found
f:LMI := dictionary() -- number which could not be factored
insert_!(n, a)
while not empty? a repeat
n := inspect a; c := count(n, a); remove_!(n, a)
prime?(n)$IntegerPrimesPackage(I) => insert_!(n, b, c)
-- test for a perfect power
(s := perfectNthRoot n).exponent > 1 =>
insert_!(s.base, a, c * s.exponent)
-- test for a difference of square
x:=approxSqrt n
if (x**2<n) then x:=x+1
(y:=perfectSqrt (x**2-n)) case I =>
insert_!(x+y,a,c)
insert_!(x-y,a,c)
(d := PollardSmallFactor n) case I =>
for m in 0.. while zero?(n rem d) repeat n := n quo d
insert_!(d, a, m * c)
if n > 1 then insert_!(n, a, c)
-- an elliptic curve factorization attempt should be made here
insert_!(n, f, c)
-- insert prime factors found
while not empty? b repeat
n := inspect b; c := count(n, b); remove_!(n, b)
flb := concat_!(flb, ["prime",n,convert c]$FFE)
-- insert non-prime factors found
while not empty? f repeat
n := inspect f; c := count(n, f); remove_!(n, f)
flb := concat_!(flb, ["nil",n,convert c]$FFE)
makeFR(u, flb)
|