Clebsch-Gordan Coefficients & Wigner 3j-6j-9j-symbols for the HP-41

Overview
 

 1°)  Clebsch-Gordan Coefficients & Wigner 3j-symbol  ( Program#1 )
 2°)  Triangular Coefficients
 3°)  Clebsch-Gordan Coefficients & Wigner 3j-symbol  ( Program#2 )
 4°)  Wigner 6j-symbol
 5°)  Wigner 9j-symbol
 

-These numbers are useful to describe the angular momentum in quantum physics.
-All the parameters j_i , m_i must be integers or half-integers.

-The standard notations for Wigner symbols use parentheses and brace brackets beside 2x3 or 3x3 arrays. This is not easy in html.
-So, I've used non-standard notations like W6j ( j₁ j₂ j₃ j₄ j₅ j₆ ) for the Wigner 6j- symbols...
 

1°)  Clebsch-Gordan Coefficients & Wigner 3j-symbol  ( Program#1 )
 

-This program returns  <  j₁ j₂  m₁ m₂ |   j₁ j₂ j m  >  in X-register
                         and      W3J (  j₁ j₂ j  m₁ m₂ -m )    in Y-register

Formulae:
 

   •  Clebsch-Gordan Coefficients  = <  j₁ j₂  m₁ m₂ |   j₁ j₂ j m  >
                                                    = [ ( 2j +1) D( j₁ j₂ j  ) ]^1/2  [ ( j₁ + m₁ )! ( j₁ - m₁ )! ( j₂ + m₂ )! ( j₂ - m₂ )! ( j + m )! ( j - m )! ]
                                   x Sum_k  (-1)^k / [ k! ( j₁ + j₂ - j - k )! ( j₁ - m₁ - k )! ( j₂ + m₂ - k )! ( j - j₂ + m₁ + k )! ( j - j₁ - m₂ + k )! ]         if   m₁ + m₂ = m

      ( The sum is to be performed for all integers k that produce non-negative integer arguments for the factorials ).

   •  <  j₁ j₂  m₁ m₂ |   j₁ j₂ j m  > = 0            if   m₁ + m₂ # m

   •  Wigner 3j-symbol  = W3J (  j₁ j₂ j  m₁ m₂ -m ) = (-1)^( j₁ - j₂ + m )  ( 2j +1) ^-1/2  <  j₁ j₂  m₁ m₂ |   j₁ j₂ j m  >

  •   Triangular coefficients  D( a b c )  are defined in paragraph 2
 

Data Registers:              R00 & R07 thru R13: temp                      ( Registers R01 thru R06 are to be initialized before executing "CGW3J" )

                                      •  R01 = j₁     •  R02 = j₂    •  R03 = j
                                      •  R04 = m₁   •  R05 = m₂   •  R06 = m
Flags: /
Subroutines: /

-If  m₁ + m₂ # m , the results = 0
-Otherwise, they are less easy to compute...
 
 

01  LBL "CGW3J"   02  RCL 04   03  RCL 05   04  +   05  RCL 06   06  X=Y?   07  GTO 00   08  CLST   09  RTN   10  LBL 00   11  RCL 02   12  RCL 03   13  -   14  STO 00   15  RCL 04   16  -   17  RCL 01   18  RCL 03   19  -   20  STO 09   21  RCL 05   22  +   23  X<Y?   24  X<>Y   25  X<0?   26  CLX   27  STO 07   28  RCL 01   29  RCL 02   30  +   31  STO 10   32  RCL 03   33  -   34  STO 11

35  RCL 01   36  RCL 04   37  -   38  STO 12            39  X>Y?   40  X<>Y   41  RCL 02   42  RCL 05   43  +   44  STO 13   45  X>Y?   46  X<>Y   47  STO 08   48  CLX   49  LBL 01   50  RCL 07   51  RCL 08   52  X<Y?   53  GTO 02   54  RDN   55  FACT   56  RCL 11   57  RCL 07   58  -   59  FACT   60  *   61  RCL 12   62  RCL 07   63  -   64  FACT   65  *   66  RCL 13   67  RCL 07   68  -

69  FACT   70  *   71  RCL 04   72  RCL 07            73  +   74  RCL 00   75  -   76  FACT   77  *   78  RCL 07   79  RCL 05   80  -   81  RCL 09   82  -   83  FACT   84  *   85  1   86  ST+ 07   87  CHS   88  RCL 07   89  Y^X   90  X<>Y   91  /   92  -   93  GTO 01   94  LBL 02   95  X<> Z   96  RCL 11   97  FACT   98  SQRT   99  * 100  RCL 01 101  RCL 00 102  -

103  FACT 104  SQRT 105  * 106  RCL 02          107  RCL 09 108  - 109  FACT 110  SQRT 111  * 112  RCL 01 113  RCL 04 114  + 115  FACT 116  SQRT 117  * 118  RCL 12 119  FACT 120  SQRT 121  * 122  RCL 13 123  FACT 124  SQRT 125  * 126  RCL 02 127  RCL 05 128  - 129  FACT 130  SQRT 131  * 132  RCL 03 133  RCL 06 134  + 135  FACT 136  SQRT

137  * 138  RCL 03 139  RCL 06          140  - 141  FACT 142  SQRT 143  * 144  RCL 03 145  RCL 10 146  + 147  1 148  + 149  FACT 150  SQRT 151  / 152  1 153  CHS 154  RCL 01 155  RCL 02 156  - 157  RCL 06 158  + 159  Y^X 160  * 161  X<>Y 162  RCL 03 163  ST+ X 164  1 165  + 166  SQRT 167  * 168  END

 
   ( 184 bytes / SIZE 014 )
 
 

STACK	INPUTS	OUTPUTS
Y	/	W3j
X	/	ClG

   where    W3j = W3j (  j₁ j₂ j  m₁ m₂ -m )     and   ClG = <  j₁ j₂  m₁ m₂ |   j₁ j₂ j m  >

Examples:

-If   j₁ = 12  ,   j₂ = 24  ,  j  = 31
      m₁ = 1   ,  m₂ = 16  , m = 17

 XEQ "CLGW3J"  returns   ClG =  0.206754081                 ---Execution time = 15s---
                             X<>Y   W3j = -0.026048566

-With   j₁ = 1  ,   j₂ = 3/2  ,  j  = 5/2
           m₁ = 0 ,  m₂ = 3/2  , m = 3/2

  XEQ "CLGW3J"  yields   ClG =  0.632455533                 ---Execution time = 5s---
                             X<>Y   W3j = -0.258198890
 

Notes:

-All the j's & m's must be integers or half-integers
-They must satisfy several inequalities so that  FACT is applied to non-negative integers only.
-Otherwise, DATA ERROR will be displayed.
-If you want to get  W3j (  j₁ j₂ j  m₁ m₂ m )  store (-m) into R06

-You can save a few bytes if you perform SQRT after the different products and quotient, but there will be more risks of overflows.
-Furthermore, the HP-41 actually computes n! up to n = 99 internally, so we can write a M-Code function to calculate sqrt(n!) and avoid an "OUT OF RANGE"

-Replace for instance  FACT  SQRT   by the M-Code function  SQRTF  listed below ( it does not check for alpha data ):
 

086   "F"
014   "T"
012   "R"
011   "Q"
013   "S"
2A0  SETDEC
0F8   READ 3(X)
006   A=0 S&X
166   A=A+1 S&X
306   ?A<C S&X
0B5   ?C GO
0A3   282D            GTO DATA ERROR if x is not an integer or if x > 99
128   WRIT 4(L)
3B1   ?NCXQ         C=
060   18EC            fact(C)
2F9   ?NCQ            C=
060   18BE            sqrt(C)
0E8   WRIT 3(X)
3E0   RTN

-After these modifications,

with   j₁ = 14  ,   j₂ = 29  ,  j  = 41
        m₁ = 1   ,  m₂ = 16  , m = 17

 XEQ "CLGW3J"   gives    ClG =  0.361681799                 ---Execution time = 10s---
                             X<>Y   W3j =  0.039699735

 ( sqrt ( 85! ) has been computed )
 

2°)  Triangular coefficients
 

