C+
C NAME:
C	CvR4
C PURPOSE:
C	Converts real*4 numbers between VMS, DOS and Unix
C CATEGORY:
C	gen/for/lib
C CALLING SEQUENCE:
	subroutine CvR4(cOp,N,XX)
C INPUTS:
C	cOp	character*3	operating system where the numbers originated
C	N	integer		# elements in XX
C	XX(N)	real*4		array
C OUTPUTS:
C	XX(N)	real*4		converted array
C CALLS:
C	CvSwap, MVBITS
C RESTRICTIONS:
C	cOp should correspond to an entry in the include file 'dirspec.h'.
C	Supported are 'DOS' (DOS/WIN/NT), 'UNX' (unix), 'VMS' and 'LNX' (Linux).
C INCLUDE:
	include		'dirspec.h'
	include		'math.h'
C PROCEDURE:
C	Usually the numbers are converted to the native OS (i.e. the OS specified
C	in 'dirspec.h'), and only the OS of origin needs to specified.
C	CvR4S and CvR8S allow specification of both origin and destination OS
C	for the conversion.
C	(This setup is not that great. Number representation is more a matter of
C	hardware architecture than OS. E.g. Linux and DOS/Windows running on Intel
C	or AMD use the same number representations. So we should be testing for
C	CPU architecture rather than OS. This is left as an excersize to the reader,
C	I suppose).
C
C 	INTEGER CONVERSIONS:
C
C	VMS and DOS use the same representation for integers.
C	Unix has the bytes in the opposite order relative to VMS and DOS (I think)
C	Linux is the same as DOS.
C
C 	REAL CONVERSIONS:

