#include <perms.h>

//          Edinburgh IMP77 Compiler - Version 8.4
//
//     1  %SYSTEMREALFN SQRT(%REAL ARG)

extern float SQRT(float ARG) {
  //     2  %INTEGER I,J

  int I;
  int J;

  //     3  %REAL Y,OLD

  float Y;
  float OLD;

  //     4  %REALARRAY A(0:15)

  float A[16];
  struct _imp_dv A_dv = {A, sizeof(float) /*, dims, <bounds> */};

  //     5  !      %IF ARG<0.0 %THEN MLIBERR(50) %AND ARG=-ARG
  //     6        %IF ARG<0.0 %THEN PRINTSTRING("SQRT ARG NEG

  if (ARG >= 0.0) goto L_0001;
  _imp_PRINTSTRING(_imp_str_literal("SQRT ARG NEG\n"));
  return -1;
L_0001:;

//     7+ ") %AND %RESULT=-1
//     8        %IF ARG=0 %THEN %RESULT=0.0

  if (ARG != 0) goto L_0002;
  return 0.0;
L_0002:;

//     9        %CYCLE I=0,1,15

  // %for I = 0, 1, 15 %cycle
  {
    I = -1;
  L_0003:
    if (I == 15) goto L_0004;
    I += 1;

    //    10           Y=ARG/2^I

    Y = (double)(ARG / (_imp_REXP(2, I)));

    //    11           %IF Y<1 %THEN ->OUT

    if (Y >= 1) goto L_0006;
    goto U_1000;
  L_0006:;

  //    12        %REPEAT

    goto L_0003;
  }

L_0004:;

//    13  OUT:  %IF Y>1 %THEN I=17 %AND Y=Y/4 %ELSE I=I+1 %AND Y=Y/2

U_1000:;
  if (Y <= 1) goto L_0007;
  I = 17;
  Y = (double)(Y / 4);
  goto L_0008;
L_0007:;
  I = I + 1;
  Y = (double)(Y / 2);
L_0008:;

//    14        %IF Y>0.34375 %THEN A(0)=Y/2 +0.4368 %ELSE A(0)=Y/2+0.381

  if (Y <= 0.34375) goto L_0009;
  A[0] = (double)(Y / 2) + 0.4368;
  goto L_000A;
L_0009:;
  A[0] = (double)(Y / 2) + 0.381;
L_000A:;

//    15        %CYCLE J=0,1,14

  // %for J = 0, 1, 14 %cycle
  {
    J = -1;
  L_000B:
    if (J == 14) goto L_000C;
    J += 1;

    //    16           OLD=A(J)

    OLD = A[J];

    //    17           A(J+1)=(OLD+Y/OLD)/2

    A[J + 1] = (double)((OLD + (double)(Y / OLD)) / 2);

    //    18           %IF MOD(A(J+1)-OLD)<0.000000000000005 %THEN %EXIT

    if (_imp_MOD(A[J + 1] - OLD) >= 0.000000000000005) goto L_000D;
    goto L_000C;
  L_000D:;

  //    19        %REPEAT

    goto L_000B;
  }

L_000C:;

//    20        %RESULT=A(J+1)*2^(I//2)

  return A[J + 1] * (_imp_REXP(2, ((long)I / (long)2)));

  //    21  %END

} /* end of SQRT */

//    22  %ENDOFFILE