/*-------------------------------------------------------------------- * The MB-system: mb_surface.c 5/2/94 * $Id: mb_surface.c 1770 2009-10-19 17:16:39Z caress $ * * Inclusion in MB-System: * Copyright (c) 1994-2009 by * David W. Caress (caress@mbari.org) * Monterey Bay Aquarium Research Institute * Moss Landing, CA 95039 * and Dale N. Chayes (dale@ldeo.columbia.edu) * Lamont-Doherty Earth Observatory * Palisades, NY 10964 * * Algorithm and original code: * Copyright (c) 1991 by P. Wessel and W. H. F. Smith * * See README file for copying and redistribution conditions. *--------------------------------------------------------------------*/ /* * SURFUNC is a function for gridding data using a minimum curvature * algorithm developed by W.H.F. Smith and P. Wessel. The source * code below is almost entirely taken from the source code * "surface.c" distributed as part of the GMT-system by Wessel * and Smith. The MB-System Copyright notice above applies only * to the inclusion of this code in MB-System and changes made to * the code as part of that inclusion. The algorithm and the * bulk of the code remains Copyright (c) by P. Wessel and W. H. F. Smith. * *--------------------------------------------------------------------* * * The original Smith and Wessel comments follow: * * surface.c: a gridding program. * reads xyz triples and fits a surface to the data. * surface satisfies (1 - T) D4 z - T D2 z = 0, * where D4 is the biharmonic operator, * D2 is the Laplacian, * and T is a "tension factor" between 0 and 1. * End member T = 0 is the classical minimum curvature * surface. T = 1 gives a harmonic surface. Use T = 0.25 * or so for potential data; something more for topography. * * Program includes overrelaxation for fast convergence and * automatic optimal grid factorization. * * See Smith & Wessel (Geophysics, 3, 293-305, 1990) for details. * * Authors: Walter H. F. Smith and Paul Wessel * Date: April, 1988. * *--------------------------------------------------------------------* * * Particulars regarding turning the program "surface" version 4.3 * (revision of 26 February, 1992) into a function "surfunc" follow: * * Author: D. W. Caress * Date: May 2, 1994 * * $Log: mb_surface.c,v $ * Revision 5.4 2008/07/10 06:43:40 caress * Preparing for 5.1.1beta20 * * Revision 5.3 2006/08/09 22:41:27 caress * Fixed programs that read or write grids so that they do not use the GMT_begin() function; these programs will now work when GMT is built in the default fashion, when GMT is built in the default fashion, with "advisory file locking" enabled. * * Revision 5.2 2005/03/25 04:09:53 caress * Problems with global variables in mb_surface.c stomping on similarly named global variables in some programs has been fixed. * * Revision 5.1 2003/04/29 20:27:48 caress * Fixed multiple definitions of "error". * * Revision 5.0 2003/03/22 03:10:36 caress * Reinserting this code into MB-System for first time in years. * * Revision 4.2 1994/10/21 13:02:31 caress * Release V4.0 * * Revision 4.1 1994/06/04 02:02:01 caress * Fixed several bugs and made some stylistic changes to * the output. Changed the data input bounds to be much * larger than the working grid bounds. * * Revision 4.0 1994/05/05 20:30:06 caress * First cut. This code derived from GMT program surface by * Walter Smith and Paul Wessel. * * */ /* standard include files */ #include #include #include /* mbio include files */ #include "../../include/mb_status.h" #include "../../include/mb_define.h" #include "../../include/mb_aux.h" #define GMT_CHUNK 2000 #ifndef FALSE #define FALSE 0 #endif #ifndef TRUE #define TRUE 1 #endif #define CNULL ((char *)NULL) #define OUTSIDE 2000000000 /* Index number indicating data is outside useable area */ struct MB_SURFACE_DATA { float x; float y; float z; int index; }; struct MB_SURFACE_BRIGGS { double b[6]; }; static int npoints=0; /* Number of data points */ static int nx=0; /* Number of nodes in x-dir. */ static int ny=0; /* Number of nodes in y-dir. (Final grid) */ static int mx = 0; static int my = 0; static int ij_sw_corner, ij_se_corner, ij_nw_corner, ij_ne_corner; static int block_nx; /* Number of nodes in x-dir for a given grid factor */ static int block_ny; /* Number of nodes in y-dir for a given grid factor */ static int max_iterations=250; /* Max iter per call to iterate */ static int total_iterations = 0; static int grid, old_grid; /* Node spacings */ static int grid_east; static int n_fact = 0; /* Number of factors in common (ny-1, nx-1) */ static int factors[32]; /* Array of common factors */ static int local_verbose = FALSE; static int local_error = MB_ERROR_NO_ERROR; static int status = MB_SUCCESS; static int n_empty; /* No of unconstrained nodes at initialization */ static int set_low = 0; /* 0 unconstrained,1 = by min data value, 2 = by user value */ static int set_high = 0; /* 0 unconstrained,1 = by max data value, 2 = by user value */ static int constrained = FALSE; /* TRUE if set_low or set_high is TRUE */ static double low_limit, high_limit; /* Constrains on range of solution */ static double xmin, xmax, ymin, ymax; /* minmax coordinates */ static float *lower, *upper; /* arrays for minmax values, if set */ static double xinc, yinc; /* Size of each grid cell (final size) */ static double grid_xinc, grid_yinc; /* size of each grid cell for a given grid factor */ static double r_xinc, r_yinc, r_grid_xinc, r_grid_yinc; /* Reciprocals */ static double converge_limit = 0.0; /* Convergence limit */ static double radius = 0.0; /* Search radius for initializing grid */ static double tension = 0.00; /* Tension parameter on the surface */ static double boundary_tension = 0.0; static double interior_tension = 0.0; static double a0_const_1, a0_const_2; /* Constants for off grid point equation */ static double e_2, e_m2, one_plus_e2; static double eps_p2, eps_m2, two_plus_ep2, two_plus_em2; static double x_edge_const, y_edge_const; static double epsilon = 1.0; static double z_mean; static double z_scale = 1.0; /* Root mean square range of z after removing planar trend */ static double r_z_scale = 1.0; /* reciprocal of z_scale */ static double plane_c0, plane_c1, plane_c2; /* Coefficients of best fitting plane to data */ static double small; /* Let data point coincide with node if distance < small */ static float *u; /* Pointer to grid array */ static char *iu; /* Pointer to grid info array */ static char mode_type[2] = {'I','D'}; /* D means include data points when iterating * I means just interpolate from larger grid */ static int offset[25][12]; /* Indices of 12 nearby points in 25 cases of edge conditions */ static double coeff[2][12]; /* Coefficients for 12 nearby points, constrained and unconstrained */ static double relax_old, relax_new = 1.4; /* Coefficients for relaxation factor to speed up convergence */ static int compare_points(); static struct MB_SURFACE_DATA *data; /* Data point and index to node it currently constrains */ static struct MB_SURFACE_BRIGGS *briggs; /* Coefficients in Taylor series for Laplacian(z) a la I. C. Briggs (1974) */ /* function prototypes */ int mb_surface(int verbose, int ndat, float *xdat, float *ydat, float *zdat, double xxmin, double xxmax, double yymin, double yymax, double xxinc, double yyinc, double ttension, float *sgrid); void set_coefficients(); void set_offset(); void fill_in_forecast (); int compare_points (struct MB_SURFACE_DATA *point_1, struct MB_SURFACE_DATA *point_2); void smart_divide (); void set_index (); void find_nearest_point(); void set_grid_parameters(); void initialize_grid(); void new_initialize_grid(); void read_data(int ndat, float *xdat, float *ydat, float *zdat); void get_output(float *sgrid); int iterate(int mode); void check_errors (); int remove_planar_trend(); int replace_planar_trend(); int throw_away_unusables(); int rescale_z_values(); void load_constraints (char *low, char *high); double guess_surface_time(int nx, int ny); int get_prime_factors(int n, int f[]); int gcd_euclid(int a,int b); static char rcs_id[] = "$Id: mb_surface.c 1770 2009-10-19 17:16:39Z caress $"; int mb_surface(int verbose, int ndat, float *xdat, float *ydat, float *zdat, double xxmin, double xxmax, double yymin, double yymax, double xxinc, double yyinc, double ttension, float *sgrid) { int i, serror = FALSE; char low[100], high[100]; char *function_name = "mb_surface"; /* print input debug statements */ if (verbose >= 2) { fprintf(stderr,"\ndbg2 MBBA function <%s> called\n",function_name); fprintf(stderr,"dbg2 Revision id: %s\n",rcs_id); fprintf(stderr,"dbg2 Input arguments:\n"); fprintf(stderr,"dbg2 verbose: %d\n",verbose); fprintf(stderr,"dbg2 xxmin: %f\n",xxmin); fprintf(stderr,"dbg2 xxmax: %f\n",xxmax); fprintf(stderr,"dbg2 yymin: %f\n",yymin); fprintf(stderr,"dbg2 yymax: %f\n",yymax); fprintf(stderr,"dbg2 xxinc: %f\n",xxinc); fprintf(stderr,"dbg2 xxinc: %f\n",xxinc); fprintf(stderr,"dbg2 ttension: %f\n",ttension); fprintf(stderr,"dbg2 ndat: %d\n",ndat); for (i=0;i 0) local_verbose = TRUE; else local_verbose = FALSE; /* New in v4.3: Default to unconstrained: */ set_low = set_high = 0; if (xmin >= xmax || ymin >= ymax) serror = TRUE; if (xinc <= 0.0 || yinc <= 0.0) serror = TRUE; if (tension != 0.0) { boundary_tension = tension; interior_tension = tension; } relax_old = 1.0 - relax_new; nx = rint((xmax - xmin)/xinc) + 1; ny = rint((ymax - ymin)/yinc) + 1; mx = nx + 4; my = ny + 4; r_xinc = 1.0 / xinc; r_yinc = 1.0 / yinc; /* New stuff here for v4.3: Check out the grid dimensions: */ grid = gcd_euclid(nx-1, ny-1); /* if (local_verbose || size_query || grid == 1) fprintf (stderr, "W: %.3lf E: %.3lf S: %.3lf N: %.3lf nx: %d ny: %d\n", xmin, xmax, ymin, ymax, nx, ny); if (grid == 1) fprintf(stderr,"surface: WARNING: Your grid dimensions are mutually prime.\n"); if (grid == 1 || size_query) suggest_sizes_for_surface(nx-1, ny-1); if (size_query) exit(0); */ /* New idea: set grid = 1, read data, setting index. Then throw away data that can't be used in end game, constraining size of briggs->b[6] structure. */ grid = 1; set_grid_parameters(); read_data(ndat,xdat,ydat,zdat); throw_away_unusables(); remove_planar_trend(); rescale_z_values(); load_constraints(low, high); /* Set up factors and reset grid to first value */ grid = gcd_euclid(nx-1, ny-1); n_fact = get_prime_factors(grid, factors); set_grid_parameters(); while ( block_nx < 4 || block_ny < 4 ) { smart_divide(); set_grid_parameters(); } set_offset(); set_index(); /* Now the data are ready to go for the first iteration. */ /* Allocate more space */ status = mb_mallocd(local_verbose, __FILE__, __LINE__, npoints * sizeof(struct MB_SURFACE_BRIGGS), (void **)&briggs, &local_error); status = mb_mallocd(local_verbose, __FILE__, __LINE__, mx * my * sizeof(char), (void **)&iu, &local_error); status = mb_mallocd(local_verbose, __FILE__, __LINE__, mx * my * sizeof(float), (void **)&u, &local_error); if (radius > 0) initialize_grid(); /* Fill in nodes with a weighted avg in a search radius */ /* if (local_verbose) fprintf(stderr,"Grid\tMode\tIteration\tMax Change\tConv Limit\tTotal Iterations\n"); */ set_coefficients(); old_grid = grid; find_nearest_point (); iterate (1); while (grid > 1) { smart_divide (); set_grid_parameters(); set_offset(); set_index (); fill_in_forecast (); iterate(0); old_grid = grid; find_nearest_point (); iterate (1); } if (local_verbose) check_errors (); replace_planar_trend(); get_output(sgrid); status = mb_freed(verbose,__FILE__, __LINE__, (void **)&data, &local_error); status = mb_freed(verbose,__FILE__, __LINE__, (void **)&briggs, &local_error); status = mb_freed(verbose,__FILE__, __LINE__, (void **)&iu, &local_error); status = mb_freed(verbose,__FILE__, __LINE__, (void **)&u, &local_error); if (set_low) status = mb_freed(verbose,__FILE__, __LINE__, (void **)&lower, &local_error); if (set_high) status = mb_freed(verbose,__FILE__, __LINE__, (void **)&upper, &local_error); /* print output debug statements */ if (verbose >= 2) { fprintf(stderr,"\ndbg2 MBIO function <%s> completed\n",function_name); fprintf(stderr,"dbg2 Return values:\n"); fprintf(stderr,"dbg2 error: %d\n",local_error); for (i=0;iindex; index_2 = point_2->index; if (index_1 < index_2) return (-1); else if (index_1 > index_2) return (1); else if (index_1 == OUTSIDE) return (0); else { /* Points are in same grid cell, find the one who is nearest to grid point */ block_i = point_1->index/block_ny; block_j = point_1->index%block_ny; x0 = xmin + block_i * grid_xinc; y0 = ymin + block_j * grid_yinc; dist_1 = (point_1->x - x0) * (point_1->x - x0) + (point_1->y - y0) * (point_1->y - y0); dist_2 = (point_2->x - x0) * (point_2->x - x0) + (point_2->y - y0) * (point_2->y - y0); if (dist_1 < dist_2) return (-1); if (dist_1 > dist_2) return (1); else return (0); } } void smart_divide () { /* Divide grid by its largest prime factor */ grid /= factors[n_fact - 1]; n_fact--; } void set_index () { /* recomputes data[k].index for new value of grid, sorts data on index and radii, and throws away data which are now outside the useable limits. */ int i, j, k, k_skipped = 0; for (k = 0; k < npoints; k++) { i = floor(((data[k].x-xmin)*r_grid_xinc) + 0.5); j = floor(((data[k].y-ymin)*r_grid_yinc) + 0.5); if (i < 0 || i >= block_nx || j < 0 || j >= block_ny) { data[k].index = OUTSIDE; k_skipped++; } else data[k].index = i * block_ny + j; } qsort ((char *)data, npoints, sizeof (struct MB_SURFACE_DATA), (void *)compare_points); npoints -= k_skipped; } void find_nearest_point() { int i, j, k, last_index, block_i, block_j, iu_index, briggs_index; double x0, y0, dx, dy, xys, xy1, btemp; double b0, b1, b2, b3, b4, b5; last_index = -1; small = 0.05 * ((grid_xinc < grid_yinc) ? grid_xinc : grid_yinc); for (i = 0; i < nx; i += grid) /* Reset grid info */ for (j = 0; j < ny; j += grid) iu[ij_sw_corner + i*my + j] = 0; briggs_index = 0; for (k = 0; k < npoints; k++) { /* Find constraining value */ if (data[k].index != last_index) { block_i = data[k].index/block_ny; block_j = data[k].index%block_ny; last_index = data[k].index; iu_index = ij_sw_corner + (block_i * my + block_j) * grid; x0 = xmin + block_i*grid_xinc; y0 = ymin + block_j*grid_yinc; dx = (data[k].x - x0)*r_grid_xinc; dy = (data[k].y - y0)*r_grid_yinc; if (fabs(dx) < small && fabs(dy) < small) { iu[iu_index] = 5; u[iu_index] = data[k].z; } else { if (dx >= 0.0) { if (dy >= 0.0) iu[iu_index] = 1; else iu[iu_index] = 4; } else { if (dy >= 0.0) iu[iu_index] = 2; else iu[iu_index] = 3; } dx = fabs(dx); dy = fabs(dy); btemp = 2 * one_plus_e2 / ( (dx + dy) * (1.0 + dx + dy) ); b0 = 1.0 - 0.5 * (dx + (dx * dx)) * btemp; b3 = 0.5 * (e_2 - (dy + (dy * dy)) * btemp); xys = 1.0 + dx + dy; xy1 = 1.0 / xys; b1 = (e_2 * xys - 4 * dy) * xy1; b2 = 2 * (dy - dx + 1.0) * xy1; b4 = b0 + b1 + b2 + b3 + btemp; b5 = btemp * data[k].z; briggs[briggs_index].b[0] = b0; briggs[briggs_index].b[1] = b1; briggs[briggs_index].b[2] = b2; briggs[briggs_index].b[3] = b3; briggs[briggs_index].b[4] = b4; briggs[briggs_index].b[5] = b5; briggs_index++; } } } } void set_grid_parameters() { block_ny = (ny - 1) / grid + 1; block_nx = (nx - 1) / grid + 1; grid_xinc = grid * xinc; grid_yinc = grid * yinc; grid_east = grid * my; r_grid_xinc = 1.0 / grid_xinc; r_grid_yinc = 1.0 / grid_yinc; } void initialize_grid() { /* * For the initial gridsize, compute weighted averages of data inside the search radius * and assign the values to u[i,j] where i,j are multiples of gridsize. */ int irad, jrad, i, j, imin, imax, jmin, jmax, index_1, index_2, k, ki, kj, k_index; double r, rfact, sum_w, sum_zw, weight, x0, y0; irad = ceil(radius/grid_xinc); jrad = ceil(radius/grid_yinc); rfact = -4.5/(radius*radius); for (i = 0; i < block_nx; i ++ ) { x0 = xmin + i*grid_xinc; for (j = 0; j < block_ny; j ++ ) { y0 = ymin + j*grid_yinc; imin = i - irad; if (imin < 0) imin = 0; imax = i + irad; if (imax >= block_nx) imax = block_nx - 1; jmin = j - jrad; if (jmin < 0) jmin = 0; jmax = j + jrad; if (jmax >= block_ny) jmax = block_ny - 1; index_1 = imin*block_ny + jmin; index_2 = imax*block_ny + jmax + 1; sum_w = sum_zw = 0.0; k = 0; while (k < npoints && data[k].index < index_1) k++; for (ki = imin; k < npoints && ki <= imax && data[k].index < index_2; ki++) { for (kj = jmin; k < npoints && kj <= jmax && data[k].index < index_2; kj++) { k_index = ki*block_ny + kj; while (k < npoints && data[k].index < k_index) k++; while (k < npoints && data[k].index == k_index) { r = (data[k].x-x0)*(data[k].x-x0) + (data[k].y-y0)*(data[k].y-y0); weight = exp (rfact*r); sum_w += weight; sum_zw += weight*data[k].z; k++; } } } if (sum_w == 0.0) { /* fprintf (stderr, "surface: Warning: no data inside search radius at: %.8lg %.8lg\n", x0, y0); */ u[ij_sw_corner + (i * my + j) * grid] = z_mean; } else { u[ij_sw_corner + (i*my+j)*grid] = sum_zw/sum_w; } } } } void new_initialize_grid() { /* * For the initial gridsize, load constrained nodes with weighted avg of their data; * and then do something with the unconstrained ones. */ int k, k_index, u_index, block_i, block_j; double sum_w, sum_zw, weight, x0, y0, dx, dy, dx_scale, dy_scale; dx_scale = 4.0 / grid_xinc; dy_scale = 4.0 / grid_yinc; n_empty = block_ny * block_nx; k = 0; while (k < npoints) { block_i = data[k].index / block_ny; block_j = data[k].index % block_ny; x0 = xmin + block_i*grid_xinc; y0 = ymin + block_j*grid_yinc; u_index = ij_sw_corner + (block_i*my + block_j) * grid; k_index = data[k].index; dy = (data[k].y - y0) * dy_scale; dx = (data[k].x - x0) * dx_scale; sum_w = 1.0 / (1.0 + dx*dx + dy*dy); sum_zw = data[k].z * sum_w; k++; while (k < npoints && data[k].index == k_index) { dy = (data[k].y - y0) * dy_scale; dx = (data[k].x - x0) * dx_scale; weight = 1.0 / (1.0 + dx*dx + dy*dy); sum_zw += data[k].z * weight; sum_w += weight; sum_zw += weight*data[k].z; k++; } u[u_index] = sum_zw/sum_w; iu[u_index] = 5; n_empty--; } } /* This function rewritten by D.W. Caress 5/3/94 */ void read_data(int ndat, float *xdat, float *ydat, float *zdat) { int i, j, k, kmax, kmin, idat; double zmin = 1.0e38, zmax = -1.0e38; status = mb_mallocd(local_verbose, __FILE__, __LINE__, ndat * sizeof(struct MB_SURFACE_DATA), (void **)&data, &local_error); /* Read in xyz data and computes index no and store it in a structure */ k = 0; z_mean = 0; for (idat=0;idat= 0 && i < block_nx && j >= 0 && j < block_ny) { data[k].index = i * block_ny + j; data[k].x = xdat[idat]; data[k].y = ydat[idat]; data[k].z = zdat[idat]; if (zmin > zdat[idat]) { zmin = zdat[idat]; kmin = k; } if (zmax < zdat[idat]) { zmax = zdat[idat]; kmax = k; } k++; z_mean += zdat[idat]; } } npoints = k; z_mean /= k; if( converge_limit == 0.0 ) { converge_limit = 0.001 * z_scale; /* c_l = 1 ppt of L2 scale */ } /* if (local_verbose) { fprintf(stderr, "surface: Minimum value of your dataset x,y,z at: %g %g %g\n", data[kmin].x, data[kmin].y, data[kmin].z); fprintf(stderr, "surface: Maximum value of your dataset x,y,z at: %g %g %g\n", data[kmax].x, data[kmax].y, data[kmax].z); } */ if (set_low == 1) low_limit = data[kmin].z; else if (set_low == 2 && low_limit > data[kmin].z) { /* low_limit = data[kmin].z; */ /* fprintf (stderr, "surface: Warning: Your lower value is > than min data value.\n"); */ } if (set_high == 1) high_limit = data[kmax].z; else if (set_high == 2 && high_limit < data[kmax].z) { /* high_limit = data[kmax].z; */ /* fprintf (stderr, "surface: Warning: Your upper value is < than max data value.\n"); */ } } /* this function rewritten from write_output() by D.W. Caress 5/3/94 */ void get_output(float *sgrid) { int index, i, j; index = ij_sw_corner; for(i = 0; i < nx; i++, index += my) for (j = 0; j < ny; j++) { sgrid[j*nx+i] = u[index + ny - j - 1]; } } int iterate(int mode) { int i, j, k, ij, kase, briggs_index, ij_v2; int x_case, y_case, x_w_case, x_e_case, y_s_case, y_n_case; int iteration_count = 0; double current_limit = converge_limit / grid; double change, max_change = 0.0, busum, sum_ij; double b0, b1, b2, b3, b4, b5; double x_0_const = 4.0 * (1.0 - boundary_tension) / (2.0 - boundary_tension); double x_1_const = (3 * boundary_tension - 2.0) / (2.0 - boundary_tension); double y_denom = 2 * epsilon * (1.0 - boundary_tension) + boundary_tension; double y_0_const = 4 * epsilon * (1.0 - boundary_tension) / y_denom; double y_1_const = (boundary_tension - 2 * epsilon * (1.0 - boundary_tension) ) / y_denom; do { briggs_index = 0; /* Reset the constraint table stack pointer */ max_change = -1.0; /* Fill in auxiliary boundary values (in new way) */ /* First set d2[]/dn2 = 0 along edges: */ /* New experiment : (1-T)d2[]/dn2 + Td[]/dn = 0 */ for (i = 0; i < nx; i += grid) { /* set d2[]/dy2 = 0 on south side: */ ij = ij_sw_corner + i * my; /* u[ij - 1] = 2 * u[ij] - u[ij + grid]; */ u[ij - 1] = y_0_const * u[ij] + y_1_const * u[ij + grid]; /* set d2[]/dy2 = 0 on north side: */ ij = ij_nw_corner + i * my; /* u[ij + 1] = 2 * u[ij] - u[ij - grid]; */ u[ij + 1] = y_0_const * u[ij] + y_1_const * u[ij - grid]; } for (j = 0; j < ny; j += grid) { /* set d2[]/dx2 = 0 on west side: */ ij = ij_sw_corner + j; /* u[ij - my] = 2 * u[ij] - u[ij + grid_east]; */ u[ij - my] = x_1_const * u[ij + grid_east] + x_0_const * u[ij]; /* set d2[]/dx2 = 0 on east side: */ ij = ij_se_corner + j; /* u[ij + my] = 2 * u[ij] - u[ij - grid_east]; */ u[ij + my] = x_1_const * u[ij - grid_east] + x_0_const * u[ij]; } /* Now set d2[]/dxdy = 0 at each corner: */ ij = ij_sw_corner; u[ij - my - 1] = u[ij + grid_east - 1] + u[ij - my + grid] - u[ij + grid_east + grid]; ij = ij_nw_corner; u[ij - my + 1] = u[ij + grid_east + 1] + u[ij - my - grid] - u[ij + grid_east - grid]; ij = ij_se_corner; u[ij + my - 1] = u[ij - grid_east - 1] + u[ij + my + grid] - u[ij - grid_east + grid]; ij = ij_ne_corner; u[ij + my + 1] = u[ij - grid_east + 1] + u[ij + my - grid] - u[ij - grid_east - grid]; /* Now set (1-T)dC/dn + Tdu/dn = 0 at each edge : */ /* New experiment: only dC/dn = 0 */ x_w_case = 0; x_e_case = block_nx - 1; for (i = 0; i < nx; i += grid, x_w_case++, x_e_case--) { if(x_w_case < 2) x_case = x_w_case; else if(x_e_case < 2) x_case = 4 - x_e_case; else x_case = 2; /* South side : */ kase = x_case * 5; ij = ij_sw_corner + i * my; u[ij + offset[kase][11]] = (u[ij + offset[kase][0]] + eps_m2*(u[ij + offset[kase][1]] + u[ij + offset[kase][3]] - u[ij + offset[kase][8]] - u[ij + offset[kase][10]]) + two_plus_em2 * (u[ij + offset[kase][9]] - u[ij + offset[kase][2]]) ); /* + tense * eps_m2 * (u[ij + offset[kase][2]] - u[ij + offset[kase][9]]) / (1.0 - tense); */ /* North side : */ kase = x_case * 5 + 4; ij = ij_nw_corner + i * my; u[ij + offset[kase][0]] = -(-u[ij + offset[kase][11]] + eps_m2 * (u[ij + offset[kase][1]] + u[ij + offset[kase][3]] - u[ij + offset[kase][8]] - u[ij + offset[kase][10]]) + two_plus_em2 * (u[ij + offset[kase][9]] - u[ij + offset[kase][2]]) ); /* - tense * eps_m2 * (u[ij + offset[kase][2]] - u[ij + offset[kase][9]]) / (1.0 - tense); */ } y_s_case = 0; y_n_case = block_ny - 1; for (j = 0; j < ny; j += grid, y_s_case++, y_n_case--) { if(y_s_case < 2) y_case = y_s_case; else if(y_n_case < 2) y_case = 4 - y_n_case; else y_case = 2; /* West side : */ kase = y_case; ij = ij_sw_corner + j; u[ij+offset[kase][4]] = u[ij + offset[kase][7]] + eps_p2 * (u[ij + offset[kase][3]] + u[ij + offset[kase][10]] -u[ij + offset[kase][1]] - u[ij + offset[kase][8]]) + two_plus_ep2 * (u[ij + offset[kase][5]] - u[ij + offset[kase][6]]); /* + tense * (u[ij + offset[kase][6]] - u[ij + offset[kase][5]]) / (1.0 - tense); */ /* East side : */ kase = 20 + y_case; ij = ij_se_corner + j; u[ij + offset[kase][7]] = - (-u[ij + offset[kase][4]] + eps_p2 * (u[ij + offset[kase][3]] + u[ij + offset[kase][10]] - u[ij + offset[kase][1]] - u[ij + offset[kase][8]]) + two_plus_ep2 * (u[ij + offset[kase][5]] - u[ij + offset[kase][6]]) ); /* - tense * (u[ij + offset[kase][6]] - u[ij + offset[kase][5]]) / (1.0 - tense); */ } /* That's it for the boundary points. Now loop over all data */ x_w_case = 0; x_e_case = block_nx - 1; for (i = 0; i < nx; i += grid, x_w_case++, x_e_case--) { if(x_w_case < 2) x_case = x_w_case; else if(x_e_case < 2) x_case = 4 - x_e_case; else x_case = 2; y_s_case = 0; y_n_case = block_ny - 1; ij = ij_sw_corner + i * my; for (j = 0; j < ny; j += grid, ij += grid, y_s_case++, y_n_case--) { if (iu[ij] == 5) continue; /* Point is fixed */ if(y_s_case < 2) y_case = y_s_case; else if(y_n_case < 2) y_case = 4 - y_n_case; else y_case = 2; kase = x_case * 5 + y_case; sum_ij = 0.0; if (iu[ij] == 0) { /* Point is unconstrained */ for (k = 0; k < 12; k++) { sum_ij += (u[ij + offset[kase][k]] * coeff[0][k]); } } else { /* Point is constrained */ b0 = briggs[briggs_index].b[0]; b1 = briggs[briggs_index].b[1]; b2 = briggs[briggs_index].b[2]; b3 = briggs[briggs_index].b[3]; b4 = briggs[briggs_index].b[4]; b5 = briggs[briggs_index].b[5]; briggs_index++; if (iu[ij] < 3) { if (iu[ij] == 1) { /* Point is in quadrant 1 */ busum = b0 * u[ij + offset[kase][10]] + b1 * u[ij + offset[kase][9]] + b2 * u[ij + offset[kase][5]] + b3 * u[ij + offset[kase][1]]; } else { /* Point is in quadrant 2 */ busum = b0 * u[ij + offset[kase][8]] + b1 * u[ij + offset[kase][9]] + b2 * u[ij + offset[kase][6]] + b3 * u[ij + offset[kase][3]]; } } else { if (iu[ij] == 3) { /* Point is in quadrant 3 */ busum = b0 * u[ij + offset[kase][1]] + b1 * u[ij + offset[kase][2]] + b2 * u[ij + offset[kase][6]] + b3 * u[ij + offset[kase][10]]; } else { /* Point is in quadrant 4 */ busum = b0 * u[ij + offset[kase][3]] + b1 * u[ij + offset[kase][2]] + b2 * u[ij + offset[kase][5]] + b3 * u[ij + offset[kase][8]]; } } for (k = 0; k < 12; k++) { sum_ij += (u[ij + offset[kase][k]] * coeff[1][k]); } sum_ij = (sum_ij + a0_const_2 * (busum + b5)) / (a0_const_1 + a0_const_2 * b4); } /* New relaxation here */ sum_ij = u[ij] * relax_old + sum_ij * relax_new; if (constrained) { /* Must check limits. Note lower/upper is v2 format and need ij_v2! */ ij_v2 = (ny - j - 1) * nx + i; if (set_low /*&& !GMT_is_fnan((double)lower[ij_v2])*/ && sum_ij < lower[ij_v2]) sum_ij = lower[ij_v2]; else if (set_high /*&& !GMT_is_fnan((double)upper[ij_v2])*/ && sum_ij > upper[ij_v2]) sum_ij = upper[ij_v2]; } change = fabs(sum_ij - u[ij]); u[ij] = sum_ij; if (change > max_change) max_change = change; } } iteration_count++; total_iterations++; max_change *= z_scale; /* Put max_change into z units */ if (local_verbose > 1) fprintf(stderr,"%4d\t%c\t%8d\t%10lg\t%10lg\t%10d\n", grid, mode_type[mode], iteration_count, max_change, current_limit, total_iterations); } while (max_change > current_limit && iteration_count < max_iterations); if (local_verbose) fprintf(stderr,"%4d\t%c\t%8d\t%10lg\t%10lg\t%10d\n", grid, mode_type[mode], iteration_count, max_change, current_limit, total_iterations); return(iteration_count); } void check_errors () { int i, j, k, ij, n_nodes, move_over[12]; /* move_over = offset[kase][12], but grid = 1 so move_over is easy */ double x0, y0, dx, dy, mean_error, mean_squared_error, z_est, z_err, curvature, c; double du_dx, du_dy, d2u_dx2, d2u_dxdy, d2u_dy2, d3u_dx3, d3u_dx2dy, d3u_dxdy2, d3u_dy3; double x_0_const = 4.0 * (1.0 - boundary_tension) / (2.0 - boundary_tension); double x_1_const = (3 * boundary_tension - 2.0) / (2.0 - boundary_tension); double y_denom = 2 * epsilon * (1.0 - boundary_tension) + boundary_tension; double y_0_const = 4 * epsilon * (1.0 - boundary_tension) / y_denom; double y_1_const = (boundary_tension - 2 * epsilon * (1.0 - boundary_tension) ) / y_denom; move_over[0] = 2; move_over[1] = 1 - my; move_over[2] = 1; move_over[3] = 1 + my; move_over[4] = -2 * my; move_over[5] = -my; move_over[6] = my; move_over[7] = 2 * my; move_over[8] = -1 - my; move_over[9] = -1; move_over[10] = -1 + my; move_over[11] = -2; mean_error = 0; mean_squared_error = 0; /* First update the boundary values */ for (i = 0; i < nx; i ++) { ij = ij_sw_corner + i * my; u[ij - 1] = y_0_const * u[ij] + y_1_const * u[ij + 1]; ij = ij_nw_corner + i * my; u[ij + 1] = y_0_const * u[ij] + y_1_const * u[ij - 1]; } for (j = 0; j < ny; j ++) { ij = ij_sw_corner + j; u[ij - my] = x_1_const * u[ij + my] + x_0_const * u[ij]; ij = ij_se_corner + j; u[ij + my] = x_1_const * u[ij - my] + x_0_const * u[ij]; } ij = ij_sw_corner; u[ij - my - 1] = u[ij + my - 1] + u[ij - my + 1] - u[ij + my + 1]; ij = ij_nw_corner; u[ij - my + 1] = u[ij + my + 1] + u[ij - my - 1] - u[ij + my - 1]; ij = ij_se_corner; u[ij + my - 1] = u[ij - my - 1] + u[ij + my + 1] - u[ij - my + 1]; ij = ij_ne_corner; u[ij + my + 1] = u[ij - my + 1] + u[ij + my - 1] - u[ij - my - 1]; for (i = 0; i < nx; i ++) { ij = ij_sw_corner + i * my; u[ij + move_over[11]] = (u[ij + move_over[0]] + eps_m2*(u[ij + move_over[1]] + u[ij + move_over[3]] - u[ij + move_over[8]] - u[ij + move_over[10]]) + two_plus_em2 * (u[ij + move_over[9]] - u[ij + move_over[2]]) ); ij = ij_nw_corner + i * my; u[ij + move_over[0]] = -(-u[ij + move_over[11]] + eps_m2 * (u[ij + move_over[1]] + u[ij + move_over[3]] - u[ij + move_over[8]] - u[ij + move_over[10]]) + two_plus_em2 * (u[ij + move_over[9]] - u[ij + move_over[2]]) ); } for (j = 0; j < ny; j ++) { ij = ij_sw_corner + j; u[ij+move_over[4]] = u[ij + move_over[7]] + eps_p2 * (u[ij + move_over[3]] + u[ij + move_over[10]] -u[ij + move_over[1]] - u[ij + move_over[8]]) + two_plus_ep2 * (u[ij + move_over[5]] - u[ij + move_over[6]]); ij = ij_se_corner + j; u[ij + move_over[7]] = - (-u[ij + move_over[4]] + eps_p2 * (u[ij + move_over[3]] + u[ij + move_over[10]] - u[ij + move_over[1]] - u[ij + move_over[8]]) + two_plus_ep2 * (u[ij + move_over[5]] - u[ij + move_over[6]]) ); } /* That resets the boundary values. Now we can test all data. Note that this loop checks all values, even though only nearest were used. */ for (k = 0; k < npoints; k++) { i = data[k].index/ny; j = data[k].index%ny; ij = ij_sw_corner + i * my + j; if ( iu[ij] == 5 ) continue; x0 = xmin + i*xinc; y0 = ymin + j*yinc; dx = (data[k].x - x0)*r_xinc; dy = (data[k].y - y0)*r_yinc; du_dx = 0.5 * (u[ij + move_over[6]] - u[ij + move_over[5]]); du_dy = 0.5 * (u[ij + move_over[2]] - u[ij + move_over[9]]); d2u_dx2 = u[ij + move_over[6]] + u[ij + move_over[5]] - 2 * u[ij]; d2u_dy2 = u[ij + move_over[2]] + u[ij + move_over[9]] - 2 * u[ij]; d2u_dxdy = 0.25 * (u[ij + move_over[3]] - u[ij + move_over[1]] - u[ij + move_over[10]] + u[ij + move_over[8]]); d3u_dx3 = 0.5 * ( u[ij + move_over[7]] - 2 * u[ij + move_over[6]] + 2 * u[ij + move_over[5]] - u[ij + move_over[4]]); d3u_dy3 = 0.5 * ( u[ij + move_over[0]] - 2 * u[ij + move_over[2]] + 2 * u[ij + move_over[9]] - u[ij + move_over[11]]); d3u_dx2dy = 0.5 * ( ( u[ij + move_over[3]] + u[ij + move_over[1]] - 2 * u[ij + move_over[2]] ) - ( u[ij + move_over[10]] + u[ij + move_over[8]] - 2 * u[ij + move_over[9]] ) ); d3u_dxdy2 = 0.5 * ( ( u[ij + move_over[3]] + u[ij + move_over[10]] - 2 * u[ij + move_over[6]] ) - ( u[ij + move_over[1]] + u[ij + move_over[8]] - 2 * u[ij + move_over[5]] ) ); /* 3rd order Taylor approx: */ z_est = u[ij] + dx * (du_dx + dx * ( (0.5 * d2u_dx2) + dx * (d3u_dx3 / 6.0) ) ) + dy * (du_dy + dy * ( (0.5 * d2u_dy2) + dy * (d3u_dy3 / 6.0) ) ) + dx * dy * (d2u_dxdy) + (0.5 * dx * d3u_dx2dy) + (0.5 * dy * d3u_dxdy2); z_err = z_est - data[k].z; mean_error += z_err; mean_squared_error += (z_err * z_err); } mean_error /= npoints; mean_squared_error = sqrt( mean_squared_error / npoints); curvature = 0.0; n_nodes = nx * ny; for (i = 0; i < nx; i++) { for (j = 0; j < ny; j++) { ij = ij_sw_corner + i * my + j; c = u[ij + move_over[6]] + u[ij + move_over[5]] + u[ij + move_over[2]] + u[ij + move_over[9]] - 4.0 * u[ij + move_over[6]]; curvature += (c * c); } } /* fprintf(stderr, "Fit info: N data points N nodes\tmean error\trms error\tcurvature\n"); fprintf (stderr,"\t%8d\t%8d\t%.8lg\t%.8lg\t%.8lg\n", npoints, n_nodes, mean_error, mean_squared_error, curvature); */ if (local_verbose) { fprintf(stderr,"\nSpline interpolation fit information:\n"); fprintf(stderr,"Data points nodes mean error rms error curvature\n"); fprintf(stderr,"%9d %9d %10g %10g %10g\n", npoints, n_nodes, mean_error, mean_squared_error, curvature); } } int remove_planar_trend() { int i; double a, b, c, d, xx, yy, zz; double sx, sy, sz, sxx, sxy, sxz, syy, syz; sx = sy = sz = sxx = sxy = sxz = syy = syz = 0.0; for (i = 0; i < npoints; i++) { xx = (data[i].x - xmin) * r_xinc; yy = (data[i].y - ymin) * r_yinc; zz = data[i].z; sx += xx; sy += yy; sz += zz; sxx +=(xx * xx); sxy +=(xx * yy); sxz +=(xx * zz); syy +=(yy * yy); syz +=(yy * zz); } d = npoints*sxx*syy + 2*sx*sy*sxy - npoints*sxy*sxy - sx*sx*syy - sy*sy*sxx; if (d == 0.0) { plane_c0 = plane_c1 = plane_c2 = 0.0; return(0); } a = sz*sxx*syy + sx*sxy*syz + sy*sxy*sxz - sz*sxy*sxy - sx*sxz*syy - sy*syz*sxx; b = npoints*sxz*syy + sz*sy*sxy + sy*sx*syz - npoints*sxy*syz - sz*sx*syy - sy*sy*sxz; c = npoints*sxx*syz + sx*sy*sxz + sz*sx*sxy - npoints*sxy*sxz - sx*sx*syz - sz*sy*sxx; plane_c0 = a / d; plane_c1 = b / d; plane_c2 = c / d; for (i = 0; i < npoints; i++) { xx = (data[i].x - xmin) * r_xinc; yy = (data[i].y - ymin) * r_yinc; data[i].z -=(plane_c0 + plane_c1 * xx + plane_c2 * yy); } return(0); } int replace_planar_trend() { int i, j, ij; for (i = 0; i < nx; i++) { for (j = 0; j < ny; j++) { ij = ij_sw_corner + i * my + j; u[ij] = (u[ij] * z_scale) + (plane_c0 + plane_c1 * i + plane_c2 * j); } } return(0); } int throw_away_unusables() { /* This is a new routine to eliminate data which will become unusable on the final iteration, when grid = 1. It assumes grid = 1 and set_grid_parameters has been called. We sort, mark redundant data as OUTSIDE, and sort again, chopping off the excess. Experimental modification 5 Dec 1988 by Smith, as part of a new implementation using core memory for b[6] coefficients, eliminating calls to temp file. */ int last_index, n_outside, k; /* Sort the data */ qsort ((char *)data, npoints, sizeof (struct MB_SURFACE_DATA), (void *)compare_points); /* If more than one datum is indexed to same node, only the first should be kept. Mark the additional ones as OUTSIDE */ last_index = -1; n_outside = 0; for (k = 0; k < npoints; k++) { if (data[k].index == last_index) { data[k].index = OUTSIDE; n_outside++; } else { last_index = data[k].index; } } /* Sort again; this time the OUTSIDE points will be thrown away */ qsort ((char *)data, npoints, sizeof (struct MB_SURFACE_DATA), (void *)compare_points); npoints -= n_outside; status = mb_reallocd(local_verbose, __FILE__, __LINE__, npoints * sizeof(struct MB_SURFACE_DATA), (void **)&data, &local_error); if (local_verbose && (n_outside)) { fprintf(stderr,"surface: %d unusable points were supplied; these will be ignored.\n", n_outside); fprintf(stderr,"\tYou should have pre-processed the data with blockmean or blockmedian.\n"); } return(0); } int rescale_z_values() { int i; double ssz = 0.0; for (i = 0; i < npoints; i++) { ssz += (data[i].z * data[i].z); } /* Set z_scale = rms(z): */ z_scale = sqrt(ssz / npoints); r_z_scale = 1.0 / z_scale; for (i = 0; i < npoints; i++) { data[i].z *= r_z_scale; } return (0); } void load_constraints (char *low, char *high) { int i, j, ij; double yy; /* struct GRD_HEADER hdr;*/ /* Load lower/upper limits, verify range, deplane, and rescale */ if (set_low > 0) { status = mb_mallocd(local_verbose, __FILE__, __LINE__, nx * ny * sizeof(float), (void **)&lower, &local_error); if (set_low < 3) for (i = 0; i < nx * ny; i++) lower[i] = low_limit; /* Comment this out: else { if (read_grd_info (low, &hdr)) { fprintf (stderr, "surface: Error opening file %s\n", low); exit (-1); } if (hdr.nx != nx || hdr.ny != ny) { fprintf (stderr, "surface: lower limit file not of proper dimension!\n"); exit (-1); } if (read_grd (low, &hdr, lower)) { fprintf (stderr, "surface: Error reading file %s\n", low); exit (-1); } n_trimmed = 0; for (i = 0; i < nx * ny; i++) if (lower[i] > low_limit) { lower[i] = low_limit; n_trimmed++; } if (n_trimmed) fprintf (stderr, "surface: %d lower limit values > min data, reset to min data!\n"); } */ for (j = ij = 0; j < ny; j++) { yy = ny - j - 1; for (i = 0; i < nx; i++, ij++) { /*if (GMT_is_fnan ((double)lower[ij])) continue;*/ lower[ij] -= (plane_c0 + plane_c1 * i + plane_c2 * yy); lower[ij] *= r_z_scale; } } constrained = TRUE; } if (set_high > 0) { status = mb_mallocd(local_verbose, __FILE__, __LINE__, nx * ny * sizeof(float), (void **)&upper, &local_error); if (set_high < 3) for (i = 0; i < nx * ny; i++) upper[i] = high_limit; /* Comment this out: else { if (read_grd_info (high, &hdr)) { fprintf (stderr, "surface: Error opening file %s\n", high); exit (-1); } if (hdr.nx != nx || hdr.ny != ny) { fprintf (stderr, "surface: upper limit file not of proper dimension!\n"); exit (-1); } if (read_grd (high, &hdr, upper)) { fprintf (stderr, "surface: Error reading file %s\n", high); exit (-1); } n_trimmed = 0; for (i = 0; i < nx * ny; i++) if (upper[i] < high_limit) { upper[i] = high_limit; n_trimmed++; } if (n_trimmed) fprintf (stderr, "surface: %d upper limit values < max data, reset to max data!\n"); } */ for (j = ij = 0; j < ny; j++) { yy = ny - j - 1; for (i = 0; i < nx; i++, ij++) { /*if (GMT_is_fnan ((double)upper[ij])) continue;*/ upper[ij] -= (plane_c0 + plane_c1 * i + plane_c2 * yy); upper[ij] *= r_z_scale; } } constrained = TRUE; } } double guess_surface_time(int nx, int ny) { /* Routine to guess a number proportional to the operations * required by surface working on a user-desired grid of * size nx by ny, where nx = (xmax - xmin)/dx, and same for * ny. (That is, one less than actually used in routine.) * * This is based on the following untested conjecture: * The operations are proportional to T = nxg*nyg*L, * where L is a measure of the distance that data * constraints must propagate, and nxg, nyg are the * current size of the grid. * For nx,ny relatively prime, we will go through only * one grid cycle, L = max(nx,ny), and T = nx*ny*L. * But for nx,ny whose greatest common divisor is a highly * composite number, we will have L equal to the division * step made at each new grid cycle, and nxg,nyg will * also be smaller than nx,ny. Thus we can hope to find * some nx,ny for which the total value of T is small. * * The above is pure speculation and has not been derived * empirically. In actual practice, the distribution of the * data, both spatially and in terms of their values, will * have a strong effect on convergence. * * W. H. F. Smith, 26 Feb 1992. */ int gcd_euclid(); /* Finds the greatest common divisor */ int get_prime_factors(); int gcd; /* Current value of the gcd */ int nxg, nyg; /* Current value of the grid dimensions */ int nfactors; /* Number of prime factors of current gcd */ int factor; /* Currently used factor */ /* Doubles are used below, even though the values will be integers, because the multiplications might reach sizes of O(n**3) */ double t_sum; /* Sum of values of T at each grid cycle */ double length; /* Current propagation distance. */ gcd = gcd_euclid(nx, ny); if (gcd > 1) { nfactors = get_prime_factors(gcd, factors); nxg = nx/gcd; nyg = ny/gcd; if (nxg < 3 || nyg < 3) { factor = factors[nfactors - 1]; nfactors--; gcd /= factor; nxg *= factor; nyg *= factor; } } else { nxg = nx; nyg = ny; } length = MAX(nxg, nyg); t_sum = nxg * (nyg * length); /* Make it double at each multiply */ /* Are there more grid cycles ? */ while (gcd > 1) { factor = factors[nfactors - 1]; nfactors--; gcd /= factor; nxg *= factor; nyg *= factor; length = factor; t_sum += nxg * (nyg * length); } return(t_sum); } int get_prime_factors(int n, int f[]) { /* Fills the integer array f with the prime factors of n. * Returns the number of locations filled in f, which is * one if n is prime. * * f[] should have been malloc'ed to enough space before * calling prime_factors(). We can be certain that f[32] * is enough space, for if n fits in a long, then n < 2**32, * and so it must have fewer than 32 prime factors. I think * that in general, ceil(log2((double)n)) is enough storage * space for f[]. * * Tries 2,3,5 explicitly; then alternately adds 2 or 4 * to the previously tried factor to obtain the next trial * factor. This is done with the variable two_four_toggle. * With this method we try 7,11,13,17,19,23,25,29,31,35,... * up to a maximum of sqrt(n). This shortened list results * in 1/3 fewer divisions than if we simply tried all integers * between 5 and sqrt(n). We can reduce the size of the list * of trials by an additional 20% by removing the multiples * of 5, which are equal to 30m +/- 5, where m >= 1. Starting * from 25, these are found by alternately adding 10 or 20. * To do this, we use the variable ten_twenty_toggle. * * W. H. F. Smith, 26 Feb 1992, after D.E. Knuth, vol. II */ int current_factor; /* The factor currently being tried */ int max_factor; /* Don't try any factors bigger than this */ int n_factors = 0; /* Returned; one if n is prime */ int two_four_toggle = 0; /* Used to add 2 or 4 to get next trial factor */ int ten_twenty_toggle = 0; /* Used to add 10 or 20 to skip_five */ int skip_five = 25; /* Used to skip multiples of 5 in the list */ int m; /* Used to keep a working copy of n */ /* Initialize m and max_factor */ m = abs(n); if (m < 2) return(0); max_factor = floor(sqrt((double)m)); /* First find the 2s */ current_factor = 2; while(!(m%current_factor)) { m /= current_factor; f[n_factors] = current_factor; n_factors++; } if (m == 1) return(n_factors); /* Next find the 3s */ current_factor = 3; while(!(m%current_factor)) { m /= current_factor; f[n_factors] = current_factor; n_factors++; } if (m == 1) return(n_factors); /* Next find the 5s */ current_factor = 5; while(!(m%current_factor)) { m /= current_factor; f[n_factors] = current_factor; n_factors++; } if (m == 1) return(n_factors); /* Now try all the rest */ while (m > 1 && current_factor <= max_factor) { /* Current factor is either 2 or 4 more than previous value */ if (two_four_toggle) { current_factor += 4; two_four_toggle = 0; } else { current_factor += 2; two_four_toggle = 1; } /* If current factor is a multiple of 5, skip it. But first, set next value of skip_five according to 10/20 toggle: */ if (current_factor == skip_five) { if (ten_twenty_toggle) { skip_five += 20; ten_twenty_toggle = 0; } else { skip_five += 10; ten_twenty_toggle = 1; } continue; } /* Get here when current_factor is not a multiple of 2,3 or 5: */ while(!(m%current_factor)) { m /= current_factor; f[n_factors] = current_factor; n_factors++; } } /* Get here when all factors up to floor(sqrt(n)) have been tried. */ if (m > 1) { /* m is an additional prime factor of n */ f[n_factors] = m; n_factors++; } return (n_factors); } /* gcd_euclid.c Greatest common divisor routine */ #define IABS(i) (((i) < 0) ? -(i) : (i)) int gcd_euclid(int a,int b) { /* Returns the greatest common divisor of u and v by Euclid's method. * I have experimented also with Stein's method, which involves only * subtraction and left/right shifting; Euclid is faster, both for * integers of size 0 - 1024 and also for random integers of a size * which fits in a long integer. Stein's algorithm might be better * when the integers are HUGE, but for our purposes, Euclid is fine. * * Walter H. F. Smith, 25 Feb 1992, after D. E. Knuth, vol. II */ int u,v,r; u = MAX(IABS(a), IABS(b)); v = MIN(IABS(a), IABS(b)); while (v > 0) { r = u%v; /* Knuth notes that u < 2v 40% of the time; */ u = v; /* thus we could have tried a subtraction */ v = r; /* followed by an if test to do r = u%v */ } return(u); }