Blog do projektu Open Source JavaHotel

środa, 28 sierpnia 2013

Clojure and polynomials

Introduction
I started reading through Clojure tutorial and found very interesting implementation of calculating polynomial and its first derivative. So I decided to stretch my math muscle and perform more polynomial arithmetic in Clojure basing on this article.
I'm assuming all the time that polynomial in Clojure is a list of coefficients from left to right (like in the written form). So list [2 8 9]  means 2x² + 8x + 9, [7 0 1 7] means: 7x³ + x + 7
Addition
Almost nothing more than (map + ...) function. The only problem is to expand shorter list (lower degree polynomial). Also improvement should be added to reduce the leading zeros from the result. Example:
(2x² + 8x + 9) + (-2x² + 2x + 4) = 0x² + 10x + 13

(defn addpol
  [pol1 pol2]
  ; degree of the sum polynomial 
  (let [nsize (max (count pol1) (count pol2) ) ]
  ; expand polynomial to the new degree by adding 0 coefficient at the frot
    (defn expandpol 
        [pol] 
            (concat ( repeat (max 0 (- nsize (count pol))) 0 ) pol )
    )
  )
  ; having both polynomials at the same degree just sum them
  (map + (expandpol pol1) (expandpol pol2))
)
Subtraction
Just use the previous function and negate the subtrahend.
(defn subtractpol
  [minuend subtrahend]
  ; negate subtrahend and run summing
  (addpol minuend ( map #(- 0 %) subtrahend))
)
Product
Multiply every coefficient of the first polynomial by every coefficient of the second polynomial. Recursion is used to iterate but this algorithm is a good candidate for loop->recur for better performance.

(defn mulpol
  [pol1 pol2]
  ( let
    ; level of the product
      [newlevel (- (+ (count pol1) (count pol2)) 1)]  
    ; multiply polynomial by digit and add zeros at the frond and end keeping it at the degree necessary
    (defn muldigit 
       [beg coeff end]
       (concat (repeat beg 0) (map #(* % coeff) pol2) (repeat end 0))
    )
    ; recursive as replacement for iteration
    (defn muladd
       ; 'beg' number of left zeros
       ; 'currentsum' list of coefficient multiplications peformed so far
       ; 'restpol1' list of coefficients not used yet
       ; 'end' number of right zeros
       [beg currentsum restpol1 end]
       (condp = restpol1
       ; end of recursion
       [] currentsum
       ; multiply and move to the next digit
       (
         map
            +
            currentsum
            (muladd 
               (+ beg 1)
               (muldigit beg (first restpol1) end)
               (rest restpol1)
               (- end 1)
            )
       )
       )
     )
  ; start of the recursion
  ( muladd 0 (repeat newlevel 0) pol1 (- newlevel (count pol2)))
  )
)
Division
Is more complicated because two results are expected: quotient and remainder. So to bring the result a map is used. Polynomial long division is implemented as easier. This function should be expanded to exclude division by zero (empty divisor) and reduce leading zeros from divisor. Example: (0 4 5)

; returns a map with two keys :quotient and :remainder
(defn divpol 
   [divident divisor]
   ; perform (lead(r)/lead(divisor)) * divisor
   (defn multerm
     [r resdiv]
       (concat (map #(* resdiv %) divisor) (repeat (- (count r) (count divisor)) 0))
   )
   ; recursive
   (defn div 
     [mapd]
   ( let [ q (get mapd :quotient)
           r (get mapd :remainder) 
           dv (/ (first r) (first divisor)) 
         ]
   ; end of recursion, pass down the result
   (if (or (< (count r) (count divisor)) (empty? r))
       mapd
   ; pull down the next digit from divident
       (div ( hash-map 
              :quotient (conj q dv)
              :remainder (rest (map - r (multerm r dv)))
            )
       )
   )
   )
   )   
   ; beginning of the recursion
   (div ( hash-map :quotient [], :remainder divident))
)   

Brak komentarzy:

Prześlij komentarz