function value = multmod ( a, s, m )

%*****************************************************************************80
%
%% MULTMOD carries out modular multiplication.
%
%  Discussion:
%
%    This procedure returns
%
%      ( A * S ) mod M
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2013
%
%  Author:
%
%    Original Pascal version by Pierre L'Ecuyer, Serge Cote.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Pierre LEcuyer, Serge Cote,
%    Implementing a Random Number Package with Splitting Facilities,
%    ACM Transactions on Mathematical Software,
%    Volume 17, Number 1, March 1991, pages 98-111.
%
%  Parameters:
%
%    Input, integer A, S, M, the arguments.
%
%    Output, integer VALUE, the value of the product of A and S,
%    modulo M.
%
  h = 32768;

  if ( a <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MULTMOD - Fatal error!\n' );
    fprintf ( 1, '  A <= 0.\n' );
    error ( 'MULTMOD - Fatal error!' );
  end

  if ( m <= a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MULTMOD - Fatal error!\n' );
    fprintf ( 1, '  M <= A.\n' );
    error ( 'MULTMOD - Fatal error!' );
  end

  if ( s <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MULTMOD - Fatal error!\n' );
    fprintf ( 1, '  S <= 0.\n' );
    error ( 'MULTMOD - Fatal error!' );
  end

  if ( m <= s )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MULTMOD - Fatal error!\n' );
    fprintf ( 1, '  M <= S.\n' );
    error ( 'MULTMOD - Fatal error!' );
  end

  if ( a < h )

    a0 = a;
    p = 0;

  else

    a1 = i4_division ( a, h );
    a0 = a - h * a1;
    qh = i4_division ( m, h );
    rh = m - h * qh;

    if ( h <= a1 )

      a1 = a1 - h;
      k = i4_division ( s, qh );
      p = h * ( s - k * qh ) - k * rh;

      while ( p < 0 )
        p = p + m;
      end

    else

      p = 0;

    end

    if ( a1 ~= 0 )

      q = i4_division ( m, a1 );
      k = i4_division ( s, q );
      p = p - k * ( m - a1 * q );

      if ( 0 < p )
        p = p - m;
      end

      p = p + a1 * ( s - k * q );

      while ( p < 0 )
        p = p + m;
      end

    end

    k = i4_division ( p, qh );
    p = h * ( p - k * qh ) - k * rh;

    while ( p < 0 )
      p = p + m;
    end

  end

  if ( a0 ~= 0 )

    q = i4_division ( m, a0 );
    k = i4_division ( s, q );
    p = p - k * ( m - a0 * q );

    if ( 0 < p )
      p = p - m;
    end

    p = p + a0 * ( s - k * q );

    while ( p < 0 )
      p = p + m;
    end

  end

  value = p;

  return
end
