hp41programs

Curvature &Torsion

Curvature & Torsion for the HP-41


Overview
 

 1°)  Derivatives

  a)  First, Second & Third Derivatives
  b)  Functions of N variables
  c)  Hessian Matrix+Gradient

 2°)  Gaussian & Mean Curvatures of a Surface

  a)  Surface
  b)  Hypersurface

 3°)  Curvature & Torsion of a Curve

  a)  3-Dimensional Problem
  b)  N-Dimensional Problem
 

-In paragraph 1 are listed several programs to calculate various derivatives
-Some of them are already listed elsewhere
 

1°)  Derivatives
 

      a) First, Second & Third Derivatives
 

-This routine uses formulae that give exact results for polynomials of degree < 11
-The formulae for  f '(x) & f ''(x)  are given in "Numerical Differentiation for the HP-41" paragraph 1-b)
-The third derivative is calculated by
 

 f '''(x) ~ [ -e f(x-5h) - d f(x-4h) - c f(x-3h) - b f(x-2h) - a f(x-h) + a f(x+h) + b f(x+2h) + c f(x+3h) + d f(x+4h) + e f(x+5h) ] / h3

  with   a = -1669/720 ,  b = 8738/5040 ,  c = -4869/10080 ,  d = 2522/30240 ,  e = -205/30240
 

Data Registers:           •  R00 = f name                               ( Register R00 is to be initialized before executing "dF" )

                                         R01 = x              R03 = f '(x)              R05 = f '''(x)        R07-R08:  temp
                                         R02 = h              R04 = f ''(x)             R06 = f(x)
Flags: /
Subroutine:   A program that takes x in X-register and returns f(x) in X-register
 
 

 01 LBL "dF"
 02 STO 01
 03 X<>Y
 04 STO 02
 05 STO 08
 06 X<>Y
 07 XEQ IND 00
 08 STO 06        
 09 6
 10 ST* 08
 11 XEQ 01
 12 8
 13 *
 14 STO 04
 15 RCL 07
 16 STO 03
 17 ST+ 03
 18 205
 19 *
 20 CHS
 21 STO 05
 22 XEQ 01
 23 125
 24 *
 25 ST- 04
 26 25
 27 RCL 07        
 28 *
 29 ST- 03
 30 2522
 31 LASTX
 32 *
 33 ST+ 05
 34 XEQ 01
 35  E3
 36 *
 37 ST+ 04
 38 150
 39 RCL 07 
 40 *
 41 ST+ 03
 42 14607
 43 LASTX
 44 *
 45 ST- 05
 46 XEQ 01
 47 6 E3
 48 *
 49 ST- 04
 50 600
 51 RCL 07        
 52 *
 53 ST- 03
 54 52428
 55 LASTX
 56 *
 57 ST+ 05
 58 XEQ 01
 59 42 E3
 60 *
 61 ST+ 04
 62 2100
 63 RCL 07 
 64 *
 65 ST+ 03
 66 70098
 67 LASTX
 68 *
 69 ST- 05
 70 RCL 02        
 71 ST/ 04
 72 2520
 73 *
 74 ST/ 03
 75 10
 76 *
 77 ST/ 04
 78 1.2
 79 *
 80 RCL 02 
 81 X^2
 82 *
 83 ST/ 05
 84 GTO 02
 85 LBL 01
 86 RCL 02
 87 ST- 08
 88 RCL 01
 89 RCL 08
 90 +
 91 XEQ IND 00
 92 STO 07        
 93 RCL 01
 94 RCL 08
 95 -
 96 XEQ IND 00
 97 RCL 07
 98 X<>Y
 99 ST- 07
