;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;                                                          ;;;;;;;;
;;;;;;  All files in this directory or any subdirectories are   ;;;;;;;;
;;;;;;  copyright 1997, 1998, 1999, 2000, 2002.                 ;;;;;;;;
;;;;;;  by Rafael D. Sorkin.  All rights reserved.              ;;;;;;;;
;;;;;;                                                          ;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; bibliotek.float       Time-stamp:<2002-Oct-23 14:25:08 15798.59780>

; Containing functions involved with floating point computations, especially
; ones which are inexact for reasons beyond truncation error of the result
; itself.

;: Roster of functions defined herein 
;  ==================================
;
;  polynomial-value
;
;  factorial           (for now just relies on log-factorial)
;  log-factorial       (fairly accurate, useful for large arguments)
;
;  solve-monotone       (solve y=f(x) for x)
;
;  The following numerical constants are defined in `bibliotek.constants'
;
;    pi, e, Euler, (log (sqrt 2 pi)),  
;    and a bunch of fractions for elisp like 1/3 or -1/2.
;
;  Infinity NaN are defined in the file preparations.el
;
;  A bunch of matrix functions are also in a separate file

;; If package info is needed here, see { ~/lisp/bibliotek.TCL.l } for sample.

;: The functions 

(defun polynomial-value (p x) 
  "\
 The first arg should be the LIST of coefficients  (a_0 a_1...a_n), 
 with 0's for absent terms, of course.
 NB: the coefficients are in order of INCREASING DEGREE"
  (if (null p) 0
    (+
     (car p)
     (* x (polynomial-value (cdr p) x)))))
 ;
 ; This is the more accurate and efficient way to evaluate a polynomial.
 ; Seems no need to localize the symbols.


(deff log-factorial (x &optional (thresh 20.0))  
  "\
 Natural log of factorial of x, where x is a real number > -1.
 Evaluation is done using Stirling series, after first promoting x to
 be greater than the threshold specified by the optional second argument
 (which defaults to 20).
  "
  (cond
  ;-----------------------
  ; check that arg is > -1
  ;-----------------------
   ((<= x -1)
    (error "`log-factorial' requires arg > -1 for now."))
  ;------------------------------------------------------------------
  ; handle x=0 and 1 specially to avoid roundoff error in those cases
  ;------------------------------------------------------------------
   ((or 
     (= x 0)
     (= x 1))
              0.0)
  ;------------------------------------------
  ; if x < thresh then augment it and recurse
  ;------------------------------------------
   ((< x thresh)
    (- (log-factorial (1+ x) thresh) (o log 1+ x)))  
  ;------------------------------------
  ; if x >= thresh then plug into the series
  ;------------------------------------
   (t (+
       (* (+ x 1/2) (log x))
       (- x)
       log_sqrt_2pi
       (polynomial-value 
         (list 0 1/12 0 -1/360 0 1/1260 0 -1/1680)
         (reciprocal x))))))
  ;
  ; The method is to make the argument large enough that a few terms in the
  ; Stirling series will be quite accurate. 
  ;
  ; It seems to work well, but behaves rather strangely.
  ; First of all, one might have expected the optimum threshold to be
  ; around 33 (because there the ratio of the smallest to biggest term is
  ; about 10^16) but empirically 20 appears to work best (see the tests).
  ; Moreover, many choices of threshshold in the range 15 to 80 (if not
  ; beyond) give precisely the same answer (see the tests).
  ;
  ; Improvement
  ;  Let it handle args < -1 as well.  One nuisance though, is that log is
  ;  complex in parts of this range, so maybe should give just its real part
  ;  (i.e. (log abs fact x))
  ;
  ; The coefficients in the series come from Dwight 851.5 (page 210)
  ;  Here is how they were computed
  ;
  ; (* 6  1 2) ; => 12
  ; (* 30 3 4) ; => 360
  ; (* 42 5 6) ; => 1260
  ; (* 30 8 7) ; => 1680
  ;
  ; (- (log (sqrt (atl 2 * pi))) log_sqrt_2pi) ; => 0.0

