!
! IAN:MATHS
!
! Maths procedures for BASIC programs.
!

%const %integer excess      = 64
%const %integer exp mask    = 2_ 0111 1111
%const %integer sign mask   = 2_ 1000 0000
%const %integer b           = 16
%const %integer b shift     = 4
%const %integer k           = b shift
%const %integer k shift     = 2

%const %real half           = 2_0.1
%const %real inv ln 2       = 1.4426 95041
%const %real ln2 high       = 8_0.543
%const %real ln2 low        = -2.1219 444 @ -4
%const %real pi             = 16_ 0.324 3f7 @ 1
%const %real pi by 2        = 16_ 0.192 1fb @ 1
%const %real two by pi      = 16_ 0.A2F 983
%const %real inv pi         = 16_ 0.517 cc1

%const %real %array twopow ( -(k-1) : (k-1) ) =
      2_0.001, 2_0.01, 2_0.1, 1, 2, 4, 8

!
! >> TC <<
!
%real %function TC ( %integer flag, %real x )
  {%const %real p0         =  1.00000 0000 }
   %const %real p1         = -0.95801 7723 @ -1
  {%const %real q0         =  1.00000 0000 }
   %const %real q1         = -0.42913 5777
   %const %real q2         =  0.97168 5835 @ -2
   %const %real pi by 2 hi = 8_1.444
   %const %real pi by 2 lo = 4.8382 67948 97 @ -4
   %const %real eps        = 0.000244141
   %const %real ymax       = 6434
   %const %real eps1       = 1.383 @ -76
   %real y = |x|
   %if flag = 1 %and y < eps1 %then %signal 10, 1, 6
   %signal 10, 1, 5 %if y > ymax
   %integer n = int(x*two by pi)
   !
   ! NB The following line must be altered on a machine
   ! with no guard digits.
   !
   %real f = (x-n*pi by 2 hi)-n*pi by 2 lo
   %real den, num, g
   %if |f| < eps %start
      num = f
      den = 1
   %finish %else %start
      g = f*f
      num = p1*g*f+f
      den = (q2*g+q1)*g+half+half
   %finish
   %if n&1 # 0 %start
      num = -num
      flag = flag !! 1
   %finish
   %if flag = 0 %start
      %result = num/den
   %finish %else %start
      %result = den/num
   %finish
%end

!
! >> TAN <<
!
%real %function tan ( %real x )
   %result = TC ( 0, x )
%end

!
! >> COT <<
!
%real %function cot ( %real x )
   %result = TC ( 1, x )
%end

!
! >> SC <<
!
! Used only by the SIN and COS procedures.
!
%real %function SC ( %integer sgn, %real x, y )
   %const %real ymax  = 13000   {arg significance lost above this}
   %const %real pi hi = 8_3.11
   %const %real pi lo = 0.96765 35897 93 @ -3
   %const %real eps   = 16_1           @ -3
   %const %real r1    = -0.16666 65668
   %const %real r2    =  0.83330 25139 @ -2
   %const %real r3    = -0.19807 41872 @ -3
   %const %real r4    =  0.26019 03036 @ -5
   %integer n = int(y/pi)
   %real xn = n
   %real f, g, rg
   %signal 10,1,5 %if y >= ymax
   sgn = -sgn %if n&1 # 0
   xn = xn-half %if |x| # y
   ! Modify following line if no guard digit
   f = (|x|-xn*pi hi)-xn*pi lo
   %if |f| >= eps %start
      g = f*f
      rg = (((r4*g+r3)*g+r2)*g+r1)*g
      f = f+f*rg
   %finish
   f = -f %if sgn # 1
   %result = f
%end

!
! >> SIN <<
!
%real %function SIN ( %real x )
   %if x < 0 %start
      %result = sc(-1,x,-x)
   %else
      %result = sc(1,x,x)
   %finish
%end

!
! >> COS <<
!
%real %function COS ( %real x )
   %result = SC(1,x,|x|+piby2)
