polynomials.scala

#!/bin/bash

exec /opt/scala-2.10/bin/scala -feature "$0" "$@"
!#
import scala.annotation.tailrec
import scala.language.postfixOps

/**
 * From https://stackoverflow.com/questions/25129721/scala-memoization-how-does-this-scala-memo-work.
 *
 * Generic way to create memoized functions (even recursive and multiple-arg ones)
 *
 * @param f the function to memoize
 * @tparam I input to f
 * @tparam K the keys we should use in cache instead of I
 * @tparam O output of f
 */
case class Memo[I <% K, K, O](f: I => O) extends (I => O) {
  import scala.collection.mutable.{Map => Dict}

  type Input = I
  type Key = K
  type Output = O

  val cache = Dict.empty[K, O]

  override def apply(x: I) = cache getOrElseUpdate(x, f(x))
}

object Memo {
  /**
   * Type of a simple memoized function e.g. when I = K
   */
  type ==>[I, O] = Memo[I, I, O]
}

class CartesianProduct(val product: Traversable[Traversable[_ <: Any]]) {
  override def toString(): String = {
    product.toString
  }

  def *(rhs: Traversable[_ <: Any]): CartesianProduct = {
    val p = product.flatMap { lhs =>
      rhs.map { r =>
        lhs.toList :+ r
      }
    }

    new CartesianProduct(p)
  }
}

object CartesianProduct {
  def apply(traversable: Traversable[_ <: Any]): CartesianProduct = {
    new CartesianProduct(
      traversable.map { t =>
        Traversable(t)
      }
    )
  }
}

implicit class _RichInt(n: Int) {
  private def fact(n: Long): BigInt = {
    lazy val f: Memo.==>[(BigInt, Long), BigInt] = Memo {
      case (a, 0) => a
      case (a, 1) => a
      case (a, n) => f(a * n, n - 1)
    }

    f(BigInt(1), n)
  }

  def ! = fact(n)
  def P(k: Long): BigInt = fact(n) / fact(n - k)
  def C(k: Long): Long = (P(k) / fact(k)).toLong
}

def gcd(lhs: Long, rhs: Long): Long = {
  if (rhs == 0) {
    lhs.abs
  } else {
    gcd(rhs, lhs % rhs)
  }
}

def lcm(lhs: Long, rhs: Long): Long = {
  (lhs * rhs).abs / gcd(lhs, rhs)
}

object Primes {
  private lazy val notDivisibleBy2: Stream[Long] = 3L #:: notDivisibleBy2.map(_ + 2)
  private lazy val notDivisibleBy2Or3: Stream[Long] = notDivisibleBy2
      .grouped(3)
      .map(_.slice(1, 3))
      .flatten
      .toStream
  private lazy val notDivisibleBy2Or3Or5: Stream[Long] = notDivisibleBy2Or3
      .grouped(10)
      .map { g =>
        g.slice(1, 7) ++ g.slice(8, 10)
      }
      .flatten
      .toStream

  lazy val primes: Stream[Long] = 2L #::
      notDivisibleBy2.head #::
      notDivisibleBy2Or3.head #::
      notDivisibleBy2Or3Or5.filter { i =>
        i < 49 || primes.takeWhile(_ <= Math.sqrt(i).toLong).forall(i % _ != 0)
      }

  def <(n: Long): Stream[Long] = primes.takeWhile(_ <= n)

  def largestPower(n: Long, p: Long): Int = {
    @tailrec
    def _largestPower(x: Int, n: Long): Int = {
      if (n % p != 0) {
        x
      } else {
        _largestPower(x + 1, n / p)
      }
    }

    _largestPower(0, n)
  }

  def factorize(n: Long): Map[Long, Int] = {
    val primeFactorsOfN = (Primes < n).filter(n % _ == 0)

    primeFactorsOfN.map { p =>
      (p, largestPower(n, p))
    }.toMap
  }
}

object Factors {
  private def dotPower(set: Iterable[(Long, Int)]): Long = {
    set.map { case (p, x) =>
      Math.pow(p, x).toLong
    }.product
  }

  private def lowerFactors(primeFactors: Map[Long, Int]): Set[Long] = {
    val primeFactorsSpanPoint = primeFactors.filter { e =>
      e._2 == primeFactors.map(_._2).max
    }.head._1

    factors(primeFactors + (primeFactorsSpanPoint -> (primeFactors(primeFactorsSpanPoint) + 1) / 2))
        .filter(_ <= Math.sqrt(dotPower(primeFactors)))
  }

