How to generate random positive integers that sum to N in Pascal

1 Answer

0 votes
program RandomSumProject;

{$mode objfpc}{$H+}{$J-}

uses
  SysUtils, Generics.Collections, Math;

type
  // Correct syntax to instantiate generic containers in Free Pascal (objfpc mode)
  TIntegerList = specialize TList<Integer>;
  TIntegerSet  = specialize THashSet<Integer>;

(*
    Generate k random positive integers that sum to n.

    Algorithm (stick-breaking / random composition):
    ------------------------------------------------
    1. We choose (k - 1) unique random "cut points" in the range [1, n - 1].
       Using unique cut points ensures that all resulting segments are 
       strictly greater than zero (positive integers).

    2. Sort the cut points in ascending order.

    3. Compute the differences between consecutive points
       (including 0 and n as boundaries). These differences
       are the k positive integers that sum to n.

    This produces a uniformly random composition of n.
*)

(*
    Paper Run (Dry Run) of the Random-Sum Program
    ---------------------------------------------

    Input:
        n = 50
        k = 7

    We need (k - 1) = 6 unique random cut points in the range [1, 49].

    Simulated Random outputs:
        12, 3, 40, 25, 7, 33

    After mapping to the [1, 49] range:
        cuts = {13, 4, 41, 26, 8, 34}

    Step 1: Sort the cut points
        cuts -> {4, 8, 13, 26, 34, 41}

    Step 2: Convert cut points into segment lengths
        prev = 0

        Result[0] = 4  - 0  = 4
        prev = 4

        Result[1] = 8  - 4  = 4
        prev = 8

        Result[2] = 13 - 8  = 5
        prev = 13

        Result[3] = 26 - 13 = 13
        prev = 26

        Result[4] = 34 - 26 = 8
        prev = 34

        Result[5] = 41 - 34 = 7
        prev = 41

        Result[6] = 50 - 41 = 9   (final segment)

    Final result:
        parts = {4, 4, 5, 13, 8, 7, 9}

    Verification:
        4 + 4 + 5 + 13 + 8 + 7 + 9 = 50

    Output:
        Random parts that sum to 50:
        4 4 5 13 8 7 9
        Sum = 50
*)

function GenerateRandomSum(N, K: Integer): TIntegerList;
var
  Cuts: TIntegerList;
  CutSet: TIntegerSet;
  Cut, Prev, NextCut: Integer;
begin
  if (K <= 0) or (N < K) then
    raise EInvalidArgument.Create('Invalid N or K: K must be positive and N >= K');

  Result := TIntegerList.Create;
  Cuts := TIntegerList.Create;
  CutSet := TIntegerSet.Create;
  try
    Cuts.Capacity := K - 1;
    Result.Capacity := K;

    // Generate (K - 1) unique random cut points
    while CutSet.Count < (K - 1) do
    begin
      // Random(Max) returns an integer in the range 0 <= X < Max.
      // We map this to 1 <= Cut <= N - 1.
      Cut := 1 + Random(N - 1);
      if CutSet.Add(Cut) then
        Cuts.Add(Cut);
    end;

    // Sort the cut points so we can compute segment lengths
    Cuts.Sort;

    Prev := 0;

    // Convert cut points into segment lengths
    for NextCut in Cuts do
    begin
      Result.Add(NextCut - Prev);
      Prev := NextCut;
    end;

    // Last segment: from last cut to N
    Result.Add(N - Prev);

  finally
    Cuts.Free;
    CutSet.Free;
  end;
end;

var
  N, K, X, TotalSum: Integer;
  Parts: TIntegerList;
begin
  // Initialize the built-in pseudo-random number generator with a random seed
  Randomize;

  N := 50;   // total sum
  K := 7;    // number of random parts

  try
    Parts := GenerateRandomSum(N, K);
    try
      WriteLn('Random parts that sum to ', N, ':');
      TotalSum := 0;
      
      for X in Parts do
      begin
        Write(X, ' ');
        TotalSum := TotalSum + X;
      end;
      WriteLn;

      // Verify sum
      WriteLn('Sum = ', TotalSum);
      
    finally
      Parts.Free;
    end;
  except
    on E: Exception do
      WriteLn('Error: ', E.Message);
  end;
end.



(*
run:

Random parts that sum to 50:
6 3 3 6 11 19 2 
Sum = 50

*)

 



answered 4 days ago by avibootz
...