'Subject: High precision natural logarithm or antilog with rational argument.
'Author : Sjoerd.J.Schaper
'Update : 07-05-2009
'Code   : QuickBasic, PDS, VBdos, FreeBasic
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DECLARE SUB RatExp (a, b, r)
'a/b:= exp(a/b), |a/b| <= .5
DECLARE SUB RatLog (a, b, r)
'a/b:= log(a/b), |a/b - 1| < .65
DECLARE SUB Expk (a, b, h, r)
'a/b:= exp(h), h >= 1, precision r
DECLARE SUB RSimpf (a, b)
'simplify ratio a / b

CONST Verb = 1, e = 100 '             set verbosity & decimal precision
CONST p = 0, q = 1, s = 2, u = 3
CONST p0 = 4, q0 = 5, pow = 6, k = 7
r = LargeInit(k, "log_exp")
DIM ms(1) AS STRING
CLS : tim# = TIMER: lin = 8

begin:
PRINT "   +------------------------------+"
PRINT "   |  input                       |"
PRINT "   |  numerator and denominator,  |"
PRINT "   |  0 to exit.                  |"
PRINT "   |                              |"
PRINT "   +------------------------------+"
LOCATE 5, 7: INPUT ; "n"; Ni$
IF INSTR(Ni$, "'") THEN
   IF INSTR(UCASE$(Ni$), "L") THEN ALOG = 0
   IF INSTR(UCASE$(Ni$), "E") THEN ALOG = -1
   GOTO begin
END IF
IF VAL(Ni$) = 0 THEN GOTO eind
INPUT ; " d"; g$
IF VAL(g$) = 0 THEN GOTO begin

LOCATE lin, 1: PRINT "working..."
Readst p, Ni$: Readst q, g$
Printn p, "  n = ", "", 0
Printn q, "/", "", 1

sg = Absf(p) XOR Absf(q)
r = 31 + INT(e * 3.322) '             scaling length
IF Verb THEN
   Letf u, 1: Lsft u, r
   Printn u, "unit = ", "", 2
END IF

IF ALOG THEN '                        Exp(p / q)
   Dup p, t0
   Divd t0, q, k: Lsft t0, 1
   IF Cmp(q, t0) < 1 THEN Inc k, 1 '  k:= round(p / q)
   IF Getl(k) > 2 THEN
      Prints "overflow", 2
      lin = CSRLIN: GOTO begin
   END IF

   kn = NOT Isf(k, 0)
   IF kn THEN
      Mult q, k, t0: Subt p, t0 '     subtract nearest integer
   END IF
   RSimpf p, q
   IF Verb THEN
      Printn k, "   k = ", ", ", 0
      Printn p, "p/q = ", "", 0
      Printn q, "/", "", 1
   END IF
   RatExp p, q, r

   IF kn THEN
      Expk p0, q0, k, r
      Lmul p, p0: Lmul q, q0 '        exp(n) = exp(n - k) * exp(k)
      Rsft p, r: Rsft q, r
   END IF
   IF sg THEN Swp p, q
   ms(0) = "Rex = ": ms(1) = "exp = "

ELSE '                                Log(p / q)
   IF sg OR Isf(p, 0) THEN
      Prints "illegal argument", 2
      lin = CSRLIN: GOTO begin
   END IF
   sw = Cmp(q, p) = 1
   IF sw THEN Swp p, q '              q > p
   c = INT(Logf(p) - Logf(q) + .5) '  k = round(Log(p / q))
   Letf k, c

   kn = c > 0
   IF kn THEN
      Expk p0, q0, k, r
      Lmul p, q0: Lmul q, p0 '        scale with factor exp(k)
   END IF
   RSimpf p, q
   IF Verb THEN
      Printn k, "   k = ", ", ", 0
      Printn p, "p/q = ", "", 0
      Printn q, "/", "", 1
   END IF
   RatLog p, q, r

   IF kn THEN
      Mult q, k, t0: Add p, t0 '      log(n) = log(n / exp(k)) + k
   END IF
   IF sw THEN Chs p
   ms(0) = "Rln = ": ms(1) = "log = "

