GPrimesD.Pas Program Listing




{ Start of file GPrimesD.PAS ***********************************************}

{$A+,B-,D+,E+,F-,G-,I+,L+,N-,O-,R-,S+,V+,X-} { Pascal defaults }
{$M 16384,0,655360}

{$N+} { Uses numeric coprocessor }

program GPrimesD; { Plots Gaussian primes }

uses Crt, Graph; { Turbo Pascal interface }

const                 
  Name    = 'GPrimesD - Plots Gaussian primes';
  Version = 'Version 1.00, last revised: 97/01/17, 1500 hours';
  Author  = 'Copyright (c) 1997 by author: Harry J Smith,';
  Address = '19628 Via Monte Dr., Saratoga CA 95070.  All rights reserved.';

{***************************************************************************}

{ Developed in TURBO Pascal 6.0 }

{ This program was inspired by the article "Complex Gaussian Integers for
  'Gaussian Graphs'" by Henry G. Baker, ACM Sigplan Notices 28, 11 (Nov 1993),
  22-27. Found on the WWW at URL ftp://ftp.netcom.com/pub/hb/hbaker/
  Gaussian.html.

  It is similar to the program PTriplesP.Pas that plotted primitive
  Pythagorean triples.
}

type
  ColorType = 0 .. MaxColors; { Turbo Pascal type for colors }
  Reals     = Double;         { Change Double to Real if no coprocessor }

const
  ESC = Char(27);

var
  a    : Reals;   { a of (a, b, c) primitive Pythagorean triples }
  b    : Reals;   { b of (a, b, c) primitive Pythagorean triples }
  i    : Integer; { Loop index }
  Ch   : Char;    { Character read from keyboard }

{ Graphics variables }
  AspectR : Reals;   { Aspect ratio, EGA Hi = 10000 / 7750 = 1.2903225806 }
  XAsp    : Word;    { Horizontal aspect }
  YAsp    : Word;    { Vertical aspect }
  S       : Reals;   { Length of bottom side of rectangle }
  IS      : LongInt; { S as an integer }
  h       : Reals;   { Height of rectangle }
  SlopeX  : Reals;   { Slope of transform of X to pixels }
  SlopeY  : Reals;   { Slope of transform of Y to pixels }
  X0      : Reals;   { X  of point in transform of X to pixels }
  Xp0     : Reals;   { Xp of point in transform of X to pixels }
  Y0      : Reals;   { Y  of point in transform of Y to pixels }
  Yp0     : Reals;   { Yp of point in transform of Y to pixels }
  Xp      : Integer; { X in pixels }
  Yp      : Integer; { Y in pixels }
  Border  : Reals;   { Size of border on screen in X pixels }
  Mess    : String;  { Character string for message to screen }
  Device  : Integer; { Turbo Pascal graphics device code }
  Mode    : Integer; { Turbo Pascal graphics device mode code }
  MaxXPix : Word;    { Max pixel number in X direction, [0, MaxXPix] }
  MaxYPix : Word;    { Max pixel number in Y direction, [0, MaxYPix] }
  Color   : ColorType;   { Turbo Pascal color code }
  Palette : PaletteType; { Turbo Pascal color palette }
  PrevX   : Integer; { Previous X in pixels }
  PrevY   : Integer; { Previous Y in pixels }

{--------------------------------------}
procedure ExitProc; { Exit the program }
begin
  RestoreCrtMode;
  CloseGraph;
  Halt(0);
end; { ExitProc }

{--------------------------------------}
function ForceColor( Color : ColorType) : ColorType;
begin
  Color:= Color mod Palette.Size;
  IF Color = 0 then Color:= 1;
  SetColor( Color);
  ForceColor:= Color;
end; { ForceColor }

{--------------------------------------}
procedure ReadLInt( Mess : String; Min, Max, Nom : LongInt;
                  var II : LongInt); { Read in a long integer from keyboard }
var
  St   : String[255];
  Stat : Integer;
  LI   : LongInt;
  I    : Integer;
begin
  repeat
    repeat
      WriteLn( Mess);
      Write( '  [', Min, ', ', Max, '] (ENTER => ', Nom, '): ');
      ReadLn( St);
      WriteLn( St);
    until IOResult = 0;
    Val( St, LI, Stat);
  until ((Stat = 0) and (LI >= Min) and (LI <= Max)) or (Length(St) = 0);
  if Length( St) = 0 then  LI:= Nom;
  II:= LI;
  WriteLn( 'Input = ', II);
  WriteLn;
