>> Rodney Roberts IS & Education Professional Homepage   >> Programming Tutorials And Downloads

Science makes it known,
Engineering makes it work,
Art makes it beautiful.


Silverfrost Fortran 95 subroutine to compute Nth Root of X

NROOT is a REAL*8 function in sttstcs.dll (called by subroutine GEOMEAN in the same .dll).  It has been tested for X up to 9536, and is callable by D, Pascal, and FORTRAN.

Iteratively computes and returns   N  X   (alternatively, X1/N) using Newton's method, converging on root.  Most difficult part was programmatically making initial guess for a[1].  Solved by trial and error with the SQRT() function. Newton's method equation

Data Table
Algorithm term FORTRAN variable Description
A X (REAL*8); passed by reference parameter; if calling from a D procedure, define as a double in the calling D procedure Variable to take the root of
a[i] A1 (REAL*8) Previous guess of root
a[i+1] A2 (REAL*8) (returned to calling procedure through ROOTN) Current guess of root
n N (INTEGER*2); passed by reference parameter Root Index; root to be taken (inverse exponent)

Source Code1

C     X**(1/N)
      REAL*8    X, ROOTN, T2, A1, A2, D
      REAL      T1, EPS
      INTEGER*2 N, N2, I
      PARAMETER (EPS=0.000005)

      ROOTN = X
      IF ((N.EQ.1).OR.(N.LE.-1)) GO TO 100
      IF (N.EQ.0) THEN
        ROOTN = 1
        GO TO 100
      I = 0
      A1 = SQRT(SQRT(X))                !Begin Seed the nth root
      IF (X.GT.448000) THEN             !to prevent overflow for large X
        A1 = SQRT(A1)                   !and reduce number of iterations
      IF (X.GT.(1.0D70)) THEN
        A1 = SQRT(A1)
      IF (X.GT.(1.0D99)) THEN
        A1 = SQRT(A1)
      IF (X.GT.(1.0D199)) THEN
        A1 = SQRT(A1)
      IF (N.GT.25) THEN
        A1 = SQRT(A1)
      ENDIF                             !End Seed the nth root
      N2 = N - 1
      T1 = 1.0/FLOAT(N)
C      PRINT  causes Invalid Floating Point Operation when 
C                   called from D command line program
      I = I + 1
      IF (I.GT.500) GO TO 100     !Avoid Infinite Loop
      T2 = X/(A1**N2)
      A2 = T1 * (T2 + (FLOAT(N2)*A1))
      D = A2 - A1
        A1 = A2
        GO TO 10
      ROOTN = A2

FORTRAN Statements
FUNCTION Subprogram type (FUNCTION or SUBROUTINE) declaration; FUNCTION statement prefixed with data type (REAL, INTEGER, etc) of returned variable
NROOT Function name
RESULT specifies variable whose value is to be returned to calling procedure
REAL*8 8 byte floating point variable declaration - same as double in D
REAL 4 byte floating point variable - same as float in D
INTEGER*2 2 byte Integer, same as short in D
PARAMETER Attribute Declaration - constant part of expression
.EQ. logical = (or ==)
.LE. logical <= (or !>)
.OR. logical or (same as ||)
.GT. logical >
1.0D70 Double Precision (or REAL*8) constant,
1.0 * 1070 (1 followed by 70 zeroes)
SQRT() Square root library function - returns the square root of its passed argument
FLOAT() casts an Integer argument into a Real
GO TO transfer statement; jump to statement with label
statement label
CONTINUE “Dummy” statement
pass control to next statement in normal control flow
used in connection with DO loops as the final statement of loop range
RETURN return control to calling program
end of subroutine or function

1. Download includes fixes to take odd roots of negative numbers, exit with zero for even root of negative number.
Source file also includes REAL*8 FUNCTION RPOWER (X, P) RESULT (R).  Using natural logarithms, raises X to real power P.  Initial tests (called by Object Pascal) have been successful.

Thanks to Scientific Psychic for the square root HTML code,
gmplib.org for the algorithm

Any and all © copyrights, ® ™ trademarks, or other intellectual property (IP) mentioned here are the property of their respective owners.  No infringement is intended.

Feel free to use any of the above in your project (without violating any intellectual property rights); please give credit (same idea as Copyleft).
The source code files available for download on these pages are enhanced and tested periodically - check for revised versions.

Page best viewed with Mozilla FireFox 3.6.13 (or higher), Q4OS Konqueror R14.0.4.

Web hosting provided by

Award Space Web Hosting

,   Free Web Hosting Free Web Hosting. Hostinger Cheap Web Hosting. , &  GigaRocket Web Hosting

>> Rodney Roberts IS & Education Professional Homepage   >> Programming Tutorials And Downloads