'===========================================================================
' Subject: COMPUTE NATURAL LOG TABLE         Date: 01-03-97 (09:16)       
' Author:  Hans Lunsing                      Code: QB, QBasic, PDS        
' Origin:  FidoNet QUIK_BAS Echo           Packet: ALGOR.ABC
'===========================================================================
' > Well, I declared each variable as DOUBLE separately, and now I
' > get correct values for ln 2, 2.5, 3, 3.5 and 4. All the rest are
' > off. Try it and see. Any suggestions? How do I find out what
' > algorithm QuickBASIC is using, short of disassembling their
' > libraries?

'This is what I got:

'2.00  .6931471805599453  .6931471805599453
'2.10  .7419373447293773  .7419373447293773
'2.20  .7884573603642703  .7884573603642703
'2.30  .8329091229351041  .8329091229351041
'2.40  .8754687373539001  .8754687373539001
'2.50  .9162907318741552  .9162907318741552
'2.60  .9555114450274366  .9555114450274366
'2.70  .9932517730102836  .9932517730102837
'2.80  1.029619417181159  1.029619417181159
'2.90  1.064710736992429  1.064710736992429
'3.00  1.09861228866811  1.09861228866811
'3.10  1.131402111491101  1.131402111491101
'3.20  1.163150809805681  1.163150809805681
'3.30  1.193922468472435  1.193922468472435
'3.40  1.223775431622116  1.223775431622116
'3.50  1.252762968495368  1.252762968495368
'3.60  1.280933845462065  1.280933845462065
'3.70  1.308332819650179  1.308332819650179
'3.80  1.33500106673234  1.33500106673234
'3.90  1.360976553135601  1.360976553135601
'4.00  1.386294361119891  1.386294361119891

'The approximation function does a really good job, isn't it? If you get
'other results you must have done something wrong. This is the source I
'used to ptoduce this table:

DEFINT A-Z
DECLARE FUNCTION Dln# (x#)
DECLARE FUNCTION Expo% (x#)

' Approx tries out an approximation routine for ln (x).

CONST TRUE = -1
CONST FALSE = 0

DIM SHARED Aerr                         ' Code for approximation error.

' Main code:

        CLS
        FOR x# = 2# TO 4.1# STEP .1#
          y# = Dln#(x#)
          z# = LOG(x#)
          PRINT USING "#.## "; x#;
          PRINT y#; z#
        NEXT x#
END

FUNCTION Dln# (x#)
' Dln finds the natural log of the passed parameter:

CONST a0 = -64.12494342374558# - 1.147E-15
CONST a1 = 16.38394356302153# + 4.222E-15
CONST a2 = -.7895611288749125# - 7.267E-17

CONST b0 = -769.4993210849487# - 9.777E-14
CONST b1 = 312.0322209192453# + 2.844E-14
CONST b2 = -35.66797773903464# - 6.171E-15

CONST c0 = .7071067811865475# + 2.44E-17
CONST c1 = .693359375#
CONST c2 = -2.121944400546906D-04

DIM F AS DOUBLE                         ' Fractional part.
DIM aw AS DOUBLE, bw AS DOUBLE          ' Temporary variables.
DIM rz AS DOUBLE, rz2 AS DOUBLE         ' Temporary variables.
DIM w AS DOUBLE, z AS DOUBLE            ' Temporary variables.
DIM Znum AS DOUBLE, Zden AS DOUBLE      ' Z numerator and denominator.
DIM N AS INTEGER                        ' Exponent of 2 in passed parameter.

' Main code:

        IF x# <= 0# THEN
                Aerr = TRUE             ' Signal an error.
                Dln# = 0#               ' Pass back zero.
                EXIT FUNCTION           ' Quit early.
        END IF

        N = Expo(x#) + 1                ' Use DEC-style exponents.
        Power# = 1#                     ' This will hold 2 ^ N.

        FOR I% = 1 TO N
                Power# = Power# * 2#    ' "Integer exponentiation."
        NEXT I%

        F = x# / Power#                 ' Find DEC-style mantissa.

        IF F > c0 THEN
                Znum = F - 1#
                Zden = .5# * F + .5#
        ELSE
                N = N - 1
                Znum = F - .5#
                Zden = .5# * Znum + .5#
        END IF

        z = Znum / Zden                 ' By definition.
        w = z * z
        aw = a0 + w * (a1 + w * a2)             ' Numerator polynomial in w.
        bw = b0 + w * (b1 + w * (b2 + w))       ' Denominator polynomial in w.

        rz2 = w * (aw / bw)
        rz = z * (1# + rz2)
        Dln# = (c2 * N + rz) + c1 * N
END FUNCTION ' Dln#

FUNCTION Expo% (x#)
' Expo returns the exponent of 2 from a double-precision number.
        DEF SEG = VARSEG(x#)            ' Find data area.
        Start% = VARPTR(x#)             ' Find x's place therein.

        a& = PEEK(Start% + 7)           ' Most significant byte.
        b& = PEEK(Start% + 6)           ' Next most significant.

        Word& = 256 * a& + b&
        IF Word& > 32767 THEN Word& = Word& - 32768     ' Erase sign bit.

        Word& = Word& \ 16              ' Erase mantissa bits.
        Expo = Word& - 1023             ' Remove bias.
END FUNCTION ' Expo