Programs that manipulate infinities

Programs that manipulate infinities

Written by

Dylan Holmes

You might know that programming languages allow you to manipulate genuinely infinite lists by deferring calculations until they're needed. So in Python, you can create a list of all primes using: 

def primes() :
    yield 2
    k = 3
    while 1 :
        for p in primes() :
            if p*p > k :
                yield k
                break
            if k % p == 0 :
                break
        k += 1

g = primes()
print([next(g) for _ in range(10)])

And now you can pull out as many primes as you need. You can do amazing, seemingly impossible, manipulations with a finite computer operating on deferred infinite calculations.

1 Generating functions

A generating function is a mathematical trick that replaces an infinite sequence of numbers \(a_0, a_1,a_2,\ldots\) with an infinitely long polynomial \(a_0x^0 + a_1x^1 + a_2x^2 + \ldots\). It turns out that the properties of the polynomial (its derivatives, its roots, its closed form) can reveal important properties of the original sequence.

Sometimes a complicated sequence has a compact closed form. For example, the fibonacci sequence [0,1,1,2,3,5,8,13,25,...] becomes the polynomial \(0 + x + x^2 + 2x^3 + 3x^4 + 5x^5 + 8x^6 + 13x^7 + \ldots\), which has a compact closed form:

$$f(x) = \frac{x}{1-x-x^2}$$

So sequences can be described in many equivalent ways:

  1. As an infinite sequence of numbers, with an exact formula for any term. For the Fibonacci sequence, that's $$a_n = \frac{1}{\sqrt{5}}\left(\frac{1+\sqrt{5}}{2}\right)^n-\frac{1}{\sqrt{5}}\left(\frac{1-\sqrt{5}}{2}\right)^n.$$
  2. As a power series, an infinitely long polynomial whose coefficients come from the sequence of numbers.
  3. As a recurrence relation, telling you how to make the next term as a function of some previous terms. For the Fibonacci sequence, that's \(a_{n+1} = a_n + a_{n-1}\), with \(a_0 = 0\) and \(a_1 = 1\).
  4. As a closed form ratio of polynomials. The infinite power series for the Fibonacci sequence is equivalent to the function \(f(x) = x/(1-x-x^2).\)

I wrote Clojure software for converting between sequences, recurrence relations, and closed form rational functions like these.

2 The cool stuff

In my software, a polynomial is a vector of numbers \([a_0, a_1, \ldots, a_n]\). This represents the polynomial \(a_0x^0 + \ldots + a_nx^n\). A power series is an infinitely long polynomial, represented by a lazy infinite sequence of coefficients \((a_0, a_1, \ldots)\).

2.1 From generating functions to infinite sequences

You can plug in the generating function \(x/(1-x-x^2)\) for the fibonacci sequence and it will compute the coefficients for you: 

(take 10 (expand-geometric [0 1] [1 -1 -1]))
;; (0 1 1 2 3 5 8 13 21 34)

You can expand other generating functions, like the one for powers of two: 

(take 10 (expand-geometric [1] [1 -2]))
;; (1 2 4 8 16 32 64 128 256 512)

2.2 From recurrences to generating functions

Suppose you know the linear recurrence relation that your sequence satisfies, such as: \(a_{n+1} = a_n + a_{n-1}.\) Put all the terms on the same side, and you have \(1\cdot a_n + 1\cdot a_{n-1} + (-1)\cdot a_{n+1} = 0.\)

You can represent the recurrence as a list of numbers [1, 1, -1]. The initial conditions \(a_0 = 0, a_1 = 1\) are also a list of numbers [0,1]. To be well-defined, the initial conditions should have exactly one less member than the recurrence list.

Given this linear recurrence [1,1,-1] with initial conditions [0,1], you can convert it into a generating function: 

(linear-recurrence-generator [-1 -1 1] [0 1])
;; returns: [[0 1] [1 -1 -1]]

The return value is a pair of polynomials \(P,Q\) so that \(P/Q\) is the generating function for the sequence that satisfies the recurrence. And indeed, that's the correct generating function for the Fibonacci sequence starting from 0.

The coolest part is that the programming language and mathematics conspire to make the definition of linear-recurrence-generator very simple. Behold, a lucid one-liner: 

(defn linear-recurrence-generator
  [coeffs init-values]
  [init-values (vec (reverse (monic coeffs)))])

The only innovation is monic, which simply divides a polynomial by its highest-term coefficient so that the highest coefficient is 1.

2.3 Finishing the sequence with the … operator.

The … operator takes a finite sequence of numbers and completes it as an infinite sequence. Here are some examples: 

(take 15 (… 1 2 3 4 5))  ;; counting numbers
;; returns: (1 2 3 4 5 6 7 8 9 10 11 12 13 14 15)

(take 15 (… 2 4 6 8 10))  ;; even numbers
;; returns: (2 4 6 8 10 12 14 16 18 20 22 24 26 28 30)

(take 15 (… 1 2 4 8)) ;; powers of two
;; returns: (1 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384)

(take 15 (… 2 6 18 54)) ;; doubled powers of three
;; returns: (2 6 18 54 162 486 1458 4374 13122 39366 118098 354294 1062882 3188646 9565938)

(take 15 (… 1 1 2 3 5)) ;; fibonnaci
;; returns: (1 1 2 3 5 8 13 21 34 55 89 144 233 377 610)

(take 15 (… 1 1 2 4 7 13)) ;; "tribonacci", each term sum of previous three
;; returns: (1 1 2 4 7 13 24 44 81 149 274 504 927 1705 3136)

(take 15 (… 1 2 3 1 2 3)) ;; cycle of 1,2,3
;; returns: (1 2 3 1 2 3 1 2 3 1 2 3 1 2 3)

