;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; LUSH Lisp Universal Shell
;;;   Copyright (C) 2002 Leon Bottou, Yann Le Cun, AT&T Corp, NECI.
;;; Includes parts of TL3:
;;;   Copyright (C) 1987-1999 Leon Bottou and Neuristique.
;;; Includes selected parts of SN3.2:
;;;   Copyright (C) 1991-2001 AT&T Corp.
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License
;;; along with this program; if not, write to the Free Software
;;; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; $Id: idx-sort.lsh,v 1.4 2005/05/16 16:17:54 profshadoko Exp $
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? *** IDX sorting and binary search functions
;; sorting and bsearch functions for vectors of doubles, floats, and ints.
;; These functions only work properly with contiguous vectors.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#? (idx-f1bsearch <am> <v>)
;; returns index of element in idx1 of float <am>
;; whose value is less than or equal to <v>
;; using a binary search.
;; this assumes the elements in <am> are sorted
;; in ascending order.
(de idx-f1bsearch (am v)
    ((-float-) v)
    ((-idx1- (-float-)) am)
    (when (<> 1 (idx-modulo am 0)) (error "vector is not contiguous"))
    (let ((m (idx-dim am 0))
	  (lo 0))
      ((-int-) m lo)
      #{{
      int hi, n;
      float *a;
      a = IDX_PTR($am, float);
      hi = $m;
      while (hi- $lo>1) {
        n = ( $lo + hi )/2;
        if ( $v >= a[n] ) $lo = n; else hi = n;
      }
      } #}
      lo))

#? (idx-d1bsearch <am> <v>)
;; returns index of element in idx1 of double <am>
;; whose value is less than or equal to <v>
;; using a binary search.
;; this assumes the elements in <am> are sorted
;; in ascending order.
(de idx-d1bsearch (am v)
    ((-double-) v)
    ((-idx1- (-double-)) am)
    (when (<> 1 (idx-modulo am 0)) (error "vector is not contiguous"))
    (let ((m (idx-dim am 0))
	  (lo 0))
      ((-int-) m lo)
      #{{
      int hi, n;
      double *a;
      a = IDX_PTR($am, double);
      hi = $m;
      while (hi- $lo>1) {
        n = ( $lo + hi )/2;
        if ( $v >= a[n] ) $lo = n; else hi = n;
      }
      } #}
      lo))

#? (idx-i1bsearch <am> <v>)
;; returns index of element in idx1 of int <am>
;; whose value is less than or equal to <v>
;; using a binary search.
;; this assumes the elements in <am> are sorted
;; in ascending order.
(de idx-i1bsearch (am v)
    ((-int-) v)
    ((-idx1- (-int-)) am)
    (when (<> 1 (idx-modulo am 0)) (error "vector is not contiguous"))
    (let ((m (idx-dim am 0))
	  (lo 0))
      ((-int-) m lo)
      #{{
      int hi, n;
      int *a;
      a = IDX_PTR($am, int);
      hi = $m;
      while (hi- $lo>1) {
        n = ( $lo + hi )/2;
        if ( $v >= a[n] ) $lo = n; else hi = n;
      }
      } #}
      lo))

#? (idx-f1sortup <m>)
;; sorts element in float vector <m> in ascending order
(de idx-f1sortup (m)
    ((-idx1- (-float-)) m)
    (when (<> 1 (idx-modulo m 0)) (error "vector is not contiguous"))
    (let ((n (idx-dim m 0)))
      ((-int-) n)
      (when (> n 1)
      #{ {
      int l,j,ir,i;
      float *ra, rra;

      ra = (float*)(IDX_PTR($m, float)) -1;
      l=($n >> 1)+1;
      ir= $n;
      for (;;) {
	   if (l > 1)
	      rra=ra[--l];
	    else {
	      rra=ra[ir];
	      ra[ir]=ra[1];
	      if (--ir == 1) {
	        ra[1]=rra;
		goto end; }}
        i=l;
	j=l << 1;
	while (j <= ir) {
	  if (j < ir && ra[j] < ra[j+1]) ++j;
	  if (rra < ra[j]) {
	    ra[i]=ra[j];
	    j += (i=j);
	    } else j=ir+1; }
        ra[i]=rra; }}
	end:  #}
	)
      ()))

