(define-module (language python module math)
  #:use-module (language python for)
  #:use-module (language python try)
  #:use-module (language python exceptions)
  #:use-module (language python number)
  #:use-module (system foreign)
  #:use-module (rnrs bytevectors)
  
  #:export (ceil copysign fabs factorial floor fmod frexp fsum gcd isclose
		 isfinite isinf isnan modf trunc exp expm1 log log1p log2 log10
		 pow sqrt acos asin atan cos sin tan atan2 hypot pi
		 degrees radians acosh asinh atanh cosh sinh tanh
		 erf erfc gamma lgamma e tau inf nan))

(define (real! s x)
  (if (complex? x)
      (if (= (imag-part x) 0)
          (real-part x)
          (raise ValueError "real math fkn result in complex number" s))
      x))

(define ceil
  (lambda (x)
    ((@ (guile) inexac->exact)
     ((@ (guile) ceiling)
      x))))
  
(define (copysign x y)
  (let ((x ((@ (guile) abs) x)))
    (if (< y 0) (- x) x)))

(define (fabs x) ((@ (guile) abs) x))

(define (factorial x)
  (if (not (and (number? x) (integer? x)))
      (raise ValueError "Not an integer"))
  (if (< x 0)
      (raise ValueError "Negative integer"))

  (let lp ((x x))
    (if (= x 0)
	1
	(* x (lp (- x 1))))))

(define (floor x)
  ((@ (guile) inexac->exact)
   ((@ (guile) floor)
    x)))

(define (fmod x y) ((@ (guile) truncate-remainder) x y))

(define frexp
  (let ((f (pointer->procedure double
                               (dynamic-func "frexp" (dynamic-link))
                               (list double '*))))
    (lambda (x)
      (let* ((v  (make-bytevector 4))
	     (vp (bytevector->pointer v)))
	(list (f x) (bytevector-s32-ref v 0 (native-endianness)))))))

(define (fsum it)
  (for ((x : it)) ((s 0))
       (+ s x)
       #:final s))
       

(define gcd (@ (guile) gcd))

(define* (isclose a b #:key (rel_tol 1e-9) (abs_tol 0.0))
  (define e-abs (abs (- a b)))
  (if (< e-abs (max (* rel_tol (max (abs a) (abs b))) abs_tol))
      #t
      #f))

(define (isfinite x)
  (if (or (inf? x) (nan? x))
      #f
      #t))

(define isinf inf?)
(define isnan nan?)

(define (ldexp x i) (* x ((@ (guile) expt) 2 i)))

(define (modf x)
  (let* ((x1 (floor x))
	 (x2 (- x x1)))
    (values x2 x1)))

(define trunc py-trunc)

;; Power and logarithms
(define (exp x) (real! 'exp ((@ (guile) exp) x)))

(define expm1
  (let ((f (pointer->procedure double
                               (dynamic-func "expm1" (dynamic-link))
                               (list double))))
    (lambda (x) (f x))))

(define* (log x #:optional (base #f))
  (real! 'log
         (if (not base)
             ((@ (guile) log) x)
             (/ ((@ (guile) log) x) ((@ (guile) log) base)))))

(define log1p
  (let ((f (pointer->procedure double
                               (dynamic-func "log1p" (dynamic-link))
                               (list double))))

    (lambda (x) (f x))))

(define log2
  (let ((f (pointer->procedure double
                               (dynamic-func "log2" (dynamic-link))
                               (list double))))
    (lambda (x) (f x))))

(define (log10 x) (real! 'log10 ((@ (guile) log10) x)))
  
(define (pow x y) (real! 'pow ((@ (guile) expt) x y)))

(define (sqrt x) (real! 'sqrt ((@ (guile) sqrt) x)))

;; Trigs
(define (acos x) (real! 'acos ((@ (guile) acos) x)))
(define (asin x) (real! 'asin ((@ (guile) asin) x)))
(define (atan x) (real! 'atan ((@ (guile) atan) x)))
(define (cos  x) (real! 'cos ( (@ (guile) cos) x)))
(define (sin  x) (real! 'sin ( (@ (guile) sin) x)))
(define (tan  x) (real! 'tan ( (@ (guile) tan) x)))

(define atan2
  (let ((f (pointer->procedure double
                               (dynamic-func "atan2" (dynamic-link))
                               (list double double))))
    (lambda (x y) (f x y))))

(define hypot
  (let ((f (pointer->procedure double
                               (dynamic-func "hypot" (dynamic-link))
                               (list double double))))

    (lambda (x y) (f x y))))


;; angular conversion
(define pi (* 4 (atan 1)))

(define degrees
  (let ((f (/ 360 (* 2 pi))))
    (lambda (x) (* f x))))

(define radians
  (let ((f (/ (* 2 pi) 360)))
    (lambda (x) (* f x))))

;; Hyperbolic funcitons
(define (acosh x) (real! 'acosh ((@ (guile) acosh) x)))
(define (asinh x) (real! 'asinh ((@ (guile) asinh) x)))
(define (atanh x) (real! 'atanh ((@ (guile) atanh) x)))
(define (cosh  x) (real! 'cosh ( (@ (guile) cosh) x)))
(define (sinh  x) (real! 'sinh ( (@ (guile) sinh) x)))
(define (tanh  x) (real! 'tanh ( (@ (guile) tanh) x)))


;; Special functions
(define erf
  (let ((f (pointer->procedure double
                               (dynamic-func "erf" (dynamic-link))
                               (list double))))
    (lambda (x) (f x))))

(define erfc
  (let ((f (pointer->procedure double
                               (dynamic-func "erfc" (dynamic-link))
                               (list double))))
    (lambda (x) (f x))))

(define gamma
  (let ((f (pointer->procedure double
                               (dynamic-func "tgamma" (dynamic-link))
                               (list double))))

    (lambda (x)
      (if (integer? x)
	  (factorial (- x 1)) 
	  (f x)))))

(define lgamma
  (let ((f (pointer->procedure double
                               (dynamic-func "lgamma" (dynamic-link))
                               (list double))))
    (lambda (x)
      (f x))))

;; constants
(define e   (exp 1))
(define tau (* 2 pi))
(define inf ((@ (guile) inf)))
(define nan ((@ (guile) nan)))