  private def factors(primeFactors: Map[Long, Int]): Set[Long] = {
    def cross(lhs: List[Int], rhs: List[List[Int]]): List[List[Int]] = {
      for { l <- lhs; r <- rhs } yield l :: r
    }

    primeFactors
        .map(x => Range(0, x._2 + 1).toList)
        .foldRight(List(List[Int]()))(cross _)
        .map { x: List[Int] =>
          dotPower(primeFactors.keys.zip(x))
        }.toSet
  }

  def factorPairs(n: Long): Set[(Long, Long)] = {
    if (n == 1) {
      Set((1, 1))
    } else {
      val primeFactorsOfN = Primes.factorize(n)

      val lf = lowerFactors(primeFactorsOfN)
          .map { f =>
            (f, n / f)
          }

      lf ++ lf.filter(p => p._1 != p._2).map(_.swap)
    }
  }
}

case class Polynomial(val coefficients: Long*) {
  override def toString: String = {
    def _toString(coefficients: Seq[Long]): String = {
      def toTermString(coefficient: Long, exponent: Int): String = {
        val coefficientPart = if (coefficient == 1 && exponent != 0) {
          ""
        } else {
          coefficient.toString
        }
        val variablePart = if (exponent == 0) {
          ""
        } else if (exponent == 1) {
          "n"
        } else {
          s"n^${exponent}"
        }

        s"${coefficientPart}${variablePart}"
      }

      val coefficientExponentPairs = coefficients
          .zipWithIndex
          .reverse
          .filter(_._1 != 0)
      val head = coefficientExponentPairs.head
      val tail = coefficientExponentPairs.tail

      val firstTerm = toTermString(head._1, head._2)
      tail
          .foldLeft(firstTerm) { case (lhs, rhs) =>
            val (coefficient, exponent) = rhs
            val operator = if (coefficient < 0) {
              "-"
            } else {
              "+"
            }

            s"${lhs} ${operator} ${toTermString(coefficient.abs, exponent)}"
          }
    }

    val factors = this.factor
    if (factors.size == 1 && factors.head._2 == 1) {
      _toString(coefficients)
    } else {
      factors.toList.sortWith { (lhs, rhs) =>
        lhs._1.order > rhs._1.order ||
          lhs._1.order == rhs._1.order && lhs._2 > rhs._2 ||
          lhs._1.order == rhs._1.order && lhs._2 == rhs._2 && lhs._1.coefficients.zip(rhs._1.coefficients)
              .dropWhile(e => e._1 == e._2)
              .take(1)
              .map(e => e._1 > e._2)
              .head
      }.map { case (polynomial, exponent) =>
        val p = if (polynomial.numberOfTerms == 1) {
          _toString(polynomial.coefficients)
        } else {
          s"(${polynomial})"
        }

        val e = if (exponent == 1) {
          ""
        } else {
          s"^${exponent}"
        }

        s"${p}${e}"
      }
    }.mkString(" × ")
  }

  def apply(n: Int): BigInt = {
    coefficients
        .zipWithIndex
        .map { case (coefficient, i) =>
          coefficient * BigInt(n).pow(i)
        }.sum
  }

  def ==(that: Polynomial): Boolean = this.coefficients == that.coefficients

  def -(that: Polynomial): Polynomial = {
    val length = Math.max(coefficients.length, that.coefficients.length)
    val lhs = coefficients.padTo(length, 0L)
    val rhs = that.coefficients.padTo(length, 0L)

    Polynomial(lhs.zip(rhs).map { case (l, r) =>
      l - r
    } :_*)
  }

  def /(n: Long): Polynomial = {
    Polynomial(coefficients.map(_ / n) :_*)
  }

  def order(): Int = {
    coefficients.length - 1
  }

  def numberOfTerms(): Int = {
    coefficients
        .filter(_ > 0)
        .length
  }

  def factor(): Map[Polynomial, Int] = {
    // TODO: @tailrec
    def _factor(coefficients: Seq[Long]): Seq[Polynomial] = {
      if (coefficients(0) == 0) {
        Polynomial(0, 1) +: _factor(coefficients.tail)
      } else if (coefficients.length < 3) {
        Seq(Polynomial(coefficients :_*))
      } else {
        val highestPowerCoefficientFactorPairs = Factors.factorPairs(coefficients.last).toList.sorted.reverse
        val lowestPowerCoefficientFactorPairs = Factors.factorPairs(coefficients.head).toList.sorted

        println(s"0: ${coefficients.last}, ${highestPowerCoefficientFactorPairs}; ${coefficients.head}, ${lowestPowerCoefficientFactorPairs}")

        (CartesianProduct(highestPowerCoefficientFactorPairs) * lowestPowerCoefficientFactorPairs).product.flatMap { case List(hpcfp, lpcfp) =>
          println(s"1: ${hpcfp}, ${lpcfp}")
          val lhsCoefficients = (hpcfp.asInstanceOf[(Long, Long)]._2, lpcfp.asInstanceOf[(Long, Long)]._2)
          val rhsCoefficients = this.rhsCoefficients(
              coefficients,
              lhsCoefficients,
              (hpcfp.asInstanceOf[(Long, Long)]._1, lpcfp.asInstanceOf[(Long, Long)]._1))

          println(s"2: ${lhsCoefficients}, ${rhsCoefficients}")

          if (rhsCoefficients == None) {
            None
          } else {
            Some((lhsCoefficients, rhsCoefficients))
          }
        }.flatMap { case (lhsCoefficients, rhsCoefficients) =>
          Polynomial(List(lhsCoefficients._1, lhsCoefficients._2) :_*) +: _factor(rhsCoefficients.get)
        }.toSeq
      }
    }

    val f = _factor(coefficients)
        .groupBy(p => p)
        .mapValues(_.length)
        .map { case (p, e) =>
          if (p == Polynomial(0L, 1L)) {
            (Polynomial(List.fill(e)(0L) :+ 1L :_*), 1)
          } else {
            (p, e)
          }
        }

    if (f.size == 1) {
      f
    } else {
      f.filter { case (p, e) =>
        p != Polynomial(1)
      }
    }
  }

  private def rhsCoefficients(
      coefficients: Seq[Long],
      lhsHeadLastCoefficients: (Long, Long),
      rhsHeadLastCoefficients: (Long, Long)): Option[Seq[Long]] = {
    @tailrec
    def _rhsCoefficients(
        accum: Seq[Long],
        reverseProductCoefficients: Seq[Long]): Option[Seq[Long]] = {
      if (reverseProductCoefficients.length == 2) {
        println(s"5: ${reverseProductCoefficients}, ${accum}")
        Some(accum)
      } else {
        val numerator = reverseProductCoefficients.head - accum.head * lhsHeadLastCoefficients._1
        val denominator = lhsHeadLastCoefficients._2

        println(s"4: ${reverseProductCoefficients.head}, ${accum.head}, ${lhsHeadLastCoefficients._1}, ${numerator}, ${denominator}")

        if (numerator % denominator != 0) {
          None
        } else {
          _rhsCoefficients(numerator / denominator +: accum, reverseProductCoefficients.tail)
        }
      }
    }

    println(s"3: ${coefficients}, ${lhsHeadLastCoefficients}, ${rhsHeadLastCoefficients}")

    _rhsCoefficients(Seq(rhsHeadLastCoefficients._2), coefficients.reverse.tail)
        .filter { rhsC =>
          println(s"6: ${coefficients.tail.head} == ${rhsC.head} * ${lhsHeadLastCoefficients._2} + ${rhsC.tail.headOption.getOrElse(1L)} * ${lhsHeadLastCoefficients._1}")
          println(s"7: ${rhsC}, ${lhsHeadLastCoefficients}\n")
          coefficients.tail.head == rhsC.head * lhsHeadLastCoefficients._2 + rhsC.tail.headOption.getOrElse(1L) * lhsHeadLastCoefficients._1
        }.map { rhsC =>
          rhsHeadLastCoefficients._1 +: rhsC
        }
  }
}

implicit class _RichLongWithPolynomial(n: Long) {
  def *(polynomial: Polynomial): Polynomial = {
    Polynomial(polynomial.coefficients.map(n * _) :_*)
  }
}

//println(Polynomial(1, 3, 3, 1).factor) // Map(n + 1 -> 3)
//println(Polynomial(1, 3, 3, 1))
//println(Polynomial(0, 0, 0, 1).factor) // Map(n^3 -> 1)
//println(Polynomial(0, 0, 1, 3, 3, 1))
println(Polynomial(24, 31, 10).factor) // FIXME: should be Map(2n + 3 -> 1, 5n + 8 -> 1) but is Map()
//1: (2,5), (3,8)
//3: WrappedArray(24, 31, 10), (5,8), (2,3)
//5: WrappedArray(31, 24), List(3)
//6: 31 == 3 * 8 + 1 * 5
//7: List(3), (5,8)
//2: (5,8), None