#? (idx-d1sortup <m>)
;; sorts element in float vector <m> in ascending order
(de idx-d1sortup (m)
    ((-idx1- (-double-)) m)
    (when (<> 1 (idx-modulo m 0)) (error "vector is not contiguous"))
    (let ((n (idx-dim m 0)))
      ((-int-) n)
      (when (> n 1)
      #{ {
      int l,j,ir,i;
      double *ra, rra;

      ra = (double*)(IDX_PTR($m, double)) -1;
      l=($n >> 1)+1;
      ir= $n;
      for (;;) {
	   if (l > 1)
	      rra=ra[--l];
	    else {
	      rra=ra[ir];
	      ra[ir]=ra[1];
	      if (--ir == 1) {
	        ra[1]=rra;
		goto end; }}
        i=l;
	j=l << 1;
	while (j <= ir) {
	  if (j < ir && ra[j] < ra[j+1]) ++j;
	  if (rra < ra[j]) {
	    ra[i]=ra[j];
	    j += (i=j);
	    } else j=ir+1; }
        ra[i]=rra; }}
	end:  #}
	)
      ()))

#? (idx-f1sortdown <m>)
;; sorts element in float vector <m> in descending order
(de idx-f1sortdown (m)
    ((-idx1- (-float-)) m)
    (when (<> 1 (idx-modulo m 0)) (error "vector is not contiguous"))
    (let ((n (idx-dim m 0)))
      ((-int-) n)
      (when (> n 1)
      #{ {
      int l,j,ir,i;
      float *ra, rra;

      ra = (float*)(IDX_PTR($m, float)) -1;
      l=($n >> 1)+1;
      ir= $n;
      for (;;) {
	   if (l > 1)
	      rra=ra[--l];
	    else {
	      rra=ra[ir];
	      ra[ir]=ra[1];
	      if (--ir == 1) {
	        ra[1]=rra;
		goto end; }}
        i=l;
	j=l << 1;
	while (j <= ir) {
	  if (j < ir && ra[j] > ra[j+1]) ++j;
	  if (rra > ra[j]) {
	    ra[i]=ra[j];
	    j += (i=j);
	    } else j=ir+1; }
        ra[i]=rra; }}
	end:  #}
	)
      ()))

#? (idx-d1sortdown <m>)
;; sorts element in double vector <m> in descending order
(de idx-d1sortdown (m)
    ((-idx1- (-double-)) m)
    (when (<> 1 (idx-modulo m 0)) (error "vector is not contiguous"))
    (let ((n (idx-dim m 0)))
      ((-int-) n)
      (when (> n 1)
      #{ {
      int l,j,ir,i;
      double *ra, rra;

      ra = (double*)(IDX_PTR($m, double)) -1;
      l=($n >> 1)+1;
      ir= $n;
      for (;;) {
	   if (l > 1)
	      rra=ra[--l];
	    else {
	      rra=ra[ir];
	      ra[ir]=ra[1];
	      if (--ir == 1) {
	        ra[1]=rra;
		goto end; }}
        i=l;
	j=l << 1;
	while (j <= ir) {
	  if (j < ir && ra[j] > ra[j+1]) ++j;
	  if (rra > ra[j]) {
	    ra[i]=ra[j];
	    j += (i=j);
	    } else j=ir+1; }
        ra[i]=rra; }}
	end:  #}
	)
      ()))


#? (idx-i1sortup <m>)
;; sorts element in integer vector <m> in ascending order
(de idx-i1sortup (m)
    ((-idx1- (-int-)) m)
    (when (<> 1 (idx-modulo m 0)) (error "vector is not contiguous"))
    (let ((n (idx-dim m 0)))
      ((-int-) n)
      (when (> n 1)
      #{ {
      int l,j,ir,i;
      int *ra, rra;

      ra = (int*)(IDX_PTR($m, int)) -1;
      l=($n >> 1)+1;
      ir= $n;
      for (;;) {
	   if (l > 1)
	      rra=ra[--l];
	    else {
	      rra=ra[ir];
	      ra[ir]=ra[1];
	      if (--ir == 1) {
	        ra[1]=rra;
		goto end; }}
        i=l;
	j=l << 1;
	while (j <= ir) {
	  if (j < ir && ra[j] < ra[j+1]) ++j;
	  if (rra < ra[j]) {
	    ra[i]=ra[j];
	    j += (i=j);
	    } else j=ir+1; }
        ra[i]=rra; }}
	end:  #}
	)
      ()))


