subroutine bis(f, a, b, eps, iter, rad, kod, ktip)
!
!	Subrutina calculeaza solutia ecuatiei F(x)=0, pe [a,b], prin bisectie.
!	'f': Functia f(x). Se declara 'external' in programul apelant.
!	'eps': Toleranta.
!	'lnit':Nr.maxim de iteratii/Nr.efectiv de iteratii (intrare/iesire)
!	'rad': Radacina calculata.
!	'kod': Cod de terminare a iteratiei:
!			0: s-a atins toleranta 'eps'
!			1: s-a depasit nr.max. de iteratii 'lnit'
!		   -1: pe intervalul specificat f(a)*f(b) > 0
!			2: radacina in 'a' sau 'b'
! 'ktip' =cod de tiparire iteratii: 0/#0

ULP(x) =1./2**(24 -exponent(x))

! Verifica a < b. Daca nu, interverteste:
if(.not. a <b) then
c = b; b =a; a =c
endif

! Initializeaza:
lnit =iter

! Initializeaza contor:
iter =0

if(ktip /=0) then
write(*,*); write(1,'(a)') '  Iter     c';
endif


! Initializeaza kod - pentru a, b = radacina:
kod =2
! Testeaza f(a), f(b):
fa =f(a); fb =f(b)
if(fa ==0.) then
rad =a
goto 500
elseif(fb ==0.) then
rad =b
goto 500
elseif(sign(1., fa)*sign(1., fb) >0) then
kod =-1
goto 500
endif

kod =0
! Iteratii:
1 iter =iter+1
c =(a +b)/2.

if(ktip /=0) then
print '(I6, 2x, 1pg15.7)', iter, c
write(1, '(I6, 2x, 1pg15.7)') iter, c
endif

if(iter > lnit) then
rad =c
kod =1
goto 500
endif

! Iteratie:
fb =f(b); fc =f(c)
! Ajusteaza eps
if(eps < ULP(c)) then
eps = ULP(c)

print *,' EPS ajustat la ULP:', eps
write(1,*) '* EPS ajustat la ULP:', eps
endif

IF (fc ==0 .or. b-c <= eps) THEN
rad =c
goto 500

ELSEIF(sign(1., fb)*sign(1., fc) <0) THEN	! Selecteaza intervalul:
a =c
ELSE
b =c
ENDIF
goto 1

500 lnit =iter
END

