Tan lambert's formula

[ELNS]

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!

[ELNS]

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

[vmagnin]

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.

[vmagnin]

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

[FortranFan]

@ELNS,

Comments by @kargl re: author attribution (copyright, link to this thread, etc.) are rather pertinent and should apply in principle to other responses on this thread as well your various other inquiries thus far on this forum.

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
      x = degx * DEG_TO_RAD
      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>

[ELNS]

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?

[simong]

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.

[milancurcic]

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

License of code posted to the Fortran Discourse

A recent thread raised an important question: How to license the code posted to this forum and can it be reused? It’s important that we discuss this, agree on some solution, and post it as a guide in the sticky “Welcome / Read this First” post. If you post your code to the forum, you should have the right to claim copyright on it and license it how you want. However, most people just post code and don’t worry about how it’s used. We should have some common (default) baseline on what happens to …

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.

[FortranFan]

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.

[certik]

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.

[certik]

@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

[certik]

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.

[certik]

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