# Tan lambert's formula

Hi all
I am going to implement this exercise in fortran:

to do so, i wrote code below:

``````module continued_fractions
implicit none
public

contains
recursive function tangente(number, max_iter) result(res)
real, intent(in) :: number
integer, intent(in) :: max_iter
integer :: n = 1
real :: res
real :: total = 1.0
res = total
if (max_iter == 1) then
return
else
res = res + (-1 * number ** 2) / ((2 * n + 1) + tangente(number, max_iter - 1))
endif
end function tangente
end module continued_fractions

program main
use continued_fractions, only: tangente
implicit none
real :: x = 0.6977, final_result
integer :: iter_num = 7
final_result = tangente(x, iter_num)
print *, "tan(40):", final_result
end program main
``````

But i got a wrong answer!

Thank you very much
Can i use you code in my repository?

Does things put on internet not automatically protected by copyright ? In the US: https://www.thebalancesmb.com/what-is-automatic-copyright-protection-3514945
In France, itâ€™s called â€śauthorâ€™s rightâ€ť and it is also automatic, as soon as the work is created.
To keep it in public domain, you should put for example a CC0 license I think.

1 Like

An URL to the thread seems to be a good practice!

2 Likes

By the way, this thread and your other ones appear quite academic in their nature, perhaps part of a university course on numerical analysis using programming languages of which Fortran might be a choice? If so, it will help readers if you share that as well.

As to your issue here with the calculation, you may have noted a loop construct is also a possibility in addition to a recursive function. The important aspect is structure and setup of the code that captures the essence of the algorithm accompanied by appropriate numerical analysis e.g., how and when to terminate a computation involving an infinite series. For the immediate case at hand, one might use the RECURSIVE attribute of procedures in Fortran to calculate the tangent using Lambertâ€™s formula like so:

``````! Program #  : Example07032020.1
! Author     : FortranFan
! Reference  : https://en.wikipedia.org/wiki/Gauss%27s_continued_fraction#cite_note-8
!
! Description:
! An example implementation that illustrates how to employ a
! RECURIVE function in Fortran to compute the tangent of x using
! Lambert's continued fraction dating back to 1768 which gives
! tan(x) = x/(1-x**2/(3-x**2/(5-x**2/(7-x**2/..))))
!

module kinds_m
integer, parameter :: WP = selected_real_kind( p=12 ) ! Select suitable precision
end module
module trig_m
use kinds_m, only : WP
! Named constants
real(WP), parameter :: ONE = 1.0_wp
real(WP), parameter :: TWO = 2.0_wp
real(WP), parameter :: PI = 3.14159265358979323846264338327950288_wp
real(WP), parameter :: DEG_TO_RAD = PI/180.0_wp
real(WP), parameter :: TOL = 1e-3_wp !<-- Suitable tolerance for continued fraction series
real(WP), parameter :: UPPER_LIMIT = ONE/TOL
contains
elemental function tand( degx ) result( tanx )
!  Calculate tangent of x in degrees using Lambert's formula
! Argument list
real(WP), intent(in) :: degx ! x in degrees
! Function result
real(WP) :: tanx
! Local variables
real(WP) :: x
tanx = x / ( ONE + CalcFracLambert(x, n=1))
return
end function
pure recursive function CalcFracLambert(x, n) result(Frac)
! Argument list
real(WP), intent(in) :: x ! x in radians
integer, intent(in)  :: n
! Function result
real(WP) :: Frac
! Local variables
real(WP) :: Term
Term = TWO*n + ONE
if ( Term > UPPER_LIMIT ) then
Frac = - x**2 / Term
else
Frac = - x**2 / ( Term + CalcFracLambert(x, n+1) )
end if
return
end function
end module
program CalcTanx
use kinds_m, only : WP
use trig_m, only : tand, DEG_TO_RAD
real(WP) :: degx, tanx
degx = 40.0_wp
tanx = tand( degx )
print *, "x(degrees): ", degx
print *, "tan(x) using Lambert's formula: ", tanx
print *, "% diff with intrinsic tan: ", (tanx/tan(degx*DEG_TO_RAD)-1.0_wp)*100.0_wp
stop
end program CalcTanx
``````