-Wigner 3j- 6j- and 9j- symbols use the calculation of

    Delta ( a , b , c ) = [ ( a + b - c )! ( a - b + c )! ( -a + b + c )! ] / ( a + b + c + 1 )!

-The routine hereunder computes  SQRT [ Delta ( a , b , c ) ]
-FACT  SQRT  may be replaced by the M-Code routine  SQRTF listed above...
   a , b , c  are stored into  R11 , R12 , R13  respectively
   we can also use synthetic registers M , N , O
 
 

01  LBL "DABC" 02  STO 13         03  X<>Y 04  STO 12 05  - 06  X<>Y 07  STO 11 08  +

09  FACT 10  SQRT 11  RCL 11 12  RCL 12         13  + 14  RCL 13 15  - 16  FACT

17  SQRT 18  * 19  RCL 12         20  RCL 11 21  - 22  RCL 13 23  + 24  FACT

25  SQRT 26  * 27  RCL 11 28  RCL 12         29  + 30  RCL 13 31  + 32  1

33  + 34  FACT 35  SQRT           36  / 37  END

 
   ( 46 bytes )
 
 

STACK	INPUTS	OUTPUTS
T	T	T
Z	a	T
Y	b	T
X	c	sqrt(delta(a,b,c))

 
Example:

   12  ENTER^
   24  ENTER^
   31  XEQ "DABC"  returns  5.963240385 E-13
 
 

