{ 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 *************************************************}
times since October 20, 2004.