1;;; -*- *cold-load:t; 1Mode:Common-Lisp; Package:SYSTEM-INTERNALS; Base:8; Fonts:(CPTFONT CPTFONTB) -*-

;;;                           RESTRICTED RIGHTS LEGEND

;;;Use, duplication, or disclosure by the Government is subject to
;;;restrictions as set forth in subdivision (c)(1)(ii) of the Rights in
;;;Technical Data and Computer Software clause at 52.227-7013.
;;;
;;;                     TEXAS INSTRUMENTS INCORPORATED.
;;;                              P.O. BOX 2909
;;;                           AUSTIN, TEXAS 78769
;;;                                 MS 2151
;;;
;;; Copyright (C) 1985-1989 Texas Instruments Incorporated. All rights reserved.
;*   1** (c) Copyright 1980 Massachusetts Institute of Technology ***

(DEFSUBST random-state-p (object)
  1"T if OBJECT is a random-state -- a seed for use by random."*
  (TYPEP object 'random-state)) 

(DEFSUBST logtest (bits word)
  1"T if the bits specified by BITS in WORD are not all zero.
   BITS is a mask in which the bits to be tested are ones."*
  (NOT (ZEROP (LOGAND bits word)))) 

(DEFF zlc:bit-test #'LOGTEST) 

(DEFSUBST ldb-test (ppss word)
  1"T if the field specified by PPSS in WORD is not zero.
   PPSS is a position (from the right) times 64., plus a size."*
  (NOT (ZEROP (LDB ppss word)))) 

(DEFSUBST %logldb-test (ppss word)
  1"T if the field specified by PPSS in WORD is not zero.
PPSS is a position (from the right) times 64., plus a size.*
Uses %LOGLDB instead of LDB, so it can access the sign bit but only works on fixnums.
1Like LDB-TEST except that when SETF'd it does a %LOGDPB rather than a dpb."*
  (NOT (ZEROP (%LOGLDB ppss word)))) 

(DEFSUBST byte (width position)
  1"Return a byte specifier for a byte WIDTH bits long starting POSITION bits from the lsb."*
  (+ width (DPB position %%byte-specifier-position 0))) 

(DEFSUBST byte-position (byte-specifier)
  1"Return the byte position specified by BYTE-SPECIFIER.
This is the index of the least significant bit included in the byte."*
  (LDB %%byte-specifier-position byte-specifier)) 

(DEFSUBST byte-size (byte-specifier)
  1"Return the byte size specified by byte-specifier."*
  (LDB %%byte-specifier-size byte-specifier)) 

(DEFSUBST logbitp (index integer)
  1"T if INTEGER's binary representation contains the 2^INDEX bit as a 1."*
  (LDB-TEST (BYTE 1 index) integer)) 

(DEFUN ldb-hard (ppss word)
   (let ((width (byte-size ppss))
	 (position (byte-position ppss))
	 (result-position 0)
	 (result 0))
     (loop until (< width 24.)
	   doing (setq result (dpb (ldb (byte 24. position) word)
				   (byte 24. result-position)
				   result)
		       width	       (- width 24.)
		       result-position (+ result-position 24.)
		       position	       (+ position 24.))
	   finally (return (dpb (ldb (byte width position) word)
				(byte width result-position)
				result)))))

(DEFUN dpb-hard (value ppss word)
   (let ((width (byte-size ppss))
	 (position (byte-position ppss))
	 (value-position 0)
	 (result word))
     (loop until (< width 24.)
	   doing (setq result (dpb (ldb (byte 24. value-position) value)
				   (byte 24. position)
				   result)
		       width            (- width 24.)
		       position		(+ position 24.)
		       value-position	(+ value-position 24.))
	   finally (return (dpb (ldb (byte width value-position) value)
				(byte width position)
				result)))))

(defun mask-field-hard (ppss integer)
   (check-type integer integer)
   (logand integer (byte-mask ppss)))		1; Taken from CLtL.*

(defun deposit-field-hard (value ppss integer)
   (check-type integer integer)
   (check-type value integer)
   (logxor integer (logand (byte-mask ppss) (logxor integer value))))	1; Taken from CLtL.*

(DEFUN floor (dividend &optional (divisor 1))
  "Return DIVIDEND divided by DIVISOR, rounded down, and the remainder."
  (DECLARE (VALUES quotient remainder))
  (FLOOR dividend divisor)) 

(proclaim '(inline sqr poly poly1))

(DEFUN sqr (x)
  (* x x))

(DEFUN poly (w cl)
1  ;;w is a floating-point number. cl is a list of
  ;;the coefficients of a polynomial in w.  The coefficients are in
  ;;descending order of the power of w.  Result is the value of the
  ;;polynomial*
  (LET ((s (CAR cl)))
    (DOLIST (c (CDR cl) s)
      (SETQ s (+ (* s w) c)))))

(DEFUN poly1 (w cl)
1  ;;w is a floating-point number.  cl is a list of
  ;;the coefficients of a polynomial in w, except that a leading one is
  ;;not in the list.  Result is the value of the polynomial.*
  (LET ((s (+ w (CAR cl))))
     (DOLIST (c (CDR cl) s)
       (SETQ s (+ (* s w) c)))))

(DEFSUBST copy-float (n)
1  "Make a copy of a floating point number."*
  (+ n (DONT-OPTIMIZE 0.0)))			;1/31/85 DNG - Original.


1;; Same as BYTE, but for use below, when BYTE is not loaded yet.*

(DEFUN xbyte (size position)
  (DPB position %%byte-specifier-position size)) 

(DEFCONSTANT most-negative-fixnum (%LOGDPB 1 %%q-boxed-sign-bit 0)
1  "Any integer smaller than this must be a bignum.") *
(DEFCONSTANT most-positive-fixnum (%LOGDPB 0 %%q-boxed-sign-bit -1)
1  "Any integer larger than this must be a bignum.") *

;;; The following constants can be used in a native IEEE environment
;;;(DEFCONSTANT most-positive-short-float 3.4027977s38
;;;1  "No short float can be greater than this number."*) 
;;;(DEFCONSTANT most-negative-short-float -3.4027977s38
;;;  1"No short float can be less than this number."*) 
;;;(DEFCONSTANT least-positive-short-float 1.1754944s-38
;;;  1"No positive short float can be closer to zero than this number."*) 
;;;(DEFCONSTANT least-negative-short-float -1.1754944s-38
;;;  1"No negative short float can be closer to zero than this number."*) 
;;;(DEFCONSTANT short-float-epsilon 7.629511s-6
;;;  1"Smallest positive short float which can be added to 1.0s0 and make a difference."*) 
;;;(DEFCONSTANT short-float-negative-epsilon 3.8147555s-6
;;;  1"Smallest positive single float which can be subtracted from 1.0s0 and make a difference."*) 

(DEFCONSTANT most-positive-short-float (%MAKE-POINTER DTP-Short-Float #xFEFFFF)
  1"No short float can be greater than this number."*) 
(DEFCONSTANT most-negative-short-float (%MAKE-POINTER DTP-Short-Float #x-10001)
  1"No short float can be less than this number."*) 
(DEFCONSTANT least-positive-short-float (%MAKE-POINTER DTP-Short-Float #x010000)
  1"No positive short float can be closer to zero than this number."*) 
(DEFCONSTANT least-negative-short-float (%MAKE-POINTER DTP-Short-Float #x-FF0000)
  1"No negative short float can be closer to zero than this number."*) 
(DEFCONSTANT short-float-epsilon (%MAKE-POINTER DTP-Short-Float #x6E0001)
  1"Smallest positive short float which can be added to 1.0s0 and make a difference."*) 
(DEFCONSTANT short-float-negative-epsilon (%MAKE-POINTER DTP-Short-Float #x6D0001)
  1"Smallest positive single float which can be subtracted from 1.0s0 and make a difference."*) 

(DEFCONSTANT most-positive-single-float 3.4028235f38
  1"No single float can be greater than this number."*) 
(DEFCONSTANT most-negative-single-float -3.4028235f38
  1"No single float can be less than this number."*) 
(DEFCONSTANT least-positive-single-float 1.1754944f-38
  1"No positive single float can be between zero and this number."*) 
(DEFCONSTANT least-negative-single-float -1.1754944f-38
  1"No negative single float can be between zero and this number."*) 
(DEFCONSTANT single-float-epsilon 5.960465f-8
  1"Smallest positive float which can be added to 1.0f0 and make a difference."*) 
(DEFCONSTANT single-float-negative-epsilon 2.9802326e-8
  1"Smallest positive single float which can be subtracted from 1.0f0 and make a difference."*) 

;;; The following constants can be used in a native IEEE environment
(DEFCONSTANT most-positive-double-float 1.7976931348623157d308
  1"No double float can be greater than this number."*) 
(DEFCONSTANT most-negative-double-float -1.7976931348623157d308 
  1"No double float can be less than this number."*) 
(DEFCONSTANT least-positive-double-float 2.2250738585072014d-308 
  1"No positive double float can be between zero and this number."*) 
(DEFCONSTANT least-negative-double-float -2.2250738585072014d-308 
  1"No negative double float can be between zero and this number."*) 
(DEFCONSTANT double-float-epsilon 1.1102230246251568d-16
  1"Smallest positive double float which can be added to 1.0d0 and make a difference."*) 
(DEFCONSTANT double-float-negative-epsilon 5.551115123125784d-17 
  1"Smallest positive double float which can be subtracted from 1.0d0 and make a difference."*) 

;;;(DEFCONSTANT most-positive-double-float most-positive-single-float
;;;  1"No double float can be greater than this number."*) 
;;;(DEFCONSTANT most-negative-double-float most-negative-single-float
;;;  1"No double float can be less than this number."*) 
;;;(DEFCONSTANT least-positive-double-float least-positive-single-float
;;;  1"No positive double float can be between zero and this number."*) 
;;;(DEFCONSTANT least-negative-double-float least-negative-single-float
;;;  1"No negative double float can be between zero and this number."*) 
;;;(DEFCONSTANT double-float-epsilon single-float-epsilon
;;;  1"Smallest positive double float which can be added to 1.0d0 and make a difference."*) 
;;;(DEFCONSTANT double-float-negative-epsilon single-float-negative-epsilon
;;;  1"Smallest positive double float which can be subtracted from 1.0d0 and make a difference."*) 

(DEFCONSTANT most-positive-long-float most-positive-double-float
  1"No long float can be greater than this number."*) 
(DEFCONSTANT most-negative-long-float most-negative-double-float
  1"No long float can be less than this number."*) 
(DEFCONSTANT least-positive-long-float least-positive-double-float
  1"No positive long float can be between zero and this number."*) 
(DEFCONSTANT least-negative-long-float least-negative-double-float
  1"No negative long float can be between zero and this number."*) 
(DEFCONSTANT long-float-epsilon double-float-epsilon
  1"Smallest positive long float which can be added to 1.02l*0 and make a difference."*) 
(DEFCONSTANT long-float-negative-epsilon double-float-negative-epsilon
  1"Smallest positive long float which can be subtracted from 1.02l*0 and make a difference."*) 

(DEFCONSTANT pi 3.14159265358979323846d0 1"pi"*) 

(DEFUN ceiling (dividend &optional (divisor 1))
  1"Return DIVIDEND divided by DIVISOR, rounded up, and the remainder."*
  (DECLARE (VALUES quotient remainder))
  (CEILING dividend divisor)) 

(DEFUN truncate (dividend &optional (divisor 1))
  1"Return DIVIDEND divided by DIVISOR, rounded toward zero, and the remainder."*
  (DECLARE (VALUES quotient remainder))
  (TRUNCATE dividend divisor)) 

(DEFUN round (dividend &optional (divisor 1))
  1"Return DIVIDEND divided by DIVISOR, rounded to nearest integer, and the remainder."*
  (DECLARE (VALUES quotient remainder))
  (ROUND dividend divisor)) 

1;; Number comparison functions.

;MAX and MIN must be written with a single REST arg, otherwise,
; the hack of (APPLY #'MAX xx) can lose because it will try to
; copy the arglist to the stack.*
(DEFUN max (&rest numbers)
  1"Return the largest of the arguments."*
  (LET ((arg0 (CAR numbers)))
    (CHECK-TYPE arg0 number)
    (DO ((REST (CDR numbers) (CDR rest)))
	((NULL rest) arg0)
      (SETQ arg0 (MAX arg0 (CAR rest)))))) 

(DEFUN min (&rest numbers)
  1"Return the smallest of the arguments."*
  (LET ((arg0 (CAR numbers)))
    (CHECK-TYPE arg0 number)
    (DO ((REST (CDR numbers) (CDR rest)))
	((NULL rest) arg0)
      (SETQ arg0 (MIN arg0 (CAR rest)))))) 

(DEFUN > (number &rest numbers)
  1"Return T if the arguments are in strictly decreasing numerical order."*
;  (CHECK-TYPE NUMBER NUMBER)
  (DO ((a number c)
       (b numbers (CDR b))
       (c))
      ((NULL b) t)
    (SETQ c (CAR b))
    (OR (> a c) (RETURN ())))) 

(DEFUN < (number &rest numbers)
  1"Return T if the arguments are in strictly increasing numerical order."*
;  (CHECK-TYPE NUMBER NUMBER)
  (DO ((a number c)
       (b numbers (CDR b))
       (c))
      ((NULL b) t)
    (SETQ c (CAR b))
    (OR (< a c) (RETURN nil)))) 

(DEFUN = (&rest numbers)
  1"Return T if the arguments are all numerically equal."*
  (DOLIST (n (CDR numbers) t)
    (UNLESS (= n (CAR numbers))
      (RETURN nil)))) 

(DEFUN <= (number &rest numbers)
  1"Return T if the arguments are in nondecreasing numerical order."*
;  (CHECK-TYPE NUMBER NUMBER)
  (DO ((a number c)
       (b numbers (CDR b))
       (c))
      ((NULL b) t)
    (SETQ c (CAR b))
    (IF (> a c) (RETURN nil)))) 

(DEFUN >= (number &rest numbers)
  1"Return T if the arguments are in nonincreasing numerical order."*
;  (CHECK-TYPE NUMBER NUMBER)
  (DO ((a number c)
       (b numbers (CDR b))
       (c))
      ((NULL b) t)
    (SETQ c (CAR b))
    (IF (< a c) (RETURN nil)))) 

(DEFUN /= (&rest numbers)
  1"Return T if no two arguments are equal."*
  (DO ((REST numbers (CDR rest)))
      ((NULL (CDR rest)) t)
    (WHEN (MEMBER (CAR rest) (CDR rest) :test '=)
      (RETURN nil)))) 

(DEFF zlc: #'>=) 
(DEFF zlc: #'<=) 
(DEFF zlc: #'/=) 

1;; Arithmetic functions.*

(DEFUN + (&rest numbers)
  1"Return the sum of the arguments."*
  (DO ((numbers numbers (CDR numbers))
       (ans 0))
      ((NULL numbers) ans)
    (SETQ ans (+ ans (CAR numbers)))))

(DEFUN zlc:difference (number &rest numbers)
  1"Return the first argument minus the remaining arguments."*
  (DO ((numbers numbers (CDR numbers))
       (ans number))
      ((NULL numbers) ans)
    (SETQ ans (- ans (CAR numbers))))) 

(DEFUN - (number &rest numbers)
  1"Return the negation of a single argument, or the first argument minus the rest."*
  (COND ((NULL numbers) (- number))
	((DO ((numbers numbers (CDR numbers))
	      (ans number))
	     ((NULL numbers) ans)
	   (SETQ ans (- ans (CAR numbers))))))) 

(DEFUN * (&rest numbers)
  1"Return the product of the arguments."*
  (DO ((numbers numbers (CDR numbers))
       (ans 1))
      ((NULL numbers) ans)
    (SETQ ans (* ans (CAR numbers))))) 

(DEFUN quotient (number &rest numbers)
  1"Return the first argument divided by the rest.
If all arguments are fixnums, the value is also a fixnum, truncated! 
If you want a correct numeric division, float one argument."*
  (DO ((numbers numbers (CDR numbers))
       (ans number))
      ((NULL numbers) ans)
    (SETQ ans (global:/ ans (CAR numbers)))))

(DEFUN mod (dividend modulus)
  1"Return DIVIDEND taken in modulus MODULUS.
It is the same as the second value (the remainder) from (FLOOR DIVIDEND MODULUS).
The result will be in the range from zero (inclusive) to MODULUS (exclusive) with the same
sign as modulus."*
  (MOD dividend modulus)) 

(DEFUN rem (dividend modulus)
  1"Return DIVIDEND taken in modulus MODULUS.
It is the same as the second value (the remainder) from (TRUNCATE DIVIDEND MODULUS).
It has the same sign as the dividend, if it is not zero."*
  (REM dividend modulus))

(DEFF zlc:remainder #'REM) 

(DEFUN global:/ (number &rest numbers)
  1"Return the reciprocal of one argument, or the first argument divided by the rest.
If all arguments are fixnums, the value is also a fixnum, truncated!
If you want a correct numeric division, float one argument."*
  (COND ((NULL numbers) (global:/ 1 number))
	((DO ((numbers numbers (CDR numbers))
	      (ans number))
	     ((NULL numbers) ans)
	   (SETQ ans (global:/ ans (CAR numbers))))))) 

(DEFUN / (number &rest numbers)
  1"Return the reciprocal of one argument, or the first argument divided by the rest.
Division of fixnums returns a rational."*
  (COND ((NULL numbers) (%DIV 1 number))
	((DO ((numbers numbers (CDR numbers))
	      (ans number))
	     ((NULL numbers) ans)
	   (SETQ ans (%DIV ans (CAR numbers))))))) 

(DEFUN gcd (&rest numbers)
  1"Return the greatest common divisor of the arguments."*
  (LOOP with answer = 0
	for z in numbers
	do (SETQ answer (gcd answer z))
	finally (RETURN answer))) 

(DEFF \\\\ #'GCD)

(DEFUN lcm (&rest numbers)
  1"Return the least common multiple of all the numbers."*
  ;;  4/29/89 DNG - Fixed to permit no arguments.  [SPR 7112]
  (if (null numbers)
      1
    (DO ((value (ABS (first numbers)))
	 (REST (rest numbers) (CDR rest)))
	((NULL rest) value)
      (SETQ value (IF (OR (ZEROP value) (ZEROP (CAR rest)))
		      (RETURN 0)
		    (TRUNCATE (ABS (* value (CAR rest)))
			      (GCD value (CAR rest))))))))

1;; this should be made faster for negative numbers.*
(DEFUN integer-length (integer)
  1"Number of bits in field needed to store INTEGER without truncating it.
For nonnegative integers, this gives the exact number of bits
needed in an unsigned field.
For all integers, it gives the number of bits aside from the sign bit
needed in a signed field."*
  (CHECK-TYPE integer integer)
  (HAULONG
    (IF (MINUSP integer)
	(1+ integer)
	integer))) 

1;; Bitwise boolean operations on numbers.*

(DEFCONSTANT boole-clr   0 1"As first arg to BOOLE, makes all output bits be zero."*) 

(DEFCONSTANT boole-and   1 1"As first arg to BOOLE, makes output bits be one if both input bits are 1."*) 

(DEFCONSTANT boole-andc1 2 1"As first arg to BOOLE, makes output bits be one if first input is 0 and second is 1."*) 

(DEFCONSTANT boole-2     3 1"As first arg to BOOLE, makes output be the second input."*) 

(DEFCONSTANT boole-andc2 4 1"As first arg to BOOLE, makes output bits be one if first input is 1 and second is 0."*) 

(DEFCONSTANT boole-1     5 1"As first arg to BOOLE, makes output be the first input."*) 

(DEFCONSTANT boole-xor   6 1"As first arg to BOOLE, makes output bits be one if an odd number of input bits are one."*) 

(DEFCONSTANT boole-ior   7 1"As first arg to BOOLE, makes output bits be one if either input bit is 1."*) 

(DEFCONSTANT boole-nor  10 1"As first arg to BOOLE, makes output bits be one if both input bits are zero."*) 

(DEFCONSTANT boole-eqv  11 1"As first arg to BOOLE, makes output bits be one if an even number of input bits are one."*) 

(DEFCONSTANT boole-c1   12 1"As first arg to BOOLE, makes output be complement of first input."*) 

(DEFCONSTANT boole-orc1 13 1"As first arg to BOOLE, makes output bits be one if first input is 0 or second is 1."*) 

(DEFCONSTANT boole-c2   14 1"As first arg to BOOLE, makes output be complement of second input."*) 

(DEFCONSTANT boole-orc2 15 1"As first arg to BOOLE, makes output bits be one if first input is 1 or second is 0."*) 

(DEFCONSTANT boole-nand 16 1"As first arg to BOOLE, makes output bits be one unless both input bits are one."*) 

(DEFCONSTANT boole-set  17 1"As first arg to BOOLE, makes all output bits be one."*) 

(DEFUN logand (&rest integers)
  1"Bitwise-AND all the arguments."*
  (DO ((ans -1 (LOGAND ans (CAR l)))
       (l integers (CDR l)))
      ((NULL l) ans)))

(DEFUN logior (&rest integers)
  1"Bitwise-OR all the arguments."*
  (DO ((ans 0 (LOGIOR ans (CAR l)))
       (l integers (CDR l)))
      ((NULL l) ans)))

(DEFUN logxor (&rest integers)
  1"Bitwise-XOR all the arguments."*
  (DO ((ans 0 (LOGXOR ans (CAR l)))
       (l integers (CDR l)))
      ((NULL l) ans)))

(DEFUN logeqv (&rest integers)
  1"Bitwise-EQV all the arguments."*
  (DO ((ans -1 (*boole boole-eqv ans (CAR l)))
       (l integers (CDR l)))
      ((NULL l) ans)))

(DEFSUBST lognand (integer1 integer2)
  1"Bitwise-NAND the arguments.  Result bit is 1 if either INTEGER1 bit or INTEGER2 bit is 0."*
  (*boole boole-nand integer1 integer2)) 

(DEFSUBST lognor (integer1 integer2)
  1"Bitwise-NOR the arguments.  Result bit is 1 if INTEGER1 bit and INTEGER2 bit are both 0."*
  (*boole boole-nor integer1 integer2)) 

(DEFSUBST logorc1 (integer1 integer2)
  1"Bitwise-ORC1 the arguments.  Result bit is 1 if INTEGER1 bit is 0 or INTEGER2 bit is 1."*
  (*boole boole-orc1 integer1 integer2)) 

(DEFSUBST logorc2 (integer1 integer2)
  1"Bitwise-ORC2 the arguments.  Result bit is 1 if INTEGER1 bit is 1 or INTEGER2 bit is 0."*
  (*boole boole-orc2 integer1 integer2)) 

(DEFSUBST logandc1 (integer1 integer2)
  1"Bitwise-ANDC1 the arguments.  Result bit is 1 if INTEGER1 bit is 0 and INTEGER2 bit is 1."*
  (*boole boole-andc1 integer1 integer2)) 

(DEFSUBST logandc2 (integer1 integer2)
  1"Bitwise-ANDC2 the arguments.  Result bit is 1 if INTEGER1 bit is 1 and INTEGER2 bit is 0."*
  (*boole boole-andc2 integer1 integer2)) 

(DEFUN boole (op arg1 &rest args)
  1"Perform any of the 16 two-operand bitwise operations on ARG1 and ARGS.
OP is a number from 0 to 15 specifying the operation to use.
If there are more than two args (aside from OP)
then the first two are combined, then the result with the next arg, etc.
OP is bit-decoded:
 the 8 bit is the result when applied to 0 and 0,
 the 4 bit is the result when applied to 0 and 1,
 the 2 bit is the result when applied to 1 and 0,
 the 1 bit is the result when applied to 1 and 1.
The constants BOOLE-AND, etc., are provided for use as op."*
  (DO ((ans arg1 (*boole op ans (CAR l)))
       (l args (CDR l)))
      ((NULL l) ans))) 

(DEFSUBST lognot (integer)
  1"Return the bitwise complement of integer."*
  (LOGXOR integer -1)) 

(proclaim '(inline count-bits))

(DEFUN count-bits (x)
  (DO ((n x (LOGAND n (1- n)))
       (cnt 0 (1+ cnt)))
      ((ZEROP n) cnt)))

(DEFUN logcount (integer)
  1"Count number of bits set in INTEGER's binary representation.
Counts number of 1's in a positive INTEGER
 or number of 0's in a negative integer."*
  (ETYPECASE integer
    (fixnum
     (count-bits (IF (MINUSP integer) (LOGNOT integer) integer)))
    (bignum
     (let* ((res 0)
	    (words (CEILING (INTEGER-LENGTH integer) 16.)))
       (cond
	 ((PLUSP integer)
	  (dotimes (i words)
	    (incf res (count-bits (LDB (BYTE 16. (ASH i 4.)) integer)))))
	 (t
	  (dotimes (i words)
	    (incf res (count-bits (logxor (LDB (BYTE 16. (ASH i 4.)) integer) #xFFFF))))))
       res))))

(DEFUN ZLC:LOAD-BYTE (from-value position width)
  1"Return a byte extracted from FROM-VALUE, of specified WIDTH and POSITION.
POSITION counts from zero at the least significant bit."*
  (LDB (BYTE width position) from-value)) 

(DEFUN zlc:deposit-byte (into-value position width byte-value)
  1"Deposit BYTE-VALUE into a byte in INTO-VALUE of specified WIDTH and POSITION.
POSITION counts from zero at the least significant bit."*
  (DPB byte-value (BYTE width position) into-value)) 

(DEFUN haipart (x n &aux tem)
  1"Return N significant bits of the absolute value of X.
N > 0 means high N bits; N < 0 means low -N bits.
If X is too small, all of it is returned."*
  ;; Get number of significant bits
  (SETQ tem (HAULONG (SETQ x (ABS x))))
  (COND						1;Positive N means get high N bits, or as many as there are*
	((> n 0) (SETQ tem (- n tem))		1;minus number of low bits to discard*
		 (COND ((< tem 0) (ASH x tem))
		       (t x)))
	((ZEROP n) 0)				1;Zero N means return no bits*
	((< (SETQ n (- n)) tem)			1;Negative N means get low -N bits, or as many as there are*
	 (REM x (ASH 1 n)))
	(t x)))

1;; Special floating arithmetic functions.*
 
(DEFUN float (number &optional (other number))
  1"Convert NUMBER to a floating point number of same precision as OTHER.
If OTHER is omitted, and NUMBER is not already a float, a single-float is returned."*
;;  (DECLARE (optimize (speed 0)(safety 3)))
  (typecase other
    (short-float (small-float number))
    (double-float (double-float number))
    (t (internal-float number))))

(DEFUN ffloor (dividend &optional (divisor 1))
  1"Like FLOOR but converts first value to a float."*
  (DECLARE (VALUES quotient remainder))
  (MULTIPLE-VALUE-BIND (quotient remainder)
      (FLOOR dividend divisor)
    (VALUES (FLOAT quotient dividend) remainder))) 

(DEFUN fceiling (dividend &optional (divisor 1))
  1"Like CEILING but converts first value to a float."*
  (DECLARE (VALUES quotient remainder))
  (MULTIPLE-VALUE-BIND (quotient remainder)
      (CEILING dividend divisor)
    (VALUES (FLOAT quotient dividend) remainder))) 

(DEFUN ftruncate (dividend &optional (divisor 1))
  1"Like TRUNCATE but converts first value to a float."*
  (DECLARE (VALUES quotient remainder))
  (MULTIPLE-VALUE-BIND (quotient remainder)
      (TRUNCATE dividend divisor)
    (VALUES (FLOAT quotient dividend) remainder))) 

(DEFUN fround (dividend &optional (divisor 1))
  1"Like ROUND but converts first value to a float."*
  (DECLARE (VALUES quotient remainder))
  (MULTIPLE-VALUE-BIND (quotient remainder)
      (ROUND dividend divisor)
    (VALUES (FLOAT quotient dividend) remainder))) 

(DEFUN float-radix (float)
  1"Returns the radix of FLOAT, which is always 2 (binary)."*
  (check-type float float)
  2) 

(DEFUN float-digits (float)
  1"Returns the number of bits of fraction part FLOAT has, including hidden bit.
This depends only on the data type of FLOAT (short-float, single-float, or double-float)."*
  (ETYPECASE float
    (short-float 17.)
    (single-float 24.)
    (double-float 53.))) 

(DEFUN float-precision (FLOAT)
  1"Returns the number of significant bits used in the representation of FLOAT."*
  (IF (ZEROP float) 0
      (ETYPECASE float
	(short-float 17.)
	(single-float 24.)
	(double-float 53.))))

(DEFUN decode-float (float)
  1"Returns three values describing the significand, exponent, and sign of FLOAT.
The first is a number of the same type as FLOAT between .5 and 1 (but zero if FLOAT is zero).
This value, times 2 to a suitable power, equals FLOAT except in sign.
The second value is an integer, the exponent of 2 needed for that calculation.
The third value is a float whose sign and type match FLOAT's, and whose magnitude is 1."*
  (DECLARE (VALUES fraction-float exponent sign-float))
  (VALUES (ABS (FLOAT-FRACTION float))
	  (FLOAT-EXPONENT float)
	  (FLOAT-SIGN float))) 

(DEFUN integer-decode-float (float)
  1"Returns three values describing the fraction part, exponent, and sign of FLOAT.
The first is an integer representing the fraction part of FLOAT.
This value floated, times two to a suitable power, equals FLOAT except in sign.
The second value is an integer, the exponent of two needed for that calculation.
The third value is either 1 or -1."*
  (DECLARE (VALUES fraction-float exponent sign-float))
  (CHECK-TYPE float float)
  (VALUES (flonum-mantissa (ABS float))
	  (flonum-exponent float)
	  (IF (ZEROP
                (ETYPECASE float
                  (short-float (LDB (Byte 1 24.) (%POINTER float)))
                  (single-float (%P-LDB-OFFSET (Byte 1 31.) float 1))
                  (double-float (%P-LDB-OFFSET (Byte 1 31.) float 2))))
              1 -1) ))

(DEFUN flonum-exponent (float)
  1"Return the exponent of FLOAT to go with (FLONUM-MANTISSA FLOAT).
This exponent, if used to scale the integer which FLONUM-MANTISSA returns,
will produce the original argument.  This is not the same as FLOAT-EXPONENT!"*
  (ETYPECASE float
    (short-float (- (LDB (byte 8. 16.) (%POINTER float)) 143.))
    (single-float (- (%P-LDB-OFFSET (Byte 8. 23.) float 1) 150.))
    (double-float (- (%P-LDB-OFFSET (Byte 11. 20.) float 2) 1075.)))) 

(DEFUN flonum-mantissa (float)
  1"Return the mantissa of FLOAT as an integer."*
  (IF (ZEROP float) 0
      (LET ((fraction) (sign-bit))
	(ETYPECASE float
	  (short-float
	   (SETQ fraction (+ (LDB (BYTE 16. 0) (%POINTER float)) #x10000)
		 sign-bit (LDB (BYTE 1 24.) (%POINTER float))))
	  (single-float
	   (SETQ fraction (+ (%P-LDB-OFFSET (BYTE 23. 0) float 1) #x800000)
		 sign-bit (%P-LDB-OFFSET (BYTE 1 31.) float 1)))
	  (double-float
	   (SETQ fraction (+ (DPB (%P-LDB-OFFSET (BYTE 20. 0) float 2)
				  (BYTE 20. 32.)
				  (DPB (%P-LDB-OFFSET (BYTE 16. 16.) float 1)
				       (BYTE 16. 16.)
				       (%P-LDB-OFFSET (BYTE 16. 0) float 1)))
			     #x10000000000000)
		 sign-bit (%p-ldb-offset (byte 1 31.) float 2))))
	(IF (ZEROP sign-bit) fraction (- fraction)))))

;;CLM for Schuurman 12/10/87 - fixes problem reported in SPR 7020,
;;function was returning the wrong value if the second parameter was negative.
(DEFUN float-sign (float1 &optional float2)
  "Returns a float whose sign matches FLOAT1 and magnitude matches FLOAT2.
If FLOAT2 is omitted, it defaults to (float 1 FLOAT1).
The type of float returned matches the type of FLOAT2."
  (SETQ float2
	(IF float2 (ABS float2) (FLOAT 1 float1)))
  (FLOAT (IF (ZEROP
	       (ETYPECASE float1
		 (short-float (LDB (Byte 1 24.) (%POINTER float1)))
		 (single-float (%P-LDB-OFFSET (Byte 1 31.) float1 1))
		 (double-float (%P-LDB-OFFSET (Byte 1 31.) float1 2))))
	     float2 (- float2)))) 

(DEFSUBST short-float (number)
  1"Convert NUMBER to a short float."*
  (SMALL-FLOAT number)) 

(DEFSUBST oddp (number)
  1"T if NUMBER is odd."*
  (LOGTEST 1 number)) 

(DEFSUBST evenp (number)
  1"T if NUMBER is even."*
  (not (LOGTEST 1 number))) 

1;;; Integer square-root*

(DEFUN isqrt (n)
  "Square root of an integer, as an integer."
  (CHECK-TYPE n (integer 0 *))
  (IF (ZEROP n) 0					;Prevent problems with 0.
      (DO ((g (ASH 1 (ASH (1- (HAULONG n)) -1))
	      (+ g e))
	   (e))
	  ((ZEROP (SETQ e (TRUNCATE (- n (* g g))	;Endtest
				    (ASH g 1))))
	   (IF (> (* g g) n) (1- g) g)))))		;Result may be 1 too high

1;;; Floating point square root.  Standard range reduction, followed by Newton
;;; iterations (either 2 or 3).  Notice that the range may be expanded before the
;;; iterations, if we're dealing with an odd power of 2 for the exponent.  This
;;; was done before the iterations to keep the loss of significance to a minimum.*

(DEFUN sqrt (number)
  1"Principal square root of a number, returned as a float or complex value."*
  (ETYPECASE number
    (float
     (COND
       ((PLUSP number)
	(MULTIPLE-VALUE-BIND (f n) (DECODE-FLOAT number)
	  (LET ((y (+ 0.41731s0 (* 0.59016s0 f))))
	    (WHEN (ODDP n)
	      (SETQ y (* y (FLOAT 0.70710678118654752440d0 y))	;(sqrt .5)
		    n (1+ n)
		    f (SCALE-FLOAT f -1)))
	    (SETQ y (+ y (/ f y))
		  y (+ (SCALE-FLOAT y -2) (/ f y)))
	    (WHEN (TYPEP y 'double-float)
	      (SETQ y (ASH (+ y (/ f y)) -1)))		1;One more iteration*
	    (SCALE-FLOAT y (ASH n -1)))))
       ((MINUSP number)
	(COMPLEX (float 0 number) (SQRT (- number))))
       (t number)))					1;Must be zerop*
    (complex
     (LET ((a (complex-real-part number))
	   (b (complex-imag-part number))
	   (r (ABS number)))
       (COND ((PLUSP a)
	      (LET ((x (SQRT (ASH (+ r a) -1))))
		(WHEN (MINUSP b) (SETQ x (- x)))
		(COMPLEX x (/ b (ASH x 1)))))
	     ((MINUSP a)
	      (LET ((y (SQRT (ASH (- r a) -1))))
		(COMPLEX (/ b (ASH y 1)) y)))
	     (t
	      (LET ((y (SQRT (ASH r -1))))
		  (COMPLEX (IF (MINUSP b)
			       (- y)
			       y)
			   y))))
       ))
    (number (SQRT (float number)))))

(DEFUN log (number &optional *base)
  1"Return the logarithm of NUMBER in the base *BASE, which defaults to e,
the base of natural logarithms."*
  (IF *base
      (LET ((common-type (MIN (REALPART number)
			      (REALPART *base))))
	(/ (LOG (FLOAT number common-type)) (LOG (FLOAT *base common-type))))
      (ETYPECASE number
        (FLOAT
         (COND
	   ((PLUSP number)
	    (LET (a b c2)
	      (ETYPECASE number
		(single-float
		 (SETQ a '(-0.5527074855F0)
		       b '(-0.6632718214F1)
		       c2 -2.121944400546905827679f-4))
		(short-float
		 (SETQ a '(-0.5527074855S0)
		       b '(-0.6632718214S1)
		       c2 -2.121944400546905827679s-4))
		(double-float
		 (SETQ a '(-0.78956112887491257267d+0 0.16383943563021534222d+2
						      -0.64124943423745581147d+2)
		       b '(-0.35667977739034646171d+2 0.31203222091924532844d+3
						      -0.76949932108494879777d+3)
		       c2 -2.121944400546905827679d-4)))
	      
	      (MULTIPLE-VALUE-BIND (f n) (DECODE-FLOAT number)
		(WHEN (< f 0.70710678118654752440d0)
		  (SETF f (SCALE-FLOAT f 1))
		  (DECF n))
		(LET* ((number (SCALE-FLOAT (/ (1- f) (1+ f)) 1))
		       (w (sqr number)))
		  (IF (ZEROP n)
		      (+ number (* number (* w (/ (poly w a)(poly1 w b)))))
		      (+ (+ number (* number (* w (/ (poly w a)(poly1 w b)))))
			 (* n c2)
			 (SCALE-FLOAT (FLOAT (* n 355.) f) -9.)))))))
	   ((MINUSP number)
		 (COMPLEX (LOG (ABS number)) (FLOAT pi number)))
	   (t (ERROR 'LOG "Arg of log is zero"))))
        (COMPLEX (COMPLEX (LOG (ABS number)) (PHASE number)))
        (number (LOG (FLOAT number))))))

(DEFUN exp (number)
  1"Returns the exponential of NUMBER as a floating-point or complex number.  "*
  (ETYPECASE number
    (float
     (LET (xn c1 c2 p q)
       (ETYPECASE number				1;Setup the correct constants*
         (single-float
          (SETQ xn (ROUND (* number 1.4426950408889634074f0))
		c1 0.693359375f0
                c2 2.1219444005469058277f-4
                p '(0.41602886268f-2 0.24999999950f0)
                q '(0.49987178778f-1 0.5f0)))
         (short-float
          (SETQ xn (ROUND (* number 1.4426950408889634074s0))
		c1 0.693359375s0
		c2 2.1219444005469058277s-4
                p '(0.41602886268s-2 0.24999999950s0)
                q '(0.49987178778s-1 0.5s0)))
	 (double-float
	  (SETQ xn (ROUND (* number 1.4426950408889634074s0))
		c1 0.693359375d0
		c2 2.1219444005469058277d-4
		p '(0.165203300268279130d-4 0.694360001511792852d-2
					    0.249999999999999993d+0)
		q '(0.495862884905441294d-3 0.555538666969001188d-1
					    0.500000000000000000d+0))))
       (LET* ((g (MULTIPLE-VALUE-BIND (x1 x2) (FLOOR number)
		   (+ (- x1 (* xn c1)) x2 (* xn c2))))
	      (w (sqr g))
	      (gp (* g (poly w p)))
	      (q (poly w q)))
	 (SCALE-FLOAT (+ .5s0 (/ gp (- q gp))) (1+ xn)))))
    (complex (* (EXP (si:complex-real-part number))
                (CIS (si:complex-imag-part number))))
    (number (EXP (FLOAT number)))))

(DEFUN sin (radians)
  1"Sine of an angle measured in radians."*
  (ETYPECASE radians
    (float
     (SIN-COS-FLOAT-AUX radians))
    (complex
     (LET ((x (FLOAT (si:complex-real-part radians)))
           (y (FLOAT (si:complex-imag-part radians))))
       (COMPLEX
         (* (SIN x) (COSH y))
         (* (COS x) (SINH y)))))
    (number (SIN-COS-FLOAT-AUX (FLOAT radians)))))

(DEFUN cos (radians)					
  1"Cosine of an angle measured in radians."*
  (ETYPECASE radians
   (FLOAT
    (SIN-COS-FLOAT-AUX (ABS radians) t))
   (complex
    (LET ((x (FLOAT (si:complex-real-part radians)))
          (y (FLOAT (si:complex-imag-part radians))))
      (COMPLEX    
	(* (COS x) (COSH y))
	(- (* (SIN x) (SINH y))))))
   (number (SIN-COS-FLOAT-AUX (FLOAT (ABS radians)) t))))

;; 5/27/87 DNG for HS - changed 1/2 to (float 0.5s0 angle) for speed. [SPR 5515]
(DEFUN sin-cos-float-aux (angle &optional cosp)
  (LET (eps c1 c2 rc pi-inv)
    (ETYPECASE angle
      (single-float
       (SETQ eps 2.44f-4				1;~ 2**(-t/2) with t=24 bits*
	     c1	3.140625f0				1; 201/64*
	     c2 -9.676535898f-4
	     rc '(0.2601903036f-5 -0.1980741872f-3
				  0.8333025139f-2
				  -0.1666665668f+0)
	     pi-inv 0.31830988618379067154f0))
      (short-float
       (SETQ eps 0.0027s0				1;~ 2**(-t/2) with t=17 bits*
	     c1 3.140625s0
	     c2 -9.676535898s-4
	     rc '(0.2601903036S-5 -0.1980741872S-3
				  0.8333025139S-2
				  -0.1666665668S+0)
	     pi-inv 0.31830988618379067154s0))
      (double-float
       (SETQ eps 1.05d-8				1;~ 2**(-t/2) with t=53 bits*
	     c1 3.1416015625d0				1; 3217/1024*
	     c2 8.908910206761537356617d-6
	     rc '(0.27204790957888846175d-14 -0.76429178068910467734d-12
					     0.16058936490371589114d-9
					     -0.25052106798274584544d-7
					     0.27557319210152756119d-5
					     -0.19841269841201840457d-3
					     0.83333333333331650314d-2
					     -0.16666666666666665052d0)
	     pi-inv 0.31830988618379067154d0)))

    (multiple-value-bind (n r) (ROUND (* angle pi-inv))
      (when (and cosp (not (minusp r)))
	(incf n))
      (let* ((xn (if cosp (- n (float 0.5s0 angle)) n))
	     (f (MULTIPLE-VALUE-BIND (x1 x2) (FLOOR angle)
		  1;; The following preserves precision by using a multiple-step*
		  1;; argument reduction scheme.  Notice that the C2 we're using*
		  1;; is equivalent to -C2 as used in Cody and Waite.*
		  (+ (- x1 (* xn c1)) x2 (* xn c2))))
	     (res
	       (IF (< (ABS f) eps)			1; sin(f) == f to machine precision*
		   f
		   (LET* ((g (sqr f)))
		     (+ f (* f g (poly g rc)))))))
	(IF (ODDP n) (- res) res)))))

(DEFUN tan (radians)
  1"Tangent of RADIANS as a floating-point or complex number.  "*
  (ETYPECASE radians
    (float
     (LET (one 2OverPi c1 c2 eps p q)
       (ETYPECASE radians
         (single-float
          (SETQ one 1.0F0
		2OverPi 0.63661977236758134308f0
		c1 1.5703125f0
		c2 -4.838267948F-4
                eps 2.44f-4				1;~ 2**(-t/2) with t=24 bits*
                p '(-0.958017723f-1)
                q '(0.971685835f-2 -0.429135777f+0)))
         (short-float
          (SETQ one 1.0S0
		2OverPi 0.63661977236758134308f0
		c1 1.5703125f0
		c2 -4.838267948F-4
                eps 0.0027				1;~ 2**(-t/2) with t=17 bits*
                p '(-0.958017723S-1)
                q '(0.971685835S-2 -0.429135777S+0)))
	 (double-float
	  (SETQ one 1.0d0
		2OverPi 0.63661977236758134308d0
		c1 1.570800781
		c2 4.454455103380768678308d-6
		eps 1.05d-8				1;~ 2**(-t/2) with t=53 bits*
		p '(-0.17861707342254426711d-4 0.34248878235890589960d-2
					       -0.13338350006421960681d0)
		q '(0.49819433993786512270d-6 -0.31181531907010027307d-3
					      0.25663832289440112864d-1
					      -0.46671683339755294240d0))))

       (LET* ((n (ROUND (* radians 2OverPi)))
	      (f (MULTIPLE-VALUE-BIND (x1 x2) (FLOOR radians)
		   (+ (- x1 (* n c1)) x2 (* n c2))))
	      (res
		(IF (< (ABS f) eps)
		    one
		    (LET* ((g (sqr f)))
		      (INCF f (* f g (poly g p)))
		      (+ (* g (poly g q)) one)))))
	 (IF (ODDP n) (/ (- res) f) (/ f res)))))
    (complex
     (LET* ((x2 (SCALE-FLOAT (FLOAT (si:complex-real-part radians)) 1))
	    (y2 (SCALE-FLOAT (FLOAT (si:complex-imag-part radians)) 1)))
       (/ (COMPLEX (SIN x2) (SINH y2))
	  (+ (COS x2) (COSH y2)))))
    (number (TAN (FLOAT radians)))))

(DEFUN tand (number)
  1"Tangent of an angle measured in degrees."*
  (TAN (* number (/ (FLOAT pi number) 180.))))

(DEFUN cosd (number)
  1"Cosine of an angle measured in degrees."*
  (COS (* number (/ (FLOAT pi number) 180.))))

(DEFUN sind (number)
  1"Sine of an angle measured in degrees."*
  (SIN (* number (/ (FLOAT pi number) 180.))))


(DEFUN expt-hard (*base *exp)
  (COND ((ZEROP *exp)
	 (IF (AND (ZEROP *base) (NEQ *exp 0))
	     (FERROR nil "(^ ~S ~S) is undefined" *base *exp)
	     (1+ (* *exp *base))))		1; Coerce 1 to the stronger of *exp and *base*
	((ZEROP *base)
	 (IF (PLUSP (REALPART *exp))
	     (* *base *exp)			1; Coerce 0 to the stronger of *base and *exp*
	     (FERROR nil "(^ ~S ~S) is undefined" *base *exp)))	1;eg 0^-1 or 0^#c(0 3)*
	((INTEGERP *exp)
	 (IF (MINUSP *exp)
	     (/ 1 (slow-binary-expt *base (- *exp)))
	     (slow-binary-expt *base *exp)))
	(t					1; Perform proper coercion through MIN function*
	 (EXP (* *exp (LOG (FLOAT *base (MIN (REALPART *base) (REALPART *exp)))))))))

(DEFUN slow-binary-expt (*base *exp)
  (DO ((ans (IF (ODDP *exp) *base 1)
	    (IF (ODDP *exp) (* ans *base) ans)))
      ((ZEROP (SETQ *exp (ASH *exp -1))) ans)
    (SETQ *base (* *base *base))))		1;to avoid oflo, procrastinate squaring*
	

(DEFVAR *random-state* () 1"Default random number generator data"*) 

(DEFSTRUCT (random-state :named-array (:constructor make-random-state-1)
			 (:print-function
			   (lambda (random-state
				    stream
				    depth)
			     (LET ((*print-array* t))
			       (print-named-structure 'random-state random-state depth stream
						      (which-operations-for-print stream)))))
			 (:callable-constructors nil)
			 (:conc-name nil)
			 (:predicate nil)
			 (:copier nil))
  random-seed
  random-pointer-1
  random-pointer-2
  random-vector)


(DEFUN make-random-state (&optional state)
  1"Create a new random-state object for RANDOM to use.
If STATE is such a state object, it is copied.
If STATE is NIL or omitted, the default random-state is copied.
If STATE is T, a new state object is created and initialized based on the microsecond clock."*
  (COND ((NULL state)
	 (LET ((new (copy-object *random-state*)))
	   (SETF (random-vector new)
		 (copy-object (random-vector new)))
	   new))
	((EQ state t)
	 (random-create-array 71. 35. (time:fixnum-microsecond-time)))
	(t (LET ((new (copy-object state)))
	     (SETF (random-vector new)
		   (copy-object (random-vector new)))
	     new)))) 

(DEFUN random-create-array (size offset seed &optional (area nil))
  (LET ((default-cons-area (OR area default-cons-area)))
    (LET ((ARRAY (make-random-state-1
		   :random-vector (MAKE-ARRAY size)
		   :random-seed seed
		   :random-pointer-1 0
		   :random-pointer-2 offset)))
      (random-initialize array)
      array)))

(DEFUN random-initialize (ARRAY &optional new-seed &aux size x pointer)
  (IF (NOT (NULL new-seed))
      (SETF (random-seed array) new-seed))
  (SETQ size (LENGTH (random-vector array))
	pointer (AP-1 (random-vector array) 0))
  (SETF (random-pointer-2 array)
	(REM (+ size (- (random-pointer-2 array)
			(random-pointer-1 array)))
	     size))
  (SETF (random-pointer-1 array) 0)
  (ARRAY-INITIALIZE (random-vector array) 0)
  (SETQ x (random-seed array))
  (DOLIST (byte-spec
	    (CASE %%q-pointer
		  (24. '(#o1414 #o14))
		  (25. '(#o1414 #o14 #o3001))
		  (31. '(#o1414 #o14 #o3011))
		  (t (FERROR () "BUG IN random-initialize"))))
    (DO ((i 0 (1+ i)))
	((= i size))
      (SETQ x (* x 4093.))			1;4093. is a prime number.*
      (SETQ x (%LOGDPB (LDB #o3001 x) #o3001 (LDB #o0030 x)))	1;Simulate a %pointer-times operation*
      1(%P-STORE-POINTER-OFFSET*
	 1(%LOGDPB (LDB #o1314 x)*
		    1byte-spec*
		    1(%P-POINTER-OFFSET pointer i))*
	 1pointer i)*)))
;; TGC      (%P-DPB-OFFSET (LDB #o1314 x) byte-spec pointer i))))

(DEFUN random (&optional high array &aux ptr1 ptr2 size ans vector)
  1"Returns a randomly chosen number.
With no argument, value is chosen randomly from all fixnums.
If HIGH is an integer, the value is a nonnegative integer and less than HIGH.
If HIGH is a float, the value is a nonnegative number of the same type, and less than HIGH.
ARRAY can be an array used for data by the random number generator (and updated);
 you can create one with RANDOM-CREATE-ARRAY or MAKE-RANDOM-STATE."*
  ;1; 01/20/88 CLM - For floating point numbers, use FLOAT to prevent type conversion*
  ;1;                from rational to floating point.*
  ;1; 01/27/88 CLM - A slight variation on the above fix, provided by HSch.  Slightly*
  ;1;                faster and a little less consing.*
  (WHEN high
    (CHECK-TYPE high (AND real (satisfies plusp)) "a positive real number"))
  (COND ((NULL array)
	 (OR (AND (VARIABLE-BOUNDP *random-state*) *random-state*)
	     (SETQ *random-state* (random-create-array 71. 35. 69.)))
	 (SETQ array *random-state*)))		;Initialization as opt arg loses on BOUNDP.
  (WITHOUT-INTERRUPTS
   (SETQ ptr1 (random-pointer-1 array)
	 ptr2 (random-pointer-2 array)
	 vector (random-vector array)
	 size (LENGTH vector))
   (OR (< (SETQ ptr1 (1+ ptr1)) size) (SETQ ptr1 0))
   (OR (< (SETQ ptr2 (1+ ptr2)) size) (SETQ ptr2 0))
   (SETF (random-pointer-1 array) ptr1)
   (SETF (random-pointer-2 array) ptr2)
   (SETQ ans (%MAKE-POINTER-OFFSET dtp-fix (AREF vector ptr1) (AREF vector ptr2)))
   (SETF (AREF vector ptr2) ans))
  (COND ((FLOATP high)
	 (* (%LOGDPB 0 %%q-boxed-sign-bit ans)         1; Multiply by a ratio 0  x < 1*
	    (/ high
	       (- (%LOGDPB 1 %%q-boxed-sign-bit 0))))) 1; and coerce to proper float*
	((NULL high) ans)
	(t
	 (DO ((bits 14. (+ bits %%q-pointer))		1; Generate as many fixnums as needed*
	      (number (%LOGDPB 0 %%q-boxed-sign-bit ans) 1; Remove sign bit from ANS*
		      (+ (%LOGDPB 0 %%q-boxed-sign-bit (RANDOM))
			 (ASH ans (1- %%q-pointer)))))	1;   and put as many fixnums together*
	     ((> bits (HAULONG high))			1;   as needed to generate a long enough bignum*
	      (MOD number high))))))			1;     then generate the final result.

;; Return a randomly chosen number at least LOW and less than HIGH.*
(DEFUN random-in-range (low high)
  1"Randomly chosen float not less than LOW and less than high."*
  (PROG* ((r (RANDOM))
	  (rnorm (/ (LOGAND r #x3FFFF ) (FLOAT #x40000))))
    (RETURN (+ low (* rnorm (- high low)))))) 

1;Force *RANDOM-STATE* to get a value.*

(EVAL-WHEN (COMPILE) (RANDOM most-positive-fixnum)) 

1;A rational is of type DTP-EXTENDED-NUMBER and occupies three words.

;The second word is the numerator, and the third is the denominator.
;Rationals with denominator 0 or 1 are never created.

;A complex looks like a rational except it contains %HEADER-TYPE-COMPLEX.
;The second word is the real part and the third is the imaginary part.
;Complexes with imaginary part 0 are not normally created
;but are created by coercion in mixed-mode arithmetic and can be made with COMPLEX.*


(DEFSUBST rational-numerator (number)
  (%P-CONTENTS-OFFSET number 1)) 

(DEFSUBST rational-denominator (number)
  (%P-CONTENTS-OFFSET number 2)) 

(DEFSUBST complex-real-part (number)
  (%P-CONTENTS-OFFSET number 1)) 

(DEFSUBST complex-imag-part (number)
  (%P-CONTENTS-OFFSET number 2)) 

(DEFSUBST %complex-cons (realpart imagpart)
  (LET ((object
	  (%ALLOCATE-AND-INITIALIZE dtp-extended-number dtp-header
				    (DPB %header-type-complex
					 %%header-type-field
					 0)
				    0 number-cons-area 3)))
    (SETF (complex-real-part object) (+ realpart (* (DONT-OPTIMIZE 0) imagpart)))
    (SETF (complex-imag-part object) (+ imagpart (* (DONT-OPTIMIZE 0) realpart)))
    object))

(DEFUN %complex (realpart imagpart)
  (LET ((object
	  (%ALLOCATE-AND-INITIALIZE dtp-extended-number dtp-header
				    (DPB %header-type-complex
					 %%header-type-field
					 0)
				    0 number-cons-area 3)))
    (SETF (complex-real-part object) (+ realpart (* (DONT-OPTIMIZE 0) imagpart)))
    (SETF (complex-imag-part object) (+ imagpart (* (DONT-OPTIMIZE 0) realpart)))
    object))

(DEFUN complex (realpart &optional (imagpart 0))
  "Return a complex number with specified real and imaginary parts."
  (CHECK-TYPE realpart real)
  (CHECK-TYPE imagpart real)
  (IF (AND (ZEROP imagpart)
	   (RATIONALP realpart)
	   (RATIONALP imagpart))
      realpart
    (LET ((object
	    (%ALLOCATE-AND-INITIALIZE dtp-extended-number dtp-header
				      (DPB %header-type-complex
					   %%header-type-field
					   0)
				      0 number-cons-area 3)))
      (SETF (complex-real-part object) (+ realpart (* (DONT-OPTIMIZE 0) imagpart)))
      (SETF (complex-imag-part object) (+ imagpart (* (DONT-OPTIMIZE 0) realpart)))
      object)))

(DEFUN make-rational (numerator denominator)
  1"Return a rational number with specified numerator and denominator.
This can be used to construct rationals not in lowest terms,
but should not normally be so used."*
  (LET ((object
	  (%ALLOCATE-AND-INITIALIZE dtp-extended-number dtp-header
				    (DPB %header-type-rational
					 %%header-type-field
					 0)
				    0 number-cons-area 3)))
    (SETF (rational-numerator object) numerator)
    (SETF (rational-denominator object) denominator)
    object))

1;Slightly higher level creation and access functions.
;NORMALIZED-mumble makes a mumble in standard form
;and knows that some complex numbers are real, and some rationals are integers.*

(DEFSUBST normalized-rational (num den) (%div num den))

(defun numerator (rational)
  1"Return the numerator of RATIONAL.  On integers, this is the identity function."*
  (check-type rational rational)
  (if (integerp rational) rational
    (rational-numerator rational)))

(DEFUN denominator (RATIONAL)
  1"Return the denominator of RATIONAL.  On integers, this returns 1."*
  (CHECK-TYPE rational rational)
  (IF (INTEGERP rational) 1
    (rational-denominator rational)))

(DEFSUBST realpart (x)
  1"Return the real part of a complex number.  The real part of a real number is itself."*
  (CHECK-TYPE x number)
  (IF (COMPLEXP x)
      (complex-real-part x)
    x))

(DEFSUBST imagpart (x)
  1"Return the imaginary part of a complex number, or 0 if given a real number."*
  (CHECK-TYPE x number)				;FORD  8-30-84
  (IF (COMPLEXP x)
      (complex-imag-part x)
    (- x x)))

(DEFUN realp (x)
  1"T if X is a real number--any number that is not COMPLEXP."*
  (and (numberp x) (not (complexp x))))

;Higher level conversion functions

(DEFUN rational (x)
  1"Convert X to a rational number.
If X is a floating point number, it is regarded as completely exact."*
  (ETYPECASE x
    (float
     (IF (ZEROP x)
	 0
	 (MULTIPLE-VALUE-BIND (m e) (INTEGER-DECODE-FLOAT x)
	   (LET* ((z (1- (INTEGER-LENGTH (LOGAND m (- m)))))
		  (d (+ e z)))
	     (IF (MINUSP d)
		 (si:make-rational
		   (IF (MINUSP x)
		       (- (ASH m (- z)))
		       (ASH m (- z)))
		   (ASH 1 (- d)))
		 (ASH (IF (MINUSP x) (- m) m) e))))))
    (rational x)))

(DEFUN rationalize (x)
  1"Generates a good rational aproximation for X, assuming floats
are only accurate to their precision."*
  (TYPECASE x
    (rational x)
    (short-float (RATIONALIZE-FLOAT x short-float-epsilon))
    (single-float (RATIONALIZE-FLOAT x single-float-epsilon))
    (double-float (RATIONALIZE-FLOAT x double-float-epsilon))
    (otherwise (FERROR () "Cannot coerce ~S into a rational" x))))

(DEFUN rationalize-float (x &optional (eps long-float-epsilon))
  (COND ((MINUSP x) (- (RATIONALIZE-FLOAT (- x) eps)))
	((ZEROP x) 0)
	(t (LET ((y ())
		 (a ()))
	     (DO ((xx x y)
		  (num (SETQ a (TRUNCATE x)) (+ (* a num) onum))
		  (den 1 (+ (* a den) oden))
		  (onum 1 num)
		  (oden 0 den))
		 ((OR (= xx a)
		      (<= (ABS (/ (- x (/ num den))
				  x))
			  eps))
		  (/ num den))			1;Generate final result*
	       (SETQ y (/ (- xx a))
		     a (TRUNCATE y)))))))

(DEFUN conjugate (number)
  1"Return the complex conjugate of NUMBER.  If NUMBER is real, NUMBER is returned."*
  (IF (TYPEP number 'complex)
      (COMPLEX (complex-real-part number) (- (complex-imag-part number)))
    number))

(DEFUN phase (number)
  1"Return the phase of NUMBER, in radians.
This is the angle in the complex plane from the positive real axis
to the ray from the origin through NUMBER.
It is between - (exclusive) and  (inclusive).
For a positive real, this is 0; for a negative real, this is . For 0, it is zero."*
  (ETYPECASE number
    (float (IF (MINUSP number) (FLOAT pi number) (FLOAT 0 number)))
    (complex (IF (ZEROP number)			1;Avoid errors with #c(0.0 0.0)*
		 (complex-real-part number)	1;Take care of proper coercion*
		 (ATAN (complex-imag-part number) (complex-real-part number))))
    (number (IF (MINUSP number) (FLOAT pi) 0.0))))

(DEFUN cis (angle)
  1"Return the value of e**(ANGLE*i).  The inverse of the function PHASE."*
  (COMPLEX (cos angle) (sin angle)))
 
(DEFUN signum (z)
  1"Returns a number with the same phase as NUMBER but with unit magnitude."*
  (IF (ZEROP z) z
      (ETYPECASE z
        (float (FLOAT (IF (MINUSP z) -1.0f0 1.0f0) z))
        (complex (/ z (ABS z)))
        (number (IF (MINUSP z) -1 1)))))



(DEFUN print-rational (number stream ignore)
  (AND *print-radix*
       (NOT (pttbl-rational-radix *readtable*))
       (NUMBERP *print-base*)
       (COND ((EQ *print-base* 8.)
	      (SEND stream :string-out "#o"))
	     ((EQ *print-base* 2)
	      (SEND stream :string-out "#b"))
	     ((EQ *print-base* 16.)
	      (SEND stream :string-out "#x"))
	     (t
	      (SEND stream :tyo #\#)
	      (LET ((*print-base* 10.) (*print-radix* nil) (*nopoint t) (tem *print-base*))
		(print-fixnum tem stream))
	      (SEND stream :tyo #\r))))       
  (LET ((*print-base* (OR (pttbl-rational-radix *readtable*) *print-base*))
	(*nopoint t)
	(*print-radix* nil))
    (print-fixnum (rational-numerator number) stream)
    (FUNCALL stream :tyo (pttbl-rational-infix *readtable*))
    (print-fixnum (rational-denominator number) stream)))

(DEFUN print-complex (cplx stream ignore)
  (SEND stream :string-out (FIRST (pttbl-complex *readtable*)))
  (PRINC (complex-real-part cplx) stream)
  (IF (SECOND (pttbl-complex *readtable*))
      (SEND stream :string-out (SECOND (pttbl-complex *readtable*)))
    (UNLESS (MINUSP (complex-imag-part cplx))
      (SEND stream :tyo #\+)))
  (PRINC (complex-imag-part cplx) stream)
  (SEND stream :string-out (THIRD (pttbl-complex *readtable*))))

;;;(DEFUN (:property rational standard-read-function) (STREAM string &aux num i (*read-base* 10.)
;;;						    (len (LENGTH string)))
;;;  stream
;;;  (MULTIPLE-VALUE-SETQ (num i)
;;;		       (xr-read-fixnum-internal string 0 len))
;;;  (VALUES
;;;    (normalized-rational
;;;      num
;;;      (xr-read-fixnum-internal string (1+ i) len))
;;;    'RATIONAL))

;;;(DEFUN (:property cl-rational standard-read-function) (STREAM string &aux num i
;;;						       (len (LENGTH string)))
;;;  stream
;;;  (MULTIPLE-VALUE-SETQ (num i)
;;;		       (xr-read-fixnum-internal string 0 len))
;;;  (VALUES
;;;    (normalized-rational
;;;      num
;;;      (xr-read-fixnum-internal string (1+ i) len))
;;;    'RATIONAL))

;;;(DEFUN (:property complex standard-read-function) (STREAM string &aux complex-start (zero 0))
;;;  stream
;;;  (DO ((i 1 (1+ i)))
;;;      ((= i (LENGTH string)))
;;;    (WHEN (AND (MEMBER (AREF string i) '(#\+ #\-) :test '=)
;;;	       (NOT (ALPHA-CHAR-P (AREF string (1- i)))))
;;;      (RETURN (SETQ complex-start i))))
;;;  (VALUES
;;;    (COMPLEX
;;;      (COND (complex-start (WITH-INPUT-FROM-STRING (strm string zero complex-start)
;;;			     (read-preserving-whitespace  strm t nil t)))
;;;	    (t (SETQ complex-start 0)))
;;;      (WITH-INPUT-FROM-STRING (strm string complex-start (1- (LENGTH string)))
;;;	(read-preserving-whitespace  strm t nil t)))
;;;    'COMPLEX))

1;Standard arithmetic functions.*

(DEFUN numeric-one-argument (code number)
  (COND ((RATIONALP number)
	 (LET ((num (rational-numerator number))
	       (den (rational-denominator number)))
	   (CASE (LOGAND #x3F code)
		 (0 (IF (OR (= num 0)			1;ABS*
			    (AND (> num 0) (> den 0))
			    (AND (< den 0) (< num 0)))
			number
			(make-rational (ABS num) (ABS den))))
		 (1 (make-rational (- num) den))	1;MINUS*
		 (2 (= num 0))			;ZEROP
		 (3 (OR (AND (> num 0) (> den 0))	1;PLUSP*
			(AND (< num 0) (< den 0))))
		 (4 (OR (AND (> num 0) (< den 0))	1;MINUSP*
			(AND (< num 0) (> den 0))))
		 (5 (make-rational (+ num den) den))	1;ADD1*
		 (6 (make-rational (- num den) den))	1;SUB1*
		 (7					1;FIX*
		  (CASE (LDB #o0603 code)
			(0 (IF (PLUSP num) (TRUNCATE num den)
			       (TRUNCATE (- num den -1) den)))
			(1 (IF (MINUSP num) (TRUNCATE num den)
			       (TRUNCATE (+ num den -1) den)))
			(2 (TRUNCATE num den))
			(3
			 (LET* ((FLOOR (IF (PLUSP num) (TRUNCATE num den)
					   (TRUNCATE (- num den -1) den)))
				(fraction-num (- num (* floor den)))
				(half-indicator (- (+ fraction-num fraction-num) den)))
			   (IF (OR (PLUSP half-indicator)
				   (AND (ZEROP half-indicator)
					(ODDP floor)))
			       (INCF floor))
			   floor))))
		 ((#o10 #o11)				1;Single-float and Short-float*
		  (LET* ((other (IF (ODDP code) 1.0s0 1.0f0))
			 (iln (INTEGER-LENGTH num))
			 (ild (INTEGER-LENGTH den))
			 (fdo (FLOAT-DIGITS other))
			 (fld (MIN (+ fdo 2) ild))	1;Provide 2 additional guard bits*
			 (fln (+ fld fdo)))
		    (WHEN (> (ASH num (- ild iln)) den)
		      (INCF iln))
		    (SCALE-FLOAT
		      (FLOAT (ROUND (ASH num (- fln iln)) (ASH den (- fld ild))) other)
		      (- iln ild fdo))))
		 (#o12 (FERROR nil "Illegal operation \"HAULONG\" on ~s" number))
		 (#o13 (FERROR nil "Illegal operation \"LDB\" on ~s" number))
		 (#o14 (FERROR nil "Illegal operation \"DPB\" on ~s" number))
		 (#o15 (FERROR nil "Illegal operation \"ASH\" on ~s" number))
		 (#o16					1;Double-float*
		  (LET* ((iln (INTEGER-LENGTH num))
			 (ild (INTEGER-LENGTH den))
;1;;*			 (fdo (FLOAT-DIGITS 1.0d0))	1;Only do this natively!!!*
			 (fdo 53.)
			 (fld (MIN (+ fdo 2) ild))	1;Provide 2 additional guard bits*
			 (fln (+ fld fdo)))
		    (WHEN (> (ASH num (- ild iln)) den)
		      (INCF iln))
		    (SCALE-FLOAT
		      (double-float (ROUND (ASH num (- fln iln)) (ASH den (- fld ild))))
		      (- iln ild fdo))))
		 (t (FERROR nil "Arith one-arg op code ~S on ~s" code number)))))
	((COMPLEXP number)
	 (LET ((real (complex-real-part number))
	       (imag (complex-imag-part number)))
	   (CASE code
		 (0 (LET ((y (ABS real))
			  (x (ABS imag)))
		      (WHEN (< y x) (PSETQ y x x y))
		      (IF (ZEROP x)			1;X is always the smaller of X and Y*
			  y				1;Already properly coerced*
			  (* y (SQRT (1+ (sqr (/ x y))))))))	1;ABS*
		 (1 (COMPLEX (- real) (- imag)))	1;MINUS*
		 (2 (AND (ZEROP real) (ZEROP imag)))	1;ZEROP*
		 (3 (FERROR nil "PLUSP applied to the complex number ~s" number))	1;PLUSP*
		 (4 (FERROR nil "MINUSP applied to the complex number ~s" number))	1;MINUSP*
		 (5 (COMPLEX (+ real 1) imag))		1;ADD1*
		 (6 (COMPLEX (- real 1) imag))		1;SUB1*
		 (7 (COMPLEX (FLOOR real) (FLOOR imag)))	1;FIX*
		 (#o10 (COMPLEX				1;Single-float*
			 (internal-float real)
			 (internal-float imag)))
		 (#o11 (COMPLEX				1;Short-float*
			 (SMALL-FLOAT real)
			 (SMALL-FLOAT imag)))
		 (#o12 (FERROR nil "Illegal operation \"HAULONG\" on ~s" number))
		 (#o13 (FERROR nil "Illegal operation \"LDB\" on ~s" number))
		 (#o14 (FERROR nil "Illegal operation \"DPB\" on ~s" number))
		 (#o15 (FERROR nil "Illegal operation \"ASH\" on ~s" number))
		 (#o16 (COMPLEX				1;Double-float*
			 (double-float real)
			 (double-float imag)))
		 (t (FERROR nil "Arith one-arg op code ~S on ~s" code number)))))
	(t (FERROR nil "Trap to macrocode for arithmetic on number ~s" number))))

1;;;
;;; Operations on 2 arguments of the type integer or float are completely
;;; handled in microcode.  Only if one of the operands is a complex or a
;;; ratio will we wind up here...
;;;*

(DEFUN numeric-two-arguments (code number1 number2 &aux function)
  (SETQ function (NTH code '(+ - * global:/ = > < min max boole /)))
  (COND ((AND (COMPLEXP number1) (COMPLEXP number2))
	 (complex-two-arguments code number1 number2))
	((COMPLEXP number1)
	 (complex-two-arguments code number1 (%complex number2 0)))
	((COMPLEXP number2)
	 (complex-two-arguments code (%complex number1 0) number2))
	((FLOATP number1)
	 (FUNCALL function number1 (FLOAT number2 number1)))
	((FLOATP number2)
	 (FUNCALL function (FLOAT number1 number2) number2))
	((AND (RATIONALP number1) (RATIONALP number2))
	 (rational-two-arguments code number1 number2))
	(t
	 (FERROR nil "Arith two-arg op code ~S on ~S and ~s" code number1 number2))))

(DEFUN rational-two-arguments (code number1 number2)
  (LET (num1 den1 num2 den2)
    (IF (INTEGERP number1) (SETQ num1 number1 den1 1)
	(SETQ num1 (rational-numerator number1)
	      den1 (rational-denominator number1)))
    (IF (INTEGERP number2) (SETQ num2 number2 den2 1)
	(SETQ num2 (rational-numerator number2)
	      den2 (rational-denominator number2)))
    (CASE code
	  (0 (normalized-rational (+ (* num1 den2) (* num2 den1)) (* den1 den2)))	1;ADD*
	  (1 (normalized-rational (- (* num1 den2) (* num2 den1)) (* den1 den2)))	1;SUB*
	  (2 (normalized-rational (* num1 num2) (* den1 den2)))	1;MUL*
	  ((3 #o12)
	   (normalized-rational (* num1 den2) (* den1 num2)))	1;DIV and IDIV*
	  ((4 #o13)
	   (AND (= num1 num2) (OR (ZEROP num1) (= den1 den2))))	1;= and EQUAL/EQL*
	  (5 (> (* num1 den2) (* num2 den1)))		1;GREATERP*
	  (6 (< (* num1 den2) (* num2 den1)))		1;LESSP*
	  (7 (IF (> (* num1 den2) (* num2 den1)) number2 number1))	1;MIN*
	  (#o10 (IF (> (* num1 den2) (* num2 den1)) number1 number2))	1;MAX*
	  (#o11 (FERROR nil "BOOLE applied to rational numbers ~S and ~s." number1 number2))	1;BOOLE*
	  (otherwise (FERROR nil "Rational two arg op code ~S on ~S and ~s."
			     code number1 number2)))))	1;Unused codes.*

(DEFUN complex-two-arguments (code number1 number2)
  (LET ((real1 (complex-real-part number1))
	(imag1 (complex-imag-part number1))
	(real2 (complex-real-part number2))
	(imag2 (complex-imag-part number2)))
    (CASE code
	  (0 (COMPLEX (+ real1 real2) (+ imag1 imag2)))	1;ADD*
	  (1 (COMPLEX (- real1 real2) (- imag1 imag2)))	1;SUB*
	  (2 (COMPLEX (- (* real1 real2) (* imag1 imag2))
		      (+ (* real1 imag2) (* imag1 real2)))) 1;MUL*
	  ((3 #o12)
	   (LET ((norm2 (+ (* real2 real2) (* imag2 imag2))))	1;DIV and IDIV*
	     (COMPLEX (%DIV (+ (* real1 real2) (* imag1 imag2)) norm2)
		      (%DIV (- (* imag1 real2) (* real1 imag2)) norm2))))
	  (4 (AND (= real1 real2) (= imag1 imag2)))	1;=*
	  (5 (FERROR nil "GREATERP applied to complex numbers ~S and ~s." number1 number2))	1;GREATERP*
	  (6 (FERROR nil "LESSP applied to complex numbers ~S and ~s." number1 number2))	1;LESSP*
	  (7 (FERROR nil "MIN applied to complex numbers ~S and ~s." number1 number2))	1;MIN*
	  (#o10 (FERROR nil "MAX applied to complex numbers ~S and ~s." number1 number2))	1;MAX*
	  (#o11 (FERROR nil "BOOLE applied to complex numbers ~S and ~s." number1 number2))	1;BOOLE*
	  (#o13 (AND (EQL real1 real2) (EQL imag1 imag2)))	1;EQUAL/EQL*
	  (otherwise (FERROR nil "Rational two arg op code ~S on ~S and ~s."
			     code number1 number2)))))	1;Unused codes.*

(DEFUN asin (number)
  1"Returns the arc sine of NUMBER, as a floating-point or complex number.  "*
  (ETYPECASE number
    (float
     (LET (local-half-pi eps p q)
       (ETYPECASE number
         (single-float
          (SETQ local-half-pi 1.57079632679489661923f0
                eps 2.44f-4				1;~2**(-24/2)*
                p '(-0.504400557f+0 0.933935835f+0)
                q '(-0.554846723f+1 0.560363004f+1)))
         (short-float
          (SETQ local-half-pi 1.57079632679489661923s0
                eps 0.0027s0				1;~2**(-17/2)*
                p '(-0.504400557S+0 0.933935835S+0)
                q '(-0.554846723S+1 0.560363004S+1)))
	 (double-float
	  (SETQ local-half-pi 1.57079632679489661923d0
		eps 1.05d-8				1;~2**(-53/2)*
		p '(-0.69674573447350646411d+0 0.10152522233806463645d+2
					       -0.39688862997504877339d+2
					       0.57208227877891731407d+2
					       -0.27368494524164255994d+2)
		q '(-0.23823859153670238830d+2 0.15095270841030604719d+3
					       -0.38186303361750149284d+3
					       0.41714430248260412556d+3
					       -0.16421096714498560795d+3))))
       (IF (> number 1.0s0)
           (COMPLEX local-half-pi (ACOSH number))
           (IF (< number -1.0s0)
               (COMPLEX (- local-half-pi) (ACOSH (- number)))
               (LET* ((y (ABS number))
                      g
                      (i nil)
                      (res
                        (COND ((< y eps) y)
                              (t
                               (IF (> y 0.5s0)
				   (SETQ i t
					 g (SCALE-FLOAT (- 1.0s0 y) -1)
					 y (SCALE-FLOAT (- (SQRT g)) 1))
				   (SETQ g (sqr y)))
                               (+ y (* y (/ (* g (poly g p))(poly1 g q))))))))
                 (WHEN i (INCF res local-half-pi))
                 (IF (MINUSP number) (- res) res))))))
    (COMPLEX
     (LET* ((x (si:complex-real-part number))
            (y (si:complex-imag-part number))
            (alpha (SCALE-FLOAT (+ (SQRT (+ (sqr (1+ x)) (sqr y)))
                                   (SQRT (+ (sqr (1- x)) (sqr y)))) -1))
            (beta (ACOSH alpha)))
       (COMPLEX (ASIN (/ x alpha))
                (IF (MINUSP y) (- beta) beta))))
    (number (ASIN (FLOAT number)))))

(DEFUN acos (number)
  1"number is a number.  Returns the angle whose cosine is number as a floating-point 
  or complex number.  "*
  (ETYPECASE number                                        
    (float
     (LET (local-pi local-half-pi)
       (ETYPECASE
         number
         (single-float
          (SETQ local-pi 3.141592653589793238462643f0
                local-half-pi 1.57079632679489661923f0))
         (short-float
          (SETQ local-pi 3.141592653589793238462643s0
                local-half-pi 1.57079632679489661923s0))
	 (double-float
	  (SETQ local-pi 3.141592653589793238462643d0
		local-half-pi 1.57079632679489661923d0)))
       (IF (> number 1.0s0)
           (COMPLEX 0 (ACOSH number))
           (IF (< number -1.0s0)
               (COMPLEX local-pi (- (ACOSH (- number))))
               (- local-half-pi (ASIN number))))))      
    (complex
     (LET* ((x (si:complex-real-part number))
            (y (si:complex-imag-part number))
            (alpha (SCALE-FLOAT (+ (SQRT (+ (sqr (1+ x)) (sqr y)))
                                   (SQRT (+ (sqr (1- x)) (sqr y)))) -1))
            (beta (ACOSH alpha)))
       (COMPLEX (ACOS (/ x alpha))
                (IF (MINUSP y) beta (- beta)))))
    (number (ACOS (FLOAT number)))))

(DEFUN atan (y &optional (x (FLOAT 1 y)))
  1"Arctangent in radians of y/x, between - and ."*
  (ETYPECASE y
    (float
     (ETYPECASE x
       (float)
       (rational (SETQ x (FLOAT x y))))
     (LET* ((n 0)
	    (f (COND
		 ((> (ABS y) (ABS x)) (SETQ n 2) (ABS (/ x y)))
		 (t (ABS (/ y x)))))
	    local-pi local-sqrt3 eps p q a-constants)
         (ETYPECASE f
           (single-float
            (SETQ local-pi 3.141592653589793238462643f0
		  local-sqrt3 1.73205080756887729353f0
                  eps 2.44f-4
                  p '(-0.5090958253F-1 -0.4708325141F+0)
                  q '(0.1412500740F+1)
                  a-constants '(0.0f0
			 0.52359877559829887308f0
			 1.57079632679489661923f0
			 1.04719755119659774615f0)))
           (short-float
            (SETQ local-pi 3.141592653589793238462643s0
		  local-sqrt3 1.73205080756887729353s0
                  eps 0.0027s0
                  p '(-0.5090958253S-1 -0.4708325141S+0)
                  q '(0.1412500740S+1)
                  a-constants '(0.0s0
			 0.52359877559829887308s0
			 1.57079632679489661923s0
			 1.04719755119659774615s0)))
	   (double-float
	    (SETQ local-pi 3.141592653589793238462643d0
		  local-sqrt3 1.73205080756887729353d0
                  eps 1.05e-8
		  p '(-0.83758299368150059274d+0 -0.84946240351320683534d+1
						 -0.20505855195861651981d+2
						 -0.13688768894191926929d+2)
		  q '(0.15024001160028576121d+2 0.59578436142597344465d+2
						0.86157349597130242515d+2
						0.41066306682575781263d+2)
		  a-constants '(0.0d0
			 0.52359877559829887308d0
			 1.57079632679489661923d0
			 1.04719755119659774615d0))))
         (WHEN (> f 0.26794919243112270647d0)		; (- 2.0 (sqrt 3.0))
           (SETQ f (/ (1- (* local-sqrt3 f)) (+ local-sqrt3 f))
                 n (1+ n)))
         (WHEN (>= (ABS f) eps)
           (LET* ((g (sqr f))
                  (r (/ (* g (poly g p)) (poly1 g q))))
             (INCF f (* f r))))
	 (WHEN (> n 1) (SETQ f (- f)))
	 (INCF f (NTH n a-constants))
	 (WHEN (MINUSP x) (SETQ f (- local-pi f)))
	 (IF (MINUSP y) (- f) f)))
    (complex
     (WHEN x (ERROR 'ATAN "When the first arg is complex, the second must be null"))
     (LET ((r (FLOAT (si:complex-real-part y)))
           (i (FLOAT (si:complex-imag-part y))))
       (COMPLEX (SCALE-FLOAT
                  (ATAN (SCALE-FLOAT r 1) (- 1 (sqr r) (sqr i)))
                  -1)
                (SCALE-FLOAT
                  (LOG (/ (+ (sqr r) (sqr (1+ i)))
                          (+ (sqr r) (sqr (1- i)))))
                  -2))))
    (number (ATAN (FLOAT y) x))))

(DEFF global:atan2 'ATAN) 

(DEFUN global:atan (y &optional x)
  1"Arctangent in radians of y/x, between 0 and 2."*
  (LET ((result (ATAN y x)))
    (IF (MINUSP result)
	(+ result (* 2 pi))
	result)))

(DEFUN sinh (number)
  1"Returns the hyperbolic sine of NUMBER as a floating-point or complex number."*
  (ETYPECASE number
    (float
     (LET ((y (ABS number))
	   max-number ln-v half-v-1 eps p q)
       (ETYPECASE
         number
         (single-float
	  (SETQ max-number 88.72283906f0
                ln-v 0.69316101074218750000f+0
                half-v-1 0.13830277879601902638F-4
                eps 2.2f-5
                p '(-0.190333399F+0 -0.713793159F+1)
                q '(-0.428277109F+2)))
         (short-float
          (SETQ max-number 88.722815s0
                ln-v 0.69316101074218750000s+0
                half-v-1 0.13830277879601902638S-4
                eps 2.8s-3
                p '(-0.190333399S+0 -0.713793159S+1)
                q '(-0.428277109S+2)))
	 (double-float
	  (SETQ max-number 709.782712893384d+0
		ln-v 0.69316101074218750000d+0
                half-v-1 0.13830277879601902638d-4
		eps 1.05e-8
		p '(-0.78966127417357099479d+0 -0.16375798202630751372d+3
					       -0.11563521196851768270d+5
					       -0.35181283430177117881d+6)
		q '(-0.27773523119650701667d+3 0.36162723109421836460d+5
					       -0.21108770058106271242d+7))))
       (COND
         ((> y 1.0s0)
          (IF (> y max-number)
              (LET ((w (EXP (- y ln-v))))
                (SETQ y (+ w (* w half-v-1))))
              (LET ((w (EXP y)))
                (SETQ y (/ (- w (/ w)) 2))))
          (IF (MINUSP number) (- y) y))
         (t
          (IF (< y eps)
              number
              (LET ((f (sqr number)))
                (+ (* (/ (poly f p) (poly1 f q)) f number) number)))))))
    (complex
     (LET ((x (FLOAT (si:complex-real-part number)))
           (y (FLOAT (si:complex-imag-part number))))
       (COMPLEX
         (* (SINH x) (COS y))
         (* (COSH x) (SIN y)))))
    (number (SINH (FLOAT number)))))

(DEFUN cosh (number)
  1"Returns the hyperbolic cosine of number as a floating-point or complex number."*
  (ETYPECASE number
    (float
     (LET ((y (ABS number)) max-number ln-v half-v-1)
       (ETYPECASE number
         (single-float
          (SETQ max-number 88.72283906f0
                ln-v 0.69316101074218750000f+0
                half-v-1 0.13830277879601902638f-4))
         (short-float
          (SETQ max-number 88.722815s0
                ln-v 0.69316101074218750000s+0
                half-v-1 0.13830277879601902638s-4))
	 (double-float
	  (SETQ max-number 709.782712893384d+0
		ln-v 0.69316101074218750000d+0
		half-v-1 0.13830277879601902638d-4)))
       (IF (> y max-number)
           (LET ((w (EXP (- y ln-v))))
             (+ w (* w  half-v-1)))
           (LET ((w (EXP y)))
             (/ (+ w (/ w)) 2)))))
    (complex
     (LET ((x (FLOAT (si:complex-real-part number)))
           (y (FLOAT (si:complex-imag-part number))))
       (COMPLEX
         (* (COSH x) (COS y))
         (* (SINH x) (SIN y)))))
    (number (COSH (FLOAT number)))))

(DEFUN tanh (number)
  1"Returns the hyperbolic tangent of number as a floating-point or complex number."*
  (ETYPECASE number
    (float
     (LET ((f (ABS number)) xbig half-ln3 eps p q)
       (ETYPECASE number
         (single-float
          (SETQ xbig 9.010914f0
                half-ln3 0.54930614433405484570f0
                eps 2.44e-4
                p '(-0.3831010665F-2 -0.8237728127F+0)
                q '(0.2471319654F+1)))
         (short-float
          (SETQ xbig 6.5849s0
                half-ln3 0.54930614433405484570s0
                eps 0.0027s0
                p '( -0.3831010665S-2 -0.8237728127S+0)
                q '( 0.2471319654S+1)))
	 (double-float
	  (SETQ xbig 19.061547465398494d+0
		half-ln3 0.54930614433405484570d0
                eps 1.05e-8
                p '(-0.96437492777225469787d+0 -0.99225929672236083313d+2
					       -0.16134119023996228053d+4)
		q '(0.11274474380534949335d+3 0.22337720718962312926d+4
					      0.48402357071988688686d+4))))
       (SETQ f
             (COND
               ((> f xbig) (FLOAT 1.0s0 number))
               ((> f half-ln3)
                (SCALE-FLOAT (- 0.5s0 (/ (1+ (EXP (SCALE-FLOAT f 1))))) 1))
               ((< f eps) f)
               (t
                (LET ((g (sqr f)))
                  (+ (* (/ (poly g p)(poly1 g q)) g f) f)))))
       (IF (MINUSP number) (- f) f)))
    (complex
     (LET ((x2 (SCALE-FLOAT (FLOAT (si:complex-real-part number)) 1))
           (y2 (SCALE-FLOAT (FLOAT (si:complex-imag-part number)) 1)))
       (/ (COMPLEX (SINH x2) (SIN y2))
          (+ (COSH x2) (COS y2)))))
    (number (TANH (FLOAT number)))))

(DEFUN asinh (number)
  1"Returns the hyperbolic sine of NUMBER as a floating-point or complex number."*
  (ETYPECASE number
    (float
     (LET* ((w (ABS number))
            (y (LOG (+ w (SQRT (1+ (sqr w)))))))
       (IF (MINUSP number) (- y) y)))
    (complex
     (LET* ((x (si:complex-real-part number))
            (y (si:complex-imag-part number))
            (w (SCALE-FLOAT (+ (SQRT (+ (sqr (1+ y)) (sqr x)))
                               (SQRT (+ (sqr (1- y)) (sqr x)))) -1))
            (ch (ACOSH w)))
       (COMPLEX
         (IF (MINUSP x) (- ch) ch)
         (ASIN (/ y w)))))
    (number (ASINH (FLOAT number)))))

(DEFUN acosh (number)
  1"Returns the hyperbolic cosine of NUMBER as a floating-point or complex number."*
  (ETYPECASE number
    (float
     (IF (>= number 1.0s0)
         (LOG (+ (SQRT (1- (sqr number))) number))
         (IF (<= number -1.0s0)
             (COMPLEX (LOG (- (SQRT (1- (sqr number))) number)) (FLOAT pi number))
             (COMPLEX 0 (ACOS number)))))
    (complex
     (LET* ((x (si:complex-real-part number))
            (y (si:complex-imag-part number))
            (w (SCALE-FLOAT (+ (SQRT (+ (sqr (1+ x)) (sqr y)))
                               (SQRT (+ (sqr (1- x)) (sqr y)))) -1))
            (cs (ACOS (/ x w))))
       (COMPLEX
         (ACOSH w)
         (IF (MINUSP y) (- cs) cs))))
    (number (ACOSH (FLOAT number)))))

(DEFUN atanh (number)
  1"Returns the hyperbolic tangent of NUMBER as a floating-point or complex number."*
  (ETYPECASE number
    (float
     (/ (LOG (/ (+ 1 number) (- 1 number))) 2))
    (complex
     (LET* ((x (si:complex-real-part number))
	    (y (si:complex-imag-part number)))
       (COMPLEX
	 (SCALE-FLOAT (LOG (/ (+ (sqr (1+ x)) (sqr y))
			      (+ (sqr (1- x)) (sqr y))))
		      -2)
	 (SCALE-FLOAT (ATAN (* 2 y) (- 1 (sqr x) (sqr y)))
		      -1))))
    (number (ATANH (FLOAT number)))))

;;PHD New function for new CLtL
(defun signed-ldb (ppss n)
  "Same as ldb but sign extends the extracted byte"
  (let ((extracted-byte (ldb ppss n)))
    (if (ldb-test (byte 1 (1- (byte-size ppss))) extracted-byte)
	(dpb extracted-byte
	     (byte (byte-size ppss) 0)
	     -1)
	extracted-byte)))

(setf (get 'signed-ldb 'setf-method) (get 'ldb 'setf-method))