3°)  Clebsch-Gordan Coefficients & Wigner 3j-symbol  ( Program#2 )
 

-We can now rewrite "CGW3J" to use "DABC"
-This will save bytes if you plan to use also "W6J" and/or "W9J".
-Otherwise, it's wasteful!
 

Data Registers:              R00 & R07 thru R13: temp                      ( Registers R01 thru R06 are to be initialized before executing "CGW3J" )

                                      •  R01 = j₁     •  R02 = j₂    •  R03 = j
                                      •  R04 = m₁   •  R05 = m₂   •  R06 = m
Flags: /
Subroutine:  "DABC"
 
 

01  LBL "CGW3J"   02  RCL 04   03  RCL 05   04  +   05  RCL 06   06  X=Y?   07  GTO 00   08  CLST   09  RTN   10  LBL 00   11  RCL 02   12  RCL 03   13  -   14  RCL 04   15  -   16  STO 00   17  RCL 01   18  RCL 03   19  -   20  RCL 05   21  +   22  STO 08   23  X<Y?   24  X<>Y   25  X<0?   26  CLX   27  STO 07   28  RCL 01   29  RCL 02

30  +   31  RCL 03   32  -   33  STO 09   34  RCL 01   35  RCL 04            36  -   37  STO 10   38  X>Y?   39  X<>Y   40  RCL 02   41  RCL 05   42  +   43  STO 11   44  X>Y?   45  X<>Y   46  STO 12   47  CLX   48  LBL 01   49  RCL 07   50  RCL 12   51  X<Y?   52  GTO 02   53  RDN   54  FACT   55  RCL 09   56  RCL 07   57  -   58  FACT

59  *   60  RCL 10   61  RCL 07   62  -   63  FACT   64  *   65  RCL 11            66  RCL 07   67  -   68  FACT   69  *   70  RCL 07   71  RCL 00   72  -   73  FACT   74  *   75  RCL 07   76  RCL 08   77  -   78  FACT   79  *   80  1   81  ST+ 07   82  CHS   83  RCL 07   84  Y^X   85  X<>Y   86  /   87  -

88  GTO 01   89  LBL 02   90  X<> Z   91  RCL 01   92  RCL 04   93  +   94  FACT   95  SQRT   96  *   97  RCL 10            98  FACT   99  SQRT 100  * 101  RCL 11 102  FACT 103  SQRT 104  * 105  RCL 02 106  RCL 05 107  - 108  FACT 109  SQRT 110  * 111  RCL 03 112  RCL 06 113  + 114  FACT 115  SQRT 116  *

117  RCL 03 118  RCL 06 119  - 120  FACT 121  SQRT 122  * 123  RCL 01 124  RCL 02 125  RCL 03 126  XEQ "DABC" 127  * 128  1 129  CHS 130  RCL 01 131  RCL 02 132  - 133  RCL 06 134  + 135  Y^X 136  * 137  X<>Y 138  RCL 03 139  ST+ X 140  1 141  + 142  SQRT 143  * 144  END

 
    ( 165 bytes / SIZE 014 )
 
 

STACK	INPUTS	OUTPUTS
Y	/	W3j
X	/	ClG

   where    W3j = W3j (  j₁ j₂ j  m₁ m₂ -m )     and   ClG = <  j₁ j₂  m₁ m₂ |   j₁ j₂ j m  >

Example:

-If   j₁ = 12  ,   j₂ = 24  ,  j  = 31
      m₁ = 1   ,  m₂ = 16  , m = 17

 XEQ "CLGW3J"  returns   ClG =  0.206754081                 ---Execution time = 15s---
                             X<>Y   W3j = -0.026048566
 

4°)  Wigner 6j-symbol
 

