#Crib: http://www.math.uni-wuppertal.de/~xsc/literatur/PXSCENGL.pdf #

INT n = 4; # x, xd, y, yd #
[]STRING cmm = ("x ", "xd ", "y ", "yd ");
MODE R = LONG LONG REAL;
MODE RF = PROC(R)R;
RF rsqrt = long long sqrt;
RF rsin = long long sin, rcos = long long cos;
MODE SV = [n]R;
PROC hlf = (SV x)SV: (SV a; FOR j TO n DO a[j] := 0.5*x[j] OD; a);
PROC prs = (STRING m, SV x)VOID: (print((m, newline));
   FOR j TO n DO print((cmm[j], x[j], newline)) OD);
OP * = (R s, SV v)SV: (SV a; FOR j TO n DO a[j] := s*v[j] OD; a);
OP * = (REAL s, SV v)SV: (SV a; FOR j TO n DO a[j] := s*v[j] OD; a);
OP + = (SV x, y)SV: (SV a; FOR j TO n DO a[j] := x[j] + y[j] OD; a);
# This code is to solve the ordinary differential equation dy/dx = f(x, y) #
# in particular a Newtonian circular orbit #
PROC f = (R tau, SV y)SV: (R d2 = y[1]*y[1] + y[3]*y[3], d = rsqrt(d2), dm3 = 1.0/(d*d2);
   (SV a; FOR i BY 2 TO n DO a[i] := y[i+1]; a[i+1] := -y[i]*dm3 OD; a));
PROC step = (R dx, x, SV y)SV: (
  SV k1 = dx*f(x,y);
  SV k2 = dx*f(x+dx/2, y + hlf(k1));
  SV k3 = dx*f(x+dx/2, y + hlf(k2));
  SV k4 = dx*f(x+dx, y+k3);
  y + 1.0/R(6)*(k1+2.0*k2 + 2.0*k3 + k4));

(SV x := (1, 0, 0, 1);
PROC rep = (INT j)VOID: (print((j, newline));
    prs("dvs", x); prs("psh", f(0, x)));
INT c = 10000;
rep(0);
FOR k TO c DO x := step (1/c, 0, x);
  (k=1 OR k=c | rep(k)) OD;
print((x[1] - rcos(1), newline, x[3] - rsin(1))))

#
         +0
dvs
x +1.0000000000000000000000000000000000000000000000000000000e  +0
xd +0.0000000000000000000000000000000000000000000000000000000e  +0
y +0.0000000000000000000000000000000000000000000000000000000e  +0
yd +1.0000000000000000000000000000000000000000000000000000000e  +0
psh
x +0.0000000000000000000000000000000000000000000000000000000e  +0
xd -1.0000000000000000000000000000000000000000000000000000000e  +0
y +1.0000000000000000000000000000000000000000000000000000000e  +0
yd +0.0000000000000000000000000000000000000000000000000000000e  +0
         +1
dvs
x +9.9999999500000000416666665885416677408854154968261735916e  -1
xd -9.9999999833333333489583333138020825093587253392537595797e  -5
y +9.9999999833333333333333332161458338216145832570393906911e  -5
yd +9.9999999500000000416666665885416659830729221700032631048e  -1
psh
x -9.9999999833333333489583333138020825093587253392537595797e  -5
xd -9.9999999500000000416666670312499977213541731296115441865e  -1
y +9.9999999500000000416666665885416659830729221700032631048e  -1
yd -9.9999999833333333333333336588541652777777796919080979424e  -5
     +10000
dvs
x +5.4030230586813971729713847590873279222936617149931327057e  -1
xd -8.4147098480789650778474645471989997252818746298909865834e  -1
y +8.4147098480789650586543092183827677525190449567089409277e  -1
yd +5.4030230586813971696690300006023594663226673250267114350e  -1
psh
x -8.4147098480789650778474645471989997252818746298909865834e  -1
xd -5.4030230586813971846156577493233868054549194547538784052e  -1
y +5.4030230586813971696690300006023594663226673250267114350e  -1
yd -8.4147098480789650767891899213668427809323381818154931309e  -1
-1.0379813153424381150294424911860895709611314300000000000e -19
-7.8707139979202222437065856512747697290381702000000000000e -19
#