(take 15 (… 1)) ;; constantly 1
;; returns: (1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)

(take 15 (… 1 3 6 10 15)) ;; sum of the first n numbers
;; returns: (1 3 6 10 15 21 28 36 45 55 66 78 91 105 120)

(take 15 (… 1 4 9 16)) ;; perfect squares
;; returns: (1 4 9 16 25 36 49 64 81 100 121 144 169 196 225)

(take 15 (… 1 1 2 2 3 3))  ;; double-counting
;; returns: (1 1 2 2 3 3 4 4 5 5 6 6 7 7 8)

(take 15 (… 1 0 2 0 3 0 4 0)) ;; alternate counting with zeroes
;; returns: (1 0 2 0 3 0 4 0 5 0 6 0 7 0 8)

Its behavior is wonderfully intuitive on a broad range of cases, and pretty magical to see. Under the hood, these results are accomplished by some straightforward linear algebra. It works by searching for a linear relationship between the first few terms and the next term, which you can do by solving a system of linear equations. The smallest recurrence that fits will complete the sequence. If a linear recurrence doesn't fit, it tries a polynomial fit of the terms—this captures sequences like the triangular numbers and the perfect squares.

2.4 And onward, to infinity

There are many other cool bits to the code, including a partial fractions algorithm and a program that finds the gcd for two polynomials.

I'll end with one final highlight—in this representation for polynomials, the definition of derivative is especially simple: 

(defn derivative
  [coeffs]
  (rest (map-indexed * coeffs)))

3 SOURCE

Code is GPL3 with attribution to me.

3.1 Version A

;; a polynomial is a vector of coefficients [a0,a1,a2,...,a_n]

(ns series)


;;;;;  POWER SERIES UTILITIES

(def finite?
  "Return true if the polynomial has finite degree."
  vector?)

(defn scale
  "Multiply the polynomial by a number."
  ([a coeffs]
   (if (finite? coeffs)
     (mapv (partial * a) coeffs)
     (map (partial * a) coeffs))))

(def constant?
  "Return true if the finite polynomial is constant."
  (comp not pos? dec count))

(defn term
  [coeffs n]
  "Return the nth term of the polynomial."
  (nth coeffs n 0))

