// Last edited on 2019-11-28 10:01:14 by jstolfi

#macro compute_drawer_heights(N, K, tot_ht, min_ht, bot_K_ht)

  // Computes heights of drawers in a stack of {N} drawers spanning
  // total height {tot_ht}.  
  
  // The heights will be increasing from top to bottom.  The bottom-most
  // {K} drawers will have combined height {bot_K_ht}.  
  // All heights will be at least {min_ht}. The heights will increase in 
  // geometric progression, except that one or more drawers at the 
  // top may be set at {min_ht} if they would otherwise be less than that. 
  
  // May be unable to honor all those constraints; in which case will
  // pick some compromise.
  
  // The result should be assigned to an array of {N} elements,
  // which are the heights of the drawers, indexed {0..N-1} from bottom to top.
  
  #debug concat("\n!! --- compute_drawer_heights( ", str(N,1,0), ", " str(K,1,0))
  #debug concat(", " str(tot_ht,1,2), ", " str(min_ht,1,2), ", " str(bot_K_ht,1,2))
  #debug concat(" ) ----------------------------------------\n")
  
  #local ht = array[N];
  #local min_tot_ht = (N-K)*min_ht + bot_K_ht; // No solution in less than this space.
  #local max_tot_ht = N*(bot_K_ht/K); // No solution in more than this space.
  #if (tot_ht <= min_tot_ht)
    // Desperation answer -- {K} drawers {bot_K_ht/K}, rest {min_ht}, all shrunk by the same ratio:
    #debug concat("\n!! warning: {compute_drawer_heights} failed - total height too small\n")
    #local shrink = tot_ht/min_tot_ht;
    #for( kk, 0, K-1, 1 )
      #local ht[kk] = bot_K_ht/K * shrink;
    #end
    #for( kk, K, N-1, 1 )
      #local ht[kk] = min_ht * shrink;
    #end
  #elseif (tot_ht >= max_tot_ht)
    // Desperation answer -- equal heights
    #debug concat("\n!! warning: {compute_drawer_heights} failed - total height too big\n")
    #local htc = tot_ht/N;
    #for (kk, 0, N-1, 1)
      #local ht[kk] = htc;
    #end
  #else
    // We should be able to size at least {K+1} drawers geometrically:
    #local M = N;           // Number of heights that are in geometric progression.
    #local geo_ht = tot_ht; // Total height of those drawers.
    #local ok = false;
    #while ((M > K) & (! ok))
      #debug concat("\n!! trying with ", str(M,1,0), " drawers\n")
      #local htM = compute_drawer_heights_no_min(M, K, geo_ht, bot_K_ht);
      #debug concat("\n!! htM[", str(M-1,1,0), "] = ", str(htM[M-1],6,2))
      #if (htM[M-1] < min_ht)
        #debug concat(" too small, setting ht[", str(M-1,1,0), "] to ", str(min_ht,6,2), "\n")
        // Top drawer too shallow; set it to min height and retry with the others:
        #local ht[M-1] = min_ht;
        #local geo_ht = geo_ht - min_ht;
        #local M = M - 1;
      #else
        #debug concat("\n")
        #local ok = true;
      #end
    #end
    #if (M <= K)
      // Maybe can't happen, but just in case:
      #debug concat("\n!! warning: no drawers in geometric progression\n")
      #for (kk, 0, M-1, 1)
        #local ht[kk] = geo_ht/M;
      #end
    #else
      #if (M < N)
        #debug concat("\n!! warning: only ", str(M,1,0), " drawers in geometric progression\n")
      #end
      #for( kk, 0, M-1, 1 )
        #debug concat("\n!! setting ht[", str(kk,1,0), "] = ", str(htM[kk],6,2))
        #local ht[kk] = htM[kk];
      #end
    #end
  #end
  // Round heights to half-millimeters:
  #local sum_ht = 0;
  #for( kk, 0, N-1, 1)
    #local ht[kk] = int(2*ht[kk] + 0.5)/2.0;
    #local sum_ht = sum_ht + ht[kk];
  #end
  #debug concat("\n!! drawer height error = ", str(sum_ht - tot_ht, 0, 2), "\n")
  #debug concat("\n!! ------------------------------------------------------------\n")
  ht
#end

#macro compute_drawer_heights_no_min(N, K, tot_ht, bot_K_ht)
  // Like compute_drawer_heights}, but without a min height constraint.
  
  #local ht = array[N] // Heights from bottom to top.
  
  #if (N <= K)
    // Should not happen, but just in case:
    #for (kk, 0, N - 1, 1)
      #local ht[kk] = tot_ht/N;
    #end
  #else
    // Geometric progression:
    #local rtot = tot_ht / bot_K_ht;
    #local rat = compute_drawer_heights_find_ratio(N, K, rtot);
    #local bot_ht = bot_K_ht/compute_drawer_heights_geom_sum(rat, K);

    #local next_ht = bot_ht; 
    #for (kk, 0, N - 1, 1)
      #local ht[kk] = int(10*next_ht)/10.0;
      #local next_ht = next_ht * rat;
    #end
  #end
  ht
#end

#macro compute_drawer_heights_geom_sum(rr, N)
  // Computes {1 + r + r^2 + ... + r^{N-1} = (r^N - 1)/(r - 1)}.
  #if (rr = 1)
    #local Sr = N;
  #else
    #local Sr = (pow(rr, N) - 1)/(rr - 1);
  #end
  Sr
#end

#macro compute_drawer_heights_geom_sum_quot(rr, N, K)
  // Computes {(1 + r + r^2 + ... + r^{N-1})/(1 + r + r^2 + ... + r^{K-1})}.
  #if (rr = 1)
    #local res = N/K;
  #else
    #local res = (pow(rr, N) - 1)/(pow(rr, K) - 1);
  #end
  res
#end

#macro compute_drawer_heights_find_ratio(N, K, S)
  // Finds the hrat {r < 1} such that 
  // {(1 + r + r^2 + ... + r^{N-1}) = S(1 + r + r^2 + ... + r^{K-1})}.
  // Uses bissection. Requires {N > K} and {S < N/K};
  // otherwise returns 1.

  #debug concat("\n!! S = ", str(S,0,6), "\n")
  
  #if ((N < K) | (S > N/K))
    #local res = 1;
  #else
    // Bisection search:

    #local rlo = 0.01;
    #local Srlo = compute_drawer_heights_geom_sum_quot(rlo, N, K);
    #debug concat("\n!! rlo = ", str(rlo,0,6), " Srlo = ", str(Srlo,0,6), "\n")

    #local rhi = 1.0;
    #local Srhi = compute_drawer_heights_geom_sum_quot(rhi, N, K);
    #debug concat("\n!! rhi = ", str(rhi,0,6), " Srhi = ", str(Srhi,0,6), "\n")

    #while (rhi - rlo > 0.0001)
      #local rm = (rlo + rhi)/2;
      #local Sm = compute_drawer_heights_geom_sum_quot(rm, N, K);
      #if (Sm >= S)
        #local rhi = rm;
      #end
      #if (Sm <= S)
        #local rlo = rm;
      #end
    #end
    #local res = (rlo + rhi)/2;
  #end
  #local Sres = compute_drawer_heights_geom_sum_quot(res, N, K);
  #debug concat("\n!! res = ", str(res,0,6), " Sres = ", str(Sres,0,6), "\n")
  res
#end