;;; Above (log-factorial) is ridiculously slow in cmucl, WHY??
;;; See ~/lisp/developing/factorial.el for attempts to improve it! 

(defun factorial (x) 
  "\
 A poor excuse for this function scrabbled together from wherever.
 Among other problems, it rejects arguments < -1"
  (if (and (integerp x) (< x 171)) (! x)
    (o exp log-factorial x)))


;;; modifiying following to make gcl accept it, LOGIC MIGHT HAVE BEEN CHANGED!!

(deff solve-monotone (f y &key 
			   ((:ii ii) (list -1.0 1.0))
			   ((:tol-x tol_x) (* 5 float-epsilon)) 
			   ((:tol-y tol_y) (* 50 float-epsilon))
			   ((:maxit maxit) 512)
			   ((:lbound lbound))
			   ((:ubound ubound)))
  "\
 Solves the equation y=f(x) for x, where f is a monotonically 
 increasing or decreasing function.  Arguments are:
   1. the function f (or a symbol for it)
   2. the desired value y
   3. :ii = (a b) is the `initial interval' in which to begin the search
            (defaults to [-1 1]) 
   4. :tol-x = (absolute) tolerance for x (defaults to (* 5 float-epsilon))
   5. :tol-y = (absolute) tolerance for y (defaults to (* 50 float-epsilon))
   6. :maxit = maximum number of iterations before quitting (defaults to 512)
   7. :lbound and :ubound = limits of search "
 ;-------------------------------
 ;/localizations and declarations
 ;-------------------------------
  (&localize 
     f y ii tol_x tol_y maxit lbound ubound
     a b c sgn@a sgn@b sgn@c 
     g x fx u v w basta j plan)

  (declare 
   (optimize (speed 3))
   (type fixnum maxit)
   (type realfloat y tol_x tol_y))

  (varbind 
   c 0.0
   u 0.0
   v 0.0
   a 0.0
   b 0.0
   w 0.0
   sgn@a 5
   sgn@b 5
   sgn@c 5
   plan nil)
  (declare 
   (type real-float a b c u v w)
   (type fixnum sgn@a sgn@b sgn@c))

 ;---------------------------------------------------------------------------
 ;/make sure arguments have the types they're declared to have or should have
 ;---------------------------------------------------------------------------
  (setf
    tol_x (coerce tol_x *float*)  ; type info doesnt seem to propagate in cmucl
    tol_y (coerce tol_y *float*))
  (assert (>= tol_x 0))
 ;-------------------------------------------------
 ;/get bounds for initial interval (and cast types)
 ;-------------------------------------------------
  (setq 
    a (coerce (car  ii) *float*)
    b (coerce (cadr ii) *float*))
  (assert (< tol_x (- b a)) 
     nil "solve-monotone: wrongly specified initial interval or bad tol-x")
  ;;(princ (list a b))			;;;
 ;-----------------------------------------------------
 ;/define g(x) := f(x) - y   (We will seek a zero of g)
 ;-----------------------------------------------------
  (fbind g(x) 
    (declare (type realfloat x))
    ;; (princ (list a b))					;;;
    (vbind fx (funcall f x))
    (if *tcl* (assert (realp fx)))
    (- fx y))
 ;---------------
 ;/define `basta'
 ;---------------
  (fbind basta ()
    (setq c (/ (+ a b) 2.0))
    (if (> (- b a) tol_x)
      (error 
        (tcl-or-elisp
	  "solve-monotone: ~s might be a solution if :tol-x were ~0,1g"
	  "solve-monotone: %s might be a solution if :tol-x were %0.1g")
	c (- b a)))
    (if (> (o abs g c) tol_y)
      (error 
        (tcl-or-elisp
	  "solve-monotone: ~s might be a solution if :tol-y were ~0,1g"
	  "solve-monotone: %s might be a solution if :tol-y were %0.1g")
	c (o abs g a)))
    (return-from solve-monotone c))
 ;-----------------------
 ;/begin iterative search
 ;-----------------------
  (psetq u (g a) v (g b))
  ;; (princ (list u v))			;;;
  (loop 
   for j from 0
  ;--------------------------------------
  ;/time to give up (por too many tries)?
  ;--------------------------------------
   if (> j maxit) 
   do (error 
       (tcl-or-elisp
	 "more than ~d tries in `solve-monotone'" 
	 "more than %d tries in `solve-monotone'")
       maxit)
  ;----------------------------------
  ;/find the signs of f at a and at b
  ;----------------------------------
   do (psetq sgn@a (sgn u) sgn@b (sgn v))
  ;----------------------------------------------
  ;/return if either a or b is already a solution
  ;----------------------------------------------
   if (= 0 sgn@a) return a
   if (= 0 sgn@b) return b
  ;-----------------------------
  ;/call `basta' if b-a <= tol_x
  ;-----------------------------
   if (<= (- b a) tol_x) do (basta)
  ;-----------------------------------------------------
  ;/determine whether root is within (a b) or outside it
  ;-----------------------------------------------------
   do (if (= sgn@a sgn@b) (setq plan 'caminar) (setq plan 'narrow))
  ;-------------------------
  ;/execute plan accordingly
  ;-------------------------
   do
   (case plan
    ;---------------------------------
    ;/case where root is outside (a b)
    ;---------------------------------
     ((caminar)
      (cond
	((< u v) (setq plan (if (> sgn@a 0) 'left  'right)))
	((> u v) (setq plan (if (< sgn@a 0) 'left  'right)))
	(t (basta)))
      (case plan
	((left)
	 (psetq 
	   a (- (* 3 a) (* 2 b))
	   b a)
	 (if (and lbound (< a lbound))(error "out of bounds, can't continue"))
	 (psetq 
	   u (g a)
	   v u))
	((right)
	 (psetq 
	   a b
	   b (- (* 3 b) (* 2 a)))
	 (if (and ubound (> b ubound))(error "out of bounds, can't continue")) 
	 (psetq 
	   u v
	   v (g b)))))
    ;--------------------------------
    ;/case where root is inside  (a b)
    ;--------------------------------
     ((narrow)
      (setq c (/ (+ a b) 2.0))
      (if (or (= a c) (= c b)) (basta))
      (setq 
       w (g c)
       sgn@c (sgn w))
      (if (= 0 sgn@c) (return c))
      (cond
       ((/= sgn@a sgn@c) (setq b c v (g b)))
       (t          (setq a c u (g a)))))))) 
 ;
 ; NOTES
 ;
 ; Some variables
 ;  c = midpoint of interval [a b] 
 ;  u = f(a), v = f(b), w = g(c)
 ;  sgn@a = sgn(u)  etc
 ;
 ; The search strategy 
 ;   if root is within interval, buscalo recursvely in the half where it is;
 ;   if it seems outside, look in the adjacent interval of twice the size.
 ;
 ; The function `basta' takes over when iteration stops. 
 ; It decides whether either issues an error mesage or returns the root found.
 ;
 ; We will have trouble if f has a big flat stretch, but there is probably no
 ; point in trying to deal with this for now (we could do it by giving extra
 ; information to know whether it is monotone increasing or decreasing, or a
 ; memory of this) 
 ;
 ; A problem that was cured was that for b-a very small c can equal one of
 ; them, then you just keep repeating the same interval (a b), now we catch 
 ; that and call basta.
 ;
 ; Localization of variables can be important. In an earlier version,
 ; `x' was left unlocalized, causing trouble. 
 ;
 ; Possible Improvements (see developing/solve-monotone.l)


;: End