END IF

RSimpf p, q
IF Verb THEN
   Printn p, ms(0), "", 0
   Printn q, "/", "", 1
END IF
Printr p, q, ms(1), e '               decimal expansion
Prints "", 1
lin = CSRLIN: GOTO begin

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

SUB Expk (a, b, h, r)
   Letf a, 1: Letf b, 4
   RatExp a, b, r '                   a / b = exp(1 / 4)
   Dup a, s: Dup b, u
   Lsft h, 2
   '
   i = Getl(h) - 1
   FOR j = i TO 0 STEP -1 '          L=>R binary exponentiation
      d = Getw(h, j): d0 = MB \ 2 '   unsigned bitvector 4h
      IF j = i THEN d0 = Hp2(h) \ 2
      WHILE d0 > 0
         Squ s, t0: Squ u, t1
         Swp t0, s: Swp t1, u '       square s/u
         Rsft s, r: Rsft u, r '       truncate
         IF d AND d0 THEN '           read bit
            Lmul s, a: Lmul u, b '    s/u times base a/b
            Rsft s, r: Rsft u, r
         END IF
         d0 = d0 \ 2
      WEND
   NEXT j
   '
   Rsft h, 2
   Swp a, s: Swp b, u
END SUB

SUB RatExp (a, b, r)
   Letf u, 1: Lsft u, r '             unit 2^ r
   Dup u, pow: Dup pow, s '           pow[0] = sum[0] = unit
   Letf t1, 0
   IF Verb THEN Prints "s(0) = unit", 1
   '
   sw = Isf(a, 1)
   FOR t = 1 TO 32766 '               compute power series
      Add t1, b
      IF sw THEN
         Swp pow, t0
      ELSE
         Mult pow, a, t0
      END IF
      Divd t0, t1, pow '              pow[t]:= pow[t-1] * a / (b * t)
      '
      IF Isf(pow, 0) THEN EXIT FOR '  precision exhausted
      Add s, pow '                    update sum
   NEXT t
   IF Verb THEN
      g$ = "s(" + LTRIM$(STR$(t - 1))
      Printn s, g$ + ") = ", "/unit", 1
   END IF
   '
   Swp a, s: Swp b, u
END SUB

SUB RatLog (a, b, r)
   Letf u, 1: Lsft u, r '             unit 2^ r
   Dup a, t0: Subt a, b '             argument |(a - b)
   Add b, t0 '                                         / (a + b)| < .25
   '
   Squ a, pow: Swp pow, a '           a'= a^2
   Lsft pow, r: Ldiv pow, b '         pow[0] = sum[0] = unit * a / b
   Squ b, t0: Swp t0, b '             b'= b^2
   od = 1: Dup pow, s
   IF Verb THEN
      Printn s, "s(0) = ", "*2/unit", 1
   END IF
   '
   FOR t = 1 TO 16383 '               compute power series
      Mult pow, a, t0
      Divd t0, b, pow '               pow[t]:= pow[t-1] * a^2 / b^2
      Dup pow, t0: od = od + 2
      rm = Divint(t0, od) '           pow[t] / (2t + 1)
      '
      IF Isf(t0, 0) THEN EXIT FOR '   precision exhausted
      Add s, t0 '                     update sum
   NEXT t
   Lsft s, 1
   '
   IF Verb THEN
      g$ = "s(" + LTRIM$(STR$(t - 1))
      Printn s, g$ + ") = ", "/unit", 1
   END IF
   '
   Swp a, s: Swp b, u
END SUB

SUB RSimpf (a, b)
   Gcd a, b, t3
   IF NOT Isf(t3, 1) THEN
      Ldiv a, t3: Ldiv b, t3 '        remove common factor
   END IF
END SUB