end; { ReadLInt }

{--------------------------------------}
procedure InitHead; { Initialize program }
begin
  TextBackground( Blue);
  TextColor( Yellow);
  ClrScr;
  WriteLn;
  WriteLn;
  WriteLn( Name);
  WriteLn( Version);
  WriteLn( Author);
  WriteLn( Address);
  WriteLn;
  WriteLn('While plotting, press any key to pause, or ESC to exit');
  WriteLn;
end; { InitHead }

{--------------------------------------}
procedure Init2; { Initialize program, 2nd part }
begin
  RestoreCrtMode;
  InitHead;
  WriteLn('Device, Mode, MaxX, MaxY = ',
           Device, ' ', Mode, ' ', MaxXPix, ' ', MaxYPix);
  Write('Palette.Size: Colors = ', Palette.Size, ':');
  FOR I:= 0 TO Palette.Size-1 DO
    Write(' ', Palette.Colors[I]);       { Show hardware information }
  WriteLn;
  WriteLn('Screen pixel aspect ratio = ', AspectR:0:9);
  WriteLn;
  WriteLn('Type any key to continue...  (or Esc to quit)');
  WriteLn;
  Ch:= ReadKey;
  if Ch = Chr(27) then  ExitProc;
  WriteLn;
  Mess:= 'Input largest x or y to be displayed, 0 => 24, S = ? ';
  ReadLInt( Mess, 0, 200, 24, IS);
  S:= IS;
  if (S = 0) then  S:= 24.0;
  h:= 2 * S;
  Border:= 0.0375 * MaxXPix;
  SetGraphMode( Mode);
  IF Palette.Size > 8 then  SetBkColor( Black);
  Color:= ForceColor( White);
end; { Init2 }

{--------------------------------------}
procedure Init; { Initialize program, can only do once }
begin
  InitHead; { Init header before and after InitGraph, InitGraph may hang }
  Device:= Detect;  Mode:= 0;
  InitGraph( Device, Mode, '');          { Determine type of graphics device }
  MaxXPix:= GetMaxX;  MaxYPix:= GetMaxY; { and initialize graphics }
  GetAspectRatio( XAsp, YAsp);
  AspectR:= YAsp / XAsp;
  GetPalette( Palette);
  IF Palette.Size < 2 then  Palette.Size:= 2;
  Init2;
end; { Init }

{--------------------------------------}
function TranX( X : Reals) : Integer; { Transform X to Pixels }
begin
  PrevX:= Round((X - X0) * SlopeX + Xp0);
  TranX:= PrevX;
end; { TranX }

{--------------------------------------}
function TranY( Y : Reals) : Integer; { Transform Y to Pixels }
begin
  PrevY:= Round((Y - Y0) * SlopeY + Yp0);
  TranY:= PrevY;
end; { TranY }

{--------------------------------------}
function Prime( P : LongInt) : Boolean; { Is P a Prime }
var
  I, R, Q : LongInt;
begin
  Prime:= TRUE;
  if P = 3 then  Exit;
  I:= 1;
  repeat
    I:= I + 2;
    Q:= P div I;
    R:= P - I * Q;
    if R = 0 then begin
      Prime:= False;  Exit;
    end;
  until I > Q;
end; { Prime }

{--------------------------------------}
procedure CompTran; { Compute constants for transforms to pixels }
                    { Using point-slope form of linear transform }
begin
  SlopeX:= (1.0 + MaxXPix - 2.0 * Border) / h;
  X0:= 0;  Xp0:= (1.0 + MaxXPix) / 2.0;
  SlopeY:= -SlopeX / AspectR;
  Y0:= 0;  Yp0:= (1.0 + MaxYPix) / 2.0;
  if (h > 10.0) or (TranY( h) * AspectR < Border) then begin
    SlopeY:= -(1 + MaxYPix - 2.0 * Border / AspectR) / h;
    Y0:= 0;  Yp0:= (1.0 + MaxYPix) / 2.0;
    SlopeX:= -SlopeY * AspectR;
    X0:= 0;  Xp0:= (1.0 + MaxXPix) / 2.0;
  end;
  if ( 0.99 < SlopeX) and (SlopeX <  1.01) then  SlopeX:=  1.0;
  if (-1.01 < SlopeY) and (SlopeY < -0.99) then  SlopeY:= -1.0;
  if (SlopeX >  3) then  SlopeX:=  2 * Int( (SlopeX - 1) / 2) + 1;
  if (SlopeY < -3) then  SlopeY:= -2 * Int( Abs((SlopeY + 1) / 2)) - 1;