100 +
101 RCL 06
102 ST+ X
103 -
104 RTN
105 LBL 02
106 RCL 05
107 RCL 04
108 RCL 03
109 END

 
      ( 190 bytes / SIZE 009 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             /         f '''(x)
           Y             h         f ''(x)
           X             x         f '(x)

 
Example:      f(x) = exp(x) + ln(x)     x = 2
 
 

 01  LBL "T"
 02  E^X
 03  LASTX
 04  LN
 05  +
 06  RTN

 
  "T"  ASTO 10

-With  h = 0.1

   0.1   ENTER^
    2     XEQ "dF"  >>>>   f '(2) = 7.889056111                                            ---Execution time = 14s---
                             RDN   f "(2) = 7.139056635
                            RDN  f '''(2) = 7.639050926

  h = 0.1 is "often" ( but not always ) a good choice.
 

      b) Functions of N Variables ( N < 10 )
 

-Formulas of order 10 are again used to estimate the partial derivatives of order 1 & 2 of a function f of n variables
 

Data Registers:   •  R00 = Function name                                 ( Registers R00 thru Rnn are to be initialized before executing "FD" or "SD"  )

                              •  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn

                                 R10 = h                                            R11 & R13 to R18 :  temp
                                 R12 =  f / xi  or  2f / xixj
Flags: F10
Subroutine:   A program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn
 
 
 

 01 LBL "FD"
 02 LBL 01
 03 STO 13
 04 CLX
 05 STO 14
 06 RCL IND 13
 07 STO 11
 08 X<> Z
 09 STO 10
 10 XEQ 03
 11 14
 12 *
 13 STO 12
 14 XEQ 03
 15 4
 16 *
 17 ST- 12
 18 XEQ 03
 19 ST+ 12
 20 6
 21 ST* 12
 22 XEQ 03
 23 ST- 12
 24 25
 25 ST* 12
 26 XEQ 03
 27 ST+ X
 28 ST+ 12
 29 2520
 30 RCL 10
 31 *
 32 ST/ 12
 33 RCL 11
 34 STO IND 13
 35 RCL 12
 36 RTN
 37 GTO 01
 38 LBL 03
 39 RCL 10
 40 ST+ 14
 41 RCL 11
 42 RCL 14
 43 -
 44 STO IND 13
 45 XEQ IND 00
 46 STO 15
 47 RCL 11
 48 RCL 14
 49 +
 50 STO IND 13
 51 XEQ IND 00
 52 RCL 15
 53 -
 54 RTN
 55 LBL "SD"
 56 LBL 02
 57 CF 10
 58 X=Y?
 59 SF 10
 60 STO 16
 61 RDN
 62 STO 17
 63 CLX
 64 STO 15
 65 RCL IND T
 66 STO 13
 67 RCL IND 17
 68 STO 14
 69 R^
 70 STO 10
 71 XEQ IND 00
 72 STO 11
 73 XEQ 03
 74 42
 75 *
 76 STO 12
 77 XEQ 03
 78 6
 79 *
 80 ST- 12
 81 XEQ 03
 82 ST+ 12
 83 XEQ 03
 84 8
 85 /
 86 ST- 12
 87  E3
 88 ST* 12
 89 XEQ 03
 90 8
 91 *
 92 ST+ 12
 93 RCL 13
 94 STO IND 16
 95 25200
 96 FC? 10
 97 ST+ X
 98 FC?C 10
 99 CHS
100 RCL 10
101 X^2
102 *
103 ST/ 12
104 RCL 12
105 RTN
106 GTO 02
107 LBL 03
108 RCL 10
109 ST+ 15
110 RCL 15
111 RCL 13
112 +
113 STO IND 16
114 XEQ IND 00
115 STO 18
116 RCL 13
117 RCL 15
118 -
119 STO IND 16
120 XEQ IND 00
121 RCL 11
122 ST+ X
123 -
124 ST+ 18
125 RCL 18
126 FS? 10
127 RTN
128 RCL 14
129 RCL 15
130 -
131 STO IND 17
132 XEQ IND 00
133 ST- 18
134 RCL 13
135 STO IND 16
136 XEQ IND 00
137 ST+ 18
138 RCL 14
139 RCL 15
140 +
141 STO IND 17
142 XEQ IND 00
143 ST+ 18
144 RCL 15
145 ST+ IND 16
146 XEQ IND 00
147 ST- 18
148 RCL 14
149 STO IND 17
150 RCL 18
151 RTN
152 END

 
     ( 254 bytes / SIZE 019 )
 
 

      STACK     FD INPUTS      OUTPUTS
           Y             h             /
           X             i    f 'xi = f / xi
      STACK     SD  INPUTS      OUTPUTS
           Z             h             /
           Y             j             /
           X             i  f "xixj = 2f / xixj

 
