Test gamma function using Maxima

maxima-gamma.rkt

#lang at-exp racket
(require racket/tcp racket/string)

;;; Maxima

;; This module starts an external Maxima process.
;; The function send will send a command to Maxima.
;; The function receive will get the output from Maxima as a list of strings.
;; The various send-* and receive-* functions sends and receives to and from Maxima.
;; The various read-* and display-* functions reads and displays to Racket (DrRacket).


;;; Configuration: Change maxima paths here.

(define PORT 8089)
(define MAXIMA-PATH "/usr/local/bin/maxima")


;;; Parameters

(define out (make-parameter #f)) ; output port for sending 
(define in  (make-parameter #f)) ; input  port for receiving
(define err (make-parameter #f)) ; error  port of maxima


;;; Sending

(define (send str)
  (sync (out))
  (display str (out))
  (flush-output (out)))

(define (send-command str)
  (send str))


;;; Receiving

(define (receive-line)
  (read-line (in)))

(define (receive-error)
  (read-line (err)))

(define (receive-welcome-message)
  ; Due to the flag --very-quiet the welcome is 
  ; a single line containing the pid.
  (list (receive-line)))

(define (maybe-receive-line)
  (if (sync/timeout 0 (in))
      (receive-line)
      #f))
  
(define (receive)
  (let ([first-line (receive-line)])
    (let loop ([lines (list first-line)])
      (let ([line (maybe-receive-line)])
        (if line
            (loop (cons line lines))
            (reverse lines))))))

(define (receive-whitespace)
  (let ([c (read-char (in))])
    (when (not (char-whitespace? c))
      (error 'read-whitespace "expected to receive whitespace " c))))


;;; String utilities

(define (string-ref-last str)
  (if (string=? "" str)
      #f
      (string-ref str (sub1 (string-length str)))))

;;; List utilies

(define (remove-last xs)
  (if (empty? xs) xs (drop-right xs 1)))


;;; Displaying

(define (display-line datum)
  (display datum)
  (newline))

(define (display-prompt prompt)
  (display prompt)
  (display " "))

(define (display-output lines)
  (unless (empty? lines)
    (display-lines
     lines)))


;;; Reading

(define (read-command)
  (let loop ([lines '()])
    (let ([line (read-line)])
      (if (memv (string-ref-last line) '(#\$ #\;))
          (string-append* (reverse (cons line lines)))
          (loop (cons line lines))))))


;;; REPL

(define (writeln x) (write x) (newline))

(define (read-send-receive-loop)
  (display-prompt ">")
  (define cmd (read-command))
  (send-command cmd)
  (define out (receive))
  (display-output out)
  (read-send-receive-loop))


;;; Start Maxima and REPL

(define (start)
  (let ([listener (tcp-listen PORT 3 #t)])
    (match-let
        ([(list pin pout pid perr status)
          (process* MAXIMA-PATH "--very-quiet" "-s" (format "~a" PORT))])
      (let-values ([(lin lout) (tcp-accept listener)])
        (parameterize ([in lin] [out lout])
          (displayln (receive-welcome-message))
          (display "Enter a Maxima command. Terminate a command with either ; or $ .\n")
          (read-send-receive-loop))))))

; (start)

;;; Testing 

(define (maxima-format z)
  ; Maxima expects complex numbers of the form   a+b%i
  (if (real? z)
      (~a (* 1.0 z))
      (~a (* 1.0 (real-part z)) " + " (* 1.0 (imag-part z)) "*%i")))

(define (maxima-read-number s)
  ; read real or complex number from string
  (define fragments
    (for/list ([x (in-port read (open-input-string (string-append* s)))])
      x))
  (match fragments
    [(list x)                          x             ]
    [(list      y '%i)                    (* y +1.0i) ]
    [(list x '- y '%i)           (-    x  (* y +1.0i))]
    [(list x '+ y '%i)           (+    x  (* y +1.0i))]
    [(list y '%i '+ x)           (+    x  (* y +1.0i))]
    [(list y '%i '- x)           (+ (- x) (* y +1.0i))]
    [(list (list '- y '%i) '- x) (- (- x) (* y +1.0i))]
    [_ (if (ormap (λ (s) (or (string-contains? s "FLOATING-POINT-OVERFLOW")
                             (string-contains? s "overflow")
                             (string-contains? s "error")))
                  (map ~a fragments))
           +nan.0
           (error 'maxima-read-number (~a "missed a case, got: " fragments)))]))



(define (start-maxima)
  (define listener (tcp-listen PORT 3 #t))
  (match-define (list pin pout pid perr status)
    (process* MAXIMA-PATH "--very-quiet" "-s" (format "~a" PORT)))
  (define-values (lin lout) (tcp-accept listener))
  (in  lin)
  (out lout)
  (receive-welcome-message)
  (void))


(define (maxima-gamma z)
  (define Z (maxima-format z))
  (send-command (~a "gamma(" Z ");"))
  (sleep .001) ; a small pause is needed 
  (maxima-read-number (receive)))


(start-maxima)
(for/list ([z '(-102.07898570979688-35.51381113633343i
                  -99.20308092474957+161.601565222518i       ; maxima overflow
                  -160.5756686534125-146.34588259270964i     ; maxima overflow
                  -44.59894031206599+8.143859564774402i
                  158.08675970929818-60.875079841822526i)])
    (maxima-gamma z))

;;;
;;; Compare Maxima and Wolfram results
;;;


;; (define wolfram-gamma-ht
;;   (with-input-from-file "gamma.rktd"
;;     (λ ()
;;       (read))))

;; (define results 
;;   (for/list ([(z expected) (in-hash wolfram-gamma-ht)])
;;     (displayln z)
;;     (list z (maxima-gamma z) (string->number expected))))

;; (define sorted-results
;;   ; sort results, biggest different first
;;   (sort results
;;         (λ (r1 r2)          
;;           (match (list r1 r2)
;;             [(list (list z1 m1 w1) (list z2 m2 w2))
;;              (cond
;;                [(eqv? m1 +nan.0) #t]
;;                [(eqv? m2 +nan.0) #f]
;;                [else             (> (magnitude (- m1 w1))
;;                                     (magnitude (- m2 w2)))])]))))