You can try the above with your compiler: with gfortran, the output is as expected:

C:\Temp>gfortran -Wall -std=f2018 p.f90 -o p.exe

C:\Temp>p.exe
x(degrees): 40.000000000000000
tan(x) using Lambertâ€™s formula: 0.83909963117727993
% diff with intrinsic tan: 0.0000000000000000

C:\Temp>

1 Like

Yes, These exercises belongs to this book:
Numerical methods of mathematics implemented in fortran - springer

Is it better to create a separate repository to share these exercises?

Nice. Iâ€™ve found that in practice I donâ€™t trust myself not to propagate copy/paste errors so I always use

``````real(wp), parameter :: PI = 4*atan(1.0_wp)
``````

But then again I always work in double precision so I donâ€™t need code for which precision needs to be variable or even a rarely used parameter.

1 Like

@ELNS @kargl @FortranFan @vmagnin This thread raises an important question do discuss. I opened a thread about it here:

The way it worked out in this thread is I think exactly how it should work: @ELNS asked for permission to reuse the code, @kargl responded with his terms. Later, @FortranFan posted his own code with his name in the header. I think these are good examples of how it can work.

2 Likes

Re: @simong comment on â€śin practice I donâ€™t trust myself not to propagate copy/paste errors so I always useâ€ť the ATAN intrinsic to setup PI as a named constant: you can find references that explain better the issues with this approach particularly in â€śproductionâ€ť code but quickly one can note

• there is potential for inconsistency and impact on portability depending on compilers and their implementation of the ATAN intrinsic,
• scientific and technical computing will often require a module of named constants anyway and including PI to sufficient digits in such a module is no more difficult than working with Planckâ€™s constant or Avogadro number, etc.
• one needs to develop good unit tests for all functionality anyway and that should reveal any error with â€ścopy-pasteâ€ť also.

Separately note what the Fortran standard views as â€śDOUBLE PRECISIONâ€ť may not be the same as what one might think of colloquially as double precision in floating-point representation.

Working with defined kinds is often the path of least problems with floating-point computations in Fortran and that is why I showed the KINDS module with the use of SELECTED_REAL_KIND intrinsic in the code above.

1 Like

What @FortranFan said is exactly how I do it in my codes, here is an example: https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/constants.f90#L12

You can notice there are more digits than double precision, so that the compiler rounds it correctly.

2 Likes

@kargl I just tested it with the code below and at least with GFortran 7.5.0 the `atan` vs direct decimal number is equivalent to all digits. But I know that such functions, for example in `openlibm`, are typically implemented to return the correctly rounded nearest floating point, but sometimes (such as in 25% of cases) they return the second nearest, based on my tests. In my applications I would not care about such a small difference, but if you can get the correctly rounded number in all cases, then it cannot hurt. And if the question is, which method is more reliable, I think no doubt the direct decimal input is more reliable.

``````program main
use iso_fortran_env, only: dp=>real64, qp=>real128
implicit none
real(dp), parameter :: pi    = 3.1415926535897932384626433832795_dp
real(dp), parameter :: pi2    = 4*atan(1._dp)

print "(es26.20)", 3.1415926535897932384626433832795_qp
print "(es23.17)", pi
print "(es23.17)", pi2
end program main
``````

Output with GFortran 7.5.0:

``````3.14159265358979323846E+00
3.14159265358979312E+00
3.14159265358979312E+00
``````
1 Like

Actually you are right, since this is done at compile time, atan will give you a correctly rounded answer. What I was talking about is runtime, which is not applicable here. Good point.

Indeed, I confirmed on my Ubuntu 18.04 with GFortran 7.5.0 that the runtime `4*atan(1._dp)` is also correctly rounded.

2 Likes