Example1:     f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 1 , z = 1

  LBL "T"        E^X               +                    "T"    ASTO 00
  RCL 01        RCL 02          LN                  1     STO 01   STO 02   STO 03
  X^2              X^2               *
  CHS            RCL 03          RTN

-With  h = 0.1

  0.1   ENTER^
    1    XEQ "FD"   >>>>   f 'x1 = f / x1 = -0.509989198   ( The last decimal should be a 5 )

  0.1   ENTER^
    2    XEQ "FD"   >>>>   f 'x2 = f / x2 = 0.367879440     ( The last decimal should be a 1 )

  0.1   ENTER^
    3    XEQ "FD"   >>>>   f 'x3 = f / x3 = 0.183939720     ( The last decimal should be a 1 )
 

Example2:       f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 1 , z = 1

  LBL "T"        E^X               +                    "T"    ASTO 00
  RCL 01        RCL 02          LN                  1     STO 01   STO 02   STO 03
  X^2              X^2               *
  CHS            RCL 03          RTN

-With  h = 0.1

  0.1  ENTER^
   1    ENTER^
         XEQ "SD"  >>>>   f "xx = 2f / x2 = 0.509989206   ( exact derivative = 0.509989195 )

 0.1   ENTER^
   1    ENTER^
   2    XEQ "SD"  >>>>   f "xy = 2f / xy = -0.735758889  ( the last decimal should be a 2 )
 

      c) Hessian Matrix + Gradient  ( n < 10 )
 

 "HESSGR" calculates and sores the coefficients of the Hessian matrix H and the coefficients of the gradient G into registers R25  R26 ..... in the following way:

       W =   H  Gt      where  Gt  = transposed G              So, we get a square matrix W of order (n+1)
                 G  0

  with        H = [ hi,j ]   with  hi,j = 2f / xixj
                G = [ gk ]    with  gk = f / xi

-This form is useful to calculate the Gaussian curvature of a hypersurface.
 
 

Data Registers:   •  R00 = Function name                                        ( Registers R00 thru Rnn are to be initialized before executing "HESSGR"  )

                              •  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn

                                 R10 = h       R11 to R24 :  temp    R25 ...... = W

