Merge pull request #102 from rikardn/log

Add log_FAC1

src/basic/,stubs_math/b9ea.log_FAC1.s

@@ -1,30 +0,0 @@
-;; #LAYOUT# STD *       #TAKE
-;; #LAYOUT# X16 BASIC_0 #TAKE-OFFSET 2000
-;; #LAYOUT# *   BASIC_0 #TAKE
-;; #LAYOUT# *   *       #IGNORE
-
-;
-; Math package - natural logarigthm
-;
-;
-; See also:
-; - [CM64] Computes Mapping the Commodore 64 - page 113
-; - https://www.c64-wiki.com/wiki/Floating_point_arithmetic
-; - https://codebase64.org/doku.php?id=base:kernal_floating_point_mathematics
-;
-
-; XXX provide implementation
-
-log_FAC1:
-
-	+STUB_IMPLEMENTATION
-
-
-poly_log:
-
-    !byte $03                          ; series length - 1
-
-    +PUT_CONST_POLY_LOG_1
-    +PUT_CONST_POLY_LOG_2
-    +PUT_CONST_POLY_LOG_3
-    +PUT_CONST_POLY_LOG_4

src/basic/math/b9ea.log_FAC1.s

@@ -0,0 +1,64 @@
+;; #LAYOUT# STD *       #TAKE
+;; #LAYOUT# X16 BASIC_0 #TAKE-OFFSET 2000
+;; #LAYOUT# *   BASIC_0 #TAKE
+;; #LAYOUT# *   *       #IGNORE
+
+;
+; Math package - natural logarigthm of FAC1
+;
+;
+; See also:
+; - [CM64] Computes Mapping the Commodore 64 - page 113
+; - https://www.c64-wiki.com/wiki/Floating_point_arithmetic
+; - https://codebase64.org/doku.php?id=base:kernal_floating_point_mathematics
+; - https://www.c64-wiki.com/wiki/LOG
+;
+; Calculate log of FAC1
+; Algorithm:
+; 1. Separate exponent and significand of FAC1.  N = exponent, XF = significand
+; 2. T = 1.0 - sqrt(2.0) / (XF + sqrt(0.5))
+; 3. Apply polynomial 0.43425594189 * T7 + 0.57658454124 * T5 + 0.96180075919 * T3 + 2.8853900731 * T - 0.5  giving log2(XF)
+; 4. N + log2(XF) giving log2(X)
+; 5. Multiply by LOG(2) giving LOG(X)
+
+
+
+log_FAC1:
+    lda FAC1_exponent           ; Store exponent for later
+    pha
+
+    ldx #$80                    ; Set exponent to 0 so that FAC1 becomes XF
+    stx FAC1_exponent
+
+    ldy #>const_INV_SQR_2       ; FAC1 <- FAC1 + sqrt(0.5)
+    lda #<const_INV_SQR_2
+    jsr add_MEM_FAC1
+
+    ldy #>const_SQR_2           ; FAC1 <- sqrt(2.0) / FAC1
+    lda #<const_SQR_2
+    jsr div_MEM_FAC1
+
+    ldy #>const_ONE             ; FAC1 <- 1.0 - FAC1
+    lda #<const_ONE
+    jsr sub_MEM_FAC1
+
+    ldy #>poly_log              ; FAC1 <- log_poly(FAC1)
+    lda #<poly_log
+    jsr poly1_FAC1
+    
+    ldy #>const_NEG_HALF        ; FAC1 <- FAC1 - 0.5
+    lda #<const_NEG_HALF
+    jsr add_MEM_FAC1
+    
+    jsr mov_FAC1_FAC2           ; FAC2 <- FAC1
+
+    pla                         ; FAC1 <- N
+    sec
+    sbc #$80
+    jsr convert_A_to_FAC1
+
+    jsr add_FAC2_FAC1           ; FAC1 <- FAC1 + FAC2
+
+    ldy #>const_LOG_2           ; FAC1 <- FAC1 * log(2.0)
+    lda #<const_LOG_2
+    jmp mul_MEM_FAC1

src/basic/math/poly_log.s

@@ -0,0 +1,15 @@
+;; #LAYOUT# STD *       #TAKE
+;; #LAYOUT# X16 BASIC_0 #TAKE-OFFSET 2000
+;; #LAYOUT# *  BASIC_0  #TAKE
+;; #LAYOUT# *   *       #IGNORE
+
+; Polynomial for the logarithm function
+
+poly_log:
+
+    !byte $03                          ; series length - 1
+
+    +PUT_CONST_POLY_LOG_1
+    +PUT_CONST_POLY_LOG_2
+    +PUT_CONST_POLY_LOG_3
+    +PUT_CONST_POLY_LOG_4