'Subject: Karatsuba multiplication. Play with a recursive subroutine
'         and find that applying K-multiplication doesn't really pay
'         for factors with less than roughly 250 digits.
'Author : Sjoerd.J.Schaper
'Date   : 06-08-2005
'Code   : Quickbasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DECLARE SUB KMult (p AS INTEGER, q AS INTEGER, r AS INTEGER)
'Save p in r, and set p to product p * q
CONST a = 0, b = 1, t = 2, ix = 32 '   stack size 6 * 5
CONST NC = 30 '                        cutoff number length
CONST nmx = 500 '                      iterations in test
n = LargeInit(ix, "kratsuba")
DIM SHARED tl, tmx, tot '              counters
DIM ti(1) AS SINGLE
CLS : tim# = TIMER: lin = 8

begin:
tl = 3: tmx = 0: tot = 0
LOCATE 1, 1
PRINT "   +--------------------------------------+"
PRINT "   |  input numbers a and b,              |"
PRINT "   |  attach 't' for a performance test,  |"
PRINT "   |  0 to exit.                          |"
PRINT "   |                                      |"
PRINT "   +--------------------------------------+"
LOCATE 5, 7: INPUT ; "a"; ai$
IF INSTR(ai$, "'") THEN GOTO begin
IF VAL(ai$) < 1 THEN GOTO eind
INPUT ; " b"; bi$
IF VAL(bi$) < 1 THEN GOTO begin
LOCATE lin, 1

IF INSTR(UCASE$(ai$ + bi$), "T") THEN
   IF LEN(ai$) + LEN(bi$) < 10 THEN
      PRINT " working..."
      FOR m = 0 TO 1
         RANDOMIZE 13
         tm! = TIMER
         FOR n = 1 TO nmx
            Rndf a, VAL(ai$) * 3.322 ' input is interpreted as
            Rndf b, VAL(bi$) * 3.322 ' length of random factors
            IF m THEN
               tl = 3: tmx = 0: tot = 0
               KMult a, b, t
            ELSE
               Mult a, b, t: Swp a, t
            END IF
         NEXT n
         ti(m) = TIMER - tm!
      NEXT m
      IF tl THEN
         LOCATE lin, 1: tmx = (tmx - 3) \ 5
         Prints " multiplied " + STR$(nmx) + " times ", 0
         Prints STR$(VAL(ai$)) + " x " + STR$(VAL(bi$)) + " digits", 1
         Prints " schoolbook: " + STR$(ti(0)) + " s", 1
         Prints "  Karatsuba: " + STR$(ti(1)) + " s", 1
         g$ = " calls " + STR$(tot) + "  depth " + STR$(tmx)
         IF ti(1) > 0 THEN g$ = g$ + "   ratio " + STR$(ti(0) / ti(1))
         Prints g$, 2
      END IF
   ELSE
      Prints " a and/or b too big for testing", 2
   END IF
   lin = CSRLIN
   GOTO begin
END IF

Readst a, ai$ '                        factors a and b
Readst b, bi$
Printn a, " a = ", "", 2
Printn b, " b = ", "", 2
KMult a, b, t
IF tl THEN Printn a, " P = ", "", 2

Prints "", 1
lin = CSRLIN
GOTO begin

eind:
LOCATE lin, 1
g$ = STR$(CSNG(TIMER - tim#))
Prints "Timer:" + g$ + " s", 1
CALL Term
SYSTEM

SUB KMult (al AS INTEGER, bt AS INTEGER, at AS INTEGER)
   la = Getl(al): tb = Getl(bt) '      number lengths
   lx = la: lm = tb: mi = bt
   IF lx < lm THEN SWAP lx, lm: mi = al
   l1 = lx \ 2 '                       splitting length
   IF l1 > lm THEN l1 = lm
   lx = l1 * 2
   IF lx < NC THEN '                   base case
      IF Isf(mi, 0) THEN
         Letf at, 0
      ELSE
         Mult al, bt, at
      END IF
      Swp al, at: EXIT SUB
   END IF

   bl = tl: abL = tl + 1: ar = tl + 2 ' get five row stackframe
   br = tl + 3: abR = tl + 4: tl = tl + 5
   IF abR > ix THEN
      Prints " out of stack space", 1
      tl = 0: EXIT SUB
   END IF
   IF tl > tmx THEN tmx = tl
   tot = tot + 1

   Setl al, l1: Setl bt, l1
   Dup al, ar: Dup bt, br '            copy lower halves
   Setl al, la: Setl bt, tb
   Lftj ar, l1 - 1: Lftj br, l1 - 1

   Dup al, at: Dup bt, bl
   Rsft al, -l1: Rsft bl, -l1 '        shift-right upper halves

   KMult al, bl, abL '                 enter recursion
   KMult ar, br, abR
   Swp al, abL: Swp ar, abR
   Add al, ar: Add bl, br
   KMult al, bl, ar
   IF tl = 0 THEN EXIT SUB
   Subt al, abL: Subt al, abR '        compute the middle term

   Lsft abL, -lx '                     shift-left upper product,
   Lsft al, -l1 '                      middle term
   Add al, abR: Add al, abL '          recombination

   tl = tl - 5
END SUB