Flag:  F10
Subroutines:   "FD"  "SD"  and a program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn
 
 
 

 01 LBL "HESSGR"
 02 STO 19
 03 STO 20            
 04 1
 05 +
 06 STO 11
 07 X^2
 08 23
 09 +
 10 STO 21
 11 STO 22
 12 X<>Y
 13 STO 10
 14 CLX
 15 SIGN
 16 +
 17 RCL 11            
 18  E-5
 19 *
 20 +
 21 STO 23
 22 CLX
 23 STO IND 23
 24 LBL 01
 25 RCL 10 
 26 RCL 20
 27 XEQ "FD"
 28 STO IND 21
 29 DSE 23
 30 STO IND 23
 31 DSE 21
 32 DSE 20
 33 GTO 01
 34 DSE 21
 35 RCL 21            
 36 RCL 23 
 37 FRC
 38 +
 39 LBL 02
 40 STO 23
 41 STO 24
 42 INT
 43 STO 21
 44 RCL 19
 45 STO 20
 46 LBL 03
 47 RCL 10
 48 RCL 19
 49 RCL 20            
 50 XEQ "SD"
 51 STO IND 21
 52 STO IND 23
 53 DSE 21 
 54 DSE 23
 55 DSE 20
 56 GTO 03
 57 RCL 24
 58 DSE X
 59 1
 60 -
 61 DSE 19
 62 GTO 02 
 63 RCL 22            
 64 1
 65 +
 66  E3
 67 /
 68 25
 69 +
 70 RCL 24
 71 FRC
 72 +
 73 END

 
        ( 130 bytes / SIZE 25+(n+1)^2 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h             /
           X          n < 10      25.eeemm

   Where  m = n+1

Example:    f(x,y,z) = x4 y3 z2 - 1    x = y = z = 1
 
 

 01  LBL "T"
 02  RCL 01
 03  X^2
 04  RCL 03
 05  *
 06  X^2
 07  RCL 02
 08  ST* Y
 09  X^2
 10  *
 11  1
 12  -
 13  RTN

 
  "T"  ASTO 00

   1  STO 01  STO 02  STO 03   and  if you choose  h = 0.1

   0.1   ENTER^
    3     XEQ "HESSGR"  >>>>    25.04003                                    ---Execution time = 124s---

-And we find the matrix W in registers  R25 to R40

      12   12   8   4
      12    6    6   3
       8     6    2   2
       4     3    2   0

                                    12  12  8
-The Hessian matrix is   12   6   6  and the gradient is  [ 4  3  2 ]
                                     8    6   2
 
 

2°)  Gaussian & Mean Curvatures of a Surface
 

     a) Surface
 

-"KGM" calculates the Gaussian curvature KG and the mean curvature Km of a surface defined by a function  z = z(x,y)

-The Gaussian Curvature KG is be computed by:

   KG = [ ( 2z / x2 ) ( 2z / y2 ) - ( 2z / xy )2 ] / [ 1 + ( z / x )2 + ( z / y )2 ]2

-The mean curvature Km is given by:

    Km = (1/2) [ t ( 1+p2 ) + r ( 1+q2 ) - 2 p q s ] / ( 1+p2+q2 )3/2

   where  p = z / x , q = z / y,  r = 2z / x2,  s = 2z / xy,  t = 2z / y2
 
 

Data Registers:           •  R00 = f name                               ( Register R00 is to be initialized before executing "KGM" )

                                         R01 to R15:  temp
Flags: /
Subroutines:   "dF2"  ( cf "Differential Geometry for the HP41" )
                          and a program that takes x in X-register, y in Y-register and returns z(x,y) in X-register
 
 

 01 LBL "KGM"
 02 XEQ "dF2"
 03 X^2
 04 X<>Y         
 05 X^2
 06 +
 07 1
 08 +
 09 STO 10 
 10 X^2
 11 X<> Z         
 12 *
 13 RCL 09
 14 X^2
 15 -
 16 X<>Y         
 17 /
 18 RCL 03 
 19 X^2
 20 RCL 07
 21 ST* Y
 22 +
 23 RCL 06
 24 X^2
 25 RCL 04 
 26 ST* Y         
 27 +
 28 +
 29 RCL 03
 30 RCL 06
 31 *
 32 RCL 09 
 33 *
 34 ST+ X         
 35 -
 36 RCL 10 
 37 ST/ Y
 38 SQRT
 39 ST+ X
 40 /
 41 X<>Y         
 42 END

 
     ( 60 bytes / SIZE 016 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             h             /
           Y             y            Km
           X             x            KG

 
Example:   A surface is defined by  z(x,y) = Ln ( 1 + x2 y3 )    Calculate the Gaussian & mean curvatures at the point  (x,y) = (1,1)
 
 

 01  LBL "T"
 02  X^2
 03  X<>Y
 04  ST* Y
 05  X^2
 06  *
 07  LN1+X
 08  RTN

 
  "T"  ASTO 00   and if you choose h = 0.1

     0.1  ENTER^
      1    ENTER^    XEQ "KGM"  >>>>   KG = -0.124572624                             ---Execution time = 58s---
                                                     X<>Y  Km = -0.171206978
 

     b) Hypersurface ( n < 10 )
 

 -Here, we assume that the hypersurface is defined implicitly by   f( x1 , ........... , xn ) = 0

