Pythagorean Prime

{
    PYTHAGOREAN_PRIME.PAS computes pythagorean prime numbers and the square
    numbers that sum up to the respective pythagorean prime.
    To compile run "fpc pythagorean_prime.pas".
    Copyright (C) <February 01, 2016> Henning Polzer,
    send comments and error reports to: h underscore polzer at gmx dot de.

    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License
    as published by the Free Software Foundation; either version 2
    of the License, or (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
}

PROGRAM pythagoreische_primzahlen_summe_zweier_quadratzahlen;
TYPE Tgz = longint;                 { Standard: integer }

VAR i, og, ug: Tgz;

PROCEDURE ppz (zahl: Tgz);              { pythag. Primzahl finden und ausgeben }
VAR diff, pythag, wurzel: Tgz;

  FUNCTION prim (zahl: Tgz): boolean;       { Primzahl? }
  VAR k: Tgz;

  BEGIN { prim }
    prim := true;
    k := 3;

    IF zahl / 2 = trunc (zahl / 2) THEN prim := false;
    IF zahl / 2 <> trunc (zahl / 2) THEN
      WHILE ((zahl / k <> trunc (zahl / k)) AND (k * k < zahl)) DO
        k := k + 2;
    IF zahl / k = trunc (zahl / k) THEN prim := false;

    CASE zahl OF
      1      : prim := false;
      2, 3, 5: prim := true
    END { case }
  END; { prim }

  FUNCTION qzl (zahl: Tgz): boolean;        { Quadratzahl? }
  BEGIN { qzl }
    IF trunc (sqrt (zahl)) = sqrt (zahl) THEN qzl := true ELSE qzl := false
  END; { qzl }

BEGIN { ppz }
  pythag := 4*zahl+1;
  IF prim (pythag) THEN
  BEGIN
    wurzel := trunc (sqrt (pythag));
    REPEAT
      diff := pythag - wurzel*wurzel;
      IF qzl (diff) THEN write (diff, '+', wurzel*wurzel);
      wurzel := wurzel-1
    UNTIL (qzl (diff)) OR (wurzel < 2);     { sichere Abbruchbedingung }
    write ('=', pythag, ' ')
  END { if }
END; { ppz }

BEGIN { Hauptprogramm }
  write ('Untergrenze: '); readln (ug);
  write ('Obergrenze : '); readln (og);
  i := trunc (ug/4);

  REPEAT
    ppz (i);
    i := i + 1
  UNTIL i > trunc (og/4);
  writeln
END.