#? (idx-i1sortdown <m>)
;; sorts element in integer vector <m> in descending order
(de idx-i1sortdown (m)
    ((-idx1- (-int-)) m)
    (when (<> 1 (idx-modulo m 0)) (error "vector is not contiguous"))
    (let ((n (idx-dim m 0)))
      ((-int-) n)
      (when (> n 1)
      #{ {
      int l,j,ir,i;
      int *ra, rra;

      ra = (int*)(IDX_PTR($m, int)) -1;
      l=($n >> 1)+1;
      ir= $n;
      for (;;) {
	   if (l > 1)
	      rra=ra[--l];
	    else {
	      rra=ra[ir];
	      ra[ir]=ra[1];
	      if (--ir == 1) {
	        ra[1]=rra;
		goto end; }}
        i=l;
	j=l << 1;
	while (j <= ir) {
	  if (j < ir && ra[j] > ra[j+1]) ++j;
	  if (rra > ra[j]) {
	    ra[i]=ra[j];
	    j += (i=j);
	    } else j=ir+1; }
        ra[i]=rra; }}
	end:  #}
	)
      ()))

#? (idx-f1i1sortup <m> <p>)
;; <m> is a vector of float, <p> is a vector of integers 
;; (same dimension as <m>).
;; on output, <m> is sorted in ascending order, and <p>
;; is sorted with the same permutation table.
(de idx-f1i1sortup (m p)
    ((-idx1- (-float-)) m)
    ((-idx1- (-int-)) p)
    (when (<> 1 (idx-modulo m 0)) (error "vector is not contiguous"))
    (when (<> 1 (idx-modulo p 0)) (error "vector is not contiguous"))
    (let ((n (idx-dim m 0))
	  (z (idx-dim p 0)))
      ((-int-) n z)
      (when (<> n z) (error "idx-f1i1sortup: vectors have different sizes"))
      (when (> n 1)
      #{ {
      int l,j,ir,i;
      float *ra, rra;
      int *rb, rrb;

      ra = (float*)(IDX_PTR($m, int)) -1;
      rb = (int*)(IDX_PTR($p, int)) -1;

      l=($n >> 1)+1;
      ir=$n;
      for (;;) {
	   if (l > 1) {
	    rra=ra[--l];
	    rrb=rb[l];
	    } else {
	      rra=ra[ir];
	      rrb=rb[ir];
	      ra[ir]=ra[1];
	      rb[ir]=rb[1];
	      if (--ir == 1) {
 	        ra[1]=rra;
		rb[1]=rrb;
		goto end; } }
        i=l;
	j=l << 1;
	while (j <= ir)	{
	  if (j < ir && ra[j] < ra[j+1]) ++j;
	  if (rra < ra[j]) {
	    ra[i]=ra[j];
	    rb[i]=rb[j];
	    j += (i=j);
	    } else j=ir+1; }
	ra[i]=rra;
	rb[i]=rrb; }}
	end: #}
	)
      ()))

#? (idx-f1i1sortdown <m> <p>)
;; <m> is a vector of float, <p> is a vector of integers 
;; (same dimension as <m>).
;; on output, <m> is sorted in descending order, and <p>
;; is sorted with the same permutation table.
(de idx-f1i1sortdown (m p)
    ((-idx1- (-float-)) m)
    ((-idx1- (-int-)) p)
    (when (<> 1 (idx-modulo m 0)) (error "vector is not contiguous"))
    (when (<> 1 (idx-modulo p 0)) (error "vector is not contiguous"))
    (let ((n (idx-dim m 0))
	  (z (idx-dim p 0)))
      ((-int-) n z)
      (when (<> n z) (error "idx-f1i1sortup: vectors have different sizes"))
      (when (> n 1)
      #{ {
      int l,j,ir,i;
      float *ra, rra;
      int *rb, rrb;

      ra = (float*)(IDX_PTR($m, int)) -1;
      rb = (int*)(IDX_PTR($p, int)) -1;

      l=($n >> 1)+1;
      ir=$n;
      for (;;) {
	   if (l > 1) {
	    rra=ra[--l];
	    rrb=rb[l];
	    } else {
	      rra=ra[ir];
	      rrb=rb[ir];
	      ra[ir]=ra[1];
	      rb[ir]=rb[1];
	      if (--ir == 1) {
 	        ra[1]=rra;
		rb[1]=rrb;
		goto end; } }
        i=l;
	j=l << 1;
	while (j <= ir)	{
	  if (j < ir && ra[j] > ra[j+1]) ++j;
	  if (rra > ra[j]) {
	    ra[i]=ra[j];
	    rb[i]=rb[j];
	    j += (i=j);
	    } else j=ir+1; }
	ra[i]=rra;
	rb[i]=rrb; }}
	end: #}
	)
      ()))


