;;; -*- Mode:Lisp; Readtable:ZL; Package:USER; Base:10; Patch-File:T -*-
;;; Patch file for System version 126.125
;;; Reason:
;;;  Portability improvements to NUMER.LISP (numerical functions), particular
;;;  proper byte-spec manipulations.
;;; Written 22-Oct-88 04:17:15 by keith (Keith Corbett) at site Gigamos Cambridge
;;; while running on Johannes Brahms from band 1
;;; with Experimental System 126.123, Experimental ZWEI 126.21, Experimental ZMail 74.9, Experimental Local-File 76.0, Experimental File-Server 25.0, Experimental Lambda-Diag 18.0, Experimental Unix-Interface 15.0, Experimental Tape 26.4, Microcode 1762, SDU Boot Tape 3.14, SDU ROM 103, 10/17.



; From modified file DJ: L.SYS2; NUMER.LISP#83 at 22-Oct-88 04:20:43
#10R SYSTEM-INTERNALS#: 
(COMPILER-LET ((*PACKAGE* (PKG-FIND-PACKAGE "SYSTEM-INTERNALS")))
  (COMPILER::PATCH-SOURCE-FILE "SYS: SYS2; NUMER  "


(defun sqrt (number)
  "Square root of a number, as a float or complex.
Result is a short-float or complex short-float according to type of NUMBER"
  (let* ((n (if (complexp number)
		(if (zerop (%complex-imag-part number))
		    (float (%complex-real-part number))
		  (%complex-cons (float (%complex-real-part number))
				 (float (%complex-imag-part number))))
	      (float number)))
	 (val
	   (cond ((complexp n)
		  (let* ((abs (abs n))
			 (real (%complex-real-part n))
			 (imag (%complex-imag-part n))
			 (r (sqrt (// (+ real abs) 2))))
		    (%complex-cons r (// imag (+ r r)))))
		 ((< n 0.0)
		  (%complex-cons (- n n) (sqrt (- n))))
		 ((= n 0.0)
		  0.0)
		 (t
		  (let ((f)
			(i2)
			(exp (- (%single-float-exponent n) single-float-exponent-offset
				-2)))
		    (let ((number-cons-area working-storage-area))
		      ;; F and I2 need to be regular-heap-consed to avoid
		      ;; the extra-pdl lossage. Identity switch on stack-group switch.
		      (setq f (+ n 0.0f0)
			    i2 (%float-double 0 1)))
                    (setf (%single-float-exponent f) single-float-exponent-offset)
                    (setf (%single-float-exponent i2)                            
                          (+ single-float-exponent-offset
			     (if (oddp exp)
				 (1+ (dpb (ldb (BYTE 23. 1.) exp) (BYTE 23. 0.) exp))
			         (dpb (ldb (BYTE 23. 1.) exp) (BYTE 23. 0.) exp))))
		    (do ((i 0 (1+ i))
			 (an (* i2 (+ 0.4826004 f (if (oddp exp) -0.25 0.0)))))
			((= i 4) an)
		      (setq an (* 0.5 (+ an (// n an))))))))))
    ;number is a complex short float, coerce val's components to shorts.
    ;But when number's imag part was 0, sqrt eliminates the unnecessary work and
    ;val will be a float!  So we have to detect those cases and make the
    ;appropriate type of complex from val.
    (if (complexp number)
	(cond ((typep (%complex-real-part number) 'short-float)
	       (if (complexp val)
	       	(%complex-cons (float (%complex-real-part val) 0s0)
			       (float (%complex-imag-part val) 0s0))
		 (%complex-cons (float val 0s0) 0s0)))
	      ((complexp val) val)
	      (t (%complex-cons val 0.0)))
      (if (typep number 'short-float)
	  (short-float val)
	val))))

))

; From modified file DJ: L.SYS2; NUMER.LISP#83 at 22-Oct-88 04:21:59
#10R SYSTEM-INTERNALS#: 
(COMPILER-LET ((*PACKAGE* (PKG-FIND-PACKAGE "SYSTEM-INTERNALS")))
  (COMPILER::PATCH-SOURCE-FILE "SYS: SYS2; NUMER  "

(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 (ALOC (RANDOM-VECTOR ARRAY) 0))
   (SETF (RANDOM-POINTER-2 ARRAY) (CL: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. (list (BYTE 12. 12.) (BYTE 12. 0.)))
	       (25. (list (BYTE 12. 12.) (BYTE 12. 0.) (BYTE 1. 24.)))
	       (31. (list (BYTE 12. 12.) (BYTE 12. 0.) (BYTE 9. 24.)))
	       (T (FERROR "Internal error -- unexpected %%Q-POINTER value = ~D." %%q-pointer))))
     (DO ((I 0 (1+ I))) ((= I SIZE))
       (SETQ X (%POINTER-TIMES X 4093.))			;4093. is a prime number.
       (%P-DPB-OFFSET (LDB (BYTE 12. 11.) X) BYTE-SPEC POINTER I))))

))