Karatsuba algorithm

[cameyo]

From: https://en.wikipedia.org/wiki/Karatsuba_algorithm">//https://en.wikipedia.org/wiki/Karatsuba_algorithm

Code Select Expand

(define (potenza n m)
  (let (pot 1L) (dotimes (x m) (setq pot (* pot n))))
)

(potenza 3 6)
;-> 729L

Iterative:

(define (karatsuba num1 num2)
  (local (m m2 high1 low1 high2 low2 z0 z1 z2)
    (cond ((or (< num1 10) (< num2 10)) (* num1 num2))
          (true
            (setq m (max (length (string num1)) (length (string num2))))
            (setq m2 (/ m 2))
            (setq n1$ (string num1))
            (setq n2$ (string num2))
            (setq high1 (int (slice n1$ 0 (- (length n1$) m2))))
            (setq low1  (int (slice n1$ (- (length n1$) m2) m2)))
            (setq high2 (int (slice n2$ 0 (- (length n2$) m2))))
            (setq low2  (int (slice n2$ (- (length n2$) m2) m2)))
            ;(println high1 { } low1)
            ;(println high2 { } low2)
            (setq z0 (karatsuba low1 low2))
            (setq z1 (karatsuba (+ low1 high1) (+ low2 high2)))
            (setq z2 (karatsuba high1 high2))
            (+ (* z2 (potenza 10 (* m2 2))) (* (- z1 z2 z0) (potenza 10 m2)) z0)
          )
    )
  );local
)

(karatsuba 12 12)
;-> 144
(karatsuba 13 17)
;-> 221
(karatsuba 120 11)
;-> 1320
(karatsuba 12345 6789)
;-> 83810205
(mul 12345 6789)
;-> 83810205
(time (karatsuba 12345 6789) 10000)
;-> 359.359

Recursive:

(define (karatsuba x y)
    (karatsuba1 x y 256)  ; in generale, opportuna potenza di 2 p (x , y < 2p)
)

(define (karatsuba1 x y p)  ; x, y, p: interi non negativi, p potenza di 2
    (if (= p 1)
        (* x y)
        (let ((x1 (/ x p)) (x0 (% x p))
              (y1 (/ y p)) (y0 (% y p))
              (q (/ p 2)))
          (let ((z2 (karatsuba1 x1 y1 q))
                (z0 (karatsuba1 x0 y0 q)))
            (let ((z1 (- (karatsuba1 (+ x1 x0) (+ y1 y0) q) (+ z2 z0))))
              (+ (* z2 p p) (* z1 p) z0)
            )
          )
        )
     )
)

(karatsuba 12 12)
;-> 144
(karatsuba 12345 6789)
;-> 83810205
(time (karatsuba 12345 6789) 10000)
;-> 33347.174

big integer "L":
 
(define (karatsuba num1 num2)
  (local (len1 len2 m m2 high1 low1 high2 low2 z0 z1 z2)
    (cond ((or (< num1 10) (< num2 10)) (* num1 num2))
          (true
            (setq len1 (length (string num1)))
            (if (= (last (string num1)) "L") (-- len1))
            (setq len2 (length (string num2)))
            (if (= (last (string num1)) "L") (-- len2))            
            (setq m (max len1 len2))
            (setq m2 (/ m 2))
            (setq n1$ (string num1))
            (if (= (last n1$) "L") (setq n1$ (chop n1$)))
            (setq n2$ (string num2))
            (if (= (last n2$) "L") (setq n2$ (chop n2$)))
            (setq high1 (bigint (slice n1$ 0 (- (length n1$) m2))))
            (setq low1  (bigint (slice n1$ (- (length n1$) m2) m2)))
            (setq high2 (bigint (slice n2$ 0 (- (length n2$) m2))))
            (setq low2  (bigint (slice n2$ (- (length n2$) m2) m2)))
            ;(println high1 { } low1)
            ;(println high2 { } low2)
            (setq z0 (karatsuba low1 low2))
            (setq z1 (karatsuba (+ low1 high1) (+ low2 high2)))
            (setq z2 (karatsuba high1 high2))
            (+ (* z2 (potenza 10 (* m2 2))) (* (- z1 z2 z0) (potenza 10 m2)) z0)
          )
    )
  );local
)

(karatsuba 12345 6789)
;-> 83810205
(karatsuba 9223372036854775807 9223372036854775807)
;-> 85070591730234615847396907784232501249L
(* 9223372036854775807L 9223372036854775807L)
;-> 85070591730234615847396907784232501249L
(time (karatsuba 12345 6789) 10000)
;-> 687.468