Formula:
 

     W6J ( j₁ j₂ j₃ j₄ j₅ j₆ )  =  [ D ( j₁ j₂ j₃ ) D ( j₁ j₅ j₆ ) D ( j₂ j₃ j₄ ) D ( j₃ j₄ j₅ ) ]^1/2

                                     x Sum_{k=k' to k"}  ( -1 )^k ( k + 1 )! / [ ( k - j₁ - j₂ - j₃ )! ( k - j₁ - j₅ - j₆ )! ( k - j₂ - j₃ - j₄ )! ( k - j₃ - j₄ - j₅ )!
                                                                                           ( j₁ + j₂ + j₄ + j₅ - k )! ( j₁ + j₃ + j₄ + j₆ - k )! ( j₂ + j₃ + j₅ + j₆ - k )! ]

       with   k' = max { j₁ + j₂ + j₃ , j₁ + j₅ + j₆ , j₂ + j₃ + j₄ , j₃ + j₄ + j₅ }
       and    k" = min { j₁ + j₂ + j₄ + j₅ , j₁ + j₃ + j₄ + j₆ , j₂ + j₃ + j₅ + j₆ }
 

Data Registers:              R00 & R07 thru R14: temp                      ( Registers R01 thru R06 are to be initialized before executing "W6J" )

                                      •  R01 = j₁     •  R02 = j₂    •  R03 = j₃
                                      •  R04 = j₄     •  R05 = j₅    •  R06 = j₆
Flags: /
Subroutine:  "DABC"
 
 

01  LBL "W6J"   02  RCL 01   03  STO 12   04  RCL 02   05  STO 13   06  +   07  STO 11            08  RCL 03   09  ST+ 13   10  +   11  STO 00   12  LASTX   13  RCL 04   14  ST+ 11   15  +   16  ST+ 12   17  RCL 05   18  ST+ 11   19  +   20  STO 08   21  X<Y?   22  X<>Y   23  RCL 01   24  LASTX   25  RCL 06

26  ST+ 12   27  +   28  ST+ 13   29  +   30  STO 09            31  X<Y?   32  X<>Y   33  RCL 02   34  RCL 04   35  +   36  RCL 06   37  +   38  STO 10   39  X<Y?   40  X<>Y   41  STO 07   42  RCL 11   43  RCL 12   44  X>Y?   45  X<>Y   46  RCL 13   47  X>Y?   48  X<>Y   49  STO 14   50  CLX

51  LBL 01   52  RCL 07   53  RCL 14            54  X<Y?   55  GTO 02   56  CLX   57  RCL 00   58  -   59  FACT   60  RCL 07   61  RCL 08   62  -   63  FACT   64  *   65  RCL 07   66  RCL 09   67  -   68  FACT   69  *   70  RCL 07   71  RCL 10   72  -   73  FACT   74  *   75  RCL 11

76  RCL 07   77  -   78  FACT   79  *   80  RCL 12   81  RCL 07            82  -   83  FACT   84  *   85  RCL 13   86  RCL 07   87  -   88  FACT   89  *   90  1   91  ST+ 07   92  CHS   93  RCL 07   94  Y^X   95  LASTX   96  FACT   97  *   98  X<>Y   99  / 100  -

101  GTO 01 102  LBL 02 103  X<> Z 104  RCL 01 105  RCL 02 106  RCL 03 107  XEQ "DABC" 108  * 109  RCL 03 110  RCL 04 111  RCL 05 112  XEQ "DABC" 113  * 114  RCL 01 115  RCL 05 116  RCL 06 117  XEQ "DABC" 118  * 119  RCL 02 120  RCL 04 121  RCL 06 122  XEQ "DABC" 123  * 124  END

 
    ( 162 bytes / SIZE 015 )
 
 

STACK	INPUTS	OUTPUTS
X	/	W6j

  with  W6j = W6j ( j₁ j₂ j₃ j₄ j₅ j₆ )

Example:

-If   j₁ = 10  ,   j₂ = 16  ,  j₃  = 21    store them     R01   R02   R03
       j₄ = 24  ,   j₅ = 12  ,  j₆  = 14         into           R04   R05   R06

 XEQ "W6J"  >>>>   W6j  =  0.006251585                 ---Execution time = 21s---
 

Notes:

-All the j's must be integers or half-integers
-They must satisfy several inequalities so that FACT is applied to non-negative integers only.
-Otherwise, DATA ERROR will be displayed.
 
 

5°)  Wigner 9j-symbol
 