Formulae:

   KG = - ( det W ) / ( - | Grad f | )n+1      where  W is defined in paragraph 1°) c)

   Km = [ Grad f x H x TGrad f  -  | Grad f |2 Trace (H) ] / (n-1) / | Grad f |3     where  H = Hessian matrix
 
 

Data Registers:   •  R00 = Function name                                        ( Registers R00 thru Rnn are to be initialized before executing "KGMN"  )

                              •  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn

                                 R10 = h       R11 to R24 :  temp    R25 ...... = W

Flag:  F10
Subroutines:   "HESSGR"  "DET"  ( the latest version listed in "Determinants for the HP41" )
                          and a program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn
 
 

 01 LBL "KGMN"
 02 XEQ "HESSGR"
 03 STO 11             
 04 FRC
 05 ENTER
 06 ISG X
 07 INT
 08 STO 12 
 09 STO 13
 10 CLX
 11  E3
 12 *
 13 INT
 14 DSE X
 15 STO 15
 16 0
 17 DSE 12
 18 LBL 01
 19 RCL IND Y
 20 X^2
 21 +
 22 DSE Y
 23 DSE 12
 24 GTO 01
 25 STO 16
 26 RCL 15             
 27 RCL 15 
 28 RCL 13
 29 STO 12
 30 STO 14
 31 -
 32 CLA
 33 DSE 12
 34 DSE 14
 35 LBL 02
 36 RCL IND Y
 37 RCL IND Y
 38 *
 39 ST+ M
 40 RDN
 41 DSE Y
 42 DSE X
 43 DSE 14
 44 GTO 02
 45 0
 46 X<> M
 47 RCL IND 15
 48 *
 49 ST+ N
 50 CLX
 51 RCL 13             
 52 DSE X
 53 STO 14 
 54 ST+ Z
 55 SIGN
 56 -
 57 DSE 15
 58 DSE 12
 59 GTO 02
 60 RCL N
 61 STO 12
 62 RCL 11
 63  E-5
 64 +
 65 0
 66 LBL 03
 67 RCL IND Y
 68 +
 69 ISG Y
 70 GTO 03
 71 RCL 16             
 72 *
 73 RCL 12 
 74 X<>Y
 75 -
 76 RCL 16
 77 ST/ Y
 78 SQRT
 79 RCL 13
 80 2
 81 -
 82 *
 83 /
 84 STO 12
 85 RCL 00
 86 STO 14
 87 RCL 11
 88 XEQ "DET"
 89 RCL 00
 90 RCL 14             
 91 STO 00
 92 CLX
 93 RCL 16 
 94 SQRT
 95 CHS
 96 RCL 13
 97 Y^X
 98 /
 99 CHS
100 STO 13
101 RCL 12
102 X<>Y
103 CLA
104 END

 
       ( 160 bytes / SIZE 25+(n+1)^2 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             /        det W
           Y             h            Km
           X         n < 10            KG

 
Example:   The surface ( n = 3 )   x4 y3 z2 - 1 = 0   with  x = y = z = 1
 
 

 01  LBL "T"
 02  RCL 01
 03  X^2
 04  RCL 03
 05  *
 06  X^2
 07  RCL 02
 08  ST* Y
 09  X^2
 10  *
 11  1
 12  -
 13  RTN

 
  "T"  ASTO 00

   1  STO 01  STO 02  STO 03   and  if you choose  h = 0.1

   0.1   ENTER^
    3     XEQ "KGMN"  >>>>      KG  =  0.256837099                             ---Execution time = 151s---
                                     RDN      Km  =  0.518666290
                                    X<>Y  det W = - 216

Warning:

-This program does not check that  f(x1, ..... ,xn) = 0
-The matrix W is not saved.
 

3°)  Curvature & Torsion of a Curve
 

     a) 3-Dimensional Problem
 

