# 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