linalg.lsp

   1 ;;; -*- mode: lisp -*-
   2 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
   3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
   4 ;;; Since 1991, ANSI was finally finished.  Edited for ANSI Common Lisp.
   5 ;;;
   6 ;;; what this should do:
   7 ;;; #1 - use CFFI (and possibly Verazanno) to import C/C++.
   8 ;;; #2 - what to do for Fortran?  Possibly: C <-> bridge, or CLapack?
   9 ;;;      problem: would be better to have access to Fortran.  For
  10 ;;;      example, think of Doug Bates comment on reverse-calls (as
  11 ;;;      distinct from callbacks).  It would be difficult if we don't
  12 ;;;      -- however, has anyone run Lapack or similar through F2CL?
  13 ;;;      Answer: yes, Matlisp does this.
  14
  15
  16 ;;;; linalg -- Lisp-Stat interface to basic linear algebra routines.
  17 ;;;;
  18 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
  19 ;;;; unrestricted use.
  20
  21 ;;;;
  22 ;;;; Package Setup
  23 ;;;;
  24
  25 (defpackage :lisp-stat-linalg
  26   (:use :common-lisp
  27         :lisp-stat-matrix)
  28   (:export chol-decomp lu-decomp lu-solve determinant inverse sv-decomp
  29            qr-decomp rcondest make-rotation spline kernel-dens kernel-smooth
  30            fft make-sweep-matrix sweep-operator ax+y numgrad numhess
  31            split-list eigen))
  32
  33
  34 ;; (in-package #:lisp-stat-basics)
  35
  36
  37 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  38 ;;;;
  39 ;;;;                 Lisp to C number conversion and checking
  40 ;;;;
  41 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  42
  43 ;;;;
  44 ;;;; Lisp to/from C sequence and matrix conversion and checking
  45 ;;;;
  46
  47 (defun is-cons (a)
  48   "FIXME:AJR this is not used anywhere?"
  49   (if (consp a) 1 0))
  50
  51 (defun check-fixnum (a)
  52   (if (/= 0 (la-data-mode a)) (error "not an integer sequence - ~s" a)))
  53
  54 (defun check-real (data)
  55   (let ((data (compound-data-seq data)))
  56     (cond
  57      ((vectorp data)
  58       (let ((n (length data)))
  59         (declare (fixnum n))
  60         (dotimes (i n)
  61                  (declare (fixnum i))
  62                  (check-one-real (aref data i)))))
  63      ((consp data) (dolist (x data) (check-one-real x)))
  64      (t (error "bad sequence - ~s" data)))))
  65
  66 (defun vec-assign (a i x) (setf (aref a i) x))
  67
  68 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  69 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  70 ;;;;
  71 ;;;;             Lisp Interfaces to Linear Algebra Routines
  72 ;;;;
  73 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  74 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  75
  76 ;;; FIXME: use dpbt[f2|rf], dpbstf, dpot[f2|rf]; dpptrf, zpbstf, zpbt[f2|rf]
  77 ;;; remember: factorization = decomposition, depending on training.
  78
  79 (defun chol-decomp (a &optional (maxoffl 0.0))
  80 "Args: (a)
  81 Modified Cholesky decomposition. A should be a square, symmetric matrix.
  82 Computes lower triangular matrix L such that L L^T = A + D where D is a
  83 diagonal matrix. If A is strictly positive definite D will be zero.
  84 Otherwise, D is as small as possible to make A + D numerically strictly
  85 positive definite. Returns a list (L (max D))."
  86   (check-square-matrix a)
  87   (check-real a)
  88   (let* ((n (array-dimension a 0))
  89          (result (make-array (list n n)))
  90          (dpars (list maxoffl 0.0)))
  91     (check-real dpars)
  92     (let ((mat (la-data-to-matrix a mode-re))
  93           (dp (la-data-to-vector dpars mode-re)))
  94       (unwind-protect
  95           (progn
  96             (chol-decomp-front mat n dp)
  97             (la-matrix-to-data mat n n mode-re result)
  98             (la-vector-to-data dp 2 mode-re dpars))
  99         (la-free-matrix mat n)
 100         (la-free-vector dp)))
 101     (list result (second dpars))))
 102
 103
 104 ;;; REPLACE with
 105 ;;;         (matlisp:lu M)
 106 ;;; i.e. result use by:
 107 ;;;         (setf (values (lu-out1 lu-out2 lu-out3)) (matlisp:lu my-matrix))
 108 ;;; for solution, ...
 109 ;;; for lu-solve:
 110 ;;;         (matlisp:gesv a b &opt ipivot)
 111
 112 (defun lu-decomp (a)
 113 "Args: (a)
 114 A is a square matrix of numbers (real or complex). Computes the LU
 115 decomposition of A and returns a list of the form (LU IV D FLAG), where
 116 LU is a matrix with the L part in the lower triangle, the U part in the
 117 upper triangle (the diagonal entries of L are taken to be 1), IV is a vector
 118 describing the row permutation used, D is 1 if the number of permutations
 119 is odd, -1 if even, and FLAG is T if A is numerically singular, NIL otherwise.
 120 Used bu LU-SOLVE."
 121   (check-square-matrix a)
 122   (let* ((n (array-dimension a 0))
 123          (mode (max mode-re (la-data-mode a)))
 124          (result (list (make-array (list n n)) (make-array n) nil nil)))
 125     (let ((mat (la-data-to-matrix a mode))
 126           (iv (la-vector n mode-in))
 127           (d (la-vector 1 mode-re))
 128           (singular 0))
 129       (unwind-protect
 130           (progn
 131             (setf singular (lu-decomp-front mat n iv mode d))
 132             (la-matrix-to-data mat n n mode (first result))
 133             (la-vector-to-data iv n mode-in (second result))
 134             (setf (third result) (la-get-double d 0))
 135             (setf (fourth result) (if (= singular 0.0) nil t)))
 136         (la-free-matrix mat n)
 137         (la-free-vector iv)
 138         (la-free-vector d)))
 139     result))
 140
 141 (defun lu-solve (lu lb)
 142 "Args: (lu b)
 143 LU is the result of (LU-DECOMP A) for a square matrix A, B is a sequence.
 144 Returns the solution to the equation Ax = B. Signals an error if A is
 145 singular."
 146   (let ((la (first lu))
 147         (lidx (second lu)))
 148     (check-square-matrix la)
 149     (check-sequence lidx)
 150     (check-sequence lb)
 151     (check-fixnum lidx)
 152     (let* ((n (num-rows la))
 153            (result (make-sequence (if (consp lb) 'list 'vector) n))
 154            (a-mode (la-data-mode la))
 155            (b-mode (la-data-mode lb)))
 156       (if (/= n (length lidx)) (error "index sequence is wrong length"))
 157       (if (/= n (length lb)) (error "right hand side is wrong length"))
 158       (let* ((mode (max mode-re a-mode b-mode))
 159              (a (la-data-to-matrix la mode))
 160              (indx (la-data-to-vector lidx mode-in))
 161              (b (la-data-to-vector lb mode))
 162              (singular 0))
 163         (unwind-protect
 164             (progn
 165               (setf singular (lu-solve-front a n indx b mode))
 166               (la-vector-to-data b n mode result))
 167           (la-free-matrix a n)
 168           (la-free-vector indx)
 169           (la-free-vector b))
 170         (if (/= 0.0 singular) (error "matrix is (numerically) singular"))
 171         result))))
 172
 173 (defun determinant (a)
 174 "Args: (m)
 175 Returns the determinant of the square matrix M."
 176   (let* ((lu (lu-decomp a))
 177          (la (first lu))
 178          (n (num-rows a))
 179          (d1 (third lu))
 180          (d2 0.d0))
 181     (declare (fixnum n))
 182     (flet ((fabs (x) (float (abs x) 0.d0)))
 183       (dotimes (i n (* d1 (exp d2)))
 184         (declare (fixnum i))
 185         (let* ((x (aref la i i))
 186                (magn (fabs x)))
 187           (if (= 0.0 magn) (return 0.d0))
 188           (setf d1 (* d1 (/ x magn)))
 189           (setf d2 (+ d2 (log magn))))))))
 190
 191 (defun inverse (a)
 192 "Args: (m)
 193 Returns the inverse of the the square matrix M; signals an error if M is ill
 194 conditioned or singular"
 195   (check-square-matrix a)
 196   (let ((n (num-rows a))
 197         (mode (max mode-re (la-data-mode a))))
 198     (declare (fixnum n))
 199     (let ((result (make-array (list n n) :initial-element 0)))
 200       (dotimes (i n)
 201                (declare (fixnum i))
 202                (setf (aref result i i) 1))
 203       (let ((mat (la-data-to-matrix a mode))
 204             (inv (la-data-to-matrix result mode))
 205             (iv (la-vector n mode-in))
 206             (v (la-vector n mode))
 207             (singular 0))
 208         (unwind-protect
 209             (progn
 210               (setf singular (lu-inverse-front mat n iv v mode inv))
 211               (la-matrix-to-data inv n n mode result))
 212           (la-free-matrix mat n)
 213           (la-free-matrix inv n)
 214           (la-free-vector iv)
 215           (la-free-vector v))
 216         (if (/= singular 0) (error "matrix is (numerically) singular"))
 217         result))))
 218
 219 ;;;;
 220 ;;;; SV Decomposition
 221 ;;;;
 222
 223 (defun sv-decomp (a)
 224 "Args: (a)
 225 A is a matrix of real numbers with at least as many rows as columns.
 226 Computes the singular value decomposition of A and returns a list of the form
 227 (U W V FLAG) where U and V are matrices whose columns are the left and right
 228 singular vectors of A and W is the sequence of singular values of A. FLAG is T
 229 if the algorithm converged, NIL otherwise."
 230   (check-matrix a)
 231   (let* ((m (num-rows a))
 232          (n (num-cols a))
 233          (mode (max mode-re (la-data-mode a)))
 234          (result (list (make-array (list m n))
 235                        (make-array n)
 236                        (make-array (list n n))
 237                        nil)))
 238     (if (< m n) (error "number of rows less than number of columns"))
 239     (if (= mode mode-cx) (error "complex SVD not available yet"))
 240     (let ((mat (la-data-to-matrix a mode))
 241           (w (la-vector n mode-re))
 242           (v (la-matrix n n mode-re))
 243           (converged 0))
 244       (unwind-protect
 245           (progn
 246             (setf converged (sv-decomp-front mat m n w v))
 247             (la-matrix-to-data mat m n mode (first result))
 248             (la-vector-to-data w n mode (second result))
 249             (la-matrix-to-data v n n mode (third result))
 250             (setf (fourth result) (if (/= 0.0 converged) t nil)))
 251         (la-free-matrix mat m)
 252         (la-free-vector w)
 253         (la-free-matrix v n))
 254       result)))
 255
 256
 257 ;;;;
 258 ;;;; QR Decomposition
 259 ;;;;
 260
 261 (defun qr-decomp (a &optional pivot)
 262 "Args: (a &optional pivot)
 263 A is a matrix of real numbers with at least as many rows as columns. Computes
 264 the QR factorization of A and returns the result in a list of the form (Q R).
 265 If PIVOT is true the columns of X are first permuted to place insure the
 266 absolute values of the diagonal elements of R are nonincreasing. In this case
 267 the result includes a third element, a list of the indices of the columns in
 268 the order in which they were used."
 269   (check-matrix a)
 270   (let* ((m (num-rows a))
 271          (n (num-cols a))
 272          (mode (max mode-re (la-data-mode a)))
 273          (p (if pivot 1 0))
 274          (result (if pivot
 275                      (list (make-array (list m n))
 276                            (make-array (list n n))
 277                            (make-array n))
 278                      (list (make-array (list m n)) (make-array (list n n))))))
 279     (if (< m n) (error "number of rows less than number of columns"))
 280     (if (= mode mode-cx) (error "complex QR decomposition not available yet"))
 281     (let ((mat (la-data-to-matrix a mode))
 282           (v (la-matrix n n mode-re))
 283           (jpvt (la-vector n mode-in)))
 284       (unwind-protect
 285           (progn
 286             (qr-decomp-front mat m n v jpvt p)
 287             (la-matrix-to-data mat m n mode (first result))
 288             (la-matrix-to-data v n n mode (second result))
 289             (if pivot (la-vector-to-data jpvt n mode-in (third result))))
 290         (la-free-matrix mat m)
 291         (la-free-matrix v n)
 292         (la-free-vector jpvt))
 293       result)))
 294
 295 ;;;;
 296 ;;;; Estimate of Condition Number for Lower Triangular Matrix
 297 ;;;;
 298
 299 (defun rcondest (a)
 300 "Args: (a)
 301 Returns an estimate of the reciprocal of the L1 condition number of an upper
 302 triangular matrix a."
 303   (check-square-matrix a)
 304   (let ((mode (max mode-re (la-data-mode a)))
 305         (n (num-rows a)))
 306     (if (= mode mode-cx)
 307         (error "complex condition estimate not available yet"))
 308     (let ((mat (la-data-to-matrix a mode))
 309           (est 0.0))
 310       (unwind-protect
 311           (setf est (rcondest-front mat n))
 312         (la-free-matrix mat n))
 313       est)))
 314
 315 ;;;;
 316 ;;;; Make Rotation Matrix
 317 ;;;;
 318
 319 (defun make-rotation (x y &optional alpha)
 320 "Args: (x y &optional alpha)
 321 Returns a rotation matrix for rotating from X to Y, or from X toward Y
 322 by angle ALPHA, in radians. X and Y are sequences of the same length."
 323   (check-sequence x)
 324   (check-sequence y)
 325   (if alpha (check-one-real alpha))
 326   (let* ((n (length x))
 327          (mode (max mode-re (la-data-mode x) (la-data-mode y)))
 328          (use-angle (if alpha 1 0))
 329          (angle (if alpha (float alpha 0.0) 0.0))
 330          (result (make-array (list n n))))
 331     (if (/= n (length y)) (error "sequences not the same length"))
 332     (if (= mode mode-cx) (error "complex data not supported yet"))
 333     (let ((px (la-data-to-vector x mode-re))
 334           (py (la-data-to-vector y mode-re))
 335           (rot (la-matrix n n mode-re)))
 336       (unwind-protect
 337           (progn
 338             (make-rotation-front n rot px py use-angle angle)
 339             (la-matrix-to-data rot n n mode-re result))
 340         (la-free-vector px)
 341         (la-free-vector py)
 342         (la-free-matrix rot n))
 343       result)))
 344
 345 ;;;;
 346 ;;;; Eigenvalues and Vectors
 347 ;;;;
 348
 349 (defun eigen (a)
 350 "Args: (a)
 351 Returns list of list of eigenvalues and list of eigenvectors of square,
 352 symmetric matrix A. Third element of result is NIL if algorithm converges.
 353 If the algorithm does not converge, the third element is an integer I.
 354 In this case the eigenvalues 0, ..., I are not reliable."
 355   (check-square-matrix a)
 356   (let ((mode (max mode-re (la-data-mode a)))
 357         (n (num-rows a)))
 358     (if (= mode mode-cx) (error "matrix must be real and symmetric"))
 359     (let ((evals (make-array n))
 360           (evecs (make-list (* n n)))
 361           (pa (la-data-to-vector (compound-data-seq a) mode-re))
 362           (w (la-vector n mode-re))
 363           (z (la-vector (* n n) mode-re))
 364           (fv1 (la-vector n mode-re))
 365           (ierr 0))
 366       (unwind-protect
 367           (progn
 368             (setf ierr (eigen-front pa n w z fv1))
 369             (la-vector-to-data w n mode-re evals)
 370             (la-vector-to-data z (* n n) mode-re evecs))
 371         (la-free-vector pa)
 372         (la-free-vector z)
 373         (la-free-vector w)
 374         (la-free-vector fv1))
 375       (list (nreverse evals)
 376             (nreverse (mapcar #'(lambda (x) (coerce x 'vector))
 377                               (split-list evecs n)))
 378             (if (/= 0 ierr) (- n ierr))))))
 379
 380 ;;;;
 381 ;;;; Spline Interpolation
 382 ;;;;
 383
 384 (defun make-smoother-args (x y xvals)
 385   (check-sequence x)
 386   (check-real x)
 387   (when y
 388         (check-sequence y)
 389         (check-real y))
 390   (unless (integerp xvals)
 391           (check-sequence xvals)
 392           (check-real xvals))
 393   (let* ((n (length x))
 394          (ns (if (integerp xvals) xvals (length xvals)))
 395          (result (list (make-list ns) (make-list ns))))
 396     (if (and y (/= n (length y))) (error "sequences not the same length"))
 397     (list x y n (if (integerp xvals) 0 1) ns xvals result)))
 398
 399 (defun get-smoother-result (args) (seventh args))
 400
 401 (defmacro with-smoother-data ((x y xvals is-reg) &rest body)
 402   `(progn
 403      (check-sequence ,x)
 404      (check-real ,x)
 405      (when ,is-reg
 406            (check-sequence ,y)
 407            (check-real ,y))
 408      (unless (integerp ,xvals)
 409              (check-sequence ,xvals)
 410              (check-real ,xvals))
 411      (let* ((supplied (not (integerp ,xvals)))
 412             (n (length ,x))
 413             (ns (if supplied (length ,xvals) ,xvals))
 414             (result (list (make-list ns) (make-list ns))))
 415        (if (and ,is-reg (/= n (length ,y)))
 416            (error "sequences not the same length"))
 417        (if (and (not supplied) (< ns 2))
 418            (error "too few points for interpolation"))
 419        (let* ((px (la-data-to-vector ,x mode-re))
 420               (py (if ,is-reg (la-data-to-vector ,y mode-re)))
 421               (pxs (if supplied
 422                        (la-data-to-vector ,xvals mode-re)
 423                        (la-vector ns mode-re)))
 424               (pys (la-vector ns mode-re)))
 425          (unless supplied (la-range-to-rseq n px ns pxs))
 426          (unwind-protect
 427              (progn ,@body
 428                     (la-vector-to-data pxs ns mode-re (first result))
 429                     (la-vector-to-data pys ns mode-re (second result)))
 430            (la-free-vector px)
 431            (if ,is-reg (la-free-vector py))
 432            (la-free-vector pxs)
 433            (la-free-vector pys))
 434          result))))
 435
 436 (defun spline (x y &key (xvals 30))
 437 "Args: (x y &key xvals)
 438 Returns list of x and y values of natural cubic spline interpolation of (X,Y).
 439 X must be strictly increasing. XVALS can be an integer, the number of equally
 440 spaced points to use in the range of X, or it can be a sequence of points at
 441 which to interpolate."
 442   (with-smoother-data (x y xvals t)
 443     (let ((work (la-vector (* 2 n) mode-re))
 444           (error 0))
 445       (unwind-protect
 446           (setf error (spline-front n px py ns pxs pys work))
 447         (la-free-vector work))
 448       (if (/= error 0) (error "bad data for splines")))))
 449
 450 ;;;;
 451 ;;;; Kernel Density Estimators and Smoothers
 452 ;;;;
 453
 454 (defun kernel-type-code (type)
 455   (cond ((eq type 'u) 0)
 456         ((eq type 't) 1)
 457         ((eq type 'g) 2)
 458         (t 3)))
 459
 460 (defun kernel-dens (x &key (type 'b) (width -1.0) (xvals 30))
 461 "Args: (x &key xvals width type)
 462 Returns list of x and y values of kernel density estimate of X. XVALS can be an
 463 integer, the number of equally spaced points to use in the range of X, or it
 464 can be a sequence of points at which to interpolate. WIDTH specifies the
 465 window width. TYPE specifies the lernel and should be one of the symbols G, T,
 466 U or B for gaussian, triangular, uniform or bisquare. The default is B."
 467   (check-one-real width)
 468   (with-smoother-data (x nil xvals nil) ;; warning about deleting unreachable code is TRUE -- 2nd arg=nil!
 469     (let ((code (kernel-type-code type))
 470           (error 0))
 471       (setf error (kernel-dens-front px n width pxs pys ns code))
 472       (if (/= 0 error) (error "bad kernel density data")))))
 473
 474 (defun kernel-smooth (x y &key (type 'b) (width -1.0) (xvals 30))
 475 "Args: (x y &key xvals width type)
 476 Returns list of x and y values of kernel smooth of (X,Y). XVALS can be an
 477 integer, the number of equally spaced points to use in the range of X, or it
 478 can be a sequence of points at which to interpolate. WIDTH specifies the
 479 window width. TYPE specifies the lernel and should be one of the symbols G, T,
 480 U or B for Gaussian, triangular, uniform or bisquare. The default is B."
 481   (check-one-real width)
 482   (with-smoother-data (x y xvals t)
 483     (let ((code (kernel-type-code type))
 484           (error 0))
 485       ;;(kernel-smooth-front px py n width pxs pys ns code)
 486       (kernel-smooth-Cport px py n width pxs pys ns code)
 487       (if (/= 0 error) (error "bad kernel density data")))))
 488
 489
 490
 491 (defun kernel-smooth-Cport (px py n width ;;wts wds ;; see above for mismatch?
 492                             xs ys ns ktype)
 493   "Port of kernel_smooth (Lib/kernel.c) to Lisp.
 494 FIXME:kernel-smooth-Cport"
 495   (cond ((< n 1) 1.0)
 496         ((and (< n 2) (<= width 0)) 1.0)
 497         (t (let* ((xmin (min px))
 498                   (xmax (max px))
 499                   (width (/ (- xmax xmin) (+ 1.0 (log n)))))
 500              (dotimes  (i (- ns 1))
 501                (setf (aref ys i)
 502                      (let ((wsum 0.0)
 503                            (ysum 0.0))
 504                        (dotimes (j (- n 1))   ) ;;remove for the following
 505 ;;;                      (let* ;; FIXME!?
 506 ;;;                          ((lwidth (if wds (* width (aref wds j)) width))
 507 ;;;                           (lwt (* (kernel-Cport (aref xs i) (aref px j) lwidth ktype) ;; px?
 508 ;;;                                   (if wts (aref wts j) 1.0))))
 509 ;;;                        (setf wsum (+ wsum lwt))
 510 ;;;                        (setf ysum (if py (+ ysum (* lwt (aref py j))))))) ;; py? y?
 511                        (if py
 512                            (if (> wsum 0.0)
 513                                (/ ysum wsum)
 514                                0.0)
 515                            (/ wsum n)))))
 516              (values ys)))))
 517
 518 (defun kernel-Cport (x y w ktype)
 519   "Port of kernel() (Lib/kernel.c) to Lisp.
 520 x,y,w are doubles, type is an integer"
 521   (if (<= w 0.0)
 522       0.0
 523       (let ((z (- x y)))
 524         (cond ((eq ktype "B")
 525                (let* ((w (* w 2.0))
 526                       (z (* z 0.5)))
 527                  (if (and (> z -0.5)
 528                           (< z 0.5))
 529                      (/ (/ (* 15.0 (* (- 1.0 (* 4 z z))  ;; k/w
 530                                       (- 1.0 (* 4 z z)))) ;; k/w
 531                            8.0)
 532                         w)
 533                      0)))
 534               ((eq ktype "G")
 535                (let* ((w (* w 0.25))
 536                       (z (* z 4.0))
 537                       (k (/ (exp (* -0.5 z z))
 538                             (sqrt (* 2 PI)))))
 539                  (/ k w)))
 540               ((eq ktype "U")
 541                (let* ((w (* 1.5 w))
 542                       (z (* z 0.75))
 543                       (k (if (< (abs z) 0.5)
 544                              1.0
 545                              0.0)))
 546                  (/ k w)))
 547               ((eq ktype "T")
 548                (cond ((and (> z -1.0)
 549                            (< z 0.0))
 550                       (+ 1.0 z))  ;; k
 551                      ((and (> z 0.0)
 552                            (< z 1.0))
 553                       (- 1.0 z))  ;; k
 554                      (t 0.0)))
 555               (t (values 0.0))))))
 556
 557
 558 ;;;;
 559 ;;;; Lowess Smoother Interface
 560 ;;;;
 561
 562 (defun |base-lowess| (s1 s2 f nsteps delta)
 563   (check-sequence s1)
 564   (check-sequence s2)
 565   (check-real s1)
 566   (check-real s2)
 567   (check-one-real f)
 568   (check-one-fixnum nsteps)
 569   (check-one-real delta)
 570   (let* ((n (length s1))
 571          (result (make-list n)))
 572     (if (/= n (length s2)) (error "sequences not the same length"))
 573     (let ((x (la-data-to-vector s1 mode-re))
 574           (y (la-data-to-vector s2 mode-re))
 575           (ys (la-vector n mode-re))
 576           (rw (la-vector n mode-re))
 577           (res (la-vector n mode-re))
 578           (error 0))
 579       (unwind-protect
 580           (progn
 581             (setf error (base-lowess-front x y n f nsteps delta ys rw res))
 582             (la-vector-to-data ys n mode-re result))
 583         (la-free-vector x)
 584         (la-free-vector y)
 585         (la-free-vector ys)
 586         (la-free-vector rw)
 587         (la-free-vector res))
 588       (if (/= error 0) (error "bad data for lowess"))
 589       result)))
 590
 591 #|
 592 static LVAL add_contour_point(i, j, k, l, x, y, z, v, result)
 593         int i, j, k, l;
 594         RVector x, y;
 595         RMatrix z;
 596         double v;
 597         LVAL result;
 598 {
 599   LVAL pt;
 600   double p, q;
 601
 602   if ((z[i][j] <= v && v < z[k][l]) || (z[k][l] <= v && v < z[i][j])) {
 603     xlsave(pt);
 604     pt = mklist(2, NIL);
 605     p = (v - z[i][j]) / (z[k][l] - z[i][j]);
 606     q = 1.0 - p;
 607     rplaca(pt, cvflonum((FLOTYPE) (q * x[i] + p * x[k])));
 608     rplaca(cdr(pt), cvflonum((FLOTYPE) (q * y[j] + p * y[l])));
 609     result = cons(pt, result);
 610     xlpop();
 611   }
 612   return(result);
 613 }
 614
 615 LVAL xssurface_contour()
 616 {
 617   LVAL s1, s2, mat, result;
 618   RVector x, y;
 619   RMatrix z;
 620   double v;
 621   int i, j, n, m;
 622
 623   s1 = xsgetsequence();
 624   s2 = xsgetsequence();
 625   mat = xsgetmatrix();
 626   v = makedouble(xlgetarg());
 627   xllastarg();
 628
 629   n = seqlen(s1); m = seqlen(s2);
 630   if (n != numrows(mat) || m != numcols(mat)) xlfail("dimensions do not match");
 631   if (data_mode(s1) == CX || data_mode(s2) == CX || data_mode(mat) == CX)
 632     xlfail("data must be real");
 633
 634   x = (RVector) data_to_vector(s1, RE);
 635   y = (RVector) data_to_vector(s2, RE);
 636   z = (RMatrix) data_to_matrix(mat, RE);
 637
 638   xlsave1(result);
 639   result = NIL;
 640   for (i = 0; i < n - 1; i++) {
 641         for (j = 0; j < m - 1; j++) {
 642           result = add_contour_point(i, j, i, j+1, x, y, z, v, result);
 643           result = add_contour_point(i, j+1, i+1, j+1, x, y, z, v, result);
 644           result = add_contour_point(i+1, j+1, i+1, j, x, y, z, v, result);
 645           result = add_contour_point(i+1, j, i, j, x, y, z, v, result);
 646         }
 647   }
 648   xlpop();
 649
 650   free_vector(x);
 651   free_vector(y);
 652   free_matrix(z, n);
 653
 654   return(result);
 655 }
 656 |#
 657
 658 ;;;
 659 ;;; FFT
 660 ;;;
 661 ;;; FIXME:ajr
 662 ;;; replace with matlisp:fft  and matlisp:ifft (the latter for inverse mapping)
 663 ;;;
 664 (defun fft (x &optional inverse)
 665 "Args: (x &optional inverse)
 666 Returns unnormalized Fourier transform of X, or inverse transform if INVERSE
 667 is true."
 668   (check-sequence x)
 669   (let* ((n (length x))
 670          (mode (la-data-mode x))
 671          (isign (if inverse -1 1))
 672          (result (if (consp x) (make-list n) (make-array n))))
 673     (let ((px (la-data-to-vector x mode-cx))
 674           (work (la-vector (+ (* 4 n) 15) mode-re)))
 675       (unwind-protect
 676           (progn
 677             (fft-front n px work isign)
 678             (la-vector-to-data px n mode-cx result))
 679         (la-free-vector px)
 680         (la-free-vector work))
 681       result)))
 682
 683 ;;;
 684 ;;; SWEEP Operator: FIXME: use matlisp
 685 ;;;
 686
 687 (defun make-sweep-front (x y w n p mode has_w x_mean result)
 688   (declare (fixnum n p mode has_w))
 689   (let ((x_data nil)
 690         (result_data nil)
 691         (val 0.0)
 692         (dxi 0.0)
 693         (dyi 0.0)
 694         (dv 0.0)
 695         (dw 0.0)
 696         (sum_w 0.0)
 697         (dxik 0.0)
 698         (dxjk 0.0)
 699         (dyj 0.0)
 700         (dx_meani 0.0)
 701         (dx_meanj 0.0)
 702         (dy_mean 0.0)
 703         (has-w (if (/= 0 has_w) t nil))
 704         (RE 1))
 705     (declare (long-float val dxi dyi dv dw sum_w dxik dxjk dyj
 706                     dx_meani dx_meanj dy_mean))
 707 ;;    (declare-double val dxi dyi dv dw sum_w dxik dxjk dyj
 708 ;;                  dx_meani dx_meanj dy_mean)
 709
 710     (if (> mode RE) (error "not supported for complex data yet"))
 711
 712     (setf x_data (compound-data-seq x))
 713     (setf result_data (compound-data-seq result))
 714
 715     ;; find the mean of y
 716     (setf val 0.0)
 717     (setf sum_w 0.0)
 718     (dotimes (i n)
 719       (declare (fixnum i))
 720       (setf dyi (makedouble (aref y i)))
 721       (when has-w
 722         (setf dw (makedouble (aref w i)))
 723         (incf sum_w dw)
 724         (setf dyi (* dyi dw)))
 725       (incf val dyi))
 726     (if (not has-w) (setf sum_w (float n 0.0)))
 727     (if (<= sum_w 0.0) (error "non positive sum of weights"))
 728     (setf dy_mean (/ val sum_w))
 729
 730     ;; find the column means
 731     (dotimes (j p)
 732       (declare (fixnum j))
 733       (setf val 0.0)
 734       (dotimes (i n)
 735         (declare (fixnum i))
 736         (setf dxi (makedouble (aref x_data (+ (* p i) j))))
 737         (when has-w
 738           (setf dw (makedouble (aref w i)))
 739           (setf dxi (* dxi dw)))
 740         (incf val dxi))
 741       (setf (aref x_mean j) (/ val sum_w)))
 742
 743     ;; put 1/sum_w in topleft, means on left, minus means on top
 744     (setf (aref result_data 0) (/ 1.0 sum_w))
 745     (dotimes (i p)
 746       (declare (fixnum i))
 747       (setf dxi (makedouble (aref x_mean i)))
 748       (setf (aref result_data (+ i 1)) (- dxi))
 749       (setf (aref result_data (* (+ i 1) (+ p 2))) dxi))
 750     (setf (aref result_data (+ p 1)) (- dy_mean))
 751     (setf (aref result_data (* (+ p 1) (+ p 2))) dy_mean)
 752
 753     ;; put sums of adjusted cross products in body
 754     (dotimes (i p)
 755       (declare (fixnum i))
 756       (dotimes (j p)
 757         (declare (fixnum j))
 758         (setf val 0.0)
 759         (dotimes (k n)
 760           (declare (fixnum k))
 761           (setf dxik (makedouble (aref x_data (+ (* p k) i))))
 762           (setf dxjk (makedouble (aref x_data (+ (* p k) j))))
 763           (setf dx_meani (makedouble (aref x_mean i)))
 764           (setf dx_meanj (makedouble (aref x_mean j)))
 765           (setf dv (* (- dxik dx_meani) (- dxjk dx_meanj)))
 766           (when has-w
 767             (setf dw (makedouble (aref w k)))
 768             (setf dv (* dv dw)))
 769           (incf val dv))
 770         (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ j 1))) val)
 771         (setf (aref result_data (+ (* (+ j 1) (+ p 2)) (+ i 1))) val))
 772       (setf val 0.0)
 773       (dotimes (j n)
 774         (declare (fixnum j))
 775         (setf dxik (makedouble (aref x_data (+ (* p j) i))))
 776         (setf dyj (makedouble (aref y j)))
 777         (setf dx_meani (makedouble (aref x_mean i)))
 778         (setf dv (* (- dxik dx_meani) (- dyj dy_mean)))
 779         (when has-w
 780           (setf dw (makedouble (aref w j)))
 781           (setf dv (* dv dw)))
 782         (incf val dv))
 783       (setf (aref result_data (+ (* (+ i 1) (+ p 2)) (+ p 1))) val)
 784       (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ i 1))) val))
 785     (setf val 0.0)
 786     (dotimes (j n)
 787       (declare (fixnum j))
 788       (setf dyj (makedouble (aref y j)))
 789       (setf dv (* (- dyj dy_mean) (- dyj dy_mean)))
 790       (when has-w
 791         (setf dw (makedouble (aref w j)))
 792         (setf dv (* dv dw)))
 793       (incf val dv))
 794     (setf (aref result_data (+ (* (+ p 1) (+ p 2)) (+ p 1))) val)))
 795
 796 ;;; FIXME: use matlisp
 797 (defun sweep-in-place-front (a rows cols mode k tol)
 798   "Sweep algorithm for linear regression."
 799   (declare (long-float tol))
 800   (declare (fixnum rows cols mode k))
 801   (let ((data nil)
 802         (pivot 0.0)
 803         (aij 0.0)
 804         (aik 0.0)
 805         (akj 0.0)
 806         (akk 0.0)
 807         (RE 1))
 808     (declare (long-float pivot aij aik akj akk))
 809
 810     (if (> mode RE) (error "not supported for complex data yet"))
 811     (if (or (< k 0) (>= k rows) (>= k cols)) (error "index out of range"))
 812
 813     (setf tol (max tol machine-epsilon))
 814     (setf data (compound-data-seq a))
 815
 816     (setf pivot (makedouble (aref data (+ (* cols k) k))))
 817
 818     (cond
 819      ((or (> pivot tol) (< pivot (- tol)))
 820       (dotimes (i rows)
 821         (declare (fixnum i))
 822         (dotimes (j cols)
 823           (declare (fixnum j))
 824           (when (and (/= i k) (/= j k))
 825             (setf aij (makedouble (aref data (+ (* cols i) j))))
 826             (setf aik (makedouble (aref data (+ (* cols i) k))))
 827             (setf akj (makedouble (aref data (+ (* cols k) j))))
 828             (setf aij (- aij (/ (* aik akj) pivot)))
 829             (setf (aref data (+ (* cols i) j)) aij))))
 830
 831       (dotimes (i rows)
 832         (declare (fixnum i))
 833         (setf aik (makedouble (aref data (+ (* cols i) k))))
 834         (when (/= i k)
 835           (setf aik (/ aik pivot))
 836           (setf (aref data (+ (* cols i) k)) aik)))
 837
 838       (dotimes (j cols)
 839         (declare (fixnum j))
 840         (setf akj (makedouble (aref data (+ (* cols k) j))))
 841         (when (/= j k)
 842           (setf akj (- (/ akj pivot)))
 843           (setf (aref data (+ (* cols k) j)) akj)))
 844
 845       (setf akk (/ 1.0 pivot))
 846       (setf (aref data (+ (* cols k) k)) akk)
 847       1)
 848      (t 0))))
 849
 850 ;;; FIXME: use matlisp
 851 (defun make-sweep-matrix (x y &optional w)
 852 "Args: (x y &optional weights)
 853 X is a matrix, Y and WEIGHTS are sequences. Returns the sweep matrix for the
 854 (possibly weighted) regression of Y on X."
 855   (check-matrix x)
 856   (check-sequence y)
 857   (if w (check-sequence w))
 858   (let ((n (num-rows x))
 859         (p (num-cols x)))
 860     (if (/= n (length y)) (error "dimensions do not match"))
 861     (if (and w (/= n (length w))) (error "dimensions do not match"))
 862     (let ((mode (max (la-data-mode x)
 863                      (la-data-mode x)
 864                      (if w (la-data-mode w) 0)))
 865           (result (make-array (list (+ p 2) (+ p 2))))
 866           (x-mean (make-array p))
 867           (y (coerce y 'vector))
 868           (w (if w (coerce w 'vector)))
 869           (has-w (if w 1 0)))
 870       (make-sweep-front x y w n p mode has-w x-mean result)
 871       result)))
 872
 873 (defun sweep-in-place (a k tol)
 874   (check-matrix a)
 875   (check-one-fixnum k)
 876   (check-one-real tol)
 877   (let ((rows (num-rows a))
 878         (cols (num-cols a))
 879         (mode (la-data-mode a)))
 880     (let ((swept (sweep-in-place-front a rows cols mode k tol)))
 881       (if (/= 0 swept) t nil))))
 882
 883 (defun sweep-operator (a columns &optional tolerances)
 884 "Args: (a indices &optional tolerances)
 885 A is a matrix, INDICES a sequence of the column indices to be swept. Returns
 886 a list of the swept result and the list of the columns actually swept. (See
 887 MULTREG documentation.) If supplied, TOLERANCES should be a list of real
 888 numbers the same length as INDICES. An index will only be swept if its pivot
 889 element is larger than the corresponding element of TOLERANCES."
 890   (check-matrix a)
 891   (check-sequence columns)
 892   (if tolerances (check-sequence tolerances))
 893   (check-real a)
 894   (check-fixnum columns)
 895   (if tolerances (check-real tolerances))
 896   (do ((tol .0000001)
 897        (result (copy-array a))
 898        (swept-columns nil)
 899        (columns (coerce columns 'list) (cdr columns))
 900        (tolerances (if (consp tolerances) (coerce tolerances 'list))
 901                    (if (consp tolerances) (cdr tolerances))))
 902       ((null columns) (list result swept-columns))
 903       (let ((col (first columns))
 904             (tol (if (consp tolerances) (first tolerances) tol)))
 905         (if (sweep-in-place result col tol)
 906             (setf swept-columns (cons col swept-columns))))))
 907
 908
 909 ;;;
 910 ;;; AX+Y
 911 ;;;
 912
 913 ;;; matlisp:axpy
 914 ;;;
 915 (defun ax+y (a x y &optional lower)
 916 "Args (a x y &optional lower)
 917 Returns (+ (matmult A X) Y). If LOWER is not nil, A is taken to be lower
 918 triangular.
 919 This could probably be made more efficient."
 920   (check-square-matrix a)
 921   (check-sequence x)
 922   (check-sequence y)
 923   (check-real a)
 924   (check-real x)
 925   (check-real y)
 926   (let* ((n (num-rows a))
 927          (result (make-list n))
 928          (a (compound-data-seq a)))
 929     (declare (fixnum n))
 930     (if (or (/= n (length x)) (/= n (length y)))
 931         (error "dimensions do not match"))
 932     (do* ((tx (make-next-element x) (make-next-element x))
 933           (ty (make-next-element y))
 934           (tr (make-next-element result))
 935           (i 0 (+ i 1))
 936           (start 0 (+ start n))
 937           (end (if lower (+ i 1) n) (if lower (+ i 1) n)))
 938          ((<= n i) result)
 939        (declare (fixnum i start end))
 940        (let ((val (get-next-element ty i)))
 941          (dotimes (j end)
 942            (declare (fixnum j))
 943            (setf val (+ val (* (get-next-element tx j)
 944                                (aref a (+ start j))))))
 945          (set-next-element tr i val)))))
 946
 947 ;;;;
 948 ;;;; Maximization and Numerical Derivatives
 949 ;;;;
 950
 951 (defvar *maximize-callback-function* nil)
 952 (defvar *maximize-callback-arg* nil)
 953
 954 (defun data2double (n data ptr)
 955   (declare (fixnum n))
 956   (let* ((seq (compound-data-seq data))
 957          (elem (make-next-element seq)))
 958     (if (/= (length seq) n) (error "bad data size"))
 959     (dotimes (i n)
 960       (declare (fixnum i))
 961       (la-put-double ptr i (get-next-element elem i)))))
 962
 963 (defun maximize-callback (n px pfval pgrad phess pderivs)
 964   (la-vector-to-data px n mode-re *maximize-callback-arg*)
 965   (let* ((val (funcall *maximize-callback-function* *maximize-callback-arg*))
 966          (derivs (if (consp val) (- (length val) 1) 0)))
 967     (la-put-integer pderivs 0 derivs)
 968     (la-put-double pfval 0 (if (consp val) (first val) val))
 969     (if (<= 1 derivs) (data2double n (second val) pgrad))
 970     (if (<= 2 derivs) (data2double (* n n) (third val) phess))))
 971
 972 (defun numgrad (f x &optional scale (h -1.0))
 973 "Args: (f x &optional scale derivstep)
 974 Computes the numerical gradient of F at X."
 975   (check-sequence x)
 976   (check-real x)
 977   (when scale
 978         (check-sequence scale)
 979         (check-real scale))
 980   (check-one-real h)
 981   (let* ((n (length x))
 982          (result (make-list n)))
 983     (if (and scale (/= n (length scale)))
 984         (error "scale not the same length as x"))
 985     (let ((*maximize-callback-function* f)
 986           (*maximize-callback-arg* (make-list n)))
 987       (let ((px (la-data-to-vector x mode-re))
 988             (pgrad (la-vector n mode-re))
 989             (pscale (la-data-to-vector
 990                      (if scale scale (make-list n :initial-element 1.0))
 991                      mode-re)))
 992         (unwind-protect
 993             (progn
 994               (numgrad-front n px pgrad h pscale)
 995               (la-vector-to-data pgrad n mode-re result))
 996           (la-free-vector px)
 997           (la-free-vector pgrad)
 998           (la-free-vector pscale))))
 999       result))
1000
1001 (defun numhess (f x &optional scale (h -1.0) all)
1002 "Args: (f x &optional scale derivstep)
1003 Computes the numerical Hessian matrix of F at X."
1004   (check-sequence x)
1005   (check-real x)
1006   (when scale
1007         (check-sequence scale)
1008         (check-real scale))
1009   (check-one-real h)
1010   (let* ((n (length x))
1011          (result (if all
1012                      (list nil (make-list n) (make-array (list n n)))
1013                      (make-array (list n n)))))
1014     (if (and scale (/= n (length scale)))
1015         (error "scale not the same length as x"))
1016     (let ((*maximize-callback-function* f)
1017           (*maximize-callback-arg* (make-list n)))
1018       (let ((hess-data (compound-data-seq (if all (third result) result)))
1019             (px (la-data-to-vector x mode-re))
1020             (pf (la-vector 1 mode-re))
1021             (pgrad (la-vector n mode-re))
1022             (phess (la-vector (* n n) mode-re))
1023             (pscale (la-data-to-vector
1024                      (if scale scale (make-list n :initial-element 1.0))
1025                      mode-re)))
1026         (unwind-protect
1027             (progn
1028               (numhess-front n px pf pgrad phess h pscale)
1029               (when all
1030                     (setf (first result) (la-get-double pf 0))
1031                     (la-vector-to-data pgrad n mode-re (second result)))
1032               (la-vector-to-data phess (* n n) mode-re hess-data))
1033           (la-free-vector pf)
1034           (la-free-vector px)
1035           (la-free-vector pgrad)
1036           (la-free-vector phess)
1037           (la-free-vector pscale))))
1038     result))
1039
1040 (defun init-minfo-ipar-values (n ipars)
1041   (let* ((TRUE 1)
1042          (FALSE 0)
1043          (k 0)
1044          (m 0)
1045          (itnlimit -1)
1046          (backtrack TRUE)
1047          (verbose 0)
1048          (vals_suppl FALSE)
1049          (exptilt TRUE)
1050          (count 0)
1051          (termcode 0))
1052     (setf (aref ipars 0) n)
1053     (setf (aref ipars 1) m)
1054     (setf (aref ipars 2) k)
1055     (setf (aref ipars 3) itnlimit)
1056     (setf (aref ipars 4) backtrack)
1057     (setf (aref ipars 5) verbose)
1058     (setf (aref ipars 6) vals_suppl)
1059     (setf (aref ipars 7) exptilt)
1060     (setf (aref ipars 8) count)
1061     (setf (aref ipars 9) termcode)))
1062
1063 (defun init-minfo-dpar-values (h dpars)
1064   (let ((typf 1.0)
1065         (gradtol -1.0)
1066         (steptol -1.0)
1067         (maxstep -1.0)
1068         (dflt 0.0)
1069         (tilt 0.0)
1070         (newtilt 0.0)
1071         (hessadd 0.0))
1072     (setf (aref dpars 0) typf)
1073     (setf (aref dpars 1) h)
1074     (setf (aref dpars 2) gradtol)
1075     (setf (aref dpars 3) steptol)
1076     (setf (aref dpars 4) maxstep)
1077     (setf (aref dpars 5) dflt)
1078     (setf (aref dpars 6) tilt)
1079     (setf (aref dpars 7) newtilt)
1080     (setf (aref dpars 8) hessadd)))
1081
1082 (defun init-minfo-internals (n h internals)
1083   (let ((ipars (aref internals 8))
1084         (dpars (aref internals 9)))
1085     (init-minfo-ipar-values n ipars)
1086     (init-minfo-dpar-values h dpars)))
1087
1088 (defun new-minfo-internals (f x &key scale ((:derivstep h) -1.0))
1089   (check-sequence x)
1090   (check-real x)
1091   (check-one-real h)
1092   (let ((n (length x)))
1093     (when scale
1094           (check-sequence scale)
1095           (check-real scale)
1096           (if (/= n (length scale)) (error "scale and x not the same length")))
1097     (let ((internals (make-array 12)))
1098       (setf (aref internals 0) f)
1099       (setf (aref internals 3) (if (consp x) (copy-list x) (coerce x 'list)))
1100       (setf (aref internals 4)
1101             (if scale (copy-seq scale) (make-array n :initial-element 1.0)))
1102       (setf (aref internals 5) (make-list (+ 1 n (* n n))))
1103       (setf (aref internals 8) (make-array 10))
1104       (setf (aref internals 9) (make-array 9))
1105       (init-minfo-internals n h internals)
1106       internals)))
1107
1108 (defun minfo-maximize (internals &optional verbose)
1109   "This function does what?"
1110   (let* ((f (aref internals 0))
1111          (x (aref internals 3))
1112          (fvals (aref internals 5))
1113          (n (length x))
1114          (v (if verbose (if (integerp verbose) verbose 1) -1)))
1115     (setf (aref internals 3) (copy-list x))
1116     (setf (aref internals 5) (copy-list fvals))
1117     (let ((*maximize-callback-function* f)
1118           (*maximize-callback-arg* (make-list n)))
1119       (let* ((x (aref internals 3))
1120              (scale (aref internals 4))
1121              (fvals (aref internals 5))
1122              (ip (aref internals 8))
1123              (dp (aref internals 9))
1124              (px (la-data-to-vector x mode-re))
1125              (pscale (la-data-to-vector scale mode-re))
1126              (pfvals (la-vector (length fvals) mode-re))
1127              (pip (la-data-to-vector ip mode-in))
1128              (pdp (la-data-to-vector dp mode-re)))
1129         (unwind-protect
1130             (progn
1131               (base-minfo-maximize px pfvals pscale pip pdp v)) ;; access to C
1132           (la-vector-to-data px n mode-re x)
1133           (la-vector-to-data pfvals (+ 1 n (* n n)) mode-re fvals)
1134           (la-vector-to-data pip (length ip) mode-in ip)
1135           (la-vector-to-data pdp (length dp) mode-re dp))
1136         (get-buf)))))
1137
1138 ;;;;
1139 ;;;; Miscellaneous Routines
1140 ;;;;
1141
1142
1143 (defun split-list (x n)
1144 "Args: (list cols)
1145 Returns a list of COLS lists of equal length of the elements of LIST.
1146 Example: (split-list '(1 2 3 4 5 6) 2) returns ((1 2 3) (4 5 6))"
1147   (check-one-fixnum n)
1148   (if (/= (rem (length x) n) 0) (error "length not divisible by ~a" n))
1149   (flet ((next-split ()
1150            (let ((result nil)
1151                  (end nil))
1152              (dotimes (i n result)
1153                (declare (fixnum i))
1154                (let ((c-elem (list (first x))))
1155                  (cond ((null result)
1156                         (setf result c-elem)
1157                         (setf end result))
1158                        (t
1159                         (setf (rest end) c-elem)
1160                         (setf end (rest end)))))
1161                (setf x (rest x))))))
1162     (let ((result nil)
1163           (end nil)
1164           (k (/ (length x) n)))
1165       (declare (fixnum k))
1166       (dotimes (i k result)
1167         (declare (fixnum i))
1168         (let ((c-sub (list (next-split))))
1169           (cond ((null result)
1170                  (setf result c-sub)
1171                  (setf end result))
1172                 (t
1173                  (setf (rest end) c-sub)
1174                  (setf end (rest end)))))))))