-We assume that the curve is parameterized by 3 functions  x(t) , y(t) , z(t)

-The curvature is calculated by   khi = | r' x r'' | / | r' |3                                      where  the vector  r = [ x(t) , y(t) , z(t) ]
                   and the torsion by   tau = ( r' x r'' ) .r''' / | r' x r'' |2                        x = cross-product  and  . = dot-product
 
 

Data Registers:              R00 & R05 to R15:  temp
                                         R01 = t  R02 = h   R03 = khi   R04 = tau
Flags: /
Subroutines:    "dF"  ( cf paragraph 1 )
                   and  3 functions named  "X"  "Y"  "Z"   that compute x(t) , y(t) , z(t)  assuming t is in X-register upon entry
 
 
 

 01 LBL "KHT"
 02 "Z"
 03 ASTO 00
 04 XEQ "dF"
 05 STO 12
 06 RDN
 07 STO 13 
 08 X<>Y
 09 STO 14
 10 "Y"
 11 ASTO 00
 12 RCL 02       
 13 RCL 01
 14 XEQ "dF"
 15 STO 09
 16 RDN
 17 STO 10
 18 X<>Y
 19 STO 11
 20 "X"
 21 ASTO 00
 22 RCL 02
 23 RCL 01 
 24 XEQ "dF"
 25 RCL 12
 26 RCL 09
 27 R-P
 28 RCL 03       
 29 R-P
 30 STO 00
 31 RCL 13
 32 RCL 09
 33 *
 34 RCL 12
 35 RCL 10
 36 *
 37 -
 38 X^2
 39 STO 06 
 40 LASTX
 41 RCL 05
 42 STO 08
 43 *
 44 RCL 12       
 45 RCL 04
 46 STO 07
 47 *
 48 RCL 13
 49 RCL 03
 50 *
 51 -
 52 ENTER
 53 X^2
 54 ST+ 06
 55 CLX
 56 RCL 11 
 57 *
 58 +
 59 RCL 10
 60 RCL 03       
 61 *
 62 RCL 09
 63 RCL 07
 64 *
 65 -
 66 ENTER
 67 X^2
 68 ST+ 06
 69 CLX
 70 RCL 14
 71 *
 72 +
 73 RCL 06 
 74 ST/ Y
 75 SQRT
 76 RCL 00       
 77 3
 78 Y^X
 79 /
 80 X<>Y
 81 STO 04
 82 X<>Y
 83 STO 03
 84 END

 
      ( 110 bytes / SIZE 015 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h          tau(t)
           X             t          khi(t)

 
Example:   The curve defined by   x(t) = et cos t  ,  y(t) = et sin t  ,  z(t) = Ln(t)        ( t expressed in radians )

-Evaluate the curvature khi (t) and the torsion tau (t) for t = 1.3

-Load the following routines in main memory
 
 

 01  LBL "X"
 02  E^X
 03  LASTX
 04  COS
 05  *
 06  RTN
 07  LBL "Y"
 08  E^X
 09  LASTX
 10  SIN
 11  *
 12  RTN
 13  LBL "Z"
 14  LN
 15  RTN

 
   XEQ "RAD"   and if you choose  h = 0.1

    0.1   ENTER^
    1.3   XEQ "KHT"  >>>>  khi (1.3)  = 0.194807823                       ---Execution time = 46s---
                                 X<>Y  tau (1.3)  = 0.123665257

-The exact values are   khi (1.3) = 0.194807823  &  tau (1.3) = 0.123665529
 

     b) N-Dimensional Problem
 

-The method follows the Gram-Schmidt process ( cf references [1] & [2] ). After some calculations, we get:

-With    r = [ x1(t) , x2(t) , ........ , xn(t) ]

     •   khi = | e'1 | / | r' |

 where   e1 = r' / | r' |    whence   e'1 = r'' / | r' | - r' < r' , r'' > / | r' |3                         < , > = dot product.

     •   tau = < e'2 , ê3 > / ( | r' | | ê3 | )

  where   e'2 = e''1 / | e'1 | - e'1 < e'1 , e''1 > / | e'1 |3
             e''1 = r''' / | r' | - 2 r'' < r' , r'' > / | r' |3 + 3 r' < r' , r'' >2 / | r' |5 - r' | r'' |2 / | r' |3 - r' < r' , r''' >2 / | r' |3
   and     ê3  = r''' - r' < r' , r''' > / | r' |2 - e'1 < e'1 , r''' > / | e'1 |2
 
 

