;; -*-theme-d-*-

;; Copyright (C) 2018, 2021 Tommi Höynälänmaa
;; Distributed under GNU Lesser General Public License version 3,
;; see file doc/LGPL-3.


(define-body (standard-library extra-math)

  ;; We trust libthemedsupport so we can use unchecked-prim-proc here.

  (add-method fmod
    (unchecked-prim-proc fmod (<real> <real>) <real> pure))

  (add-method r-remainder
    (unchecked-prim-proc r-remainder (<real> <real>) <real> pure))

  (add-method fma
    (unchecked-prim-proc fma (<real> <real> <real>) <real> pure))

  (add-method fast-fma?
    (unchecked-prim-proc fast-fma? () <boolean> pure))

  (add-method fmin
    (unchecked-prim-proc fmin (<real> <real>) <real> pure))

  (add-method fmax
    (unchecked-prim-proc fmax (<real> <real>) <real> pure))

  (add-method fdim
    (unchecked-prim-proc fdim (<real> <real>) <real> pure))

  (add-method r-exp2
    (unchecked-prim-proc r-exp2 (<real>) <real> pure))

  (add-method r-expm1
    (unchecked-prim-proc r-expm1 (<real>) <real> pure))

  (add-method r-log2
    (unchecked-prim-proc r-log2 (<real>) <real> pure))

  (add-method r-log1p
    (unchecked-prim-proc r-log1p (<real>) <real> pure))

  (add-method logb
    (unchecked-prim-proc logb (<real>) <real> pure))

  (add-method ilogb
    (unchecked-prim-proc ilogb (<real>) <integer> pure))

  (add-method ilogb0
    (unchecked-prim-proc ilogb0 () <integer> pure))

  (add-method ilogbnan
    (unchecked-prim-proc ilogbnan () <integer> pure))

  (add-method r-cbrt
    (unchecked-prim-proc r-cbrt (<real>) <real> pure))

  (add-method r-hypot
    (unchecked-prim-proc r-hypot (<real> <real>) <real> pure))

  (add-method r-erf
    (unchecked-prim-proc r-erf (<real>) <real> pure))

  (add-method r-erfc
    (unchecked-prim-proc r-erfc (<real>) <real> pure))

  (add-method r-lgamma
    (unchecked-prim-proc r-lgamma (<real>) <real> pure))

  (add-method r-tgamma
    (unchecked-prim-proc r-tgamma (<real>) <real> pure))

  (add-method r-nearbyint
    (unchecked-prim-proc r-nearbyint (<real>) <real> pure))

  (add-method rint
    (unchecked-prim-proc rint (<real>) <real> pure))

  (add-method frexp
    (unchecked-prim-proc frexp (<real>) (:pair <real> <integer>) pure))

  (add-method ldexp
    (unchecked-prim-proc ldexp (<real> <integer>) <real> pure))

  (add-method modf
    (unchecked-prim-proc modf (<real>) (:pair <real> <real>) pure))

  (add-method r-nextafter
    (unchecked-prim-proc r-nextafter (<real> <real>) <real> pure))

  (add-method r-copysign
    (unchecked-prim-proc r-copysign (<real> <real>) <real> pure))

  (add-method fpclassify
    (unchecked-prim-proc fpclassify (<real>) <integer> pure))

  (add-method fpclassify-nan
    (unchecked-prim-proc fpclassify-nan () <integer> pure))

  (add-method fpclassify-infinite
    (unchecked-prim-proc fpclassify-infinite () <integer> pure))

  (add-method fpclassify-zero
    (unchecked-prim-proc fpclassify-zero () <integer> pure))

  (add-method fpclassify-subnormal
    (unchecked-prim-proc fpclassify-subnormal () <integer> pure))

  (add-method fpclassify-normal
    (unchecked-prim-proc fpclassify-normal () <integer> pure))

  (add-method r-isnormal?
    (unchecked-prim-proc r-isnormal? (<real>) <boolean> pure))

  (add-method r-signbit
    (unchecked-prim-proc r-signbit (<real>) <integer> pure))

  (define-simple-method r-log2-neg (((r <real>)) <complex> pure)
    (complex (r-log2 (r-abs r)) gl-r-pi/ln2))

  (define-simple-method i-cbrt (((i <integer>)) (:union <real> <integer>) pure)
    (let ((r-result (r-cbrt (integer->real i))))
      (if (integer-float? r-result) (real->integer r-result) r-result)))

  (define-simple-method rat-cbrt (((rat <rational>))
				  (:union <real> <rational> <integer>) pure)
    (let ((r-numer (r-cbrt (integer->real (numerator rat))))
	  (r-denom (r-cbrt (integer->real (denominator rat)))))
      (if (and (integer-float? r-numer) (integer-float? r-denom))
	  (let ((i-numer (real->integer r-numer))
		(i-denom (real->integer r-denom)))
	    (if (equal? (remainder i-numer i-denom) 0)
		(quotient i-numer i-denom)
		(rational i-numer i-denom)))
	  (/ r-numer r-denom))))

  (define-simple-virtual-method exp2 (((nr <number>)) <number> pure)
    (raise-simple 'exp2:dispatch-error))

  (define-simple-virtual-method exp2 (((nr <real-number>)) <real-number> pure)
    (raise-simple 'exp2:dispatch-error))

  (include-virtual-methods exp2 c-exp2)

  (include-virtual-methods exp2 r-exp2)

  (define-simple-virtual-method exp2 (((rat <rational>)) <real-number> pure)
    (let* ((rat-s (simplify-rational rat))
	   (i-numer (numerator rat-s))
	   (i-denom (denominator rat-s)))
      (if (equal? i-denom 1)
	  (i-expt 2 i-numer)
	  (r-exp2 (rational->real rat)))))

  (define-simple-virtual-method exp2 (((i <integer>)) <rational-number> pure)
    (i-expt 2 i))

  (define-simple-virtual-method expm1 (((nr <number>)) <number> pure)
    (raise-simple 'expm1:dispatch-error))

  (define-simple-virtual-method expm1 (((nr <real-number>)) <real-number> pure)
    (raise-simple 'expm1:dispatch-error))

  (define-simple-virtual-method expm1 (((cx <complex>)) <complex> pure)
    (- (c-exp cx) 1.0))

  (include-virtual-methods expm1 r-expm1)

  (define-simple-virtual-method expm1 (((rat <rational>))
				       (:union <real> <integer>)
				       pure)
    (if (rat-zero? rat)
	0
	(r-expm1 (rational->real rat))))

  (define-simple-virtual-method expm1 (((i <integer>)) (:union <real> <integer>)
				       pure)
    (if (equal? i 0)
	0
	(r-expm1 (integer->real i))))

  (define-simple-virtual-method log2 (((nr <number>)) <number> pure)
    (raise-simple 'log2:dispatch-error))

  (define-simple-virtual-method log2 (((cx <complex>)) <complex> pure)
    (let ((r-abs-value (c-abs cx))
	  (r-phase (c-angle cx)))
      (create <complex> (r-log2 r-abs-value) (/ r-phase gl-r-ln2))))

  (define-simple-virtual-method log2 (((r <real>)) (:union <real> <complex>)
				      pure)
    (if (>= r 0.0)
	(r-log2 r)
	(r-log2-neg r)))

  (define-simple-virtual-method log2 (((rat <rational>))
				      (:union <integer> <real> <complex>) pure)
    (cond
     ((rat-zero? rat) (raise-numerical-overflow 'log2))
     ((rat-one? rat) 0)
     ((> rat 0)
      (match-type (rat-log2-exact rat)
	((<null>) (r-log2 (rational->real rat)))
	((i2 <integer>) i2)))
     (else (r-log2-neg (rational->real rat)))))

  (define-simple-virtual-method log2 (((i <integer>))
				      (:union <integer> <real> <complex>) pure)
    (cond
     ((= i 0) (raise-numerical-overflow 'log2))
     ((= i 1) 0)
     ((> i 0)
      (match-type (i-log2-exact i)
	((<null>) (r-log2 (integer->real i)))
	((i2 <integer>) i2)))
     (else (r-log2-neg (integer->real i)))))

  (define-simple-virtual-method log1p (((nr <number>)) <number> pure)
    (raise-simple 'log1p:dispatch-error))

  (define-simple-virtual-method log1p (((cx <complex>)) <complex> pure)
    (c-log (+ 1 cx)))

  (define-simple-virtual-method log1p (((r <real>))
				       (:union <complex> <real>) pure)
    (cond
     ((> r -1.0) (r-log1p r))
     ((equal? r -1.0) (neg-inf))
     ((< r -1.0) (r-log-neg (+ 1.0 r)))
     (else (nan))))

  (define-simple-virtual-method log1p (((rat <rational>))
				       (:union <complex> <real> <integer>) pure)
    (cond
     ((> rat -1) (r-log1p (rational->real rat)))
     ((equal? rat -1) (raise-numerical-overflow 'log1p))
     ((< rat -1) (r-log-neg (rational->real (+ 1 rat))))
     (else
      ;; We should never arrive here.
      (raise-simple 'log1p:internal-error))))

  (define-simple-virtual-method log1p (((i <integer>))
				       (:union <complex> <real> <integer>) pure)
    (cond
     ((> i -1) (r-log1p (integer->real i)))
     ((equal? i -1) (raise-numerical-overflow 'log1p))
     ((< i -1) (r-log-neg (integer->real (+ 1 i))))
     (else
      ;; We should never arrive here.
      (raise-simple 'log1p:internal-error))))

  (define-simple-virtual-method cbrt (((nr <number>)) <number> pure)
    (raise-simple 'cbrt:dispatch-error))

  (define-simple-virtual-method cbrt (((r <real-number>)) <real-number> pure)
    (raise-simple 'cbrt:dispatch-error))

  (define-simple-virtual-method cbrt (((cx <complex>)) <complex> pure)
    (let ((r-abs (c-abs cx))
	  (r-angle (c-angle cx)))
      (make-polar (r-cbrt r-abs) (/ r-angle 3))))

  (include-virtual-methods cbrt r-cbrt)

  (include-virtual-methods cbrt rat-cbrt)

  (include-virtual-methods cbrt i-cbrt)

  (define-simple-virtual-method hypot (((r1 <real-number>) (r2 <real-number>))
				       <real-number> pure)
    (raise-simple 'hypot:dispatch-error))

  (include-virtual-methods hypot r-hypot)

  (define-simple-virtual-method hypot (((r1 <real>) (rat2 <rational>)) <real>
				       pure)
    (r-hypot r1 (rational->real rat2)))

  (define-simple-virtual-method hypot (((r1 <real>) (i2 <integer>)) <real> pure)
    (r-hypot r1 (integer->real i2)))

  (define-simple-virtual-method hypot (((rat1 <rational>) (r2 <real>)) <real>
				       pure)
    (r-hypot (rational->real rat1) r2))

  (define-simple-virtual-method hypot (((rat1 <rational>) (rat2 <rational>))
				       <real-number> pure)
    (rat-sqrt (+ (square rat1) (square rat2))))

  (define-simple-virtual-method hypot (((rat1 <rational>) (i2 <integer>))
				       <real-number> pure)
    (rat-sqrt (+ (square rat1) (square i2))))

  (define-simple-virtual-method hypot (((i1 <integer>) (r2 <real>)) <real> pure)
    (r-hypot (integer->real i1) r2))

  (define-simple-virtual-method hypot (((i1 <integer>) (rat2 <rational>))
				       <real-number> pure)
    (rat-sqrt (+ (square i1) (square rat2))))

  (define-simple-virtual-method hypot (((i1 <integer>) (i2 <integer>))
				       (:union <real> <integer>) pure)
    (i-sqrt (+ (square i1) (square i2))))

  (define-simple-virtual-method erf (((r <real-number>))
				     (:union <real> <integer>)
				     pure)
    (raise-simple 'erf:dispatch-error))

  (include-virtual-methods erf r-erf)

  (define-simple-virtual-method erf (((rat <rational>))
				     (:union <real> <integer>)
				     pure)
    (if (rat-zero? rat) 0 (r-erf (rational->real rat))))

  (define-simple-virtual-method erf (((i <integer>)) (:union <real> <integer>)
				     pure)
    (if (equal? i 0) 0 (r-erf (integer->real i))))

  (define-simple-virtual-method erfc (((r <real-number>))
				      (:union <real> <integer>)
				      pure)
    (raise-simple 'erfc:dispatch-error))

  (include-virtual-methods erfc r-erfc)

  (define-simple-virtual-method erfc (((rat <rational>))
				      (:union <real> <integer>)
				      pure)
    (if (rat-zero? rat) 0 (r-erfc (rational->real rat))))

  (define-simple-virtual-method erfc (((i <integer>)) (:union <real> <integer>)
				      pure)
    (if (equal? i 0) 1 (r-erfc (integer->real i))))

  (define-simple-virtual-method lgamma (((r <real-number>))
					(:union <real> <integer>)
					pure)
    (raise-simple 'lgamma:dispatch-error))

  (include-virtual-methods lgamma r-lgamma)

  (define-simple-virtual-method lgamma (((rat <rational>))
					(:union <real> <integer>)
					pure)
    (let* ((rat-s (simplify-rational rat))
	   (i-numer (numerator rat-s))
	   (i-denom (denominator rat-s)))
      (if (and (equal? i-denom 1) (or (equal? i-numer 1) (equal? i-numer 2)))
	  0
	  (r-lgamma (rational->real rat-s)))))

  (define-simple-virtual-method lgamma (((i <integer>))
					(:union <real> <integer>)
					pure)
    (if (or (equal? i 1) (equal? i 2))
	0
	(r-lgamma (integer->real i))))

  (define-simple-virtual-method tgamma (((r <real-number>))
					(:union <real> <integer>)
					pure)
    (raise-simple 'tgamma:dispatch-error))

  (include-virtual-methods tgamma r-tgamma)

  (define-simple-virtual-method tgamma (((rat <rational>))
					(:union <real> <integer>)
					pure)
    (let* ((rat-s (simplify-rational rat))
	   (i-numer (numerator rat-s))
	   (i-denom (denominator rat-s)))
      (if (and (equal? i-denom 1) (>= i-numer 1))
	  (factorial (- i-numer 1))
	  (r-tgamma (rational->real rat-s)))))

  (define-simple-virtual-method tgamma (((i <integer>))
					(:union <real> <integer>)
					pure)
    (if (>= i 1)
	(factorial (- i 1))
	(r-tgamma (integer->real i)))))