#? (idx-i1i1sortup <m> <p>)
;; <m> is a vector of int, <p> is a vector of integers 
;; (same dimension as <m>).
;; on output, <m> is sorted in ascending order, and <p>
;; is sorted with the same permutation table.
(de idx-i1i1sortup (m p)
    ((-idx1- (-int-)) m)
    ((-idx1- (-int-)) p)
    (when (<> 1 (idx-modulo m 0)) (error "vector is not contiguous"))
    (when (<> 1 (idx-modulo p 0)) (error "vector is not contiguous"))
    (let ((n (idx-dim m 0))
	  (z (idx-dim p 0)))
      ((-int-) n z)
      (when (<> n z) (error "idx-i1i1sortup: vectors have different sizes"))
      (when (> n 1)
      #{ {
      int l,j,ir,i;
      int *ra, rra;
      int *rb, rrb;

      ra = (int*)(IDX_PTR($m, int)) -1;
      rb = (int*)(IDX_PTR($p, int)) -1;

      l=($n >> 1)+1;
      ir=$n;
      for (;;) {
	   if (l > 1) {
	    rra=ra[--l];
	    rrb=rb[l];
	    } else {
	      rra=ra[ir];
	      rrb=rb[ir];
	      ra[ir]=ra[1];
	      rb[ir]=rb[1];
	      if (--ir == 1) {
 	        ra[1]=rra;
		rb[1]=rrb;
		goto end; } }
        i=l;
	j=l << 1;
	while (j <= ir)	{
	  if (j < ir && ra[j] < ra[j+1]) ++j;
	  if (rra < ra[j]) {
	    ra[i]=ra[j];
	    rb[i]=rb[j];
	    j += (i=j);
	    } else j=ir+1; }
	ra[i]=rra;
	rb[i]=rrb; }}
	end: #}
	)
      ()))

#? (idx-i1i1sortdown <m> <p>)
;; <m> is a vector of int, <p> is a vector of integers 
;; (same dimension as <m>).
;; on output, <m> is sorted in descending order, and <p>
;; is sorted with the same permutation table.
(de idx-i1i1sortdown (m p)
    ((-idx1- (-int-)) m)
    ((-idx1- (-int-)) p)
    (when (<> 1 (idx-modulo m 0)) (error "vector is not contiguous"))
    (when (<> 1 (idx-modulo p 0)) (error "vector is not contiguous"))
    (let ((n (idx-dim m 0))
	  (z (idx-dim p 0)))
      ((-int-) n z)
      (when (<> n z) (error "idx-i1i1sortup: vectors have different sizes"))
      (when (> n 1)
      #{ {
      int l,j,ir,i;
      int *ra, rra;
      int *rb, rrb;

      ra = (int*)(IDX_PTR($m, int)) -1;
      rb = (int*)(IDX_PTR($p, int)) -1;

      l=($n >> 1)+1;
      ir=$n;
      for (;;) {
	   if (l > 1) {
	    rra=ra[--l];
	    rrb=rb[l];
	    } else {
	      rra=ra[ir];
	      rrb=rb[ir];
	      ra[ir]=ra[1];
	      rb[ir]=rb[1];
	      if (--ir == 1) {
 	        ra[1]=rra;
		rb[1]=rrb;
		goto end; } }
        i=l;
	j=l << 1;
	while (j <= ir)	{
	  if (j < ir && ra[j] > ra[j+1]) ++j;
	  if (rra > ra[j]) {
	    ra[i]=ra[j];
	    rb[i]=rb[j];
	    j += (i=j);
	    } else j=ir+1; }
	ra[i]=rra;
	rb[i]=rrb; }}
	end: #}
	)
      ()))

(dhc-make ()
	  idx-f1bsearch idx-d1bsearch idx-i1bsearch
	  idx-f1sortdown idx-d1sortdown idx-i1sortdown 
	  idx-f1sortup idx-d1sortup idx-i1sortup
	  idx-i1i1sortup idx-i1i1sortdown
	  idx-f1i1sortup idx-f1i1sortdown
)          