path: root/Transcen.scm

;; Copyright (C) 1992, 1993, 1995, 1997, 2005 Free Software Foundation, Inc.
;;
;; This program is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 2, or (at your option)
;; any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this software; see the file COPYING.  If not, write to
;; the Free Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111, USA.
;;
;; As a special exception, the Free Software Foundation gives permission
;; for additional uses of the text contained in its release of SCM.
;;
;; The exception is that, if you link the SCM library with other files
;; to produce an executable, this does not by itself cause the
;; resulting executable to be covered by the GNU General Public License.
;; Your use of that executable is in no way restricted on account of
;; linking the SCM library code into it.
;;
;; This exception does not however invalidate any other reasons why
;; the executable file might be covered by the GNU General Public License.
;;
;; This exception applies only to the code released by the
;; Free Software Foundation under the name SCM.  If you copy
;; code from other Free Software Foundation releases into a copy of
;; SCM, as the General Public License permits, the exception does
;; not apply to the code that you add in this way.  To avoid misleading
;; anyone as to the status of such modified files, you must delete
;; this exception notice from them.
;;
;; If you write modifications of your own for SCM, it is your choice
;; whether to permit this exception to apply to your modifications.
;; If you do not wish that, delete this exception notice.

;;;; "Transcen.scm", Complex trancendental functions for SCM.
;;; Author: Jerry D. Hedden.
;;;; 2005-05 SRFI-70 extensions.
;;; Author: Aubrey Jaffer