Formula:
 

     W9J ( j₁ j₂ j₃ j₄ j₅ j₆ j₇ j₈ j₉ )  =  Sum_{k=k' to k"}  (-1)^2k (2k+1)  W6J ( j₁ j₉ k j₈ j₄ j₇ )  W6J ( j₈ j₄ k j₆ j₂ j₅ )  W6J ( j₆ j₂ k j₁ j₉ j₃ )
 

    where    k' = max { | j₁ - j₉ | , | j₄ - j₈ | , | j₂ - j₆ | }    and    k" =  min {  j₁ + j₉  ,  j₄ + j₈  ,  j₂ + j₆ }
 
 

Data Registers:              R00 & R10 thru R27: temp                      ( Registers R01 thru R09 are to be initialized before executing "W9J" )

                                      •  R01 = j₁     •  R02 = j₂    •  R03 = j₃
                                      •  R04 = j₄     •  R05 = j₅    •  R06 = j₆              registers R01 to R09 are saved in registers R15 to R23
                                      •  R07 = j₇     •  R08 = j₈    •  R09 = j₉              and they are restored at the end
Flags: /
Subroutine:  "W6J"

-If you don't have an X-Functions Module, lines 02-03-04  may be replaced by  1.009   15   XEQ "LCO"
-And lines 94-95  may be replaced by    15.023   1   XEQ "LCO"
  ( cf "Miscellaneous Short Routines for the HP-41" )
 
 

01  LBL "W9J" 02  1.015009 03  STO 27 04  REGMOVE 05  RCL 01       06  RCL 09 07  - 08  ABS 09  RCL 04 10  RCL 08 11  - 12  ABS 13  X<Y? 14  X<>Y 15  RCL 02 16  RCL 06 17  - 18  ABS 19  X<Y? 20  X<>Y

21  STO 03 22  RCL 01 23  RCL 09 24  + 25  RCL 04       26  RCL 08 27  + 28  X>Y? 29  X<>Y 30  RCL 02 31  RCL 06 32  + 33  X>Y? 34  X<>Y 35  STO 25 36  CLX 37  STO 24 38  LBL 01 39  RCL 03 40  RCL 25

41  X<Y? 42  GTO 02 43  RCL 15 44  STO 01       45  RCL 23 46  STO 02 47  RCL 22 48  STO 04 49  RCL 18 50  STO 05 51  RCL 21 52  STO 06 53  XEQ "W6J" 54  STO 26 55  RCL 04 56  STO 01 57  RCL 05 58  STO 02 59  RCL 20 60  STO 04

61  RCL 16 62  STO 05 63  RCL 19 64  STO 06 65  XEQ "W6J" 66  ST* 26 67  RCL 04 68  STO 01 69  RCL 05 70  STO 02 71  RCL 15 72  STO 04 73  RCL 23 74  STO 05 75  RCL 17 76  STO 06 77  XEQ "W6J" 78  RCL 26 79  * 80  1

81  CHS 82  RCL 03 83  ST+ X 84  Y^X 85  LASTX 86  1 87  ST+ 03 88  + 89  * 90  * 91  ST+ 24 92  GTO 01 93  LBL 02 94  RCL 27       95  REGSWAP 96  RCL 24 97  END

 
    ( 149 bytes / SIZE 028 )
 
 

STACK	INPUTS	OUTPUTS
X	/	W9j

    where  W9j = W9J ( j₁ j₂ j₃ j₄ j₅ j₆ j₇ j₈ j₉ )

Examples:

         j₁ = 3  ,   j₂ = 7  ,  j₃  = 5              R01   R02   R03
-If     j₄ = 6  ,   j₅ = 8  ,  j₆  = 9      =      R04   R05   R06
         j₇ = 4  ,   j₈ = 5  ,  j₉  = 7              R07   R08   R09

 XEQ "W9J"  >>>>   W9j  =  0.001032402067               ---Execution time = 4m32s---

               j₁ = 1/2  ,   j₂ =  1    ,  j₃  = 3/2
-With      j₄ =  1    ,   j₅ = 1/2  ,  j₆  = 1/2
               j₇ = 3/2  ,   j₈ = 1/2  ,  j₉  =  1

     XEQ "W9J"  >>>>   W9j  =  -0.02777777771   ( in fact -1/36 )         ---Execution time = 58s---

Notes:

-All the j's must be integers or half-integers
-They must satisfy several inequalities so that  FACT is applied to non-negative integers only.
-Otherwise, DATA ERROR will be displayed.
 

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   http://functions.wolfram.com/