(defn index-highest
  "Return the index of the highest order term, if finite. Otherwise, return inf."
  [coeffs]
  (if (finite? coeffs)
    (last (keep-indexed #(when-not (zero? %2) %1) coeffs))
    ##Inf))

  (def degree "Numerical degree of the polynomial, i.e. the index of the highest-order term. Infinite for infinite power series.

  Alias for index-highest." index-highest)

(defn index-lowest-nonzero
  "Return the index of the lowest order term, excluding the constant term. If polynomial is constant, return -inf instead."
  [coeffs]
  (or (first (map inc (keep-indexed #(when-not (zero? %2) %1) (rest coeffs)))) ##-Inf))


(defn monic
  "Divide the polynomial by its highest-order coefficient so that its
  highest-order coefficient is 1. If the polynomial is zero, return it."
  [coeffs]
  (if-let [n (index-highest coeffs)]
    (scale (/ (term coeffs n)) coeffs)))



;;;;;  ARITHMETIC WITH POLYNOMIALS
(defn poly-sum
  "Add polynomials."
  ([] [])
  ([coeffs] coeffs)
  ([coeffs-1 coeffs-2]
   (cond (empty? coeffs-1) coeffs-2
         (empty? coeffs-2) coeffs-1
         (and (finite? coeffs-1) (finite? coeffs-2))
         (remove-leading-zeroes
          (mapv +
               (concat coeffs-1 (repeat (max 0 (- (count coeffs-2) (count coeffs-1))) 0))
               (concat coeffs-2 (repeat (max 0 (- (count coeffs-1) (count coeffs-2))) 0))))
         :default
         (lazy-seq
          (cons (+ (first coeffs-1) (first coeffs-2))
                (poly-sum (rest coeffs-1) (rest coeffs-2)))))))

(defn poly-diff
  "Subtract polynomials."
  ([] [])
  ([coeffs] (scale -1 coeffs))
  ([coeffs-1 & others]
   (reduce poly-sum coeffs-1 (map (partial scale -1) others))))


(defn remove-leading-zeroes [coeffs]
  "Remove the leading zeroes from the polynomial, returning a
  polynomial whose highest order term is nonzero. If the polynomial is
  infinite, return it unmodified."
  (if (finite? coeffs)
    (vec (reverse (drop-while zero? (reverse coeffs))))
    coeffs))

(defn poly-product
  "Multiply the polynomials."
  ([] [1])
  ([coeffs] (vec coeffs))
  ([coeffs-1 coeffs-2]
   ;; let's do this in a way that handles infinite polynomials nicely.
   (if (and (finite? coeffs-1) (finite? coeffs-2))
     (vec (remove-leading-zeroes
           (reduce poly-sum
                   (for [i (range (count coeffs-1))]
                     (scale (nth coeffs-1 i) (concat (repeat i 0) coeffs-2))))))
     (map
      (fn [index]
        (reduce +
                (for [i (range (inc index))]
                  (* (term coeffs-1 i) (term coeffs-2 (- index i))))))
      (range)))))



(defn poly-power
  "Raise the polynomial to the nth power."
  [coeffs n]
  (reduce poly-product (repeat n coeffs)))




;;;;;  EXPANDING RATIONAL FUNCTIONS INTO POWER SERIES

(defn powers
  "An infinite list of poly^0 poly^1 poly^2, etc."
  [coeffs]
  (iterate (partial poly-product coeffs) [1]))   

(defn expand-geometric
  "Return the power series expansion of sum_(n=0)^infty (P1/P2)^n.
  When P2 is not supplied, use 1 as default."

  ([coeffs]
   ;;(println "attempting to expand" coeffs)
  ((fn loop [powers index]
     (if (and (not (zero? (term coeffs 0)))
              (not (zero? (term coeffs index))))
       (lazy-seq
        (cons ##Inf ;; infinitely many terms appear here.
              (loop powers (inc index))))
       (lazy-seq
        (let [;;_ (println (take 5 powers))
              powers* (drop-while #(< (index-highest %) index) powers)
              powers** (take-while #(<= (index-lowest-nonzero %) index) powers*)]
          (cons (reduce + (map #(term % index) powers**))
                (loop powers* (inc index))))))
     )
   (powers coeffs) 0))

  ([coeffs-1 coeffs-2]
   (let [const (term coeffs-2 0)]
     (->> coeffs-2
          (rest)(cons 0)(vec)     ;; zero out the constant term
          (scale (/ -1 const))    ;; scale by constant term
          (expand-geometric)      ;; use 1/(1-f(x)) = 1 + f + f^2 + f^3 + ...
          (poly-product coeffs-1) ;; reintroduce numerator
          (scale (/ const))       ;; rescale by constant term
          ))))


(def seq-fibonacci (expand-geometric [0 1] [1 -1 -1])) ;; 1/(1-x-x^2)
(def seq-powers-2 (expand-geometric [1] [1 -2])) ;; 1/(1-2x)





;;;;;  EUCLIDEAN DIVISION ALGORITHMS, SQUAREFREE FACTORIZATION, PARTIAL FRACTIONS
(defn derivative
  "Return the derivative of the power series as another power series."
  [coeffs]
  ((if (finite? coeffs) vec identity)
   (rest (map-indexed * coeffs))))

(defn monomial
  "Return a polynomial representing ax^n (by default, a=1.)"
  ([n a]
   (vec (concat (repeat n 0) [a])))
  ([n] (monomial n 1)))

(defn poly-divide
  "Divide a high degree polynomial by a low degree polynomial, producing
  a pair of polynomials [quotient remainder] as a result."
  [coeffs-1 coeffs-2]
  ;; (assert (>= (index-highest coeffs-1) (index-highest coeffs-2))
  ;;         "Polynomial division requires degree(dividend) >= degree(divisor).")

  (loop [divisor coeffs-2,
         dividend coeffs-1,
         quotient []]
    (if (every? zero? dividend) {:quotient quotient :remainder []}
    (let [i (index-highest divisor),   a (term divisor i),
          j (index-highest dividend),  b (term dividend j),

          q (monomial (- j i) (/ b a))

          ]
      (if (= i j) {:quotient (poly-sum q quotient),
                   :remainder (poly-diff dividend (poly-product q divisor))}
          (recur divisor
                 (poly-diff dividend
                            (poly-product q divisor))
                 (poly-sum quotient q)))))))


(comment :example "from wikipedia on polynomial long division"
         (poly-divide [-4 0 -2 1] [-3 1]))

(defn poly-gcd
  "Euclidean algorithm finds the gcd polynomial of the two given polynomials."
  [coeffs-1 coeffs-2]
  (cond (every? zero? coeffs-1) coeffs-2
        (every? zero? coeffs-2) coeffs-1
        (< (index-highest coeffs-1) (index-highest coeffs-2))
        (poly-gcd coeffs-2 coeffs-1)
        :default
        (let [m (poly-divide coeffs-1 coeffs-2)]
          (if (every? zero? (:remainder m))
            (monic coeffs-2)
            (recur coeffs-2 (:remainder m))))))

(defn poly-gcd*
  "Extended euclidean algorithm returns a map containing the monic
  polynomial :gcd of the two polynomials, as well as the unique
  smallest polynomials q_1 and q_2 such that poly1 * q1 + poly2 * q2 =
  gcd."
  [coeffs-1 coeffs-2]
  (if (< (index-highest coeffs-1) (index-highest coeffs-2))
    (let [subproblem (poly-gcd* coeffs-2 coeffs-1)]
      {:gcd (:gcd subproblem)
       :bezout-1 (:bezout-2 subproblem)
       :bezout-2 (:bezout-1 subproblem)})
    (loop [r1 coeffs-1, r2 coeffs-2,
           s1 [1], s2 [],
           t1 [],  t2 [1]]
      (if (every? zero? r2)
        (let [lead (term r1 (index-highest r1))]
          {:gcd (scale (/ lead) r1)
           :bezout-1 (scale (/ lead) s1)
           :bezout-2 (scale (/ lead) t1)})
          (let [q (:quotient (poly-divide r1 r2))]
            (recur r2 (poly-diff r1 (poly-product r2 q))
                   s2 (poly-diff s1 (poly-product s2 q))
                   t2 (poly-diff t1 (poly-product t2 q))))))))

(comment :example
         "https://www.youtube.com/watch?v=Y3CQn_KOAkk"
         (poly-gcd* [-2 -3 -2 0 1] [1 4 4 1]))


(defn extended-gcd
  "Extended gcd algorithm for integers."
  [m n]
  (if (< m n) (recur n m)
        (loop [r1 m, r2 n,
               s1 1, s2 0,
               t1 0, t2 1]
          ;;(println r1 r2 \tab s1 \tab t1)
          (if-not (zero? r2) (println s1 "-" (quot r1 r2) "x" s2 "=" (- s1 (* (quot r1 r2) s2))))
          (if (zero? r2) {:gcd r1 :bezout-1 s1 :bezout-2 t1}
              (let [q (quot r1 r2) r (rem r1 r2)]
                (recur r2, (- r1 (* q r2)),
                       s2, (- s1 (* q s2)),
                       t2 (- t1 (* q t2))))))))



          )))



(comment :example "from wikipedia on polynomial gcd"
         "https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor"
         (poly-gcd [6 7 1] [-6 -5 1]))


;; ref: https://people.eecs.berkeley.edu/~fateman/282/F%20Wright%20notes/week6.pdf
(defn squarefree
  "Factor the polynomial into squarefree components. Returns a map
  associating each squarefree term with the number of times it
  appears. (Yun's algorithm)"
  [coeffs]
  (let [g₁ (poly-gcd coeffs (derivative coeffs))
        h₁ (:quotient (poly-divide coeffs g₁))]
    (loop [term-number 1, terms {}, g g₁, h h₁]
      (println term-number terms g h)
      (if (= [1] h) (into {} (remove (comp constant? key) terms))
          (let [h* (poly-gcd g h)
                g* (:quotient (poly-divide g h*))]
            (recur (inc term-number)
                   (assoc terms (:quotient (poly-divide h h*)) term-number)
                   g* h*))))))

(comment :example
         (squarefree [1 2 1]))



;; ref: https://www.eecs.harvard.edu/~htk/publication/1977-siam-journal-on-computing-kung-tong.pdf
(defn partial-fractions
  "Compute the partial fraction decomposition of poly1/poly2. Return a
  list of polynomial pairs [p_i q_i] representing rational functions p_i/q_i. "
  [coeffs-1 coeffs-2]
  (if (> (degree coeffs-1) (degree coeffs-2))
    (let [m (poly-divide coeffs-1 coeffs-2)]
      (cons [(:quotient m) [1]]
            (partial-fractions (:remainder coeffs-1) coeffs-2)))
  (let [denoms (squarefree coeffs-2)
        Ri (map (fn [[k v]] (poly-power k v)) denoms)
        R
        (reduce poly-sum
        (for [i (range (count Ri))]
          (reduce poly-product
                  (for [j (range (count Ri))
                        :when (not= i j)]
                    (nth Ri j)))))
        Di (map #(:remainder (poly-divide R %)) Ri)
        Ai (map (comp :bezout-1 poly-gcd*) Di Ri)
        Ki (map (comp :remainder #(poly-divide coeffs-1 %)) Ri)

        Ci (map (comp :remainder
                      #(poly-divide (poly-product %1 %2) %3))
                      Ki Ai Ri)
        ]
    ;; Todo: solve case 2 for C_i/Q_i^n_i by recursive halving of n_i.
    (map vector Ci Ri)

    )))


;; https://www2.math.upenn.edu/~wilf/gfology2.pdf
;; p7 : (1 - 2x + 2x^2) / (1-x)^2(1-2x)
;; (partial-fractions [1 -2 2] [1 -4 5 -2])





;;;;;  LINEAR RECURRENCES

;; Next: how do you get a general /formula/ for the ith term in the
;; geometric expansion?  If the term is a monomial ax^k, it's a little
;; easier. The nth term is 0 unless n divides k, in which case it's a^(n/k).

;; Then: how do you convert a recurrence relation into a solution? (and how to represent
;; recurrence relations? Macro transformation from recurrence to solution?)
     ;; (scale (/ const)
     ;;        (poly-product coeffs-1
     ;;                      (expand-geometric (scale (/ -1 const)
     ;;                                               (vec (cons 0 (rest coeffs-2))))))))))



(defn linear-recurrence-series
  "Return a series representing the solution to the linear recurrence."
  [coeffs & init-values]
  (assert (= (index-highest coeffs)  (count init-values)) "Mismatch between recurrence size and number of initial values.")
  (let [solve-for-xn (vec (butlast (scale (/ -1 (term coeffs (index-highest coeffs))) coeffs)))]
    (lazy-cat init-values
              (map (fn [terms]
                     (reduce + (map * solve-for-xn terms)))
              (partition (count solve-for-xn) 1
                         (apply linear-recurrence-series coeffs init-values))))))

(defn linear-recurrence-generator
  "Return a pair of polynomials representing the rational generating
  function P/Q for the recurrence.
  The recurrence for a sequence F is Σ_n a_n F_n = 0, with initial conditions F_i=b_i.

  Note that (expand-geometric P Q) yields the terms of this sequence."
  [coeffs init-values]
  [init-values (vec (reverse (monic coeffs)))])

;; ! wow. such a simple form.



(comment :example
         (linear-recurrence-series [1 1 -1] 0 1) ;; fibonacci a_n + a_{n+1} - a_{n+2} = 0
         (linear-recurrence-series [2 -1] 1)) ;; powers of two: 2a_n - a_{n+1} = 0

3.2 Version B

(ns ai.logical.series
  (:gen-class)
  (:import
    [javax.swing SwingUtilities JFrame JPanel JScrollPane JLabel JSpinner JTabbedPane JTextArea SwingConstants]
    [java.awt BorderLayout Color Container]
    javax.swing.event.ChangeEvent))

(import javax.swing.event.ChangeEvent)
(import net.sourceforge.jeuclid.swing.JMathComponent)
;;(import org.clojars.chapmanb)

(defn formula-polynomial
  ([coeffs] (formula-polynomial "x" coeffs))
  ([var coeffs]
   (let [s (clojure.string/join "+"
                  (for [[coeff exponent]
                        (identity (map vector coeffs (range)))
                        :when (not (zero? coeff))]
                    (str "[" coeff "]" "X**{" exponent "}")))
         s (clojure.string/replace s "X**{1}" "X")
         s (clojure.string/replace s "X**{0}" "")
         s (clojure.string/replace s "+[-" "-[")
         s (clojure.string/replace s #"\[(\d+)\]"
                                   #(if (= "1" (% 1)) ""
                                        (str "<mn>" (% 1) "</mn>")))

         s (clojure.string/replace s #"X\*\*\{(\d+)\}" #(str "<msup><mi>" var "</mi><mn>" (% 1) "</mn></msup>"))
         s (clojure.string/replace s #"([\+\-])([X<])" #(str "<mo>" (% 1) "</mo>" (% 2)))
         s (clojure.string/replace s #"\[([\d\-]+)\]" #(str "<mn>" (% 1) "</mn>"))
         s (clojure.string/replace s "X" (str "<mi>" var "</mi>"))
         s (clojure.string/replace s #"\-$" (str "<mo>-</mo><mn>1</mn>"))
         s (clojure.string/replace s #"^<mo>" "<mn>1</mn><mo>")
         s (clojure.string/replace s #"^$" "<mn>1</mn>")
         ;;_ (println ">>>>>>" s)
         ]

  s)))


(defn mathml-wrap [& strings]
   (str "<?xml version=\"1.0\"?>"
        "<!DOCTYPE math PUBLIC \"-//W3C//DTD MathML 2.0//EN\" \"http://www.w3.org/TR/MathML2/dtd/mathml2.dtd\">"
        "<math mode=\"display\">"
        (apply str strings)
        "</math>"))

(defn show-rational
  ([coeffs]
   (show-rational coeffs 1))
  ([coeffs-1 coeffs-2]
     (let [my-frame (doto (JFrame. "My Frame")
                   (.setDefaultCloseOperation JFrame/EXIT_ON_CLOSE))

        formula (doto
                    (new
                     net.sourceforge.jeuclid.swing.JMathComponent))

        _ (.add my-frame formula)
        _ (doto formula
            (.setFontSize 56)
            (.setContent
             (if (= 1 coeffs-2)
               (mathml-wrap (formula-polynomial "x" coeffs-1))
               (mathml-wrap
                "<mfrac><mrow>"
                (formula-polynomial "x" coeffs-1)
                "</mrow><mrow>"
                (formula-polynomial "x" coeffs-2)
                "</mrow></mfrac>"
                ))))]
    (.pack my-frame)
    (.setSize my-frame 600 200)
    (.setVisible my-frame true))))



(defn show-formula
  ([formula]
   (let [my-frame (doto (JFrame. "My Frame")
                    (.setDefaultCloseOperation JFrame/EXIT_ON_CLOSE))

         formula-component (doto
                     (new
                      net.sourceforge.jeuclid.swing.JMathComponent))

         _ (.add my-frame formula-component)
         _ (doto formula-component
            (.setFontSize 56)
            (.setContent (mathml-wrap formula)))]
    (.pack my-frame)
    (.setSize my-frame 600 200)
    (.setVisible my-frame true))))


(def show-polynomial show-rational)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  POWER SERIES UTILITIES
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(derive clojure.lang.PersistentVector ::finite-polynomial)
(derive ::finite-polynomial ::power-series)
(derive clojure.lang.Seqable ::power-series)


(def finite?
  "Return true if the polynomial has finite degree."
  #(isa? (type %) ::finite-polynomial))


(defn monomial
  "Return a polynomial representing ax^n (by default, a=1.)"
  ([n a]
   (vec (concat (repeat n 0) [a])))
  ([n] (monomial n 1)))

(defn scale
  "Multiply the polynomial by a number."
  ([a coeffs]
   (if (finite? coeffs)
     (mapv (partial * a) coeffs)
     (map (partial * a) coeffs))))

(def constant?
  "Return true if the finite polynomial is constant."
  (comp not pos? dec count))

(defn term
  [coeffs n]
  "Return the nth term of the polynomial."
  (nth coeffs n 0))

(defn index-highest
  "Return the index of the highest order term, if finite. Otherwise, return inf."
  [coeffs]
  (if (finite? coeffs)
    (last (keep-indexed #(when-not (zero? %2) %1) coeffs))
    ##Inf))

(def degree "Numerical degree of the polynomial, i.e. the index of the highest-order term. Infinite for infinite power series.

  Alias for index-highest."
  index-highest)

(defn index-lowest-nonzero
  "Return the index of the lowest order term, excluding the constant term. If polynomial is constant, return -inf instead."
  [coeffs]
  (or (first (map inc (keep-indexed #(when-not (zero? %2) %1) (rest coeffs)))) ##-Inf))


(defn monic
  "Divide the polynomial by its highest-order coefficient so that its
  highest-order coefficient is 1. If the polynomial is zero, return it."
  [coeffs]
  (if-let [n (index-highest coeffs)]
    (scale (/ (term coeffs n)) coeffs)))


(defn remove-leading-zeroes [coeffs]
  "Remove the leading zeroes from the polynomial, returning a
  polynomial whose highest order term is nonzero. If the polynomial is
  infinite, return it unmodified."
  (if (finite? coeffs)
    (vec (reverse (drop-while zero? (reverse coeffs))))
    coeffs))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  ARITHMETIC WITH POLYNOMIALS
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defmulti ++ (fn [& args] (vec (map class args))))
(defmethod ++ [Number Number] [x y] (clojure.core/+ x y))
(defmethod ++
  [::finite-polynomial ::finite-polynomial]
  [coeffs-1 coeffs-2]
  (remove-leading-zeroes
   (mapv +
         (concat coeffs-1 (repeat (max 0 (- (count coeffs-2) (count coeffs-1))) 0))
         (concat coeffs-2 (repeat (max 0 (- (count coeffs-1) (count coeffs-2))) 0)))))

(defmethod ++
  [::power-series ::power-series]
  [coeffs-1 coeffs-2]
  (cond (empty? coeffs-2) coeffs-1
        (empty? coeffs-1) coeffs-2
        :default
        (cons (clojure.core/+ (first coeffs-1) (first coeffs-2))
              (++ (rest coeffs-1) (rest coeffs-2)))))

(defmethod ++
  [Number ::finite-polynomial]
  [x y]
  (++ (monomial 0 x) y))

(defmethod ++
  [::finite-polynomial Number ]
  [x y]
  (++ x (monomial 0 y)))



(defmulti * (fn [& args] (vec (map class args))))
(defmethod * [Number] [x] x)
(defmethod * [Number Number] [x y] (clojure.core/* x y))
(defmethod * [Number ::finite-polynomial] [a coeffs]
  (scale a coeffs))
(defmethod * [::finite-polynomial Number ] [coeffs a]
  (scale a coeffs))
(defmethod * [::finite-polynomial ::finite-polynomial]
  [coeffs-1 coeffs-2]
  (->> (for [i (range (count coeffs-1))]
         (scale (nth coeffs-1 i)
                (concat (repeat i 0) coeffs-2)))
       (reduce ++)
       (remove-leading-zeroes)
       (vec)))
(defmethod * [::power-series] [x] x)
(defmethod * [::power-series ::power-series]
  [coeffs-1 coeffs-2]
  (map (fn [index]
         (reduce +
                 (for [i (range (inc index))]
                   (* (term coeffs-1 i) (term coeffs-2 (- index i))))))
   (range)))

(prefer-method *
               [::finite-polynomial ::finite-polynomial]
               [::power-series ::power-series])


(defmulti - (fn [& args] (vec (map class args))))
(defmethod - [Number] [x] (clojure.core/- x))
(defmethod - [Number Number] [x y] (clojure.core/- x y))
(defmethod - [::power-series] [coeffs] (* -1 coeffs))
(defmethod - [::power-series ::power-series] [coeffs-1 coeffs-2]
  (++ coeffs-1 (* -1 coeffs-2)))

(defmethod - :default [x y]
  (++ x (* -1 y)))




(defn poly-divide
  "Divide a high degree polynomial by a low degree polynomial, producing
  a pair of polynomials {:quotient :remainder} as a result."
  [coeffs-1 coeffs-2]
  (loop [divisor coeffs-2,
         dividend coeffs-1,
         quotient []]

    (if (every? zero? dividend) {:quotient quotient :remainder []}
    (let [i (index-highest divisor),   a (term divisor i),
          j (index-highest dividend),  b (term dividend j),

          q (monomial (- j i) (/ b a))
          ]
      (if (= i j) {:quotient (++ q quotient),
                   :remainder (- dividend (* q divisor))}
          (recur divisor
                 (- dividend (* q divisor))
                 (++ quotient q)))))))


(defmulti divide
  "Divide two numbers or two polynomials, producing a map
  {:quotient :remainder} as a result."
  (fn [& args] (vec (map class args))))
(defmethod divide [Number Number] [x y] {:quotient (quot x y),
                                         :remainder (rem x y)})
(defmethod divide [::finite-polynomial ::finite-polynomial]
  [coeffs-1 coeffs-2]
  (loop [divisor coeffs-2,
         dividend coeffs-1,
         quotient []]

    (if (every? zero? dividend) {:quotient quotient :remainder []}
        (let [i (index-highest divisor),   a (term divisor i),
              j (index-highest dividend),  b (term dividend j),

              q (monomial (- j i) (/ b a))
              ]
          (if (= i j) {:quotient (++ q quotient),
                       :remainder (- dividend (* q divisor))}
              (recur divisor
                     (- dividend (* q divisor))
                     (++ quotient q)))))))

(defn poly-power
  "Raise the polynomial to the nth power."
  [coeffs n]
  (reduce * (repeat n coeffs)))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  EXPANDING RATIONAL FUNCTIONS
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  INTO POWER SERIES
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defn powers
  "An infinite list of poly^0 poly^1 poly^2, etc."
  [coeffs]
  (iterate (partial * coeffs) [1]))   

(defn expand-geometric
  "Return the power series expansion of sum_(n=0)^infty (P1/P2)^n.
  When P2 is not supplied, use 1 as default."

  ([coeffs]
  ((fn loop [powers index]
     (if (and (not (zero? (term coeffs 0)))
              (not (zero? (term coeffs index))))
       (lazy-seq
        (cons ##Inf ;; infinitely many terms appear here.
              (loop powers (inc index))))
       (lazy-seq
        (let [powers* (drop-while #(< (index-highest %) index) powers)
              powers** (take-while #(<= (index-lowest-nonzero %) index) powers*)]
          (cons (reduce + (map #(term % index) powers**))
                (loop powers* (inc index))))))
     )
   (powers coeffs) 0))

  ([coeffs-1 coeffs-2]
   (let [const (term coeffs-2 0)]
     (->> coeffs-2
          (rest)(cons 0)(vec)     ;; zero out the constant term
          (scale (/ -1 const))    ;; scale by constant term
          (expand-geometric)      ;; use 1/(1-f(x)) = 1 + f + f^2 + f^3 + ...
          (* coeffs-1)            ;; reintroduce numerator
          (scale (/ const))       ;; rescale by constant term
          ))))


(def seq-fibonacci (expand-geometric [0 1] [1 -1 -1])) ;; 1/(1-x-x^2)
(def seq-powers-2 (expand-geometric [1] [1 -2])) ;; 1/(1-2x)





;;;;;  EUCLIDEAN DIVISION ALGORITHMS, SQUAREFREE FACTORIZATION, PARTIAL FRACTIONS
(defn derivative
  "Return the derivative of the power series as another power series."
  [coeffs]
  ((if (finite? coeffs) vec identity)
   (rest (map-indexed * coeffs))))


(defn identically-zero? [x]
  (or (and (number? x) (zero? x))
      (and (finite? x) (every? zero? x))))


(defn poly-divide
  "Divide a high degree polynomial by a low degree polynomial, producing
  a pair of polynomials {:quotient :remainder} as a result."
  [coeffs-1 coeffs-2]
  ;; (assert (>= (index-highest coeffs-1) (index-highest coeffs-2))
  ;;         "Polynomial division requires degree(dividend) >= degree(divisor).")

  (loop [divisor coeffs-2,
         dividend coeffs-1,
         quotient []]
    (if (identically-zero? dividend) {:quotient quotient :remainder []}
    (let [i (index-highest divisor),   a (term divisor i),
          j (index-highest dividend),  b (term dividend j),

          q (monomial (- j i) (/ b a))

          ]
      (if (= i j) {:quotient (++ q quotient),
                   :remainder (- dividend (* q divisor))}
          (recur divisor
                 (- dividend (* q divisor))
                 (++ quotient q)))))))

(defmulti divide
  "Divide two numbers or two polynomials, producing a map
  {:quotient :remainder} as a result."
  (fn [& args] (vec (map class args))))
(defmethod divide [Number Number] [x y] {:quotient (quot x y),
                                         :remainder (rem x y)})
(defmethod divide [::finite-polynomial ::finite-polynomial]
  [coeffs-1 coeffs-2] (poly-divide coeffs-1 coeffs-2))

(defn poly-power
  "Raise the polynomial to the nth power."
  [coeffs n]
  (reduce * (repeat n coeffs)))


(defn gcd*
  "Extended euclidean algorithm returns a map containing the :gcd
  of the two arguments, as well as the unique smallest q_1 and q_2
  such that q_1 * x + q_2 * y = gcd.

  Works on numbers and polynomials. If given polynomial inputs,
  returns polynomial outputs, including a monic polynomial gcd."
  [x y]
  (if (or (and (finite? x) (finite? y)
               (< (index-highest x) (index-highest y)))
          (and (number? x) (number? y) (< x y)))
    (let [subproblem (gcd* y x)]
      {:gcd (:gcd subproblem)
       :bezout-1 (:bezout-2 subproblem)
       :bezout-2 (:bezout-1 subproblem)})

    (loop [r1 x, r2 y,
           s1 1, s2 0,
           t1 0,  t2 1]
      (if (identically-zero? r2)
        (let [lead (if (number? r1) 1 (term r1 (index-highest r1)))]
          {:gcd (* (/ lead) r1)
           :bezout-1 (* (/ lead) s1)
           :bezout-2 (* (/ lead) t1)})
          (let [q (:quotient (divide r1 r2))]
            (recur r2 (- r1 (* r2 q))
                   s2 (- s1 (* s2 q))
                   t2 (- t1 (* t2 q))))))))



(defn factorial [n]
  (reduce clojure.core/* (take n (iterate inc 1))))

(defn stirling-2
  "Count the number of ways of partitioning of n into k nonempty
  groups."
  [n k]
  (cond (= 1 k) 1
        (> k n) 0
        (neg? n) 0
        (neg? k) 0
        :default
        (+ (stirling-2 (dec n) (dec k))
           (clojure.core/* k (stirling-2 (dec n) k)))))

(defn add-rationals
  "Return the sum of the rational functions P_i/Q_i in reduced form."
  [r1 r2 & rs]
  (let [rationals (concat [r1 r2] rs)
        numerators (map first rationals)
        denominators (map second rationals)
        gcd (reduce (comp :gcd gcd*) denominators)
        lcm (:quotient (divide (reduce * denominators) gcd))

        num (reduce ++ (for [[num den] rationals]
                         (:quotient (divide (* num lcm) den))))
        gcd (:gcd (gcd* num lcm))
        num (:quotient (divide num gcd))
        den (:quotient (divide lcm gcd))
        ]

    [num den]

    ))



(def monomial-generator-coefficients
  "Return an infinite list of the polynomials in 1/(1-z)
  which comprise the generating function for a_n = n^c."
  (map (fn [c]
         (vec (cons 0
                    (for [j (range (inc c))]
                      (clojure.core/* (stirling-2 (inc c) (inc j))
                         (factorial j)
                         (if (even? (- c j)) 1 -1))))))
       (range)))



(comment
  defn monomial-generator
  "Return the generating function for the sequence a_n = n^c"
  [c]
  (vec (for [j (range (inc c))]
         (* (stirling-2 (inc c) (inc j))
            (factorial j)
            (if (even? (- c j)) 1 -1)))))


(defn n-generator
  "Return the generating function for the sequence a_n = P(n)"
  [coeffs]
  (apply add-rationals
         (map vector 
              (map vector (reduce ++ (map * coeffs monomial-generator-coefficients)))
              (powers [1 -1]))))

(defrecord QuadraticRoot [a b c])
(derive first.core.QuadraticRoot ::root)

(defn discriminant [root]
  (let [a (:a root), b (:b root), c (:c root)]
  (- (* b b) (clojure.core/* 4 a c))))

(defn formula-root [root]
  (str "<mfrac><mrow>"
       "<mn>" (- (:b root)) "</mn>"
       "<mo>±</mo>"  
       "<msqrt>"
       "<mn>" (discriminant root) "</mn>"
       "</msqrt></mrow><mrow>"
       "<mn>" (* 2 (:a root)) "</mn>"
       "</mrow></mfrac>"))


(defmethod ++ [first.core.QuadraticRoot Number]
  [root k]
   (let [a (:a root) b (:b root) c (:c root)]
     (QuadraticRoot. a (- b (clojure.core/* 2 a k))
                     (+ c (clojure.core/* -1 b k)
                          (clojure.core/* a k k)))))




(defn linear-recurrence-series
  "Return a series representing the solution to the linear recurrence."
  [coeffs & init-values]
  (assert (= (index-highest coeffs)  (count init-values)) "Mismatch between recurrence size and number of initial values.")
  (let [solve-for-xn (vec (butlast (scale (/ -1 (term coeffs (index-highest coeffs))) coeffs)))]
    (lazy-cat init-values
              (map (fn [terms]
                     (reduce + (map * solve-for-xn terms)))
              (partition (count solve-for-xn) 1
                         (apply linear-recurrence-series coeffs init-values))))))

(comment :example
         (linear-recurrence-series [-1 -1 1] 0 1))


(defn mediant [p q]
  (/
   (+ (if (integer? p) p (numerator p))
      (if (integer? q) q (numerator q)))
   (+ (if (integer? p) 1 (denominator p))
      (if (integer? q) 1 (denominator q)))))

(defn farey-iterate [coll]
  (concat (interleave coll (map (partial apply mediant) (partition 2 1 coll))) [(last coll)]))

(def nice-fractions
  (distinct
   (lazy-seq (apply concat (iterate farey-iterate [0 1])))))

(defn nice-fraction [x]
  (cond (neg? x)
        (- (nice-fraction (- x)))
        (> x 1) (+ (int (quot x 1)) (nice-fraction (rem x 1)))
        :default
    (first (filter #(< (Math/abs (- x %)) 1e-3)
                   nice-fractions))))

(defn polynomial-regression [degree terms]
  (let [X (new Jama.Matrix
               (into-array
                (map
                 #(double-array
                   (take (count terms) ;;(inc degree);;(count terms)
                         (reductions * 1 (repeat (+ 0 %)))))
                 (range (count terms)))))
        y (new Jama.Matrix
               (into-array
                (map (comp double-array vector)
                     terms)))

        least-squares (->> (.transpose X)
                           (.times X)
                           (.inverse )
                           (.times (.transpose X)))
        coeffs (.times least-squares y)
        errors (map first (.getArray (.minus (.times X coeffs) y)))
        ]
    ;; (map (comp nice-fraction first)
    ;;  (.getArray (.times
    ;;             (.times (.inverse (.times (.transpose X) X))
    ;;                     (.transpose X)) y)))
    {:errors errors,
     :coeffs (identity ;; remove-leading-zeroes
              (vec (map (comp ;;nice-fraction
                         first) (.getArray coeffs))))
     }))

(defn guess-polynomial
  "Guess the polynomial P(n) given the first few terms starting at n=0"
  [terms]
  (let [fit (polynomial-regression (count terms) terms)]
    (when (every? #(< (Math/abs %) 1e-5) (:errors fit))
      (remove-leading-zeroes (vec (map nice-fraction (:coeffs fit))))
    )))

  ;; (first
  ;;  (for [k (range (inc (count terms)))]
  ;;    (let [fit (polynomial-regression k terms)]
  ;;      (when (and ;;(= (count (:coeffs fit)) (inc k))
  ;;                 (every? #(< (Math/abs %) 1e-8)
  ;;                         (:errors fit)))
  ;;        (:coeffs fit))))))

(defn guess-recurrence
  "Given the first few terms of a sequence, infer the complete sequence.
  Find the lowest-dimensional linear recurrence that satisfies all
  terms, or nil if the sequence is underspecified."
  [terms]
  (first
  (filter some?
  (for [k (range 1 (count terms))
        :while true]
    (let [b (new Jama.Matrix (into-array (map (comp double-array vector) (take k (drop k terms)) )))
          M (new Jama.Matrix (into-array (take k (map double-array (partition k 1 terms)))))
          Q (try (.getArray (.getQ (.qr (.transpose M)))) (catch Exception e nil))

          x (try (first (.getArray (.transpose (.solve M b)))) (catch Exception e nil))
          coeffs  (when x (conj (vec (map rationalize x)) -1))
          ]

      (when (and coeffs
                 (=  terms
                     (take (count terms) (apply linear-recurrence-series coeffs
                                                (take k terms)))))
        coeffs))))))




(defn … [& terms]
  "Given the first few terms of a sequence, infer the complete sequence.
  Find the lowest-dimensional linear recurrence that satisfies all
  terms, or nil if the sequence is underspecified.

  Alias for guess-recurrence."
  (if-let [coeffs (guess-recurrence terms)]
    (map int
         (apply linear-recurrence-series  coeffs (take (dec (count coeffs)) terms)))

    (if-let [coeffs (guess-polynomial terms)]
      (map #(int (reduce +
                    (map *
                         coeffs
                         (reductions * 1 (repeat %)))))
           (range)))))


(defn -main
  "entry point for app"
  [& args]
  ;;  (SwingUtilities/invokeLater show-polynomial)
  ;; (show-rational [1 -2] [1 -1 -1])

  ;; (println (factorial 5))

  (show-rational [1] [1 -1 -1])
  )

Date: 23/May/2022

Author: Dylan Holmes

Created: 2022-05-27 Fri 09:48

Emacs 27.2 (Org mode 8.3beta)

Validate