end; { CompTran }

{--------------------------------------}
procedure PlotSq( a, b : Reals); { Plot 1 square from a, b }
var
  Sq : array[1..4] of PointType;
  d  : Reals;
begin
  d:= (S - 1) / (2 * S);
  Sq[1].x:= TranX(a - d);  Sq[1].y:= TranY(b - d);
  Sq[2].x:= TranX(a + d);  Sq[2].y:= TranY(b - d);
  Sq[3].x:= TranX(a + d);  Sq[3].y:= TranY(b + d);
  Sq[4].x:= TranX(a - d);  Sq[4].y:= TranY(b + d);
  FillPoly(4, Sq);
end; { PlotSq }

{--------------------------------------}
procedure Plot( a, b : Reals); { Plot 8 points from a, b }
begin
  PlotSq( a,  b); { Plot a Square }
  PlotSq( b,  a);
  PlotSq(-a,  b);
  PlotSq(-b,  a);
  PlotSq(-a, -b);
  PlotSq(-b, -a);
  PlotSq( a, -b);
  PlotSq( b, -a);
end; { Plot }

{--------------------------------------}
procedure SearchP; { Search for Gaussian primes }
var
  IS, P : LongInt;
begin
  MoveTo( 10, TextHeight('M') + 3);
  IS:= Trunc(S);
  Plot(1, 1);
  P:= 1;
  repeat
    repeat
      repeat          { Get next prime P}
        P:= P + 2;
      until Prime(P);
      if (P mod 4 = 3) then begin
        if (P < IS) then  Plot(P, 0);
      end;
    until (P mod 4 = 1);

    a:= 0;    { Reduce Prime integer P = 4*k + 1 to a^2 + b^2 }
    repeat
      a:= a + 1;
      b:= INT( SQRT( P - a * a));
    until (P = a * a + b * b) or (a>b);

    if (a < IS) and (b < IS) and (a ESC) then
          Ch:= ReadKey;             { Pause }
      end;
      if (Ch = ESC) then begin      { If ESC typed }
        a:= IS + 1;
        b:= IS + 1;                 { Break out of loops }
      end;
    end;
  until (a > IS) and (b > IS);
  while (Ch <> ESC) do  Ch:= ReadKey; { Pause }
end; { SearchP }

{--------------------------------------}
begin { Main program, PTripleP }
  Init; { Initialize program }
  repeat
    CompTran;                         { Compute constants for transforms }
    MoveTo( TranX(-S), TranY(-S));    { Draw an 2S by 2S square }
    LineTo( TranX(-S), TranY( S));
    LineTo( TranX( S), TranY( S));
    LineTo( TranX( S), TranY(-S));
    LineTo( TranX(-S), TranY(-S));
    a:= -S;
    repeat
      b:= -S;
      repeat
        MoveTo( TranX(a),  TranY(b)); { Plot a point }
        LineTo( PrevX,     PrevY);
        b:= b + 1.0;
      until (b >= S);
      a:= a + 1.0;

      while (Keypressed) do begin   { If key pressed }
        Ch:= ReadKey;
        if (Ch <> ESC) then
          Ch:= ReadKey;             { Pause }
      end;
      if (Ch = ESC) then begin      { If ESC typed }
        a:= S + 1;
        b:= S + 1;                  { Break out of loops }
      end;
    until (a >= S);
    MoveTo( 10, 1);
    Str(S:0:0, Mess);
    OutText('Plot of Gaussian primes, max x, y < ' + Mess + ':');
    SearchP;
    Init2; { Re-initialize program }
  until False;
end. { PTripleP }

{ End of file GPrimesD.PAS *************************************************}

Return to Gaussian Primes
Return to Harry's Home Page


This page accessed times since October 20, 2004.
Page created by: hjsmithh@sbcglobal.net
Changes last made on Saturday, 14-May-05 12:48:08 PDT

1