Convex Hull Software (v20140708) ================================ The following software implements the algorithms described in the PhD thesis "Efficient Construction of Convex Hulls from Spherical Distributions in Higher Dimensions" by Wm. Stewart in 1992. The programs are written in what was then IBM's main programming language "PL/1", which could do absolutely anything, and often did. The program constructs convex hulls in any dimension higher than 1. They were tested in up to 20 dimensions. Many of the variable names are only one letter, however note correspond to the mathematical naming from the thesis, e.g. "S" for a set of points, “f” for a facet, and "v" for a set of vertices. Rather than using pointers, a self-initializing doubly linked list is implemented, so that it grows as new nodes are requested. This approach was chosen to avoid storage overflow, and to minimize the possibility of anomalies in the cpu timings caused by dynamic memory allocation. The recursion of the Traverse algorithm is unwrapped, so it can nest as deep as required without performance issues. The Quicksort function at the end is highly optimized, and about as fast as a sort algorithm can be implemented. (The file is formatted for a fixed-width font, such as Courier New.) /*====================================================================*/ Convex : procedure ( S, n, d, S_index, P, p_size, m, proc_type, error_code ) ; /*---------------------------------------------------------------------- | | | GPL 3 Licensed. | | | | (c) Copyright W. Stewart, 5 September 1991. | | | | 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 3 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 can find a copy of the full license at | | www.gnu.org/licenses/gpl.txt. | | | ----------------------------------------------------------------------*/ /*------------------------------------------------------------------- | | | Name: | | | | Convex | | | | Description: | | | | Find conv(S). Return a linked set of facets P. | | | | Usage: | | | | call Convex ( S, n, d, S_index, P, p_size, m, proc_type, | | error_code ) | | | | S ( n, d ) - Set of n points in Rd. | | | | n - Size of S. | | | | d - Dimension. | | | | S_index - Permutation index of S if sorted. | | | | P ( * ) - Set of facets. | | | | vertices ( d ) - Indices of points of facet. | | neighbor ( d ) - Pointers to neighboring facets. | | h_y ( d ) - Coefficients of halfspace equation (h_y,h_b)| | h_b - RHS constant of halfspace equation. | | h_length - Length of halfspace vector h_y. | | non_visible_facet - Set on creation, used by Interlink. | | visible_facet - Set on creation, used by Interlink. | | in_v - Set if facet is put in visible set V. | | f_link - Forward link. | | b_link - Backward link. | | t_link - Temporary link, used to put facets in the | | visible set V and the cap C. | | | | p_size - Size of P data structure. Must be as large as | | the maximum number of facets required. | | | | m - Number of facets in the convex hull. | | | | error_code - | | 0 = No error. | | 1 = Overflow of facet data structure. | | 2 = Wrong path taken by Two_Flat. | | 3 = Failure of Create_Facet. | | | | proc_type - | | 1 = Preprocessing construction. | | 2 = Online construction. | | | | | | Logical Program Structure: | | | | Convex_initialization | | Sort1 (if proc_type is preprocessing) | | Create_simplex | | Get_new_facet | | Create_H | | Lneqs | | | | Two_Flat (if proc_type is online) | | Lneqs | | | | Traverse | | Create_and_link_new_facet | | Get_new_facet | | Find_common_subfacet | | Create_H | | Lneqs | | | | Interlink | | Rotate | | Find_common_subfacet | | | | Delete_V | | | | Convex_termination | | | | Reference: | | An Algorithm to Find the Facets of the Convex Hull in | | Higher Dimensions, SCCCG Proceedings, August 1990. | | | | Author: | | Wm. M. Stewart | | Computer Science | | University of New Brunswick | | Fredericton, New Brunswick | | Canada, E3B 5A3 | | | | Creation Date: | | 4 December 1989 | | | | Last Revision Date: | | 5 September 1991 | | | -------------------------------------------------------------------*/ /*------------------------------------------------------------------- | | | Variable dictionary: | | | | avail - Head of available list of facets in P. | | c_header - Head of list of cap facets in C. | | centroid ( d ) - Centroid of first simplex P, used to set | | the correct sign of halfspace equations. | | d2 - Constant set to d+2, for declaration of | | arrays of size d+2. | | f - A facet. | | i - Index of current point in S. | | k - Visible facet found by two_flat algorithm. | | link_full - Switch, turned on when the dynamic | | initialization of the linked list of facets| | finishes. | | link_init_counter - Counter for dynamic initialization of the | | linked list of facets. | | p ( d ) - Current point. | | stack ( p_size ) - Stack to unwrap the recursion of Traverse. | | subsubfacet_size - Constant set to d-2, for declaration of | | subsubfacet by Interlink. | | v_header - Head of list of visible facets in V. | | | |------------------------------------------------------------------*/ declare s ( *, * ) float binary ( 21 ), n fixed binary ( 31 ), d fixed binary ( 31 ), s_index ( * ) fixed binary ( 31 ), 1 P ( * ) controlled, 2 vertices ( * ) fixed binary ( 31 ), 2 neighbor ( * ) fixed binary ( 31 ), 2 h_y ( * ) float binary ( 53 ), 2 h_b float binary ( 53 ), 2 h_length float binary ( 53 ), 2 non_visible_facet fixed binary ( 31 ), 2 visible_facet fixed binary ( 31 ), 2 in_v bit ( 1 ), 2 f_link fixed binary ( 31 ), 2 b_link fixed binary ( 31 ), 2 t_link fixed binary ( 31 ), p_size fixed binary ( 31 ), m fixed binary ( 31 ), proc_type fixed binary ( 31 ), error_code fixed binary ( 31 ), avail fixed binary ( 31 ), c_header fixed binary ( 31 ), centroid ( d ) float binary ( 53 ), d2 fixed binary ( 31 ), f fixed binary ( 31 ), i fixed binary ( 31 ), k fixed binary ( 31 ), link_full bit ( 1 ), link_init_counter fixed binary ( 31 ), p ( d ) float binary ( 53 ), stack ( p_size ) fixed binary ( 31 ) ctl, subsubfacet_size fixed binary ( 31 ), v_header fixed binary ( 31 ) ; allocate stack ; /* Create initial simplex, centroid, initialize s_index, etc. */ call convex_initialization ; /* Add points i=d+1...n. */ i = d+1 ; do while ( i < n & error_code = 0 ) ; i = i + 1 ; p ( * ) = s ( i, * ) ; if proc_type = 2 then k = two_flat ; if k > 0 then do ; /* Update P with p_i. */ f = k ; k = 0 ; in_v ( f ) = '1'b ; t_link ( f ) = 0 ; v_header = f ; c_header = 0 ; call traverse ( f ) ; if error_code = 0 then do ; call interlink ; call delete_v ; end ; end ; end ; end ; if error_code ^= 0 then do ; put skip list ( '*** ERROR *** Error_code =', error_code, 'I = ', i ) ; end ; return ; /*=================================================================*/ traverse : procedure ( f ) ; /*------------------------------------------------------------------- | Traverse and identify the set of visible facets V, create the | | cap C, and link C to the nonvisible facets on the horizon. | -------------------------------------------------------------------*/ declare stack_top fixed binary ( 31 ), f fixed binary ( 31 ), g fixed binary ( 31 ), j fixed binary ( 31 ), more bit ( 1 ), sum builtin, t fixed binary ( 31 ) ; stack_top = 0 ; j = 0 ; more = '1'b ; do while ( more & error_code = 0 ) ; j = j + 1 ; g = neighbor ( f, j ) ; if ^ in_v ( g ) then if sum (h_y(g,*)*p(*)) - h_b(g) > 0.0 then do ; in_v ( g ) = '1'b ; t_link ( g ) = v_header ; v_header = g ; stack_top = stack_top + 1 ; stack ( stack_top ) = g ; end ; else do ; t = create_and_link_new_facet ; t_link ( t ) = c_header ; c_header = t ; end ; if j = d then do ; if stack_top = 0 then more = '0'b ; else do ; f = stack ( stack_top ) ; stack_top = stack_top - 1 ; j = 0 ; end ; end ; end ; return ; /*=================================================================*/ create_and_link_new_facet : procedure ; /*------------------------------------------------------------------- | Create a new facet and link to neighbor on the horizon. | -------------------------------------------------------------------*/ declare jj fixed binary ( 31 ), t fixed binary ( 31 ) ; t = get_new_facet ; if error_code = 0 then do ; vertices ( t, * ) = vertices ( f, * ) ; non_visible_facet ( t ) = g ; visible_facet ( t ) = f ; do jj = j to d - 1 ; vertices ( t, jj ) = vertices ( t, jj+1 ) ; end ; vertices ( t, d ) = i ; neighbor ( t, d ) = g ; jj = find_common_subfacet ( t, g ) ; neighbor ( g, jj ) = t ; call create_h ( t ) ; if proc_type = 1 & i < n & k = 0 then if sum (h_y(t,*)*s(i+1,*)) - h_b(t) > 0.0 then k = t ; end ; return ( t ) ; end create_and_link_new_facet ; end traverse ; /*=================================================================*/ interlink : procedure ; /*------------------------------------------------------------------- | Interlink the new facets of the cap with each other. | -------------------------------------------------------------------*/ declare c_size fixed binary ( 31 ), j fixed binary ( 31 ), t fixed binary ( 31 ), /* the variables below are used by rotate. */ a fixed binary ( 31 ), b fixed binary ( 31 ), g fixed binary ( 31 ), h fixed binary ( 31 ), l fixed binary ( 31 ), phi ( subsubfacet_size ) fixed binary ( 31 ), x fixed binary ( 31 ), w fixed binary ( 31 ) ; c_size = 0 ; t = c_header ; do while ( t ^= 0 ) ; c_size = c_size + 1 ; do j = 1 to d-1 ; if neighbor ( t, j ) = 0 then call rotate ; end ; t = t_link ( t ) ; end ; m = m + c_size ; return ; /*=================================================================*/ rotate : procedure ; /*------------------------------------------------------------------- | Rotate around a subsubfacet phi. | -------------------------------------------------------------------*/ do l = 1 to j-1 ; phi ( l ) = vertices ( t, l ) ; end ; do l = j to d-2 ; phi ( l ) = vertices ( t, l+1 ) ; end ; g = non_visible_facet ( t ) ; h = visible_facet ( t ) ; do while ( in_v ( h ) ) ; a = 1 ; b = 1 ; x = 0 ; do while ( b <= d-2 & x = 0 ) ; if vertices ( h, a ) = phi ( b ) then do ; a = a + 1 ; b = b + 1 ; end ; else if neighbor ( h, a ) ^= g then x = neighbor ( h, a ) ; else a = a + 1 ; end ; do while ( x = 0 ) ; if neighbor ( h, a ) ^= g then x = neighbor ( h, a ) ; else a = a + 1 ; end ; g = h ; h = x ; end ; l = find_common_subfacet ( g, h ) ; w = neighbor ( h, l ) ; neighbor ( t, j ) = w ; l = find_common_subfacet ( t, w ) ; neighbor ( w, l ) = t ; return ; end rotate ; end interlink ; /*=================================================================*/ find_common_subfacet : procedure ( f, g ) ; /*------------------------------------------------------------------- | Find j such that sbf^F_j is shared by F and G. | -------------------------------------------------------------------*/ declare a fixed binary ( 31 ), b fixed binary ( 31 ), f fixed binary ( 31 ), g fixed binary ( 31 ), j fixed binary ( 31 ) ; j = 0 ; a = 1 ; b = 1 ; do while ( j = 0 ) ; if vertices ( f, a ) = vertices ( g, b ) then do ; a = a + 1 ; b = b + 1 ; end ; else if vertices ( f, a ) < vertices ( g, b ) then a = a + 1 ; else j = b ; if a > d then j = d ; end ; return ( j ) ; end find_common_subfacet ; /*=================================================================*/ two_flat : procedure returns ( fixed binary ( 31 ) ) ; /*------------------------------------------------------------------- | Find a first visible facet by 2-flat intersection, defined by | | p, the centroid of the first facet in P, and the neighboring | | facet G of F closest to p. The triangle defined by these three | | points is then shifted to ensure it passes through the subfacet | | shared by F and G, and then made into an isoceles. | -------------------------------------------------------------------*/ declare a ( d2, d2 ) float binary ( 53 ), b ( d2 ) float binary ( 53 ), b_save ( d2 ) float binary ( 53 ), d_save fixed binary ( 31 ), det float binary ( 53 ), distance float binary ( 53 ), distance_t float binary ( 53 ), f fixed binary ( 31 ), float builtin, g fixed binary ( 31 ), h fixed binary ( 31 ), half_path_switch bit ( 1 ), ii fixed binary ( 31 ), intersect_switch bit ( 1 ), j fixed binary ( 31 ), jj fixed binary ( 31 ), l fixed binary ( 31 ), last_angle float binary ( 53 ), last_last_angle float binary ( 53 ), llast_last_angle float binary ( 53 ), last_f fixed binary ( 31 ), llast_f fixed binary ( 31 ), last_g fixed binary ( 31 ), llast_g fixed binary ( 31 ), length1 float binary ( 53 ), length2 float binary ( 53 ), length float binary ( 53 ), mod builtin, more bit ( 1 ), next_angle float binary ( 53 ), next_facet fixed binary ( 31 ), q1 ( d ) float binary ( 53 ), q1t ( d ) float binary ( 53 ), q2 ( d ) float binary ( 53 ), q2t ( d ) float binary ( 53 ), q3 ( d ) float binary ( 53 ), start_facet fixed binary ( 31 ), sum builtin, sqrt builtin, v1 ( d ) float binary ( 53 ), v2 ( d ) float binary ( 53 ) ; f = b_link ( 1 ) ; /* First facet in P. */ if sum (h_y(f,*)*p(*)) - h_b(f) > 0.0 then return ( f ) ; jj = 1 ; g = neighbor ( f, jj ) ; distance = (sum(h_y(g,*) * p(*)) - h_b(g)) / h_length(g) ; do j = 2 to d ; g = neighbor ( f, j ) ; distance_t = (sum(h_y(g,*) * p(*)) - h_b(g)) / h_length(g) ; if distance_t > distance then do ; distance = distance_t ; jj = j ; end ; end ; g = neighbor ( f, jj ) ; /* Find the midpoints of f, g, and their common subfacet. */ q1t ( * ) = 0.0 ; q2t ( * ) = 0.0 ; q3 ( * ) = 0.0 ; do j = 1 to d ; q1t ( * ) = q1t ( * ) + s ( vertices ( f, j ), * ) ; q2t ( * ) = q2t ( * ) + s ( vertices ( g, j ), * ) ; end ; do j = 1 to d-1 ; l = mod ( jj + j - 1, d ) + 1 ; q3 ( * ) = q3 ( * ) + s ( vertices ( f, l ), * ) ; end ; q1t ( * ) = q1t ( * ) / float ( d ) ; q2t ( * ) = q2t ( * ) / float ( d ) ; q3 ( * ) = q3 ( * ) / float ( d-1 ) ; /* Set q1 and q2. Make triangle (p,q1,q2) an isoceles. */ | q1 ( * ) = q3 ( * ) + 3 * ( q2t ( * ) - q1t ( * ) ) ; q2 ( * ) = q3 ( * ) - 3 * ( q2t ( * ) - q1t ( * ) ) ; q1t ( * ) = q1 ( * ) - p ( * ) ; q2t ( * ) = q2 ( * ) - p ( * ) ; length1 = sqrt ( sum ( q1t ( * ) * q1t ( * ) ) ) ; length2 = sqrt ( sum ( q2t ( * ) * q2t ( * ) ) ) ; if length1 > length2 then q2 ( * ) = p ( * ) + ( length1 / length2 ) * q2t ( * ) ; else do ; q1 ( * ) = p ( * ) + ( length2 / length1 ) * q1t ( * ) ; end ; do ii = 1 to d ; b_save ( ii ) = 0.0 ; end; b_save ( d+1 ) = 1.0 ; b_save ( d+2 ) = 1.0 ; last_f = f ; llast_f = f ; last_g = g ; llast_g = g ; start_facet = f ; last_angle = 2.0 ; last_last_angle = last_angle ; llast_last_angle = last_angle ; half_path_switch = '1'b; v1 ( * ) = q3 ( * ) - p ( * ) ; length = sqrt ( sum ( v1 ( * ) * v1 ( * ) ) ) ; v1 ( * ) = v1 ( * ) / length ; /* Walk around 2-flat path. */ do while ( error_code = 0 & g ^= start_facet & half_path_switch & sum (h_y(g,*)*p(*)) - h_b(g) < 0.0 ) ; next_facet = 0 ; do ii = 1 to d while ( next_facet = 0 ) ; h = neighbor ( g, ii ) ; if h ^= f then do ; do j = 1 to d-1 ; l = mod ( ii + j - 1, d ) + 1 ; do jj = 1 to d ; a ( jj, j ) = s ( vertices ( g, l ), jj ) ; end ; a ( d+1, j ) = 1.0 ; a ( d+2, j ) = 0.0 ; end ; do j = 1 to d ; a ( j, d ) = -p ( j ) ; a ( j, d+1 ) = -q1 ( j ) ; a ( j, d+2 ) = -q2 ( j ) ; end ; do j = d to d+2 ; a ( d+1, j ) = 0.0 ; a ( d+2, j ) = 1.0 ; end ; d_save = d + 2 ; b ( * ) = b_save ( * ) ; call lneqs ( b, det, a, d_save ) ; if det = 0.0 then do ; error_code = 2 ; end ; else do ; intersect_switch = '1'b ; do j = 1 to d-1 ; if b ( j ) < 0.0 then intersect_switch = '0'b ; end ; if intersect_switch then next_facet = h ; end ; end ; end ; if next_facet = 0 then error_code = 2 ; if error_code = 0 then do ; llast_f = last_f ; last_f = f ; f = g ; llast_g = last_g ; last_g = g ; g = next_facet ; /* Determine if path has gone more than half way around. */ v2 ( * ) = ( b ( d ) * p ( * ) + b ( d+1 ) * q1 ( * ) + b ( d+2 ) * q2 ( * ) ) ; v2 ( * ) = v2 ( * ) - p ( * ) ; length = sqrt ( sum ( v2 ( * ) * v2 ( * ) ) ) ; v2 ( * ) = v2 ( * ) / length ; next_angle = sum ( v1 ( * ) * v2 ( * ) ) ; if next_angle <= last_last_angle then do ; llast_last_angle = last_last_angle ; last_last_angle = last_angle ; last_angle = next_angle ; end ; else do ; half_path_switch = '0'b ; end ; end ; end ; if error_code ^= 0 | ^half_path_switch | sum (h_y(g,*)*p(*)) - h_b(g) < 0.0 then do ; g = 0 ; end ; return ( g ) ; end two_flat ; /*=================================================================*/ delete_v : procedure ; /*------------------------------------------------------------------- | Delete the set of visible facets V. | -------------------------------------------------------------------*/ declare f fixed binary ( 31 ), v_size fixed binary ( 31 ) ; v_size = 0 ; f = v_header ; do while ( f ^= 0 ) ; f_link ( b_link ( f ) ) = f_link ( f ) ; b_link ( f_link ( f ) ) = b_link ( f ) ; f_link ( f ) = avail ; avail = f ; v_size = v_size + 1 ; f = t_link ( f ) ; end ; m = m - v_size ; return ; end delete_v ; /*=================================================================*/ convex_initialization : procedure ; /*------------------------------------------------------------------- | Initialize S_index. If proc_type=1 sort S lexicographically. | | Create the initial simplex and centroid of the simplex. | -------------------------------------------------------------------*/ declare d3 fixed binary ( 31 ), float builtin, j fixed binary ( 31 ), k fixed binary ( 31 ), least_p ( d, d ) float binary ( 21 ), max_p ( d, d ) float binary ( 21 ) ; error_code = 0 ; subsubfacet_size = d - 2 ; d2 = d + 2 ; do j = 1 to n ; s_index ( j ) = j ; end ; if proc_type = 1 then call sort1 ( s, s_index, 1, n, d ) ; centroid ( * ) = 0.0 ; do j = 1 to d + 1 ; centroid ( * ) = centroid ( * ) + s ( j, * ) ; end ; centroid ( * ) = centroid ( * ) / float ( d + 1 ) ; /* Initialize linked list of facets P. */ m = 0 ; f_link ( 1 ) = 1 ; b_link ( 1 ) = 1 ; link_init_counter = 1 ; link_full = '0'b ; avail = 0 ; d3 = d + 1 ; call create_simplex ; m = d+1 ; return ; /*=================================================================*/ Create_simplex : procedure ; /*------------------------------------------------------------------- | Create initial simplex from first d+1 points of S. | -------------------------------------------------------------------*/ declare f fixed binary ( 31 ), f_index_list ( d3 ) fixed binary ( 31 ), i fixed binary ( 31 ), j fixed binary ( 31 ), l fixed binary ( 31 ), sum builtin, temp fixed binary ( 31 ) ; k = 0 ; do i = 1 to d + 1 ; f_index_list ( i ) = get_new_facet ; end ; do i = 1 to d + 1 ; l = 0 ; f = f_index_list ( i ) ; do j = 1 to d + 1 ; if i ^= j then do ; l = l + 1 ; vertices ( f, l ) = j ; neighbor ( f, l ) = f_index_list ( j ) ; end ; end ; call create_h ( f ) ; /* If preprocessing case, set first visible facet K. */ if proc_type = 1 & k = 0 & n > d+1 then if sum (h_y(f,*)*s(d+2,*)) - h_b(f) > 0.0 then k = f ; end ; return ; end create_simplex ; end convex_initialization ; /*=================================================================*/ get_new_facet : procedure returns ( fixed binary ( 31 ) ) ; /*------------------------------------------------------------------- | Get a new facet from available storage. | -------------------------------------------------------------------*/ declare last_node fixed binary ( 31 ), t fixed binary ( 31 ) ; if ^ link_full then do ; link_init_counter = link_init_counter + 1 ; t = link_init_counter ; if link_init_counter >= p_size then link_full = '1'b ; end ; else if avail ^= 0 then do ; t = avail ; avail = f_link ( avail ) ; end ; else do ; error_code = 1 ; t = 1 ; end ; last_node = b_link ( 1 ) ; f_link ( t ) = 1 ; b_link ( t ) = last_node ; f_link ( last_node ) = t ; b_link ( 1 ) = t ; neighbor ( t, * ) = 0 ; in_v ( t ) = '0'b; return ( t ) ; end get_new_facet ; /*=================================================================*/ create_h : procedure ( f ) ; /*------------------------------------------------------------------- | Create the plane equation for facet F. | -------------------------------------------------------------------*/ declare a ( d, d ) float binary ( 53 ), b ( d ) float binary ( 53 ), dd fixed binary ( 31 ), det float binary ( 53 ), f fixed binary ( 31 ), j fixed binary ( 31 ), sqrt builtin, sum builtin ; do j = 1 to d ; a ( j, * ) = s ( vertices ( f,j ), * ) ; end ; b ( * ) = 1.0 ; h_b ( f ) = 1.0 ; dd = d ; call lneqs ( b, det, a, dd ) ; if det = 0.0 then error_code = 3 ; else do ; h_y ( f, * ) = 0.0 ; h_y ( f, * ) = b ( * ) ; h_length(f) = sqrt ( sum ( h_y(f,*) ** 2 ) ) ; if sum (h_y(f,*) * centroid(*)) - h_b(f) > 0.0 then do ; h_y ( f, * ) = -h_y ( f, * ) ; h_b ( f ) = -h_b ( f ) ; end ; end ; return ; end create_h ; /*=================================================================*/ sort1 : procedure ( a, b, s, e, d ) ; /*--------------------------------------------------------------------- | Name: | | Quicksort | | Description | | This subroutine implement the quicksort algorithm to | | sort a floating point matrix A in lexicographic order. | | As well, array B (probably an index array) is sorted. | | Reference | | Sedgewick, R. (1978) "Implementing Quicksort Algorithms", | | Communications of the ACM, Volume 21, Number 10, October. | ---------------------------------------------------------------------*/ declare a ( *, * ) float binary ( 21 ), b ( * ) fixed binary ( 31 ), d fixed binary ( 31 ), done bit ( 1 ), e fixed binary ( 31 ), i fixed binary ( 31 ), j fixed binary ( 31 ), k fixed binary ( 31 ), l fixed binary ( 31 ), m fixed binary ( 31 ), lstack ( 100 ) fixed binary ( 31 ), r fixed binary ( 31 ), rstack ( 100 ) fixed binary ( 31 ), s fixed binary ( 31 ), stackp fixed binary ( 31 ), tempa ( d ) float binary ( 21 ), tempb fixed binary ( 31 ), trunc builtin , v ( d ) float binary ( 21 ), x fixed binary ( 31 ); l = s ; r = e ; m = 10 ; done = '0'b ; stackp = 0 ; if r - l + 1 <= m then done = '1'b; do while ( ^ done ) ; x = trunc ( ( l + r ) / 2 ) ; tempa = a ( x, * ) ; a ( x, * ) = a ( l + 1, * ) ; a ( l + 1, * ) = tempa ; tempb = b ( x ) ; b ( x ) = b ( l + 1 ) ; b ( l + 1 ) = tempb ; do k = 1 to d while ( a ( l + 1, k ) = a ( r, k ) ) ; end ; if a ( l + 1, k ) > a ( r, k ) then do; tempa = a ( l + 1, * ) ; a ( l + 1, * ) = a ( r, * ) ; a ( r, * ) = tempa ; tempb = b ( l + 1 ) ; b ( l + 1 ) = b ( r ) ; b ( r ) = tempb ; end ; do k = 1 to d while ( a ( l, k ) = a ( r, k ) ) ; end ; if a ( l, k ) > a ( r, k ) then do ; tempa = a ( l, * ) ; a ( l, * ) = a ( r, * ) ; a ( r, * ) = tempa ; tempb = b ( l ) ; b ( l ) = b ( r ) ; b ( r ) = tempb ; end ; do k = 1 to d while ( a ( l + 1, k ) = a ( l, k ) ) ; end ; if a ( l + 1, k ) > a ( l, k ) then do; tempa = a ( l + 1, * ) ; a ( l + 1, * ) = a ( l, * ) ; a ( l, * ) = tempa ; tempb = b ( l + 1 ) ; b ( l + 1 ) = b ( l ) ; b ( l ) = tempb ; end ; i = l + 1 ; j = r ; v = a ( l, * ) ; do while ( j >= i ) ; i = i + 1 ; do k = 1 to d while ( a ( i, k ) = v ( k ) ) ; end ; do while ( a ( i, k ) < v ( k ) ) ; i = i + 1 ; do k = 1 to d while ( a ( i, k ) = v ( k ) ) ; end ; end ; j = j - 1 ; do k = 1 to d while ( a ( j, k ) = v ( k ) ) ; end ; do while ( a ( j, k ) > v ( k ) ) ; j = j - 1 ; do k = 1 to d while ( a ( j, k ) = v ( k ) ) ; end ; end ; if j < i then goto exit_loop1 ; tempa = a ( i, * ) ; a ( i, * ) = a ( j, * ) ; a ( j, * ) = tempa ; tempb = b ( i ) ; b ( i ) = b ( j ) ; b ( j ) = tempb ; end ; exit_loop1: tempa = a ( l, * ) ; a ( l, * ) = a ( j, * ) ; a ( j, * ) = tempa ; tempb = b ( l ) ; b ( l ) = b ( j ) ; b ( j ) = tempb ; if j - l <= m & r - i + 1 <= m then do ; if stackp = 0 then do ; done = '1'b ; end ; else do ; l = lstack ( stackp ) ; r = rstack ( stackp ) ; stackp = stackp - 1 ; end ; end ; else do ; if j - l <= m | r - i + 1 <= m then do ; if j - l >= r - i + 1 then do ; r = j - 1 ; end ; else do ; l = i ; end ; end ; else do ; stackp = stackp + 1 ; if j - l >= r - i + 1 then do; lstack ( stackp ) = l ; rstack ( stackp ) = j - 1 ; l = i ; end ; else do ; lstack ( stackp ) = i ; rstack ( stackp ) = r ; r = j - 1 ; end ; end ; end ; end ; do i = 2 to e ; tempa = a ( i, * ) ; tempb = b ( i ) ; j = i - 1 ; do k = 1 to d while ( a ( j, k ) = tempa ( k ) ) ; end ; do while ( a ( j, k ) > tempa ( k ) ) ; a ( j+1, * ) = a ( j, * ) ; b ( j+1 ) = b ( j ) ; j = j - 1 ; if j <= 0 then goto exit_loop2 ; do k = 1 to d while ( a ( j, k ) = tempa ( k ) ) ; end ; end ; exit_loop2 : a ( j+1, * ) = tempa ; b ( j+1 ) = tempb ; end ; return ; end sort1; end convex ; /*=================================================================*/ /*=================================================================*/ /*=================================================================*/