/*****************************************************************************
*                                                                            *
*   ISPRIME.SPL  Copyright (C) 2011 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Returns 1 if value is prime, else 0                         *
*                                                                            *
*   Revisions:   25 Feb 2011  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#if @HELP_ISPRIME

    ISPRIME

    Purpose: Returns 1 for each series value that is prime, else 0.

    Syntax:  ISPRIME(series)

                 series - A series or scalar, the input to test.

    Returns: A series or scalar if the input is a scalar.

    Example:
             W1: 0..8
             W2: isprime(w1)

             W2 contains the series {0, 0, 1, 1, 0, 1, 0, 1, 0}
             indicating the values {2, 3, 5, 7} are prime and
             {0, 1, 4, 6, 8} are not prime.

    Example:
             W1: 0..8
             W2: W1[find(isprime(w1))]

             W2 contains the series {2, 3, 5, 7}, the prime values
             of W1.

    Remarks:
             A prime number is a natural number (positive integer) with 
             exactly two divisors, 1 and itself. Thus, 0 and 1 are not
             prime.

    See Also:
             Primes
             Sermatch
#endif


/* find primes */
ITERATE isprime(n)
{
        local p, scalar;

        if ((scalar = isscalar(n))) n = {n};

        if ((max(n) > 1000000) || (max(n) / length(n) > 40))
        {
                /* remainder method */
                p = isprime_rem(n);
        }
        else
        {
                /* lookup method */
                p = isprime_match(n);
        }

        if (scalar) p = castlogical(p);

        return(p);
}


/* find primes by matching against a list */
isprime_match(n)
{
        local p, m;

        /* generate primes */
        p = primes(ceil(max(n)));

        /* indices of matches, p is sorted  */
        m = sermatch(n, p, 1, 1);
        p = logical(zeros(length(n), 1));

        /* these are prime */
        p[m] = 1;

        return(p);
}


/* find primes by division remainder */
isprime_rem(n)
{
        local p, q, r, k;

        q = logical(zeros(length(n), 1));
        p = 2..ceil(sqrt(max(n)));

        loop (k = 1..length(n))
        {
                r    = n[k];
                q[k] = all(r % p[find(p < r)]);
        }

        /* handle degenerate case */
        q[find(n < 2)] = 0;

        return(q);
}