Data Registers:              R00 & R05 to R09 & R11 to R5n+12:  temp

                                         R01 = t    R03 = khi
                                         R02 = h   R04 = tau
Flags: /
Subroutines:    "dF"  ( cf paragraph 1 )
                   and  n functions named  "X1"  "X2" ........... "XN"   that compute x1(t) , x2(t) , ........ , xn(t)  assuming t is in X-register upon entry
 
 
 

 01 LBL "KHTN"
 02 STO 01
 03 RDN
 04 STO 02         
 05 X<>Y
 06 STO 09
 07 STO 10
 08 FIX 00
 09 CF 29
 10 LBL 01
 11 "X"
 12 ARCL 09
 13 ASTO 00
 14 RCL 02
 15 RCL 01
 16 XEQ "dF"
 17 RCL 09
 18 12
 19 +
 20 X<>Y
 21 STO IND Y 
 22 CLX
 23 RCL 10
 24 +
 25 RCL 04
 26 STO IND Y
 27 CLX
 28 RCL 10
 29 +
 30 RCL 05
 31 STO IND Y
 32 DSE 09
 33 GTO 01
 34 FIX 09
 35 SF 29
 36 RCL 10
 37 12.012
 38 +
 39 STO 05
 40 INT
 41 RCL 10
 42 +
 43 STO 06         
 44 LASTX
 45 +
 46 STO 07
 47 LASTX
 48 +
 49 STO 08
 50 LASTX
 51 +
 52 STO 09
 53 RCL 05
 54 INT
 55 XEQ 11
 56 STO 03
 57 RCL 06
 58 RCL 05
 59 INT
 60 XEQ 13
 61 STO 04
 62 RCL 06
 63 RCL 05
 64 LBL 02
 65 RCL IND X 
 66 RCL 04
 67 *
 68 RCL 03
 69 X^2
 70 /
 71 RCL IND Z
 72 X<>Y
 73 -
 74 RCL 03
 75 /
 76 STO IND 08
 77 RDN
 78 DSE 08
 79 DSE Y
 80 DSE X
 81 GTO 02
 82 RCL 10         
 83 ST+ 08
 84 RCL 08
 85 XEQ 11
 86 STO 11
 87 RCL 06
 88 XEQ 11
 89 STO 00
 90 RCL 07
 91 RCL 05
 92 INT
 93 XEQ 13
 94 STO 12
 95 LBL 03
 96 RCL 04
 97 RCL 03
 98 /
 99 X^2
