common/trace/kracea,kracem,kraced,kracef kracea = 0 kracem = 0 kraced = 0 kracef = 0 call qnormal stop end subroutine qnormal double precision zd,ad integer vec(64),avec(64) vec(1) = 64 avec(1) = 64 1 continue write(6,*) write(6,*)' Enter 1 to calculate the area under normal curve' write(6,*)' Enter 2 to calculate the z belonging to given area.' write(6,*)' Enter 0 to exit program.' read(5,*,end=13)kch if(kch.eq.0) go to 13 if(kch.eq.1)go to 2 if(kch.eq.2)go to 2 write(6,*)' please enter either 1 or 2.' go to 1 2 continue write(6,*)' Enter 1 to calculate in single precision.' write(6,*)' Enter 2 to calculate in double precision.' write(6,*)' Enter 3 to calculate in 60 digit arithmetic.' write(6,*)' Enter 0 to exit program.' read(5,*,end=13)kp if(kp.eq.0)go to 13 if(kp.eq.1)go to 3 if(kp.eq.2)go to 3 if(kp.eq.3)go to 3 write(6,*)' please enter either 1,2 or 3.' go to 2 3 continue if(kch.eq.2) go to 6 write(6,*)' Enter the z value corresponding to the area desired.' read(5,*)z write(6,*) write(6,*)' The area under the normal curve from 0 to ',z,' is' if(kp.eq.2) go to 4 if(kp.eq.3) go to 5 call area1(z,a) write(6,*)a go to 1 4 continue zd = z call area2(zd,ad) write(6,*)ad go to 1 5 continue call rset(vec,z) call areav(vec,avec) call vprint(avec) go to 1 6 continue write(6,*)' Enter the area corresponding to the desired z.' read(5,*,end=13)a if(a.lt.0)go to 7 if(a.eq.0)go to 11 if(a.lt.0.5)go to 8 if(a.eq.0.5)go to 12 7 continue write(6,*)' Please enter an area between 0. and 0.5' go to 6 8 continue write(6,*) write(6,*)' The z value corresponding to the area ',a,' is ' if(kp.eq.2) go to 9 if(kp.eq.3) go to 10 call findz(a,z) write(6,*)z go to 1 9 ad = a call findzd(ad,zd) write(6,*)zd go to 1 10 call rset(avec,a) call findvz(avec,vec) call vprint(vec) go to 1 11 write(6,*)' Exactly zero.' go to 1 12 write(6,*)' infinity.' go to 1 13 stop end subroutine testvexp integer ex(64),x(64),ex2(64),x2(64) ex(1) = 64 x(1) = 64 ex2(1) = 64 x2(1) = 64 1 write(6,*)' enter value to test in vexp routine' read(5,*)y if(y.eq.-1.0)return expy = exp(y) call rset(x,y) write(6,*)' value of exponent is ',y write(6,*)' value of power is approximately ',expy write(6,*)' x = ' call vprint(x) call vexp(x,ex) write(6,*)'e^x = ' call vprint(ex) call add(x,x,x2) write(6,*)' 2 x = ' call vprint(x2) call mpy(ex,ex,ex2) write(6,*)' square of e^x = ' call vprint(ex2) call vexp(x2,ex2) write(6,*)'e^(2x) = ' call vprint(ex2) print * go to 1 end subroutine area1(z,a) real pi,sqrt2p integer k real z,a,sum,term,oldsum,z2,term1,term2 pi = atan(1.0) * 4 sqrt2pi = sqrt(2.0 * pi) z2 = z * z sum = 0 k = 0 term1 = z 1 continue term2 = -term1 * z2 * (2*k+1) / ( 2*(k+1) * (2*k+3) ) k = k + 1 term = term1 + term2 oldsum = sum sum = sum + term if(oldsum.eq.sum)go to 2 term1 = -term2 * z2 * (2*k+1) / ( 2*(k+1) * (2*k+3) ) k = k + 1 go to 1 2 continue a = sum/sqrt2pi return end subroutine findz(a,z) real z,a real pi,sqrt2p,zt0,zt1,at,ztsq,adiff pi = atan(1.0) * 4 sqrt2pi = sqrt(2.0 * pi) zt0 = 0.0 zt1 = a 1 continue ztm = zt0 zt0 = zt1 call area1(zt0,at) ztsq = zt0 * zt0 adiff =( at - a) aexp = exp(-.5 * ztsq) aprime = aexp/sqrt2pi zt1 = zt0 - adiff / aprime if(zt1.eq.zt0) go to 2 if(zt1.eq.ztm) go to 2 go to 1 2 continue z = zt1 return end subroutine findzd(ad,zd) double precision ad,zd double precision pi,sqrt2pi,zt0,zt1,at,ztsq double precision adiff,aexp,aprime double precision one double precision ztm one = 1.0 pi = datan(one) * 4.0 sqrt2pi = dsqrt(2.0 * pi) zt0 = 0 zt1 = ad 1 continue ztm = zt0 zt0 = zt1 call area2(zt0,at) ztsq = zt0 * zt0 adiff =( at - ad) aexp = dexp(-.5 * ztsq) aprime = aexp/sqrt2pi zt1 = zt0 - adiff / aprime if(zt1.eq.zt0) go to 2 if(zt1.eq.ztm) go to 2 go to 1 2 continue zd = zt1 return end subroutine findvz(avec,vec) common/trace/kracea,kracem,kraced,kracef integer avec(64),vec(64) integer zt1(64),zt0(64),ztsq(64),at(64),eztsq(64),ztdif(64) integer adiff(64) common/pix/pi(64),twopi(64),sqrt2pi(64),vsqrt2pi(64),two(64) integer pi integer vaprime(64) integer ztm(64) zt1(1) = 64 zt0(1) = 64 ztsq(1) = 64 at(1) = 64 eztsq(1) = 64 ztdif(1) = 64 adiff(1) = 64 pi(1) = 64 two(1) = 64 twopi(1) = 64 sqrt2pi(1) = 64 vsqrt2pi(1) = 64 vaprime(1) = 64 ztm(1) = 64 call setpi(pi) if(kracef.ne.0)then print *,' pi = ' call vprint(pi) pause endif call iset(two,2,0) if(kracef.ne.0)then print *,' two = ' call vprint(two) endif call mpy(two,pi,twopi) if(kracef.ne.0)then print *,' twopi = ' call vprint(twopi) endif call msqrt(twopi,sqrt2pi) if(kracef.ne.0)then print *,' sqrt2pi = ' call vprint(sqrt2pi) pause endif call recip(sqrt2pi,vsqrt2pi) if(kracef.ne.0)then print *,' vsqrt2pi = ' call vprint(vsqrt2pi) pause endif call iset(zt0,0,0) call copy(avec,zt1) 1 continue call copy(zt0,ztm) if(kracef.ne.0)then print *,' call copy(zt0,ztm) ' endif call copy(zt1,zt0) print *,' call copy(zt1,zt0) ' call areav(zt0,at) if(kracef.ne.0)then print *,' call areav(zt0,at) ' endif call mpy(zt0,zt0,ztsq) if(kracef.ne.0)then print *,' call mpy(zt0,zt0,ztsq) ' print *, ' zt0 = ' call vprint(zt0) call vprint(ztsq) pause endif call divbytwo(ztsq) if(kracef.ne.0)then print *,' call divbytwo(ztsq) ' endif call vexp(ztsq,eztsq) if(kracef.ne.0)then print *,' call vexp(ztsq,eztsq) ' endif call mpy(eztsq,vsqrt2pi,vaprime) if(kracef.ne.0)then print *,' call mpy(eztsq,vsqrt2pi,vaprime) ' endif call sub(at,avec,adiff) if(kracef.ne.0)then print *,' call sub(at,avec,adiff) ' endif call mpy(adiff,vaprime,ztdif) if(kracef.ne.0)then print *,' call mpy(adiff,vaprime,ztdif) ' print *,' adiff = ' call vprint(adiff) print *,' vaprime = ' call vprint(adiff) print *,' ztdif = ' call vprint(ztdif) pause endif call sub(zt0,ztdif,zt1) if(kracef.ne.0)then print *,' call sub(zt0,ztdif,zt1) ' print * print *,' zt0 = ' call vprint(zt0) print *,' ztdif = ' call vprint(ztdiff) print *,' zt1 = ' call vprint(zt1) pause endif call equal(zt1,zt0,kans) if(kans.eq.0)go to 2 call equal(zt1,ztm,kans) if(kans.eq.0)go to 2 go to 1 2 continue call copy(zt1,vec) return end subroutine vexp(ztsq,eztsq) integer ztsq(64),eztsq(64) integer vz(64),vsum(64),vterm1(64),vterm2(64),vx(64),vy(64) integer voldsum(64) vz(1) = 64 vsum(1) = 64 vterm1(1) = 64 vterm2(1) = 64 vx(1) = 64 vy(1) = 64 voldsum(1) = 64 call copy(ztsq,vz) call iset(vsum,1,0) k = 1 call iset(vterm1,1,0) 1 continue call iset(vy,k,0) call mpy(vterm1,vz,vx) call divide(vx,vy,vterm2) k = k + 1 call copy(vsum,voldsum) call add(vsum,vterm2,vx) call copy(vx,vsum) call copy(vterm2,vterm1) call equal(voldsum,vsum,kans) if(kans.eq.0)go to 2 go to 1 2 continue call copy(vsum,eztsq) return end subroutine area2(zd,ad) double precision zd,ad double precision pi,sqrt2pi integer k double precision sum,term,oldsum,z2,term1,term2 double precision one one = 1.0 pi = datan(one) * 4.0 sqrt2pi = dsqrt(2.0 * pi) z2 = zd * zd sum = 0 k = 0 term1 = zd 1 continue term2 = -term1 * z2 * (2*k+1) / ( 2*(k+1) * (2*k+3) ) k = k + 1 term = term1 + term2 oldsum = sum sum = sum + term if(oldsum.eq.sum)go to 2 term1 = -term2 * z2 * (2*k+1) / ( 2*(k+1) * (2*k+3) ) k = k + 1 go to 1 2 continue ad = sum/sqrt2pi return end subroutine areav(vecz,veca) integer vecz(64),veca(64) common/space/s(64),t(64),u(64),v(64),w(64),x(64),y(64),z(64) common/space2/ss(64),tt(64),uu(64),vv(64) integer s,t,u,v,w,x,y,z integer ss,tt,uu,vv common/pix/pi(64),twopi(64),sqrt2pi(64),vsqrt2pi(64),two(64) integer pi integer iz,k integer vz(64),vsum(64),vterm(64),vz2(64),vterm1(64),vterm2(64) integer voldsum(64),vsum2(64) s(1) = 64 t(1) = 64 u(1) = 64 v(1) = 64 w(1) = 64 x(1) = 64 y(1) = 64 z(1) = 64 ss(1) = 64 tt(1) = 64 uu(1) = 64 vv(1) = 64 vz(1) = 64 vsum(1) = 64 vterm(1) = 64 vz2(1) = 64 vterm1(1) = 64 vterm2(1) = 64 voldsum(1) = 64 vsum2(1) = 64 if(sqrt2pi(2).eq.59)go to 999 call setpi(pi) twopi(1) = 64 two(1) = 64 sqrt2pi(1) = 64 call iset(two,2,0) call mpy(two,pi,twopi) call msqrt(twopi,sqrt2pi) 999 continue call copy(vecz,vz) call mpy(vz,vz,vz2) call iset(vsum,0,0) k = 0 call copy(vz,vterm1) 1 continue call iset(y,2 * (k+1) * (2 * k + 3),0) call iset(x,2*k+1,0) call recip(y,z) call mpy(z,x,y) call mpy(vz2,y,x) call mpy(x,vterm1,vterm2) vterm2(3) = - vterm2(3) k = k + 1 call add(vterm1,vterm2,vterm) call copy(vsum,voldsum) call add(vsum,vterm,x) call copy(x,vsum) call equal(voldsum,vsum,kans) if(kans.eq.0)go to 2 call iset(y,2 * (k+1) * (2 * k + 3),0) call iset(x,2*k+1,0) call recip(y,z) call mpy(z,x,y) call mpy(vz2,y,x) call mpy(x,vterm2,vterm1) vterm1(3) = - vterm1(3) k = k + 1 go to 1 2 continue call divide(vsum,sqrt2pi,veca) return end subroutine vnorm character * 80, lineout common/space/s(64),t(64),u(64),v(64),w(64),x(64),y(64),z(64) common/space2/ss(64),tt(64),uu(64),vv(64) integer s,t,u,v,w,x,y,z integer ss,tt,uu,vv common/pix/pi(64),twopi(64),sqrt2pi(64),vsqrt2pi(64),two(64) integer pi * The standardized normal density function is * f(z) = exp( - .5 z^2 )/sqrt(2 pi) * * The normal distribution function is * z * integral ( f(z) dz) * 0 * * infinity * exp(y) = sum ( [ y^k ]/ k! ) * k=0 * * * * infinity * exp( - .5 z^2 ) = sum ( [-1]^k [z^2]^k /[ 2^k k! ] ) * k = 0 * * Integrate term by term to get * * infinity * area(z) = [1/sqrt(2 pi) ] sum ( [-1]^k z^[2 k + 1] / [ 2^k (2 k +1) k!] ) * k=0 * * real pi,sqrt2pi * integer iz,k integer iz,k * real z,sum,term,oldsum,z2,term1,term2 integer vz(64),vsum(64),vterm(64),vz2(64),vterm1(64),vterm2(64) integer voldsum(64),vsum2(64) s(1) = 64 t(1) = 64 u(1) = 64 v(1) = 64 w(1) = 64 x(1) = 64 y(1) = 64 z(1) = 64 ss(1) = 64 tt(1) = 64 uu(1) = 64 vv(1) = 64 vz(1) = 64 vsum(1) = 64 vterm(1) = 64 vz2(1) = 64 vterm1(1) = 64 vterm2(1) = 64 voldsum(1) = 64 vsum2(1) = 64 * open(1,file='normtable',status='unknown') * rewind 1 open(1,file='vnormtable',status='unknown') rewind 1 * pi = atan(1.0) * 4 call setpi(pi) call vencode(pi,lineout,leno) * sqrt2pi = sqrt(2.0 * pi) twopi(1) = 64 two(1) = 64 sqrt2pi(1) = 64 call iset(two,2,0) call mpy(two,pi,twopi) call vencode(twopi,lineout,leno) call msqrt(twopi,sqrt2pi) call vencode(sqrt2pi,lineout,leno) * do 4 iz = 1,300 do 4 iz = 1,650 write(6,*)' iz = ',iz * z = float(iz)/100. call iset(vz,iz,2) * z2 = z * z call mpy(vz,vz,vz2) * sum = 0 call iset(vsum,0,0) * k = 0 k = 0 * term1 = z call copy(vz,vterm1) * 1 continue 1 continue * term2 = -term1 * z2 * (2*k+1) / ( 2*(k+1) * (2*k+3) ) call iset(y,2 * (k+1) * (2 * k + 3),0) call iset(x,2*k+1,0) call recip(y,z) call mpy(z,x,y) call mpy(vz2,y,x) call mpy(x,vterm1,vterm2) vterm2(3) = - vterm2(3) * k = k + 1 k = k + 1 * term = term1 + term2 call add(vterm1,vterm2,vterm) * oldsum = sum call copy(vsum,voldsum) * sum = sum + term call add(vsum,vterm,x) call copy(x,vsum) * if(oldsum.eq.sum)go to 2 call equal(voldsum,vsum,kans) if(kans.eq.0)go to 2 * term1 = -term2 * z2 * (2*k+1) / ( 2*(k+1) * (2*k+3) ) call iset(y,2 * (k+1) * (2 * k + 3),0) call iset(x,2*k+1,0) call recip(y,z) call mpy(z,x,y) call mpy(vz2,y,x) call mpy(x,vterm2,vterm1) vterm1(3) = - vterm1(3) * k = k + 1 k = k + 1 * go to 1 go to 1 * 2 continue 2 continue * sum = sum/sqrt2pi call divide(vsum,sqrt2pi,x) call copy(x,vsum) * sum2 = 2.0 * sum call mpy(two,vsum,vsum2) * write(1,3)z,sum,sum2 call vwrite(vz,vsum,vsum2) * 3 format(1x,f5.3,2f12.8) * 4 continue 4 continue * rewind 1 rewind 1 stop end subroutine testadd common/space/s(64),t(64),u(64),v(64),w(64),x(64),y(64),z(64) common/space2/ss(64),tt(64),uu(64),vv(64) integer s,t,u,v,w,x,y,z integer ss,tt,uu,vv 1 continue m1 = 16777216.0 * rand() n1 = 90.0 * rand() - 45.0 m2 = 16777216.0 * rand() n2 = 90 * rand() - 45.0 m2 = -m2 call iset(s,m1,n1) call iset(t,m2,n2) call add(s,t,u) go to 1 return end subroutine vwrite(vz,vsum,vsum2) integer vz(64),vsum(64),vsum2(64) character * 80,linevz,linevs,linevs2,lino call vencode(vz,linevz,lenvz) call vencode(vsum,linevs,lenvs) call vencode(vsum2,linevs2,lenvs2) limit = 80 lino = linevz lino(10:limit)=linevs write(1,1)lino 1 format(a80) lino = ' * 2 = ' lino(10:limit)=linevs2 write(1,1)lino return end subroutine vencode(vz,linevz,lenvz) integer vz(*) integer lenvz character * 80,linevz character * 1,vsign izero = ichar('0') kzsig = vz(2) if(kzsig .eq. 0) go to 14 vsign = ' ' if(vz(3).lt.0)vsign = '-' kzdec = vz(4) kzint = kzsig - kzdec linevz = ' ' linevz(1:1) = vsign write(1,*) if(kzint)1,4,6 1 continue kzp = - kzint linevz(2:2) = '.' do 2 j = 1,kzp linevz(2+j:2+j)='0' 2 continue lim = min(78 - kzp,kzsig) do 3 jj = 1,lim j = kzsig + 1 - jj jpt = 2 + kzp + jj linevz(jpt:jpt) = char(izero +vz(4+j)) 3 continue lenvz = 2 + kzp + lim return 4 continue linevz(2:2) = '.' lim = min(78,kzsig) do 5 jj = 1,lim j = kzsig + 1 - jj jpt = 2 + jj linevz(jpt:jpt) = char(izero +vz(4+j)) 5 continue lenvz = lim + 2 return 6 continue if(kzint.gt.kzsig)go to 10 do 7 jj = 1,kzint j = kzsig + 1 - jj jpt = jj + 1 linevz(jpt:jpt) = char(izero +vz(4+j)) 7 continue jpt = kzint + 2 linevz(jpt:jpt) = '.' jjs = kzint + 1 do 8 jjj = 1,kzdec jj = jjj + kzint j = kzdec + 1 - jjj jpt = jj +2 if(jpt.gt.80)go to 9 linevz(jpt:jpt) = char(izero +vz(4+j)) 8 continue 9 continue lenvz = kzsig + 2 return 10 continue do 11 jj = 1,kzsig j = kzsig + 1 - jj jpt = jj + 1 linevz(jpt:jpt) = char(izero +vz(4+j)) 11 continue if(kzdec.ge.0)go to 13 ksp = - kzdec do 12 jj = 1,ksp jpt = kzsig + jj + 1 linevz(jpt:jpt) = '0' 12 continue 13 continue jpt = kzsig + ksp + 2 linevz(jpt:jpt) = '.' lenvz = kzsig - kzdec + 2 return 14 linevz(1:3) = ' 0.' lenvz = 3 return end subroutine vprint(vz) integer vz(64) character * 80,line integer len call vencode(vz,line,len) write(6,*)line(1:len) write(6,*)' size = ',vz(1),' sig = ',vz(2),' sign = ',vz(3) write(6,*)' decimal places = ',vz(4) write(6,*) return end subroutine dnorm * The standardized normal density function is * f(z) = exp( - .5 z^2 )/sqrt(2 pi) * * The normal distribution function is * z * integral ( f(z) dz) * 0 * * infinity * exp(y) = sum ( [ y^k ]/ k! ) * k=0 * * * * infinity * exp( - .5 z^2 ) = sum ( [-1]^k [z^2]^k /[ 2^k k! ] ) * k = 0 * * Integrate term by term to get * * infinity * area(z) = [1/sqrt(2 pi) ] sum ( [-1]^k z^[2 k + 1] / [ 2^k (2 k +1) k!] ) * k=0 * double precision pi,sqrt2p integer iz,k double precision z,sum,term,oldsum,z2,term1,term2 double precision one one = 1.0 open(1,file='dnormtable',status='unknown') rewind 1 pi = datan(one) * 4.0 sqrt2pi = sqrt(2.0 * pi) do 4 iz = 1,650 z = float(iz)/100. z2 = z * z sum = 0 k = 0 term1 = z 1 continue term2 = -term1 * z2 * (2*k+1) / ( 2*(k+1) * (2*k+3) ) k = k + 1 term = term1 + term2 oldsum = sum sum = sum + term if(oldsum.eq.sum)go to 2 term1 = -term2 * z2 * (2*k+1) / ( 2*(k+1) * (2*k+3) ) k = k + 1 go to 1 2 continue sum = sum/sqrt2pi sum2 = 2.0 * sum write(1,3)z,sum,sum2 3 format(1x,d12.3,2d15.8) 4 continue rewind 1 close(1) return end subroutine setpi common/pix/pi(64),twopi(64),sqrt2pi(64),vsqrt2pi(64),two(64) integer pi character * 64, piline data piline/ .'3.1415926535897932384626433832795028841971693993751058209749'/ if(pi(2).eq.59)return pi(1) = 64 pi(2) = 59 pi(3) = 1 pi(4) = 58 pi(63) = 3 izero = ichar('0') do 1 j = 1,58 k = 61-j kd = ichar(piline(k:k)) - izero pi(4+j) = kd 1 continue return end subroutine mpytest integer a(64),b(64),c(64) open(2,file='mpytest',status='unknown') rewind 2 a(1) = 64 b(1) = 64 c(1) = 64 a(2) = 1 a(3) = 1 a(4) = 0 a(5) = 2 call mpy(a,a,b) m = 1 1 continue m = m + 1 kb2 = b(2) write(2,2)m,(b(j),j=1,5),(b(5+kb2-j),j=1,kb2) 2 format(6i5,1x,60i1) if(m.gt.144)go to 3 call mpy(a,b,c) m = m + 1 kc2 = c(2) write(2,2)m,(c(j),j=1,5),(c(5+kc2-j),j=1,kc2) call mpy(a,c,b) go to 1 3 stop end subroutine anorm * The standardized normal density function is * f(z) = exp( - .5 z^2 )/sqrt(2 pi) * * The normal distribution function is * z * integral ( f(z) dz) * 0 * * infinity * exp(y) = sum ( [ y^k ]/ k! ) * k=0 * * * * infinity * exp( - .5 z^2 ) = sum ( [-1]^k [z^2]^k /[ 2^k k! ] ) * k = 0 * * Integrate term by term to get * * infinity * area(z) = [1/sqrt(2 pi) ] sum ( [-1]^k z^[2 k + 1] / [ 2^k (2 k +1) k!] ) * k=0 * real pi,sqrt2p integer iz,k real z,sum,term,oldsum,z2,term1,term2 open(1,file='normtable',status='unknown') rewind 1 pi = atan(1.0) * 4 sqrt2pi = sqrt(2.0 * pi) do 4 iz = 1,650 z = float(iz)/100. z2 = z * z sum = 0 k = 0 term1 = z 1 continue term2 = -term1 * z2 * (2*k+1) / ( 2*(k+1) * (2*k+3) ) k = k + 1 term = term1 + term2 oldsum = sum sum = sum + term if(oldsum.eq.sum)go to 2 term1 = -term2 * z2 * (2*k+1) / ( 2*(k+1) * (2*k+3) ) k = k + 1 go to 1 2 continue sum = sum/sqrt2pi sum2 = 2.0 * sum write(1,3)z,sum,sum2 3 format(1x,f5.3,2f12.8) 4 continue rewind 1 close(1) return end subroutine mpy(a,b,c) common/trace/kracea,kracem,kraced,kracef * * c = a * b * integer a(*),b(*),c(*) * * a(1) = dimension of a * b(1) = dimension of b * c(1) = dimension of c * a(2) = number of significant digits in a * b(2) = number of significant digits in b * c(2) = number of significant digits in c * a(3) = sign of a * b(3) = sign of b * c(3) = sign of c * a(4) = number of decimal places in a * b(4) = number of decimal places in b * c(4) = number of decimal places in c if(kracem.ne.0)then print *,' Entering subroutine mpy.' print *,' a = ' call vprint(a) print *,' b = ' call vprint(b) endif ka2 = a(2) kb2 = b(2) kc2 = c(1) - 4 kab2 = ka2 + kb2 if(kc2 .gt. kab2) kc2 = kab2 kaj2 = kab2 - kc2 if(ka2 .eq. 0) go to 6 if(kb2 .eq. 0) go to 6 nc = c(1) kab4 = a(4) + b(4) do 1 j = 1,kc2 c(4+j) = 0 1 continue jem = 0 do 3 ja = 1,ka2 jbstrt = max(1,2 + kaj2 - ja) do 2 jb = jbstrt,kb2 jem = jem + 1 jab = ja + jb - kaj2 c(3+jab) = c(3+jab) + a(4+ja) * b(4+jb) if(kracem.ne.0)then print *,' ja = ',ja,' jb = ',jb,' jab = ',jab print *,' c = ',c(3+jab),'a = ',a(4+ja),' b = ',b(4+jb) if(mod(jem,10).eq.0)pause endif 2 continue if(kracea.ne.0)then pause endif 3 continue kcm1 = kc2 - 1 do 4 jc = 1,kcm1 kq = c(4+jc)/10 kr = c(4+jc) - 10 * kq c(4+jc) = kr c(5+jc) = c(5 + jc) + kq if(kracem.ne.0)then print *,' jc = ',jc,' kq = ',kq,' kr = ',kr,' c(5+) = ',c(5+jc) if(mod(jc,10).eq.0)pause endif 4 continue if(kracem.ne.0)then pause end if if(c(4+kc2).ne.0)go to 5 kc2 = kcm1 5 continue 6 continue c(2) = kc2 c(3) = a(3) * b(3) c(4) = a(4) + b(4) - kaj2 if(kracem.ne.0)then print *,' c = ' call vprint(c) print *,' leaving subroutine mpy.' endif return end subroutine recip(a,b) integer a(*),b(*) common/space/s(64),t(64),u(64),v(64),w(64),x(64),y(64),z(64) integer s,t,u,v,w,x,y,z common/space2/ss(64),tt(64),uu(64),vv(64) integer ss,tt,uu,vv ss(1) = 64 tt(1) = 64 uu(1) = 64 vv(1) = 64 * a(1) = dimension of a * b(1) = dimension of b * a(2) = number of significant digits of a * b(2) = number of significant digits of b * a(3) = sign of a * b(3) = sign of b * a(4) = number of decimal places of a * b(4) = number of decimal places of b s(2) = a(2) s(3) = a(3) s(4) = s(2) kim = a(2) - a(4) ka2 = a(2) do 1 ja = 1,ka2 s(4+ja) = a(4+ja) 1 continue mad = s(4+ka2) if(mad.eq.0) then write(6,*) ' un-normalized real input to recip.' write(6,'(4i4,60i1/(16x,60i1))')s stop 1 endif mae = 9/mad 999 continue call iset(tt,mae,0) call mpy(s,tt,uu) if(uu(2).eq.uu(4))go to 998 mae = mae - 1 if(mae.eq.0)stop 2 go to 999 998 continue call iset(t,1,0) call sub(t,uu,w) call add(t,w,s) call mpy(w,w,u) 2 continue call add(s,u,v) call equal(s,v,ktest) if(ktest.eq.0)go to 3 call mpy(u,w,t) call add(v,t,s) call equal(v,s,ktest) if(ktest.eq.0) go to 3 call mpy(t,w,u) go to 2 3 continue call mpy(s,tt,vv) b(2) = vv(2) b(3) = vv(3) b(4) = vv(4) + kim kb2 = b(2) do 4 jb = 1,kb2 b(4 + jb) = vv(4 + jb) 4 continue return end subroutine msqrt(a,b) integer a(*),b(*) common/space/s(64),t(64),u(64),v(64),w(64),x(64),y(64),z(64) integer s,t,u,v,w,x,y,z s(1) = 64 t(1) = 64 u(1) = 64 v(1) = 64 w(1) = 64 x(1) = 64 y(1) = 64 z(1) = 64 call iset(x,1,0) iterlim = 0 1 continue iterlim = iterlim + 1 if(iterlim.gt.30) go to 3 call divide(a,x,z) call add(x,z,y) call divbytwo(y) call equal(x,y,kans) if(kans.eq.0)go to 2 iterlim = iterlim + 1 if(iterlim.gt.30) go to 3 call divide(a,y,z) call add(y,z,x) call divbytwo(x) if(kans.eq.0)go to 2 go to 1 2 call copy(x,b) return 3 continue write(6,*)' msqrt failed to converge after 30 iterations.' stop 2 end subroutine copy(a,b) integer a(*),b(*) nasig = a(2) nbsig = min(nasig,b(1) - 4) do 1 jj = 1,nbsig ja = nasig + 1 - jj jb = nbsig + 1 - jj b(4+jb) = a(4+ja) 1 continue ndifsig = nasig - nbsig b(2) = nbsig b(3) = a(3) b(4) = a(4) - ndifsig return end subroutine divbytwo(a) integer a(*) nsig = a(2) nsign = a(3) ndec = a(4) kr = 0 do 1 jj = 1,nsig j = nsig + 1 - jj kd = a(4+j) kd2 = 10 * kr + kd kq = kd2/2 kr = kd2 - 2 * kq a(4+j) = kq 1 continue return end subroutine unset(a,real) integer a(*) real real,term,preal1,preal2 nsig = a(2) ksign = a(3) ndec = a(4) nint = nsig - ndec term = 10.0 ** nint preal2 = 0 do 1 jj = 1,nsig j = nsig + 1 - jj preal1 = preal2 + a(j) * term if(preal1.eq.preal2) go to 2 preal2 = preal1 term = term / 10.0 1 continue 2 continue real = preal2 return end subroutine divide(a,b,c) integer a(*),b(*),c(*) common/space/s(64),t(64),u(64),v(64),w(64),x(64),y(64),z(64) integer s,t,u,v,w,x,y,z common/space2/ss(64),tt(64),uu(64),vv(64) integer ss,tt,uu,vv * a(1) = dimension of a * b(1) = dimension of b * c(1) = dimension of c * a(2) = number of significant digits of a * b(2) = number of significant digits of b * c(2) = number of significant digits of c * a(3) = sign of a * b(3) = sign of b * c(3) = sign of c * a(4) = number of decimal places of a * b(4) = number of decimal places of b * c(4) = number of decimal places of c s(1)=64 t(1)=64 u(1)=64 v(1)=64 w(1)=64 x(1)=64 y(1)=64 z(1)=64 ss(1)=64 tt(1)=64 uu(1)=64 vv(1)=64 call recip(b,ss) call mpy(a,ss,c) return end subroutine add(a,b,c) integer a(*),b(*),c(*) common/trace/kracea,kracem,kraced,kracef integer firstt firstt = 1 999 continue * a(1) is the dimension of a * b(1) is the dimension of b * c(1) is the dimension of c * a(2) is the number of significant digits of a * b(2) is the number of significant digits of b * c(2) is the number of significant digits of c * a(3) is the sign of a * b(3) is the sign of b * c(3) is the sign of c * a(4) is the number of decimal places of a * b(4) is the number of decimal places of b * c(4) is the number of decimal places of c if(kracea.ne.0)then write(6,*)' Entering subroutine add.' write(6,*) write(6,*)' a = ' call vprint(a) write(6,*)' b = ' call vprint(b) endif kasign = a(3) kbsign = b(3) kabsign = kasign * kbsign kcsign = kasign ka2 = a(2) if(kracea.ne.0)then write(6,*)' ka2 = ',ka2 endif if(ka2 .eq. 0) go to 15 ka4 = a(4) kaint = max(0, ka2 - ka4) kb2 = b(2) if(kracea.ne.0)then write(6,*)' ka4 = ',ka4,' kaint = ',kaint,' kb2 = ',kb2 endif if(kb2 .eq. 0) go to 17 kb4 = b(4) kbint = max(0,kb2 - kb4) if(kracea.ne.0)then write(6,*)' kb4 = ',kb4,' kbint = ',kbint,' kabsign = ',kabsign endif if(kabsign .eq. 1) go to 3 if( kaint .gt. kbint ) go to 3 if( kaint .lt. kbint ) go to 2 kint = max(kaint,kbint,ka2,kb2) if(kracea.ne.0)then write(6,*)' kint = ',kint endif do 1 j = 1,kint jl = 4 + kint - j ja = jl + ka4 jb = jl + kb4 if( (ja.le.0).and.(jb.le.0)) go to 1 if( (ja.le.0).and.(jb.gt.0)) go to 2 if( (ja.gt.0).and.(jb.le.0)) go to 3 if( (ja.gt.ka2).and.(jb.gt.kb2))go to 1 if( (ja.gt.ka2).and.(jb.le.kb2)) go to 2 if( (ja.le.ka2).and.(jb.gt.kb2)) go to 3 if( a(ja) .gt. b(jb) ) go to 3 if( a(ja) .lt. b(jb) ) go to 2 1 continue kcsign = 1 go to 3 2 kcsign = kbsign 3 continue if(kracea.ne.0)then write(6,*)' kasign = ',kasign,' kbsign = ',kbsign, . ' kcsign =',kcsign write(6,*) endif c(3) = kcsign ksign = (kasign - kbsign) * kcsign / 2 if(kracea.ne.0)then write(6,*)' ksign = ',ksign endif kcint = max(kaint,kbint) kc2x = c(1) - 4 kc4y = max(ka4,kb4) kc2y = kcint + kc4y kc2 = min(kc2x,kc2y) kc2diff = kc2y - kc2 if(kracea.ne.0)then write(6,*)' kcint = ',kcint,' kc2x = ',kc2x,' kc4y = ',kc4y write(6,*)' kc2y = ',kc2y,' kc2 = ',kc2,' kc2diff = ',kc2diff endif kc4 = kc2 - kcint ja4diff = ka4 - kc4y jb4diff = kb4 - kc4y if(kracea.ne.0)then write(6,*)' kc4 = ',kc4,' ja4diff = ',ja4diff,' jb4diff = ',jb4diff pause endif do 10 jc = 1,kc2y if(kracea.ne.0)then if(mod(jc,5).eq.1)pause endif kda = 0 ja = jc + ja4diff if(ja.lt.1)go to 4 if(ja.gt.ka2) go to 4 kda = a(4+ja) 4 continue if(kracea.ne.0)then write(6,*)' jc = ',jc,' ja = ',ja,' kda = ',kda endif kdb = 0 jb = jc + jb4diff if(jb.lt.1)go to 5 if(jb.gt.kb2)go to 5 kdb = b(4+jb) 5 continue if(kracea.ne.0)then write(6,*)' jb = ',jb,' kdb = ',kdb endif if(ksign)6,7,8 6 kdc = kdb - kda go to 9 7 kdc = kda + kdb go to 9 8 kdc = kda - kdb 9 continue if(kracea.ne.0)then write(6,*)' kda = ',kda,' kdb = ',kdb,' kdc = ',kdc endif jcc = jc - kc2diff if(jcc.le.0)go to 10 if(jcc.gt.kc2)go to 10 c(4 + jcc) = kdc 10 continue if(kracea.ne.0)then write(6,*)' jc = ',jc,' kc2diff = ',kc2diff,' jcc = ',jcc write(6,*)' kc2 = ',kc2,' kdc = ',kdc endif kcsiz = kc2 + 4 kc2m = kc2 - 1 if(kracea.ne.0)then write(6,*)' kc2 = ',kc2,' kcsiz = ',kcsiz,' kc2m = ',kc2m endif do 12 jc = 1,kc2m if(kracea.ne.0)then if(mod(jc,5).eq.1)pause write(6,*)' jc = ',jc endif kdc = c(4+jc) if(kracea.ne.0)then write(6,*)' kdc = ',kdc endif if(kdc.ge.0) go to 11 kq = 1 - (kdc+1)/10 kr = kdc + 10 * kq c(4+jc)= kr c(5+jc) = c(5+jc) - kq if(kracea.ne.0)then write(6,*)' kr = ',kr,' kq = ',kq,' c(5+jc) = ',c(5+jc) endif go to 12 11 continue if(kracea.ne.0)then write(6,*)' kdc = ',kdc endif if(kdc .lt. 10) go to 12 kq = kdc/10 kr = kdc - 10 * kq c(4 + jc) = kr c(5 + jc) = c(5 + jc) + kq if(kracea.ne.0)then write(6,*)' kr = ',kr,' kq = ',kq,' c(5+jc) = ',c(5+jc) endif 12 continue 13 continue if(kracea.ne.0)then write(6,*)' kc2 = ',kc2,' c(4+kc2) = ',c(4+kc2) endif if(c(4+kc2).ne.0)go to 14 kc2 = kc2 - 1 if(kc2.gt.0)go to 13 14 continue c(2) = kc2 c(4) = kc4 go to 19 15 continue c(2) = b(2) c(3) = b(3) c(4) = b(4) kb2 = b(2) do 16 jb = 1,kb2 c(4+jb) = b(4+jb) 16 continue go to 19 17 c(2) = a(2) c(3) = a(3) c(4) = a(4) ka2 = a(2) do 18 ja = 1,ka2 c(4+ja) = a(4+ja) 18 continue 19 continue kc2 = c(2) do 20 jc = 1,kc2 if(c(4+jc).lt.0)go to 21 if(c(4+jc).gt.9)go to 21 20 continue return 21 continue kracea = 1 go to 999 end subroutine sub(a,b,c) integer a(*),b(*),c(*) * a(1) is the dimension of a * b(1) is the dimension of b * c(1) is the dimension of c * a(2) is the number of significant digits of a * b(2) is the number of significant digits of b * c(2) is the number of significant digits of c * a(3) is the sign of a * b(3) is the sign of b * c(3) is the sign of c * a(4) is the number of decimal places of a * b(4) is the number of decimal places of b * c(4) is the number of decimal places of c b(3) = - b(3) call add(a,b,c) b(3) = - b(3) return end subroutine iset(a,m,kpow) integer a(*) k = m ja = 4 lim = a(1) a(3) = 1 if( k.ge.0) go to 1 k = - k a(3) = - a(3) 1 continue 2 continue if(k.eq.0)go to 3 ja = ja + 1 if(ja.gt.lim) go to 4 kq = k/10 kr = k - 10 * kq a(ja) = kr k = kq go to 2 3 a(2) = ja - 4 a(4) = kpow return 4 continue write(6,*)' error in iset program' write(6,*)' Overflow of array of dimension ',a(1) write(6,*)' in attempting to store ',m write(6,*) stop 1 end subroutine rset(a,z) integer a(*) kdec = 0 rk = z ja = 4 lim = a(1) a(3) = 1 if( rk.ge.0) go to 1 rk = - rk a(3) = - a(3) 1 continue 2 continue if(rk.eq.0)go to 3 ja = ja + 1 if(ja.gt.lim) go to 4 if(rk.lt.1)go to 25 k = rk kq = k/10 kr = k - 10 * kq a(ja) = kr rk = rk - k + kq go to 2 25 continue kdec = kdec + 1 do 26 jj = 6,ja jm = 6 + ja - jj a(jm) = a(jm - 1) 26 continue rk = 10 * rk kr = rk a(5) = kr rk = rk - kr go to 2 29 continue ja = a(1) 3 a(2) = ja - 4 a(4) = kdec return 4 continue if(rk.lt.1)go to 29 write(6,*)' error in rset program' write(6,*)' Overflow of array of dimension ',a(1) write(6,*)' in attempting to store integer part of ',z write(6,*) stop 1 end subroutine equal(a,b,ktest) integer a(*),b(*) ktest = 1 if(a(2).ne.b(2)) go to 2 if(a(3).ne.b(3)) go to 2 if(a(4).ne.b(4)) go to 2 ka2 = a(2) do 1 ja = 1,ka2 if( a(4+ja).ne.b(4+ja))go to 2 1 continue ktest = 0 2 continue return end