bugfix for f24amean and f24geomean

f24/f24amean.z80

@@ -2,7 +2,24 @@
 #define included_f24amean
 
 f24amean:
-;(AHL+BDE) ==> AHL
+;(AHL+CDE) ==> AHL
+
+;Make sure x+y won't overflow
+;save A
+  ld b,a
+  and $7F
+  jr z,f24amean_nooverflow
+  cp $7E
+  jr z,f24amean_overflow
+
+  ld a,c
+  and $7F
+  jr z,f24amean_nooverflow
+  cp $7E
+  jr z,f24amean_overflow
+
+f24amean_nooverflow:
+  ld a,b
   call f24add
 #ifdef included_f24div2
   jp f24div2
@@ -22,4 +39,21 @@ f24div2:
   dec a
   ret
 #endif
+
+f24amean_overflow:
+;need to decrement the exponents first
+  ld a,b
+  push de
+  push bc
+  call f24div2
+  pop bc
+  ex (sp),hl
+  push af
+  ld a,c
+  call f24div2
+  pop bc
+  pop de
+  ld c,b
+  jp f24add
+
 #endif

f24/f24geomean.z80

@@ -1,17 +1,49 @@
-;!! Note !!
-; This has numerical issues in that it is directly computing the the product,
-; then the square root. This can cause overflow when it shouldn't (note that as
-; long as both inputs are finite, the output should be strictly finite). There
-; is a slower and more difficult algorithm that can and probably should be made.
-
 #ifndef included_f24geomean
 #define included_f24geomean
 
 #include "f24mul.z80"
 #include "f24sqrt.z80"
 
 f24geomean:
-;sqrt(AHL*BDE) ==> AHL
+;sqrt(AHL*CDE) ==> AHL
+
+;return NaN if the product is negative
+  xor c
+  jp p,+_
+  ld a,$7F
+  ld h,a
+  ret
+_:
+
+;the product is positive, so we can just set the inputs to positive
+  xor c
+  and $7F
+  ret z   ;may as well exit if the input is 0
+  cp $7F
+  jr z,f24geomean_sub
+  ld b,a
+
+  ld a,c
+  and $7F
+  ret z
+  cp $7F
+  jr z,f24geomean_sub
+
+;now calculate the output exponent
+  add a,b
+  rra
+  push af
+
+  ld a,63
+  ld c,a
+  adc a,0
+  call f24geomean_sub
+  pop bc
+  add a,b
+  sub 63
+  ret
+
+f24geomean_sub:
   call f24mul
   jp f24sqrt
 #endif