C	The only difference between DOS and Unix appears to be that the order of the
C	bytes is reversed (as for integers).
C	Linux is the same as DOS.
C
C	The real work is the conversion between DOS and VMS and v.v.
C
C	Real*8 structure (64 bits):
C	<------------------------ fraction -----------------> S<--- Exp --><-->	VMS
C	S<--- Exp --><-------------------- fraction -------------------------->	INTEL
C	32109876 54321098 76543210 98765432 10987654 32109876 54321098 76543210
C
C	Real*4 structure (32 bits):
C	<--- fraction --> S<- Exp -><----->	VMS
C	S<- Exp -><--- fraction ---------->	INTEL
C	10987654 32109876 54321098 76543210
C
C	For both real*4 and real*8 the position of exponent and fraction in VMS and DOS
C	are matched by reversing the order of the 2-byte words. 
C
C 	Real*8 (G-floating) numbers originating on the VAX have the following structure:
C	Bit 0..3,16..63 : normalized 52-bit fraction; the 53th bit (J-bit) is omitted, and is
C		always 1, i.e. the numbers are always normalized
C	Bit 4..14	: 11-bit biased exponent (biasing constant is 1025)
C	Bit 64		: sign bit
C
C	Real*8 for the INTEL FPU have the following structure
C	Bit 0..51  : normalized 52-bit fraction; the 53th bit (J-bit) is omitted, but contrary to
C		the VAX it can be zero (for very small, 'denormalized', numbers).
C	Bit 52..62 	: 11-bit biased exponent (biasing constant is 1023)
C	Bit 63     	: sign bit
C
C 	Real*4 (F-floating) numbers originating on the VAX have the following structure:
C	Bit 0..6,16..31 : normalized 24-bit fraction; the 24th bit (J-bit) is omitted, and is
C		always 1, i.e. the numbers are always normalized
C	Bit 7..14	: 8-bit biased exponent (biasing constant is 129)
C	Bit 15		: sign bit
C	As far as I can tell the exponent is always non-zero (1<=exp<=255), except for the
C	number zero (when all bits are zero). An exponent of zero combined with a non-zero fraction
C		are not valid real numbers.
C
C	Real*4 for the INTEL FPU have the following structure
C	Bit 0..22  : normalized 24-bit fraction; the 24th bit (J-bit) is omitted, but contrary to
C		the VAX it can be zero (for very small, 'denormalized', numbers).
C	Bit 23..30 	: 8-bit biased exponent (biasing constant is 127)
C	Bit 31     	: sign bit
C
C	VAX -> INTEL:
C	The VAX structure is mapped to the INTEL structure by reversing the order of the
C	words (2 for real*4, 4 for real*8.
C	The VAX exponent is usually two more than the INTEL exponent (the difference of the
C	biasing constants). If after subtraction of 2 the exponent is still positive, the
C	subtraction completes the conversion from a VAX to an INTEL real*4 or real*8. If the
C	corrected exponent is 0 or -1 than the INTEL real is denormalized and the fraction
C	must be updated also. This is done by shifting the bits in the fraction downward
C	(decreasing). Note that this means that the J-bit is shifted into the 'visible' 23-
C	or 52-bit part of the fraction.
C
C	INTEL -> VMS:
C	The Intel architecture uses exponent 255 (real*4) and 2047 (real*8) for special
C	purposes: quiet NaN, signal NaN, plus and minus infinity. These are not available
C	on VMS. These values will be set to the parameter values MATH__NARN or MATH__NARN8
C	and +/-MATH__PINF or +/-MATH__PINF8 defined in include file math.h.
C	Numbers with exponent 254 (real*4) and 2046 (real*8) are outside the range of
C	real values covered by VMS, and are set to +/-MATH__PINF or +/-MATH__PINF8.
C	Numbers with exponent in the range [1,253] (real*4) and [1,2045] (normalized numbers)
C	are converted to VMS by adding 2 to the exponent.
C	Numbers with exponent 0 (denormalized numbers) can be normalized on the VMS only
C	if the fraction can be made >=1 by shifting bits at most 2 positions. If it takes
C	more than 2 the number will be set to zero.
C MODIFICATION HISTORY:
C	OCT-1998, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C TEST PROGRAM
C	program test
C	include		'dirspec.h'
C	include		'openfile.h'
C	include		'math.h'
C
C	character	cOp*3
C	logical		bWrite
C
C	real*4		xx
C	real*8		dd
C	integer		i4(2)
C	equivalence	(dd,i4),(xx,i4)
C
C	real*4		xdos
C	real*8		ddos
C	integer		idos(2)
C	equivalence	(ddos,idos),(xdos,idos)
C
C	real*4		xvms
C	real*8		dvms
C	integer		ivms(2)
C	equivalence	(dvms,ivms),(xvms,ivms)
C
C	dd = MATH__NARN8
C
C	cOp = OS__VMS
C	if (cOpSys .eq. OS__VMS) cOp = OS__DOS
C
C	xx = .29E-36
C
C	I0 = I
C	do while (I-I0 .lt. 60)
C	    I = I+1
C
C	    ii = 0
C	    if (cOpSys .eq. OS__DOS) call MVBITS(i4(1),23,8,ii,0)
C	    if (cOpSys .eq. OS__VMS) call MVBITS(i4(1), 7,8,ii,0)
C
C	    xvms = xx
C	    call CvR4S(cOpSys,cOp,1,xvms)
C
C	    xdos = xvms
C	    call CvR4S(cOp,cOpSys,1,xdos)
C
C	    print *, ' exp ',ii, xx, xdos, xdos-xx, idos(1)-i4(1)
C
C	    xx = xx/1.1
C	end do
C
C	xx = MATH__PINF/100.
C	I0 = I
C
C	do while (I-I0 .lt. 10)
C	    I = I+1
C
C	    ii = 0
C	    if (cOpSys .eq. OS__DOS) call MVBITS(i4(1),23,8,ii,0)
C	    if (cOpSys .eq. OS__VMS) call MVBITS(i4(1), 7,8,ii,0)
C
C	    xvms = xx
C	    call CvR4S(cOpSys,cOp,1,xvms)
C
C	    xdos = xvms
C	    call CvR4S(cOp,cOpSys,1,xdos)
C
C	    print *, ' exp ',ii, xx, xdos, xdos-xx, idos(1)-i4(1)
C
C	    if (cOpSys .ne. OS__VMS .or. xx .lt. MATH__PINF/2.) xx = xx*2
C	end do
C
C	dd = 0.56D-306
C	I0 = I
C	do while (I-I0 .lt. 60)
C	    I = I+1
C
C	    ii = 0
C	    if (cOpSys .eq. OS__DOS) call MVBITS(i4(2),20,11,ii,0)
C	    if (cOpSys .eq. OS__VMS) call MVBITS(i4(1), 4,11,ii,0)
C
C	    dvms = dd
C	    call CvR8S(cOpSys,cOp,1,dvms)
C
C	    ddos = dvms
C	    call CvR8S(cOp,cOpSys,1,ddos)
C
C	    print *, ' exp ',ii, dd, ddos, ddos-dd, idos(1)-i4(1), idos(2)-i4(2)
C
C	    dd = dd/1.1d0
C	end do
C
C	dd = MATH__PINF8/100.
C	I0 = I
C
C	do while (I-I0 .lt. 10)
C	    I = I+1
C
C	    ii = 0
C	    if (cOpSys .eq. OS__DOS) call MVBITS(i4(2),20,11,ii,0)
C	    if (cOpSys .eq. OS__VMS) call MVBITS(i4(1), 4,11,ii,0)
C
C	    dvms = dd
C	    call CvR8S(cOpSys,cOp,1,dvms)
C
C	    ddos = dvms
C	    call CvR8S(cOp,cOpSys,1,ddos)
C
C	    print *, ' exp ',ii, dd, ddos, ddos-dd, idos(1)-i4(1), idos(2)-i4(2)
C
C	    if (cOpSys .ne. OS__VMS .or. dd .lt. MATH__PINF8/2) dd = dd*2d0
C	end do
C
C	end
C-
	    character	cOp*3
	    integer	N
	    real	XX(*)

	!~-------
	! For entry point CvR4S

	character	cOpDest*3

	!--------

	parameter	(nb = 4)

	character	cOrigin*3
	character	cDestin*3
	equivalence	(X,i4)

	cOrigin = cOp
	cDestin = cOpSys

	go to 1
C+
C NAME:
C	CvR4S
C PURPOSE:
C	Convert real*4 numbers between VMS, DOS & UNIX
C CALLING SEQUENCE:
	entry CvR4S(cOp,cOpDest,N,XX)
C INPUTS:
C	cOp	character*3	operating system where the numbers originated
C	cOpDest	character*3	operating system to be converted to
C	N	integer		# elements in XX
C	XX(N)	real*4		array
C OUTPUTS:
C	XX(N)	real*4		converted array
C RESTRICTIONS:
C	cOp and cOpDest should correspond to an entry in the include file 'dirspec.h'.
C CALLS:
C	CvSwap, MVBITS
C PROCEDURE:
C	See href=CvR4=
C-
	cOrigin = cOp
	cDestin = cOpDest

 1	continue

	!-------
	! Linux real*4 are the same as DOS real*4, so treat Linux as DOS

	if (cOrigin .eq. OS__LINUX) cOrigin = OS__DOS
	if (cDestin .eq. OS__LINUX) cDestin = OS__DOS

	!-------
	! Convert from Unix to DOS by reversing the byte order

	if (cOrigin .eq. OS__UNIX .and. cDestin .ne. OS__UNIX) then
	    call CvSwap(1,nb,N,XX)
	    cOrigin = OS__DOS				! XX is now in DOS form
	end if

	!-------
	! Unix has been changed to DOS. The origin is now DOS or VMS

	if (cOrigin .ne. OS__VMS .and. cDestin .eq. OS__VMS) then	! Change to VMS
	    do I=1,N
		X = XX(I)				! X equivalenced to i4

		ii = 0
		call MVBITS(i4,23,8,ii,0)		! Extract exponent
		if (ii .eq. 255) then
		    call MVBITS(i4,0,23,ii,0)		! Extract 23 bit fraction
		    if (ii .ne. 0) then			! Non-zero fraction
			i4 = MATH__NARN_VMS		! .. SNaN or QNaN
		    else				! Zero fraction
			ii = ISHFT(i4,-31)		! Sign bit: 0=+, 1=-
			i4 = MATH__PINF_VMS		! .. +/- Infinity
			if (ii .eq. 1) i4 = MATH__MINF_VMS
		    end if

		else if (ii .eq. 254) then		! Outside of VAX real*4 range
		    ii = ISHFT(i4,-31)			! Sign bit 0=+, 1=-
		    i4 = MATH__PINF_VMS			! +/- Infinity
		    if (ii .eq. 1) i4 = MATH__MINF_VMS

		else
		    if (ii. gt. 0) then			! Normalized number
			call MVBITS(ii+2,0,8,i4,23)	! Fix exponent (will be >= 3)
		    else				! Denormalized number
			jj = 0				! Redundant 24th J-bit is zero
			call MVBITS(i4,0,23,jj,0)	! Copy 23 bit fraction into jj
			if (jj .ge. Z'200000') then	! >=2^21: Big enough to normalize for VMS
			    ii = 1+ISHFT(jj,-22)	! New exponent = 2,1
			    call MVBITS(ii,0,8,i4,23)	! Fix exponent
			    jj = ISHFT(jj,3-ii)		! Shift by 1 or 2 positions
			    call MVBITS(jj,0,23,i4,0)	! Insert fixed 23 bit fraction into i4
			else				! Too small to normalize for VMS
			    X = 0.			! Set to zero
			end if
		    end if

		    call CvSwap(2,nb,1,i4)		! Reverse order of 16-bit words

		end if

		XX(I) = X

	    end do

	    cOrigin = OS__VMS				! XX now in VMS form

	else if (cOrigin .eq. OS__VMS .and. cDestin .ne. OS__VMS) then

	    do I=1,N
		X = XX(I)

		call CvSwap(2,nb,1,i4)			! Reverse order of 16-bit words

		ii = 0
		call MVBITS(i4,23,8,ii,0)		! Extract exponent
		if (ii .gt. 2) then
		    call MVBITS(ii-2,0,8,i4,23)		! Fix exponent (remains >= 1)
		else if (ii .gt. 0) then		! Exponent=2,1: denormalized numbers on Intel
		    call MVBITS(0,0,8,i4,23)		! Set exponent=0
		    jj = Z'800000'			! Set redundant 24th J-bit
		    call MVBITS(i4,0,23,jj,0)		! Copy 23 bit fraction into jj
		    jj = ISHFT(jj,ii-3)			! Shift by -1 or -2 positions
		    call MVBITS(jj,0,23,i4,0)		! Insert fixed 23 bit fraction into i4

		else					! Zero exp
		    ii = 0
		    call MVBITS(i4,0,23,ii,0)		! Copy 23 bit fraction into ii
		    i4 = 0
		    if (ii .ne. 0) i4 = MATH__NARN_DOS

		end if

		XX(I) = X

	    end do

	    cOrigin = OS__DOS

	end if

	if (cOrigin .eq. OS__DOS .and. cDestin .eq. OS__UNIX) call CvSwap(1,nb,N,XX)

	return
	end