--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

)abbrev package INTHEORY IntegerNumberTheoryFunctions
++ Author: Michael Monagan, Martin Brock, Robert Sutor, Timothy Daly
++ Date Created: June 1987
++ Date Last Updated:
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: number theory, integer
++ Examples:
++ References: Knuth, The Art of Computer Programming Vol.2
++ Description:
++ This package provides various number theoretic functions on the integers.
IntegerNumberTheoryFunctions(): Exports == Implementation where
 I ==> Integer
 RN ==> Fraction I
 SUP ==> SparseUnivariatePolynomial
 NNI ==> NonNegativeInteger

 Exports ==> with
  bernoulli : I -> RN
    ++ \spad{bernoulli(n)} returns the nth Bernoulli number.
    ++ this is \spad{B(n,0)}, where \spad{B(n,x)} is the \spad{n}th Bernoulli
    ++ polynomial.
  chineseRemainder: (I,I,I,I) -> I
    ++ \spad{chineseRemainder(x1,m1,x2,m2)} returns w, where w is such that
    ++ \spad{w = x1 mod m1} and \spad{w = x2 mod m2}. Note: \spad{m1} and
    ++ \spad{m2} must be relatively prime.
  divisors : I -> List I
    ++ \spad{divisors(n)} returns a list of the divisors of n.
  euler : I -> I
    ++ \spad{euler(n)} returns the \spad{n}th Euler number. This is
    ++ \spad{2^n E(n,1/2)}, where \spad{E(n,x)} is the nth Euler polynomial.
  eulerPhi : I -> I
    ++ \spad{eulerPhi(n)} returns the number of integers between 1 and n
    ++ (including 1) which are relatively prime to n. This is the Euler phi
    ++ function \spad{\phi(n)} is also called the totient function.
  fibonacci : I -> I
    ++ \spad{fibonacci(n)} returns the nth Fibonacci number. the Fibonacci
    ++ numbers \spad{F[n]} are defined by \spad{F[0] = F[1] = 1} and
    ++ \spad{F[n] = F[n-1] + F[n-2]}.
    ++ The algorithm has running time \spad{O(log(n)^3)}.
    ++ Reference: Knuth, The Art of Computer Programming
    ++ Vol 2, Semi-Numerical Algorithms.
  harmonic : I -> RN
    ++ \spad{harmonic(n)} returns the nth harmonic number. This is
    ++ \spad{H[n] = sum(1/k,k=1..n)}.
  jacobi : (I,I) -> I
    ++ \spad{jacobi(a,b)} returns the Jacobi symbol \spad{J(a/b)}.
    ++ When b is odd, \spad{J(a/b) = product(L(a/p) for p in factor b )}.
    ++ Note: by convention, 0 is returned if \spad{gcd(a,b) ~= 1}.
    ++ Iterative \spad{O(log(b)^2)} version coded by Michael Monagan June 1987.
  legendre : (I,I) -> I
    ++ \spad{legendre(a,p)} returns the Legendre symbol \spad{L(a/p)}.
    ++ \spad{L(a/p) = (-1)**((p-1)/2) mod p} (p prime), which is 0 if \spad{a}
    ++ is 0, 1 if \spad{a} is a quadratic residue \spad{mod p} and -1 otherwise.
    ++ Note: because the primality test is expensive,
    ++ if it is known that p is prime then use \spad{jacobi(a,p)}.
  moebiusMu : I -> I
    ++ \spad{moebiusMu(n)} returns the Moebius function \spad{mu(n)}.
    ++ \spad{mu(n)} is either -1,0 or 1 as follows:
    ++ \spad{mu(n) = 0} if n is divisible by a square > 1,
    ++ \spad{mu(n) = (-1)^k} if n is square-free and has k distinct
    ++ prime divisors.
  numberOfDivisors: I -> I
    ++ \spad{numberOfDivisors(n)} returns the number of integers between 1 and n
    ++ (inclusive) which divide n. The number of divisors of n is often
    ++ denoted by \spad{tau(n)}.
  sumOfDivisors : I -> I
    ++ \spad{sumOfDivisors(n)} returns the sum of the integers between 1 and n
    ++ (inclusive) which divide n. The sum of the divisors of n is often
    ++ denoted by \spad{sigma(n)}.
  sumOfKthPowerDivisors: (I,NNI) -> I
    ++ \spad{sumOfKthPowerDivisors(n,k)} returns the sum of the \spad{k}th
    ++ powers of the integers between 1 and n (inclusive) which divide n.
    ++ the sum of the \spad{k}th powers of the divisors of n is often denoted
    ++ by \spad{sigma_k(n)}.
 Implementation ==> add
  import IntegerPrimesPackage(I)

  -- we store the euler and bernoulli numbers computed so far in
  -- a Vector because they are computed from an n-term recurrence
  E: IndexedFlexibleArray(I,0)   := new(1, 1)
  B: IndexedFlexibleArray(RN,0)  := new(1, 1)
  H: Record(Hn:I,Hv:RN) := [1, 1]

  harmonic n ==
    s:I; h:RN
    negative? n => error("harmonic not defined for negative integers")
    if n >= H.Hn then (s,h) := H else (s := 0; h := 0)
    for k in s+1..n repeat h := h + 1/k
    H.Hn := n
    H.Hv := h
    h

  fibonacci n ==
    n = 0 => 0
    negative? n => (odd? n => 1; -1) * fibonacci(-n)
    f1: I := 0
    f2: I := 1
    for k in length(n)-2 .. 0 by -1 repeat
      t := f2**2
      (f1,f2) := (t+f1**2,t+2*f1*f2)
      if bit?(n,k) then (f1,f2) := (f2,f1+f2)
    f2

  euler n ==
    negative? n => error "euler not defined for negative integers"
    odd? n => 0
    l := (#E) :: I
    n < l => E(n)
    concat!(E, new((n+1-l)::NNI, 0)$IndexedFlexibleArray(I,0))
    for i in 1 .. l by 2 repeat E(i) := 0
    -- compute E(i) i = l+2,l+4,...,n given E(j) j = 0,2,...,i-2
    t,e : I
    for i in l+1 .. n by 2 repeat
      t := e := 1
      for j in 2 .. i-2 by 2 repeat
        t := (t*(i-j+1)*(i-j+2)) quo (j*(j-1))
        e := e + t*E(j)
      E(i) := -e
    E(n)

  bernoulli n ==
    negative? n => error "bernoulli not defined for negative integers"
    odd? n =>
      n = 1 => -1/2
      0
    l := (#B) :: I
    n < l => B(n)
    concat!(B, new((n+1-l)::NNI, 0)$IndexedFlexibleArray(RN,0))
    for i in 1 .. l by 2 repeat B(i) := 0
    -- compute B(i) i = l+2,l+4,...,n given B(j) j = 0,2,...,i-2
    for i in l+1 .. n by 2 repeat
      t:I := 1
      b := (1-i)/2
      for j in 2 .. i-2 by 2 repeat
        t := (t*(i-j+2)*(i-j+3)) quo (j*(j-1))
        b := b + (t::RN) * B(j)
      B(i) := -b/((i+1)::RN)
    B(n)

  inverse : (I,I) -> I
  inverse(a,b) ==
    borg:I:=b
    c1:I := 1
    d1:I := 0
    while b ~= 0 repeat
      q:I := a quo b
      r:I := a-q*b
      (a,b):=(b,r)
      (c1,d1):=(d1,c1-q*d1)
    not one? a => error("moduli are not relatively prime")
    positiveRemainder(c1,borg)
  
  chineseRemainder(x1,m1,x2,m2) ==
    negative? m1 or negative? m2 => error "moduli must be positive"
    x1 := positiveRemainder(x1,m1)
    x2 := positiveRemainder(x2,m2)
    x1 + m1 * positiveRemainder(((x2-x1) * inverse(m1,m2)),m2)


  jacobi(a,b) ==
    -- Revised by Clifton Williamson January 1989.
    -- Previous version returned incorrect answers when b was even.
    -- The formula J(a/b) = product ( L(a/p) for p in factor b) is only
    -- valid when b is odd (the Legendre symbol L(a/p) is not defined
    -- for p = 2).  When b is even, the Jacobi symbol J(a/b) is only
    -- defined for a = 0 or 1 (mod 4).  When a = 1 (mod 8),
    -- J(a/2) = +1 and when a = 5 (mod 8), we define J(a/2) = -1.
    -- Extending by multiplicativity, we have J(a/b) for even b and
    -- appropriate a.
    -- We also define J(a/1) = 1.
    -- The point of this is the following: if d is the discriminant of
    -- a quadratic field K and chi is the quadratic character for K,
    -- then J(d/n) = chi(n) for n > 0.
    -- Reference: Hecke, Vorlesungen ueber die Theorie der Algebraischen
    -- Zahlen.
    if negative? b then b := -b
    b = 0 => error "second argument of jacobi may not be 0"
    b = 1 => 1
    even? b and positiveRemainder(a,4) > 1 =>
      error "J(a/b) not defined for b even and a = 2 or 3 (mod 4)"
    even? b and even? a => 0
    k : Integer := 0
    while even? b repeat
      b := b quo 2
      k := k + 1
    j:I := (odd? k and positiveRemainder(a,8) = 5 => -1; 1)
    b = 1 => j
    a := positiveRemainder(a,b)
    -- assertion: 0 < a < b and odd? b
    while a > 1 repeat
      if odd? a then
        -- J(a/b) = J(b/a) (-1) ** (a-1)/2 (b-1)/2
        if a rem 4 = 3 and b rem 4 = 3 then j := -j
        (a,b) := (b rem a,a)
      else
        -- J(2*a/b) = J(a/b) (-1) (b**2-1)/8
        k := 0
        until odd? a repeat
          a := a quo 2
          k := k + 1
        if odd? k and (b+2) rem 8 > 4 then j := -j
    a = 0 => 0
    j

  legendre(a,p) ==
    prime? p => jacobi(a,p)
    error "characteristic of legendre must be prime"

  eulerPhi n ==
    n = 0 => 0
    r : RN := 1
    for entry in factors factor n repeat
      r := ((entry.factor - 1) /$RN entry.factor) * r
    numer(n * r)

  divisors n ==
    oldList : List Integer := [1]
    for f in factors factor n repeat
      newList := oldList
      for k in 1..f.exponent repeat
        pow := f.factor ** k
        for m in oldList repeat
          newList := concat(pow * m,newList)
      oldList := newList
    sort(#1 < #2,oldList)

  numberOfDivisors n ==
    n = 0 => 0
    */[1+entry.exponent for entry in factors factor n]

  sumOfDivisors n ==
    n = 0 => 0
    r : RN := */[(entry.factor**(entry.exponent::NNI + 1)-1)/
      (entry.factor-1) for entry in factors factor n]
    numer r

  sumOfKthPowerDivisors(n,k) ==
    n = 0 => 0
    r : RN := */[(entry.factor**(k*entry.exponent::NNI+k)-1)/
      (entry.factor**k-1) for entry in factors factor n]
    numer r

  moebiusMu n ==
    n = 1 => 1
    t := factor n
    for k in factors t repeat
      k.exponent > 1 => return 0
    odd? numberOfFactors t => -1
    1

)abbrev package PNTHEORY PolynomialNumberTheoryFunctions
++ Author: Michael Monagan, Clifton J. Williamson
++ Date Created: June 1987
++ Date Last Updated: 10 November 1996 (Claude Quitte)
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: polynomial, number theory
++ Examples:
++ References: Knuth, The Art of Computer Programming Vol.2
++ Description:
++ This package provides various polynomial number theoretic functions
++ over the integers.
PolynomialNumberTheoryFunctions(): Exports == Implementation where
 I ==> Integer
 RN ==> Fraction I
 SUP ==> SparseUnivariatePolynomial
 NNI ==> NonNegativeInteger

 Exports ==> with
  bernoulli : I -> SUP RN
    ++ bernoulli(n) returns the nth Bernoulli polynomial \spad{B[n](x)}.
    ++ Note: Bernoulli polynomials denoted \spad{B(n,x)} computed by solving the
    ++ differential equation  \spad{differentiate(B(n,x),x) = n B(n-1,x)} where
    ++ \spad{B(0,x) = 1} and initial condition comes from \spad{B(n) = B(n,0)}.
  chebyshevT: I -> SUP I
    ++ chebyshevT(n) returns the nth Chebyshev polynomial \spad{T[n](x)}.
    ++ Note: Chebyshev polynomials of the first kind, denoted \spad{T[n](x)},
    ++ computed from the two term recurrence.  The generating function
    ++ \spad{(1-t*x)/(1-2*t*x+t**2) = sum(T[n](x)*t**n, n=0..infinity)}.
  chebyshevU: I -> SUP I
    ++ chebyshevU(n) returns the nth Chebyshev polynomial \spad{U[n](x)}.
    ++ Note: Chebyshev polynomials of the second kind, denoted \spad{U[n](x)},
    ++ computed from the two term recurrence.  The generating function
    ++ \spad{1/(1-2*t*x+t**2) = sum(T[n](x)*t**n, n=0..infinity)}.
  cyclotomic: I -> SUP I
    ++ cyclotomic(n) returns the nth cyclotomic polynomial \spad{phi[n](x)}.
    ++ Note: \spad{phi[n](x)} is the factor of \spad{x**n - 1} whose roots
    ++ are the primitive nth roots of unity.
  euler     : I -> SUP RN
    ++ euler(n) returns the nth Euler polynomial \spad{E[n](x)}.
    ++ Note: Euler polynomials denoted \spad{E(n,x)} computed by solving the
    ++ differential equation  \spad{differentiate(E(n,x),x) = n E(n-1,x)} where
    ++ \spad{E(0,x) = 1} and initial condition comes from \spad{E(n) = 2**n E(n,1/2)}.
  fixedDivisor: SUP I -> I
    ++ fixedDivisor(a) for \spad{a(x)} in \spad{Z[x]} is the largest integer
    ++ f such that f divides \spad{a(x=k)} for all integers k.
    ++ Note: fixed divisor of \spad{a} is
    ++ \spad{reduce(gcd,[a(x=k) for k in 0..degree(a)])}.
  hermite   : I -> SUP I
    ++ hermite(n) returns the nth Hermite polynomial \spad{H[n](x)}.
    ++ Note: Hermite polynomials, denoted \spad{H[n](x)}, are computed from
    ++ the two term recurrence.  The generating function is:
    ++ \spad{exp(2*t*x-t**2) = sum(H[n](x)*t**n/n!, n=0..infinity)}.
  laguerre  : I -> SUP I
    ++ laguerre(n) returns the nth Laguerre polynomial \spad{L[n](x)}.
    ++ Note: Laguerre polynomials, denoted \spad{L[n](x)}, are computed from
    ++ the two term recurrence.  The generating function is:
    ++ \spad{exp(x*t/(t-1))/(1-t) = sum(L[n](x)*t**n/n!, n=0..infinity)}.
  legendre  : I -> SUP RN
    ++ legendre(n) returns the nth Legendre polynomial \spad{P[n](x)}.
    ++ Note: Legendre polynomials, denoted \spad{P[n](x)}, are computed from
    ++ the two term recurrence.  The generating function is:
    ++ \spad{1/sqrt(1-2*t*x+t**2) = sum(P[n](x)*t**n, n=0..infinity)}.
 Implementation ==> add
  import IntegerPrimesPackage(I)

  x := monomial(1,1)$SUP(I)
  y := monomial(1,1)$SUP(RN)

  -- For functions computed via a fixed term recurrence we record
  -- previous values so that the next value can be computed directly

  E : Record(En:I, Ev:SUP(RN)) := [0,1]
  B : Record( Bn:I, Bv:SUP(RN) ) := [0,1]
  H : Record( Hn:I, H1:SUP(I), H2:SUP(I) ) := [0,1,x]
  L : Record( Ln:I, L1:SUP(I), L2:SUP(I) ) := [0,1,x]
  P : Record( Pn:I, P1:SUP(RN), P2:SUP(RN) ) := [0,1,y]
  CT : Record( Tn:I, T1:SUP(I), T2:SUP(I) ) := [0,1,x]
  U : Record( Un:I, U1:SUP(I), U2:SUP(I) ) := [0,1,0]

  MonicQuotient: (SUP(I),SUP(I)) -> SUP(I)
  MonicQuotient (a,b) ==
    not one? leadingCoefficient(b) => error "divisor must be monic"
    b = 1 => a
    da := degree a
    db := degree b            -- assertion: degree b > 0
    q:SUP(I) := 0
    while da >= db repeat
      t := monomial(leadingCoefficient a, (da-db)::NNI)
      a := a - b * t
      q := q + t
      da := degree a
    q

  cyclotomic n ==
    --++ cyclotomic polynomial denoted phi[n](x)
    p:I; q:I; r:I; s:I; m:NNI; c:SUP(I); t:SUP(I)
    negative? n => error "cyclotomic not defined for negative integers"
    n = 0 => x
    k := n; s := p := 1
    c := x - 1
    while k > 1 repeat
      p := nextPrime p
      (q,r) := divide(k, p)
      if r = 0 then
        while r = 0 repeat (k := q; (q,r) := divide(k,p))
        t := multiplyExponents(c,p::NNI)
        c := MonicQuotient(t,c)
        s := s * p
    m := (n quo s) :: NNI
    multiplyExponents(c,m)

  euler n ==
    p : SUP(RN); t : SUP(RN); c : RN; s : I
    negative? n => error "euler not defined for negative integers"
    if n < E.En then (s,p) := (0$I,1$SUP(RN)) else (s,p) := E
    -- (s,p) := if n < E.En then (0,1) else E
    for i in s+1 .. n repeat
      t := (i::RN) * integrate p
      c := euler(i)$IntegerNumberTheoryFunctions / 2**(i::NNI) - t(1/2)
      p := t + c::SUP(RN)
    E.En := n
    E.Ev := p
    p

  bernoulli n ==
    p : SUP RN; t : SUP RN; c : RN; s : I
    negative? n => error "bernoulli not defined for negative integers"
    if n < B.Bn then (s,p) := (0$I,1$SUP(RN)) else (s,p) := B
    -- (s,p) := if n < B.Bn then (0,1) else B
    for i in s+1 .. n repeat
      t := (i::RN) * integrate p
      c := bernoulli(i)$IntegerNumberTheoryFunctions
      p := t + c::SUP(RN)
    B.Bn := n
    B.Bv := p
    p

  fixedDivisor a ==
    g:I; d:NNI; SUP(I)
    d := degree a
    g := coefficient(a, minimumDegree a)
    for k in 1..d while g > 1 repeat g := gcd(g,a k)
    g

  hermite n ==
    s : I; p : SUP(I); q : SUP(I)
    negative? n => error "hermite not defined for negative integers"
    -- (s,p,q) := if n < H.Hn then (0,1,x) else H
    if n < H.Hn then (s := 0; p := 1; q := x) else (s,p,q) := H
    for k in s+1 .. n repeat (p,q) := (2*x*p-2*(k-1)*q,p)
    H.Hn := n
    H.H1 := p
    H.H2 := q
    p

  legendre n ==
    s:I; t:I; p:SUP(RN); q:SUP(RN)
    negative? n => error "legendre not defined for negative integers"
    -- (s,p,q) := if n < P.Pn then (0,1,y) else P
    if n < P.Pn then (s := 0; p := 1; q := y) else (s,p,q) := P
    for k in s+1 .. n repeat
      t := k-1
      (p,q) := ((k+t)$I/k*y*p - t/k*q,p)
    P.Pn := n
    P.P1 := p
    P.P2 := q
    p

  laguerre n ==
    s:I; t:I; p:SUP(I); q:SUP(I)
    negative? n => error "laguerre not defined for negative integers"
    -- (s,p,q) := if n < L.Ln then (0,1,x) else L
    if n < L.Ln then (s := 0; p := 1; q := x) else (s,p,q) := L
    for k in s+1 .. n repeat
      t := k-1
      (p,q) := ((((k+t)$I)::SUP(I)-x)*p-t**2*q,p)
    L.Ln := n
    L.L1 := p
    L.L2 := q
    p

  chebyshevT n ==
    s : I; p : SUP(I); q : SUP(I)
    negative? n => error "chebyshevT not defined for negative integers"
    -- (s,p,q) := if n < CT.Tn then (0,1,x) else CT
    if n < CT.Tn then (s := 0; p := 1; q := x) else (s,p,q) := CT
    for k in s+1 .. n repeat (p,q) := ((2*x*p - q),p)
    CT.Tn := n
    CT.T1 := p
    CT.T2 := q
    p

  chebyshevU n ==
    s : I; p : SUP(I); q : SUP(I)
    negative? n => error "chebyshevU not defined for negative integers"
    if n < U.Un then (s := 0; p := 1; q := 0) else (s,p,q) := U
    for k in s+1 .. n repeat (p,q) := ((2*x*p - q),p)
    U.Un := n
    U.U1 := p
    U.U2 := q
    p