%end

!
! >> EXP <<
!
%real %function EXP ( %real x )
   %const %real bigx = 174.6730895
   %const %real smlx = -180.21827
   %const %real eps  = 3 @ -8
   %const %real p0   = 0.2499999995
   %const %real p1   = 0.41602 88626 @ -2
 { %const %real q0   = half }
   %const %real q1   = 0.49987 17887 @ -1
   %real g, z, gpz, qz, rg
   %integer n
   %signal 10,1,4 %if x > bigx
   %result = 0 %if x < smlx
   %result = 1 %if |x| < eps
   n = int(x*inv ln2)
   !
   ! NB the next line must be modified if the
   ! machine in use has no guard digits.
   !
   g = (x-n*ln2 high)-n*ln2 low
   z = g*g
   gpz = (p1*z+p0)*g
   qz  = (q1*z+{q0}half)
   rg  = gpz/(qz-gpz)+half
   n   = n+1
   !
   ! %result = rg * 2\n
   ! NB rg known to be positive.
   !
   rg  = rg * twopow(rem(n,k))
   byte integer(addr(rg)) = byte integer(addr(rg)) + n//k
   %result = rg
%end

!
! >> LOG <<
!
%real %function LOG ( %real f )
   %const %real a0 = -0.55270 74855
   %const %real b0 = -6.63271 8214
   %const %real c0 =  0.70710 67812
   %real znum, zden, w, z, r
   %signal 10, 1, 3 %if f <= 0
   %integer n = byteinteger(addr(f))-excess
   byteinteger(addr(f)) = excess
   n = n*k
   %while f < half %cycle
      f = f+f
      n = n-1
   %repeat
   %if f > c0 %start
      znum = f-half-half
      zden = f*half+half
   %finish %else %start
      n = n-1
      znum = f-half
      zden = znum*half+half
   %finish
   z = znum/zden
   w = z*z
   r = (w*a0/(w+b0))*z+z
   %result = (n*ln2 low+r)+n*ln2 high
%end

!
! >> LOG10 <<
!
%real %function LOG10 ( %real x )
   %const %real c3 = 8_0.33626 75425
   %result = c3 * LOG(x)
%end

!
! SQRT
!
%real %function sqrt ( %real x )
   %real y, f
   %integer n
   %result = 0 %if x = 0
   %signal 10,1,2 %if x < 0
   n = byteinteger(addr(x))-64
   byteinteger(addr(x)) = 64
   f = x
   !
   ! y <- initial estimate (linear poly)
   !
   y = 0.894427*f + 0.223607
   !
   ! 3 N-R iterations for 31 bits accuracy (24 needed)
   ! First two are Y = 0.5*(y+f/y)
   !
   y = (y + f/y)/2
   y = (y + f/y)/2
   !
   ! last is special
   !
   y = y+(f/y-y)/2
   !
   ! Now, for ODD N
   !
   %if n&1 # 0 %start
      y = y/4
      n = n+1
   %finish
   byteinteger(addr(y)) = byteinteger(addr(y)) + n//2
   %result = y
%end

!
! >> RAN <<
!
!
! Return a random real number in the range [0..1)
! Roughly 29 bits random part, so should NOT be
! used in tests of systems with > 29 bits of FP
! significand.
!
%real %function ran
   %own %integer y, called = 0
   %record(parm fm) p
   %const %integer m = 2796203
   %if called = 0 %start
      !
      ! Random initialisation from system time
      !
      called = 1
      svc(19,p)
      y = rem(p_p1 {system up time}, m)
   %finish
   y = rem(y*125,m)
   %result = y/m
%end

!
! >> RANDL <<
!
%real %function randl ( %real x )
   %external %real %function %spec exp ( %real x )
   %result = exp(x*ran)
%end

!
! >> POW <<
!
%real %function pow ( %real b, e )
   %result = exp(e*log(b))
%end

%end %of %file
