# SPSEARCH --- search for special primes (p = N^2 + 1)
  #  
  # spsearch(a,b) prints a list of the primes in the
  # interval [a,b] which are of the form "one plus a square,"
  # and it returns the number it finds.
  #
  # It is believed that there are infinitely many
  # such primes, but no one has ever proved this.

      spsearch := proc(a,b)
          local i, count;
          count := 0;
          for i from a to b do
            if isprime(i) then
              if i = (trunc(sqrt(i)))^2 + 1 then
                print(i);
                count := count + 1;
              fi;
            fi;
          od;
        RETURN (count );
      end:

  # jac, 6/10/2001