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