C C Программа INIT_CONST - получения констант C C Dushanov E. LCTA. Phone 6-48-45 C C Публикация: C Э.Б.Душанов, М.Г.Емельяненко, Г.Ю.Коновалова C О форматах представления вещественных чисел и алгоритме C автоматического определения констант вещественной C арифметики ЭВМ. ОИЯИ,Р11-2000-163,Дубна,2000. C C beta - основание системы счисления, C (l<0,u>0) - нижняя и верхняя границы порядков представления C нормализованного числа в ЭВМ, C t - число beta-ичных разрядов, отведенных для хранения C мантиссы числа в ЭВМ, C e_8 - максимальное число, представимое в данной ЭВМ, C e_0 - минимальное число, представимое в данной ЭВМ, C e_1 - относительная погрешность вычислений с нормализованными C числами данной ЭВМ, C C а также получения параметров настройки функций GTEXP,CHEXP [89] C C Общие блоки: COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN C COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ C COMMON /F499_ICHGT/ BINV,INO,ILT,MNE1,MNE2,INB,LWT,KB C Обpащение: CALL INIT_CONST, если печать вычисленных костант не требуется; C CALL INITP, если константы печатаются. C Выходные данные: вещественные константы e_8,e_0,e_1,e_2=e_1/beta C и целые константы beta,u,l,t. Выдаются в блоках F499_RCONST и C F499_ICONST в указанных последовательностях соответственно. C В блоке F499_ICHGT содержатся следующие выделенные параметры C настройки функции GTEXP, CHEXP: C BINV - (real*8) вещественное число равное 1/beta; C KB - (integer) номер половины формата, в котором хранится C младший разряд порядка нормализованного числа; C LWT - (integer) представление целого порядка e=0 в разрядах C формата, отведенных для порядка; C INO - (integer) количество beta-ичных разрядов, отведенных для C порядка с учетом его знака (в формате для хранения C вещественного числа); C INB - (integer) количество beta-ичных разрядов, отведенных для C мантиссы с учетом ее знака (в формате для хранения C вещественного числа); C MNE1 - (integer) признак кода (e>0)-положительного порядка. C При этом MNE1=[e]^k_(код)-e; C MNE2 - (integer) признак кода (e<0)-отрицательного порядка. C При этом MNE2=[e]^k_(код)-e. C SUBROUTINE INIT_CONST IMPLICIT REAL*8 (A-H, O-Z) INTEGER RADIX,J1(2),J2(2),J3(2),N1(2),N2(2),N3(2) COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ COMMON /F499_ICHGT/ BINV,INO,ILT,MNE1,MNE2,INB,LWT,KB EQUIVALENCE (J1,R1),(J2,R2),(J3,R3) EQUIVALENCE (OVFLO,N1),(UNDFLO,N2),(EPS,N3) data radix/0/ C lpconst=0 ! TSap: for GNU Fortran G77 if(radix.ne.0) return ! Sap: already calculated GOTO 1 ENTRY INITP LPCONST=1 1 R1=1.D0 RADIX=1 KB=1 C ЦИКЛ ПОИСКА (BETA) - ОСНОВАНИЯ СИСТЕМЫ СЧИСЛЕНИЯ 2 RADIX=RADIX+1 R2=RADIX R3=1D0/R2 ILT=J2(1)-J1(1) INO=J2(2)-J1(2) IF(ILT.NE.(J1(1)-J3(1)).OR.INO.NE.(J1(2)-J3(2))) GOTO 2 C ЦИКЛ ПОИСКА t - РАЗРЯДОВ ОТВЕДЕННЫХ, ДЛЯ ХРАНЕНИИ МАНТИССЫ m D=R3*(R2-R1) UNDFLO=D MANTSZ=0 BINV=R3 IF(INO.EQ.0) GOTO 3 ILT=INO KB=2 3 OVFLO=UNDFLO MANTSZ=MANTSZ+1 D=D*R3 UNDFLO=UNDFLO+D EPS=UNDFLO IF(EPS.LT.R1) GOTO 3 EPSMIN=R1-OVFLO MINEXP=2-MANTSZ EPS=EPSMIN UNDFLO=EPS*R2 MNE2=N2(KB) IF(MNE2.LT.N3(KB)) MNE2=-MNE2 MNE2=MNE2/ILT-MINEXP I=IABS(N3(KB)-N2(KB)) C ЦИКЛ ВЫЧИСЛЕНИЯ e_0 - МИНИМАЛЬНО ПРЕДСТАВИМОГО C В ЭВМ ЧИСЛА И ЕГО ПОРЯДКА 4 UNDFLO=EPS MINEXP=MINEXP-1 D=UNDFLO*R3 EPS=D IF(EPS.GT.0D0.AND.I.EQ.IABS(N3(KB)-N2(KB))) GOTO 4 C ЦИКЛ ВЫЧИСЛЕНИЯ u - ПОРЯДКА МАКСИМАЛЬНО ПРЕДСТАВИМОГО C В ЭВМ ЧИСЛА e_8 EPS=R1 D=R2 I=1 5 MAXEXP=2*I INO=IOR(N3(KB),INO) ! for G77: INO=N3(KB).OR.INO IF(MAXEXP.GE.-MINEXP) GOTO 6 EPS=EPS*D D=D*D I=MAXEXP GOTO 5 6 IF((MAXEXP+MINEXP).LE.1) MAXEXP=MAXEXP-1 C ЦИКЛ ВЫЧИСЛЕНИЯ e_8 - МАКСИМАЛЬНО ПРЕДСТАВИМОГО В ЭВМ ЧИСЛА OVFLO=OVFLO*EPS DO 7 K=I,MAXEXP OVFLO=OVFLO*R2 7 CONTINUE C БЛОК ВЫЧИСЛЕНИЯ ПАРАМЕТРОВ ВЫДЕЛЕНИЯ ПОРЯДКА И МАНТИССЫ ЧИСЛА EPS=R3*R3 I=MIN0(N3(KB),N2(KB),J3(KB)) INO=IOR(IOR((INO-I),(N2(KB)-I)),(J3(KB)-I)) ! for G77: INO=(INO-I).OR.(N2(KB)-I).OR.(J3(KB)-I) INO=IOR(IOR(INO,(N3(KB)-I)),(J1(KB)-I)) ! INO=INO.OR.(N3(KB)-I).OR.(J1(KB)-I) INB=NOT(INO) ! INB=.NOT.INO LWT=IAND(J3(KB),INO) ! LWT=J3(KB).AND.INO MNE1=LWT/ILT EPS=EPSMIN*R2 IF(LPCONST.NE.1) GOTO 8 write(*,*)'The constants of mashin`s arithmetic ...' write(*,9) MAXEXP,OVFLO,MINEXP,UNDFLO,RADIX,EPS,MANTSZ,EPSMIN 9 FORMAT(1X,'MAXEXP=',I8,' OVFLO=',1PE24.15E3,/, * 1X,'MINEXP=',I8,' UNDFLO=',1PE24.15E3,/, * 1X,' RADIX=',I8,' EPS=',1PE24.15E3,/, * 1X,'MANTSZ=',I8,' EPSMIN=',1PE24.15E3) 8 RETURN END C C Подпрограмма-функция GTEXP - выделения мантиссы и порядка вещественного C числа C Общие блоки: COMMON /F499_ICHGT/ BINV,INO,ILT,MNE1,MNE2,INB,LWT,KB C Обpащение: RM=GTEXP(R,E) C Входные данные: R - (real*8) вещественное число, порядок и мантисса C которого выделяется. Блок F499_ICHGT содержит настроечные C параметры, определенные в INIT_CONST (или заранее известные); C Выходные данные: RM - (real*8) содержит выделенную мантиссу числа R; C E - (integer) содержит выделенный порядок числа R C REAL*8 FUNCTION GTEXP(R,E) INTEGER I(2),E REAL*8 R,Z,BINV COMMON /F499_ICHGT/ BINV,INO,ILT,MNE1,MNE2,INB,LWT,KB EQUIVALENCE (Z,I) IF(R.EQ.0D0) GOTO 2 Z=R J=MNE1 E=IAND(I(KB),INO)/ILT ! for G77: E=(I(KB).AND.INO)/ILT IF(DABS(Z).GE.BINV) GOTO 3 IF(MNE2.LT.0) E=-E J=MNE2 3 E=E-J I(KB)=IOR(IAND(I(KB),INB),LWT) ! for G77: I(KB)=(I(KB).AND.INB).OR.LWT GTEXP=Z 1 RETURN 2 GTEXP=0 E=0 GOTO 1 END C C Подпрограмма-функция CHEXP - восстанавления вещественного числа по его C мантиссе и порядку C C Общие блоки: COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ C COMMON /F499_ICHGT/ BINV,INO,ILT,MNE1,MNE2,INB,LWT,KB C Обpащение: R=CHEXP(RM,E) C Входные данные: RM - (real*8) мантисса числа R; E - (integer) порядок C числа R. Блоки F499_ICONST и F499_ICHGT содержат элементы, C определенные в INIT_CONST (или заранее известные); C Выходные данные: R - (real*8) содержит вещественное число beta^E RM=R C REAL*8 FUNCTION CHEXP(R,E) INTEGER I(2),L,J,E,RADIX REAL*8 R,Z,BINV COMMON /F499_ICHGT/ BINV,INO,ILT,MNE1,MNE2,INB,LWT,KB COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ EQUIVALENCE (Z,I) IF(E.EQ.0.OR.R.EQ.0D0) GOTO 2 Z=R K=MNE1 J=IAND(I(KB),INO)/ILT ! for G77: J=(I(KB).AND.INO)/ILT IF(DABS(Z).GE.BINV) GOTO 4 IF(MNE2.LT.0) J=-J K=MNE2 4 L=J-K+E IF(L.GT.MAXEXP) L=MAXEXP K=MNE1 IF(L.GT.0) GOTO 5 IF(L.LT.MINEXP) GOTO 3 K=IABS(MNE2) IF(MNE2.LT.0) L=-L 5 I(KB)=IOR(IAND(I(KB),INB),((L+K)*ILT)) ! for G77: I(KB)=(I(KB).AND.INB).OR.((L+K)*ILT) CHEXP=Z 1 RETURN 2 CHEXP=R GOTO 1 3 CHEXP=0 GOTO 1 END