100 3
101 *
102 RCL 00
103 X^2
104 -
105 RCL 12
106 -
107 RCL IND 05
108 *
109 RCL 04
110 ST+ X
111 RCL IND 06
112 *
113 -
114 RCL 03
115 X^2
116 /
117 RCL IND 07
118 +
119 RCL 03         
120 /
121 STO IND 09
122 DSE 09
123 DSE 07
124 DSE 06
125 DSE 05
126 GTO 03
127 RCL 10
128 ST+ 05
129 ST+ 06
130 ST+ 07
131 ST+ 09
132 RCL 08
133 RCL 09
134 XEQ 13
135 STO 00
136 RCL 10
137 STO M
138 LBL 04
139 RCL IND 09
140 RCL IND 08
141 RCL 00
142 *
143 RCL 11
144 X^2
145 /
146 -
147 RCL 11
148 /
149 STO IND 09
150 DSE 09
151 DSE 08
152 DSE M
153 GTO 04
154 RCL 10
155 ST+ 08
156 ST+ 09
157 RCL 07
158 RCL 08         
159 XEQ 13
160 STO 00
161 LBL 05
162 RCL IND 07
163 RCL IND 05
164 RCL 12
165 *
166 RCL 03
167 X^2
168 /
169 -
170 RCL IND 08
171 RCL 00
172 *
173 RCL 11
174 X^2
175 /
176 -
177 STO IND 08
178 DSE 08
179 DSE 07
180 DSE 05
181 GTO 05
182 RCL 10
183 ST+ 05
184 ST+ 07
185 ST+ 08
186 RCL 08
187 XEQ 11
188 STO 00
189 RCL 08
190 RCL 09
191 XEQ 13
192 GTO 00
193 LBL 11
194 RCL 10
195 STO M
196 CLX
197 LBL 12
198 RCL IND Y
199 X^2
200 +
201 DSE Y
202 DSE M
203 GTO 12
204 SQRT
205 RTN
206 LBL 13
207 CLA
208 RCL 10         
209 STO N
210 RDN
211 LBL 14
212 RCL IND Y
213 RCL IND Y
214 *
215 ST+ M
216 RDN
217 DSE Y
218 DSE X
219 DSE N
220 GTO 14
221 CLX
222 X<> M
223 RTN
224 LBL 00
225 RCL 03
226 RCL 00
227 *
228 /
229 STO 04
230 RCL 11
231 RCL 03
232 /
233 STO 03
234 RCL 10
235 3
236 X#Y?
237 GTO 00
238 RCL 14 
239 RCL 18
240 *
241 RCL 15
242 RCL 17
243 *
244 -
245 RCL 19         
246 *
247 RCL 13
248 RCL 18
249 *
250 RCL 15
251 RCL 16
252 *
253 -
254 RCL 20
255 *
256 -
257 RCL 13
258 RCL 17
259 *
260 RCL 14
261 RCL 16
262 *
263 -
264 RCL 21
265 *
266 +
267 SIGN
268 ST* 04
269 LBL 00
270 RCL 04
271 RCL 03
272 END

 
      ( 389 bytes / SIZE 5n+13 )
 
 

      STACK        INPUTS      OUTPUTS
           Z          n < 62             /
           Y             h          tau(t)
           X             t          khi(t)

 
Example:   In a 4-dimensional space, the curve defined by

     x1(t) = t4       x3(t) = t3      at the point  t = 1
     x2(t) = t2       x4(t) = et
 
 

 01  LBL "X1"
 02  X^2
 03  X^2
 04  RTN
 05  LBL "X2"
 06  X^2
 07  RTN
 08  LBL "X3"
 09  ENTER^
 10  X^2
 11  *
 12  RTN
 13  LBL "X4"
 14  E^X
 15  RTN

 
-If you choose  h = 0.1

     4   ENTER^
   0.1  ENTER^
     1   XEQ "KHTN"   >>>>   khi (1)  = 0.142277239                       ---Execution time = 71s---
                                   X<>Y  tau (1)  = 0.121469324

Note:

-In the Gram-Schmidt process, the orthogonal basis does not always satisfy the right-hand rule.
-So, if n = 3 , the torsion given by the formulas above is multiplied by sign [ det ( r' , r'' , r''' ) ] ( lines 238 to 268 )
-Otherwise, the sign could be different from the result obtained by "KHT"
 
 
 

References:

[1]  https://en.wikipedia.org/w/index.php?title=Differential_geometry_of_curves&oldid=740874421
[2]  https://en.wikipedia.org/w/index.php?title=Frenet%E2%80%93Serret_formulas&oldid=731811089
[3]  Ron Goldman - "Curvature formulas for implicit curves and surfaces"