(define compile-allnumbers #t)		;for HOBBIT compiler

(define $pi (* 4 ($atan 1)))
(define pi $pi)
(define (pi* z) (* $pi z))
(define (pi/ z) (/ $pi z))

(define (exp z)
  (if (real? z) ($exp z)
      (make-polar ($exp (real-part z)) (imag-part z))))

(define (log z)
  (if (and (real? z) (>= z 0))
      ($log z)
      (make-rectangular ($log (magnitude z)) (angle z))))

(define (sqrt z)
  (if (real? z)
      (if (negative? z) (make-rectangular 0 ($sqrt (- z)))
	  ($sqrt z))
      (make-polar ($sqrt (magnitude z)) (/ (angle z) 2))))

(define (sinh z)
  (if (real? z) ($sinh z)
      (let ((x (real-part z)) (y (imag-part z)))
	(make-rectangular (* ($sinh x) ($cos y))
			  (* ($cosh x) ($sin y))))))
(define (cosh z)
  (if (real? z) ($cosh z)
      (let ((x (real-part z)) (y (imag-part z)))
	(make-rectangular (* ($cosh x) ($cos y))
			  (* ($sinh x) ($sin y))))))
(define (tanh z)
  (if (real? z) ($tanh z)
      (let* ((x (* 2 (real-part z)))
	     (y (* 2 (imag-part z)))
	     (w (+ ($cosh x) ($cos y))))
	(make-rectangular (/ ($sinh x) w) (/ ($sin y) w)))))

(define (asinh z)
  (if (real? z) ($asinh z)
      (log (+ z (sqrt (+ (* z z) 1))))))

(define (acosh z)
  (if (and (real? z) (>= z 1))
      ($acosh z)
      (log (+ z (sqrt (- (* z z) 1))))))

(define (atanh z)
  (if (and (real? z) (> z -1) (< z 1))
      ($atanh z)
      (/ (log (/ (+ 1 z) (- 1 z))) 2)))

(define (sin z)
  (if (real? z) ($sin z)
      (let ((x (real-part z)) (y (imag-part z)))
	(make-rectangular (* ($sin x) ($cosh y))
			  (* ($cos x) ($sinh y))))))
(define (cos z)
  (if (real? z) ($cos z)
      (let ((x (real-part z)) (y (imag-part z)))
	(make-rectangular (* ($cos x) ($cosh y))
			  (- (* ($sin x) ($sinh y)))))))
(define (tan z)
  (if (real? z) ($tan z)
      (let* ((x (* 2 (real-part z)))
	     (y (* 2 (imag-part z)))
	     (w (+ ($cos x) ($cosh y))))
	(make-rectangular (/ ($sin x) w) (/ ($sinh y) w)))))

(define (asin z)
  (if (and (real? z) (>= z -1) (<= z 1))
      ($asin z)
      (* -i (asinh (* +i z)))))

(define (acos z)
  (if (and (real? z) (>= z -1) (<= z 1))
      ($acos z)
      (+ (/ (angle -1) 2) (* +i (asinh (* +i z))))))

(define (atan z . y)
  (if (null? y)
      (if (real? z) ($atan z)
	  (/ (log (/ (- +i z) (+ +i z))) +2i))
      ($atan2 z (car y))))

;;;; SRFI-70
(define expt
  (let ((integer-expt integer-expt))
    (lambda (z1 z2)
      (cond ((and (exact? z2) (not (and (zero? z1) (not (positive? z2)))))
	     (integer-expt z1 z2))
	    ((and (real? z2) (real? z1) (positive? z1))
	     ($expt z1 z2))
	    (else
	     (exp (* (if (zero? z1) (real-part z2) z2) (log z1))))))))

(set! quotient
      (let ((integer-quotient quotient))
	(lambda (x1 x2)
	  (if (and (exact? x1) (exact? x2))
	      (integer-quotient x1 x2)
	      (truncate (/ x1 x2))))))

(set! remainder
      (let ((integer-remainder remainder))
	(lambda (x1 x2)
	  (if (and (exact? x1) (exact? x2))
	      (integer-remainder x1 x2)
	      (- x1 (* x2 (quotient x1 x2)))))))

(set! modulo
      (let ((integer-modulo modulo))
	(lambda (x1 x2)
	  (if (and (exact? x1) (exact? x2))
	      (integer-modulo x1 x2)
	      (- x1 (* x2 (floor (/ x1 x2))))))))

(define (infinite? z) (and (= z (* 2 z)) (not (zero? z))))
(define (finite? z) (not (infinite? z)))

(define (invintp f1 f2 f3)
  (define f1^2 (* f1 f1))
  (define f2^2 (* f2 f2))
  (define f3^2 (expt f3 2))
  (let ((c (+ (* -3 f1^2 f2)
	      (* 3 f1 f2^2)
	      (* (- (* 2 f1^2) f2^2) f3)
	      (* (- f2 (* 2 f1)) f3^2)))
	(b (+ (- f1^2 (* 2 f2^2)) f3^2))
	(a (- (* 2 f2) f1 f3)))
    (define disc (- (* b b) (* 4 a c)))
    (if (negative? (real-part disc))
	(/ b -2 a)
	(let ((sqrt-disc (sqrt disc)))
	  (define root+ (/ (- sqrt-disc b) 2 a))
	  (define root- (/ (+ sqrt-disc b) -2 a))
	  (if (< (magnitude (- root+ f1)) (magnitude (- root- f1)))
	      root+
	      root-)))))

(define (extrapolate-0 fs)
  (define n (length fs))
  (define (choose n k)
    (do ((kdx 1 (+ 1 kdx))
	 (prd 1 (/ (* (- n kdx -1) prd) kdx)))
	((> kdx k) prd)))
  (do ((k 1 (+ 1 k))
       (lst fs (cdr lst))
       (L 0 (+ (* -1 (expt -1 k) (choose n k) (car lst)) L)))
      ((null? lst) L)))

(define (sequence->limit proc sequence)
  (define lval (proc (car sequence)))
  (if (finite? lval)
      (let ((val (proc (cadr sequence))))
	(define h_n*nsamps (* (length sequence) (magnitude (- val lval))))
	(if (finite? val)
	    (let loop ((sequence (cddr sequence))
		       (fxs (list val lval))
		       (trend #f)
		       (ldelta (- val lval))
		       (jdx (+ -1 (length sequence))))
	      (cond ((null? sequence)
		     (case trend
		       ((diverging) (and (real? val) (* ldelta 1/0)))
		       ((bounded) (invintp val lval (caddr fxs)))
		       (else (cond ((zero? ldelta) val)
				   ((not (real? val)) #f)
				   (else (extrapolate-0 fxs))))))
		    (else
		     (set! lval val)
		     (set! val (proc (car sequence)))
		     (if (finite? val)
			 (let ((delta (- val lval)))
			   (define h_j (/ h_n*nsamps jdx))
			   (cond ((case trend
				    ((converging) (<= (magnitude delta) h_j))
				    ((bounded)    (<= (magnitude ldelta) (magnitude delta)))
				    ((diverging)  (>= (magnitude delta) h_j))
				    (else #f))
				  (loop (cdr sequence) (cons val fxs) trend delta (+ -1 jdx)))
				 (trend #f)
				 (else
				  (loop (cdr sequence) (cons val fxs)
					(cond ((> (magnitude delta) h_j) 'diverging)
					      ((< (magnitude ldelta) (magnitude delta)) 'bounded)
					      (else 'converging))
					delta (+ -1 jdx)))))
			 (and (eq? trend 'diverging) val)))))
	    (and (real? val) val)))
      (and (real? lval) lval)))

(define (limit proc x1 x2 . k)
  (set! k (if (null? k) 8 (car k)))
  (cond ((not (finite? x2)) (slib:error 'limit 'infinite 'x2 x2))
	((not (finite? x1))
	 (or (positive? (* x1 x2)) (slib:error 'limit 'start 'mismatch x1 x2))
	 (limit (lambda (x) (proc (/ x))) 0.0 (/ x2) k))
	((= x1 (+ x1 x2)) (slib:error 'limit 'null 'range x1 (+ x1 x2)))
	(else (let ((dec (/ x2 k)))
		(do ((x (+ x1 x2 0.0) (- x dec))
		     (cnt (+ -1 k) (+ -1 cnt))
		     (lst '() (cons x lst)))
		    ((negative? cnt)
		     (sequence->limit proc (reverse lst))))))))