code: def test_couple_3_states_numerical():
      # Default coupling
      # j1=1/2,j2=1/2,j3=1/2
      assert couple(
          TensorProduct(
              JzKet(S.Half, S.Half),
              JzKet(S.Half, S.Half),
              JzKet(S.Half, S.Half),
          )
      ) == JzKetCoupled(
          Rational(3, 2),
          S(3) / 2,
          (S.Half, S.Half, S.Half),
          ((1, 2, 1), (1, 3, Rational(3, 2))),
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
              )
          )
          == sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, S.Half),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 3
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, S.One / 2),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, S.Half),
              ((1, 2, 0), (1, 3, S.Half)),
          )
          / 2
          - sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, S.Half),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, S.One / 2),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, S.Half),
              ((1, 2, 0), (1, 3, S.Half)),
          )
          / 2
          + sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, S.Half),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, S.One / 2),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
              )
          )
          == -sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, S.Half),
              ((1, 2, 0), (1, 3, S.Half)),
          )
          / 2
          - sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, S.Half),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, S.One / 2),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
              )
          )
          == -sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, S.Half),
              ((1, 2, 0), (1, 3, S.Half)),
          )
          / 2
          + sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, S.Half),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, S.One / 2),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
              )
          )
          == -sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, S.Half),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 3
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, S.One / 2),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 3
      )
      assert couple(
          TensorProduct(
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(S.Half, Rational(-1, 2)),
          )
      ) == JzKetCoupled(
          Rational(3, 2),
          -S(3) / 2,
          (S.Half, S.Half, S.Half),
          ((1, 2, 1), (1, 3, Rational(3, 2))),
      )
      # j1=S.Half, j2=S.Half, j3=1
      assert couple(
          TensorProduct(
              JzKet(S.Half, S.Half), JzKet(S.Half, S.Half), JzKet(1, 1)
          )
      ) == JzKetCoupled(
          2, 2, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(1, 0),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              1, 1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 1))
          )
          / 2
          + sqrt(2)
          * JzKetCoupled(
              2, 1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
          )
          / 2
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(1, -1),
              )
          )
          == sqrt(3)
          * JzKetCoupled(
              0, 0, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 0))
          )
          / 3
          + sqrt(2)
          * JzKetCoupled(
              1, 0, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 1))
          )
          / 2
          + sqrt(6)
          * JzKetCoupled(
              2, 0, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
          )
          / 6
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 1),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              1, 1, (S.Half, S.Half, 1), ((1, 2, 0), (1, 3, 1))
          )
          / 2
          - JzKetCoupled(
              1, 1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 1))
          )
          / 2
          + JzKetCoupled(
              2, 1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
          )
          / 2
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 0),
              )
          )
          == -sqrt(6)
          * JzKetCoupled(
              0, 0, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 0))
          )
          / 6
          + sqrt(2)
          * JzKetCoupled(
              1, 0, (S.Half, S.Half, 1), ((1, 2, 0), (1, 3, 1))
          )
          / 2
          + sqrt(3)
          * JzKetCoupled(
              2, 0, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, -1),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              1, -1, (S.Half, S.Half, 1), ((1, 2, 0), (1, 3, 1))
          )
          / 2
          + JzKetCoupled(
              1, -1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 1))
          )
          / 2
          + JzKetCoupled(
              2, -1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
          )
          / 2
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(1, 1),
              )
          )
          == -sqrt(2)
          * JzKetCoupled(
              1, 1, (S.Half, S.Half, 1), ((1, 2, 0), (1, 3, 1))
          )
          / 2
          - JzKetCoupled(
              1, 1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 1))
          )
          / 2
          + JzKetCoupled(
              2, 1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
          )
          / 2
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(1, 0),
              )
          )
          == -sqrt(6)
          * JzKetCoupled(
              0, 0, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 0))
          )
          / 6
          - sqrt(2)
          * JzKetCoupled(
              1, 0, (S.Half, S.Half, 1), ((1, 2, 0), (1, 3, 1))
          )
          / 2
          + sqrt(3)
          * JzKetCoupled(
              2, 0, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(1, -1),
              )
          )
          == -sqrt(2)
          * JzKetCoupled(
              1, -1, (S.Half, S.Half, 1), ((1, 2, 0), (1, 3, 1))
          )
          / 2
          + JzKetCoupled(
              1, -1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 1))
          )
          / 2
          + JzKetCoupled(
              2, -1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
          )
          / 2
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 1),
              )
          )
          == sqrt(3)
          * JzKetCoupled(
              0, 0, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 0))
          )
          / 3
          - sqrt(2)
          * JzKetCoupled(
              1, 0, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 1))
          )
          / 2
          + sqrt(6)
          * JzKetCoupled(
              2, 0, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
          )
          / 6
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 0),
              )
          )
          == -sqrt(2)
          * JzKetCoupled(
              1, -1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 1))
          )
          / 2
          + sqrt(2)
          * JzKetCoupled(
              2, -1, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
          )
          / 2
      )
      assert couple(
          TensorProduct(
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(1, -1),
          )
      ) == JzKetCoupled(
          2, -2, (S.Half, S.Half, 1), ((1, 2, 1), (1, 3, 2))
      )
      # j1=S.Half, j2=1, j3=1
      assert couple(
          TensorProduct(JzKet(S.Half, S.Half), JzKet(1, 1), JzKet(1, 1))
      ) == JzKetCoupled(
          Rational(5, 2),
          Rational(5, 2),
          (S.Half, 1, 1),
          ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, 1), JzKet(1, 0)
              )
          )
          == sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 5
          + sqrt(10)
          * JzKetCoupled(
              S(5) / 2,
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, 1), JzKet(1, -1)
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, S.Half)),
          )
          / 2
          + sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 5
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, 0), JzKet(1, 1)
              )
          )
          == sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          - 2
          * sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              S(5) / 2,
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, 0), JzKet(1, 0)
              )
          )
          == JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, S.Half)),
          )
          / 3
          - sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, S.Half)),
          )
          / 3
          + sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          + sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              S(5) / 2,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, 0), JzKet(1, -1)
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, S.Half)),
          )
          / 3
          + JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, S.Half)),
          )
          / 3
          + JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          + 4
          * sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, -1), JzKet(1, 1)
              )
          )
          == -2
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, S.Half)),
          )
          / 3
          + sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, S.Half)),
          )
          / 6
          + sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          - 2
          * sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, -1), JzKet(1, 0)
              )
          )
          == -sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, S.Half)),
          )
          / 3
          - JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, S.Half)),
          )
          / 3
          + 2
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          - sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, -1), JzKet(1, -1)
              )
          )
          == sqrt(6)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          + sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 1),
                  JzKet(1, 1),
              )
          )
          == -sqrt(6)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          - sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              S(5) / 2,
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 1),
                  JzKet(1, 0),
              )
          )
          == -sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, S.Half)),
          )
          / 3
          - JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, S.Half)),
          )
          / 3
          - 2
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          + sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              S(5) / 2,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 1),
                  JzKet(1, -1),
              )
          )
          == -2
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, S.Half)),
          )
          / 3
          + sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, S.Half)),
          )
          / 6
          - sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          + 2
          * sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 0),
                  JzKet(1, 1),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, S.Half)),
          )
          / 3
          + JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, S.Half)),
          )
          / 3
          - JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          - 4
          * sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              S(5) / 2,
              S.Half,
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 0),
                  JzKet(1, 0),
              )
          )
          == JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, S.Half)),
          )
          / 3
          - sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, S.Half)),
          )
          / 3
          - sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          - sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 0),
                  JzKet(1, -1),
              )
          )
          == -sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 2, S.Half), (1, 3, Rational(3, 2))),
          )
          / 3
          + 2
          * sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, -1),
                  JzKet(1, 1),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, S.Half)),
          )
          / 2
          - sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 5
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, -1),
                  JzKet(1, 0),
              )
          )
          == -sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(3, 2))),
          )
          / 5
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert couple(
          TensorProduct(
              JzKet(S.Half, Rational(-1, 2)), JzKet(1, -1), JzKet(1, -1)
          )
      ) == JzKetCoupled(
          S(5) / 2,
          Rational(-5, 2),
          (S.Half, 1, 1),
          ((1, 2, Rational(3, 2)), (1, 3, Rational(5, 2))),
      )
      #  j1=1, j2=1, j3=1
      assert couple(
          TensorProduct(JzKet(1, 1), JzKet(1, 1), JzKet(1, 1))
      ) == JzKetCoupled(3, 3, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
      assert (
          couple(TensorProduct(JzKet(1, 1), JzKet(1, 1), JzKet(1, 0)))
          == sqrt(6)
          * JzKetCoupled(2, 2, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 3
          + sqrt(3)
          * JzKetCoupled(3, 2, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 3
      )
      assert (
          couple(TensorProduct(JzKet(1, 1), JzKet(1, 1), JzKet(1, -1)))
          == sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 5
          + sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 3
          + sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, 1), JzKet(1, 0), JzKet(1, 1)))
          == sqrt(2)
          * JzKetCoupled(2, 2, (1, 1, 1), ((1, 2, 1), (1, 3, 2)))
          / 2
          - sqrt(6)
          * JzKetCoupled(2, 2, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + sqrt(3)
          * JzKetCoupled(3, 2, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 3
      )
      assert (
          couple(TensorProduct(JzKet(1, 1), JzKet(1, 0), JzKet(1, 0)))
          == JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          - sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 10
          + JzKetCoupled(2, 1, (1, 1, 1), ((1, 2, 1), (1, 3, 2))) / 2
          + sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + 2
          * sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, 1), JzKet(1, 0), JzKet(1, -1)))
          == sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 0)))
          / 6
          + JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 10
          + sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 2)))
          / 6
          + JzKetCoupled(2, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 2))) / 2
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 10
      )
      assert (
          couple(TensorProduct(JzKet(1, 1), JzKet(1, -1), JzKet(1, 1)))
          == sqrt(3)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 0), (1, 3, 1)))
          / 3
          - JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 30
          + JzKetCoupled(2, 1, (1, 1, 1), ((1, 2, 1), (1, 3, 2))) / 2
          - sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, 1), JzKet(1, -1), JzKet(1, 0)))
          == -sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 0)))
          / 6
          + sqrt(3)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 0), (1, 3, 1)))
          / 3
          - sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 15
          + sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 2)))
          / 3
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 10
      )
      assert (
          couple(TensorProduct(JzKet(1, 1), JzKet(1, -1), JzKet(1, -1)))
          == sqrt(3)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 0), (1, 3, 1)))
          / 3
          + JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 30
          + JzKetCoupled(2, -1, (1, 1, 1), ((1, 2, 1), (1, 3, 2))) / 2
          + sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, 0), JzKet(1, 1), JzKet(1, 1)))
          == -sqrt(2)
          * JzKetCoupled(2, 2, (1, 1, 1), ((1, 2, 1), (1, 3, 2)))
          / 2
          - sqrt(6)
          * JzKetCoupled(2, 2, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + sqrt(3)
          * JzKetCoupled(3, 2, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 3
      )
      assert (
          couple(TensorProduct(JzKet(1, 0), JzKet(1, 1), JzKet(1, 0)))
          == -JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          - sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 10
          - JzKetCoupled(2, 1, (1, 1, 1), ((1, 2, 1), (1, 3, 2))) / 2
          + sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + 2
          * sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, 0), JzKet(1, 1), JzKet(1, -1)))
          == -sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 0)))
          / 6
          - JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 10
          - sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 2)))
          / 6
          + JzKetCoupled(2, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 2))) / 2
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 10
      )
      assert (
          couple(TensorProduct(JzKet(1, 0), JzKet(1, 0), JzKet(1, 1)))
          == -sqrt(3)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 0), (1, 3, 1)))
          / 3
          + sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 15
          - sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 3
          + 2
          * sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, 0), JzKet(1, 0), JzKet(1, 0)))
          == -sqrt(3)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 0), (1, 3, 1)))
          / 3
          - 2
          * sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 15
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 5
      )
      assert (
          couple(TensorProduct(JzKet(1, 0), JzKet(1, 0), JzKet(1, -1)))
          == -sqrt(3)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 0), (1, 3, 1)))
          / 3
          + sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 15
          + sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 3
          + 2
          * sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, 0), JzKet(1, -1), JzKet(1, 1)))
          == sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 0)))
          / 6
          - JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 10
          + sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 2)))
          / 6
          - JzKetCoupled(2, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 2))) / 2
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 10
      )
      assert (
          couple(TensorProduct(JzKet(1, 0), JzKet(1, -1), JzKet(1, 0)))
          == -JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          - sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 10
          + JzKetCoupled(2, -1, (1, 1, 1), ((1, 2, 1), (1, 3, 2))) / 2
          - sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + 2
          * sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, 0), JzKet(1, -1), JzKet(1, -1)))
          == sqrt(2)
          * JzKetCoupled(2, -2, (1, 1, 1), ((1, 2, 1), (1, 3, 2)))
          / 2
          + sqrt(6)
          * JzKetCoupled(2, -2, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + sqrt(3)
          * JzKetCoupled(3, -2, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 3
      )
      assert (
          couple(TensorProduct(JzKet(1, -1), JzKet(1, 1), JzKet(1, 1)))
          == sqrt(3)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 0), (1, 3, 1)))
          / 3
          + JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 30
          - JzKetCoupled(2, 1, (1, 1, 1), ((1, 2, 1), (1, 3, 2))) / 2
          - sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, -1), JzKet(1, 1), JzKet(1, 0)))
          == sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 0)))
          / 6
          + sqrt(3)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 0), (1, 3, 1)))
          / 3
          - sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 15
          - sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 2)))
          / 3
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 10
      )
      assert (
          couple(TensorProduct(JzKet(1, -1), JzKet(1, 1), JzKet(1, -1)))
          == sqrt(3)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 0), (1, 3, 1)))
          / 3
          - JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 30
          - JzKetCoupled(2, -1, (1, 1, 1), ((1, 2, 1), (1, 3, 2))) / 2
          + sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, -1), JzKet(1, 0), JzKet(1, 1)))
          == -sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 0)))
          / 6
          + JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 10
          - sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 2, 1), (1, 3, 2)))
          / 6
          - JzKetCoupled(2, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 2))) / 2
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 10
      )
      assert (
          couple(TensorProduct(JzKet(1, -1), JzKet(1, 0), JzKet(1, 0)))
          == JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 1), (1, 3, 1))) / 2
          - sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 10
          - JzKetCoupled(2, -1, (1, 1, 1), ((1, 2, 1), (1, 3, 2))) / 2
          - sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + 2
          * sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, -1), JzKet(1, 0), JzKet(1, -1)))
          == -sqrt(2)
          * JzKetCoupled(2, -2, (1, 1, 1), ((1, 2, 1), (1, 3, 2)))
          / 2
          + sqrt(6)
          * JzKetCoupled(2, -2, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 6
          + sqrt(3)
          * JzKetCoupled(3, -2, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 3
      )
      assert (
          couple(TensorProduct(JzKet(1, -1), JzKet(1, -1), JzKet(1, 1)))
          == sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 1)))
          / 5
          - sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 3
          + sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 15
      )
      assert (
          couple(TensorProduct(JzKet(1, -1), JzKet(1, -1), JzKet(1, 0)))
          == -sqrt(6)
          * JzKetCoupled(2, -2, (1, 1, 1), ((1, 2, 2), (1, 3, 2)))
          / 3
          + sqrt(3)
          * JzKetCoupled(3, -2, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
          / 3
      )
      assert couple(
          TensorProduct(JzKet(1, -1), JzKet(1, -1), JzKet(1, -1))
      ) == JzKetCoupled(3, -3, (1, 1, 1), ((1, 2, 2), (1, 3, 3)))
      # j1=S.Half, j2=S.Half, j3=Rational(3, 2)
      assert couple(
          TensorProduct(
              JzKet(S.Half, S.Half),
              JzKet(S.Half, S.Half),
              JzKet(Rational(3, 2), Rational(3, 2)),
          )
      ) == JzKetCoupled(
          Rational(5, 2),
          S(5) / 2,
          (S.Half, S.Half, Rational(3, 2)),
          ((1, 2, 1), (1, 3, Rational(5, 2))),
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), S.Half),
              )
          )
          == sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 5
          + sqrt(15)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), Rational(-1, 2)),
              )
          )
          == sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 6
          + 2
          * sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), Rational(-3, 2)),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 2
          + sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 5
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              -S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), Rational(3, 2)),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 0), (1, 3, Rational(3, 2))),
          )
          / 2
          - sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 10
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), S.Half),
              )
          )
          == -sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 6
          + sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 0), (1, 3, Rational(3, 2))),
          )
          / 2
          - sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 30
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), Rational(-1, 2)),
              )
          )
          == -sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 6
          + sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 0), (1, 3, Rational(3, 2))),
          )
          / 2
          + sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 30
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              -S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), Rational(-3, 2)),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 0), (1, 3, Rational(3, 2))),
          )
          / 2
          + sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 10
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), Rational(3, 2)),
              )
          )
          == -sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 0), (1, 3, Rational(3, 2))),
          )
          / 2
          - sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 10
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), S.Half),
              )
          )
          == -sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 6
          - sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 0), (1, 3, Rational(3, 2))),
          )
          / 2
          - sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 30
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), Rational(-1, 2)),
              )
          )
          == -sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 6
          - sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 0), (1, 3, Rational(3, 2))),
          )
          / 2
          + sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 30
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              -S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), Rational(-3, 2)),
              )
          )
          == -sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 0), (1, 3, Rational(3, 2))),
          )
          / 2
          + sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 10
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), Rational(3, 2)),
              )
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 2
          - sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 5
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), S.Half),
              )
          )
          == sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, S.Half)),
          )
          / 6
          - 2
          * sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 15
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              -S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), Rational(-1, 2)),
              )
          )
          == -sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 2, 1), (1, 3, Rational(3, 2))),
          )
          / 5
          + sqrt(15)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 2, 1), (1, 3, Rational(5, 2))),
          )
          / 5
      )
      assert couple(
          TensorProduct(
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(Rational(3, 2), Rational(-3, 2)),
          )
      ) == JzKetCoupled(
          Rational(5, 2),
          -S(5) / 2,
          (S.Half, S.Half, Rational(3, 2)),
          ((1, 2, 1), (1, 3, Rational(5, 2))),
      )
      # Couple j1 to j3
      # j1=1/2, j2=1/2, j3=1/2
      assert couple(
          TensorProduct(
              JzKet(S.Half, S.Half),
              JzKet(S.Half, S.Half),
              JzKet(S.Half, S.Half),
          ),
          ((1, 3), (1, 2)),
      ) == JzKetCoupled(
          Rational(3, 2),
          S(3) / 2,
          (S.Half, S.Half, S.Half),
          ((1, 3, 1), (1, 2, Rational(3, 2))),
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, S.Half),
              ((1, 3, 0), (1, 2, S.Half)),
          )
          / 2
          - sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, S.Half),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, S.One / 2),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, S.Half),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 3
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, S.One / 2),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, S.Half),
              ((1, 3, 0), (1, 2, S.Half)),
          )
          / 2
          + sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, S.Half),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, S.One / 2),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, S.Half),
              ((1, 3, 0), (1, 2, S.Half)),
          )
          / 2
          - sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, S.Half),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, S.One / 2),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, S.Half),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 3
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, S.One / 2),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, S.Half),
              ((1, 3, 0), (1, 2, S.Half)),
          )
          / 2
          + sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, S.Half),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, S.One / 2),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 3
      )
      assert couple(
          TensorProduct(
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(S.Half, Rational(-1, 2)),
          ),
          ((1, 3), (1, 2)),
      ) == JzKetCoupled(
          Rational(3, 2),
          -S(3) / 2,
          (S.Half, S.Half, S.Half),
          ((1, 3, 1), (1, 2, Rational(3, 2))),
      )
      # j1=1/2, j2=1/2, j3=1
      assert couple(
          TensorProduct(
              JzKet(S.Half, S.Half), JzKet(S.Half, S.Half), JzKet(1, 1)
          ),
          ((1, 3), (1, 2)),
      ) == JzKetCoupled(
          2, 2, (S.Half, S.Half, 1), ((1, 3, Rational(3, 2)), (1, 2, 2))
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(1, 0),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(3)
          * JzKetCoupled(
              1, 1, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 1))
          )
          / 3
          - sqrt(6)
          * JzKetCoupled(
              1,
              1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 1)),
          )
          / 6
          + sqrt(2)
          * JzKetCoupled(
              2,
              1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 2)),
          )
          / 2
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(1, -1),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(3)
          * JzKetCoupled(
              0, 0, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 0))
          )
          / 3
          + sqrt(3)
          * JzKetCoupled(
              1, 0, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 1))
          )
          / 3
          - sqrt(6)
          * JzKetCoupled(
              1,
              0,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 1)),
          )
          / 6
          + sqrt(6)
          * JzKetCoupled(
              2,
              0,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 2)),
          )
          / 6
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 1),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(3)
          * JzKetCoupled(
              1,
              1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 1)),
          )
          / 2
          + JzKetCoupled(
              2,
              1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 2)),
          )
          / 2
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 0),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(
              0, 0, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 0))
          )
          / 6
          + sqrt(6)
          * JzKetCoupled(
              1, 0, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 1))
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              1,
              0,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 1)),
          )
          / 3
          + sqrt(3)
          * JzKetCoupled(
              2,
              0,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 2)),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, -1),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(
              1, -1, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 1))
          )
          / 3
          + sqrt(3)
          * JzKetCoupled(
              1,
              -1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 1)),
          )
          / 6
          + JzKetCoupled(
              2,
              -1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 2)),
          )
          / 2
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(1, 1),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(6)
          * JzKetCoupled(
              1, 1, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 1))
          )
          / 3
          - sqrt(3)
          * JzKetCoupled(
              1,
              1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 1)),
          )
          / 6
          + JzKetCoupled(
              2,
              1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 2)),
          )
          / 2
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(1, 0),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(
              0, 0, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 0))
          )
          / 6
          - sqrt(6)
          * JzKetCoupled(
              1, 0, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 1))
          )
          / 6
          - sqrt(3)
          * JzKetCoupled(
              1,
              0,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 1)),
          )
          / 3
          + sqrt(3)
          * JzKetCoupled(
              2,
              0,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 2)),
          )
          / 3
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(1, -1),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(3)
          * JzKetCoupled(
              1,
              -1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 1)),
          )
          / 2
          + JzKetCoupled(
              2,
              -1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 2)),
          )
          / 2
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 1),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(3)
          * JzKetCoupled(
              0, 0, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 0))
          )
          / 3
          - sqrt(3)
          * JzKetCoupled(
              1, 0, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 1))
          )
          / 3
          + sqrt(6)
          * JzKetCoupled(
              1,
              0,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 1)),
          )
          / 6
          + sqrt(6)
          * JzKetCoupled(
              2,
              0,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 2)),
          )
          / 6
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 0),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(3)
          * JzKetCoupled(
              1, -1, (S.Half, S.Half, 1), ((1, 3, S.Half), (1, 2, 1))
          )
          / 3
          + sqrt(6)
          * JzKetCoupled(
              1,
              -1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 1)),
          )
          / 6
          + sqrt(2)
          * JzKetCoupled(
              2,
              -1,
              (S.Half, S.Half, 1),
              ((1, 3, Rational(3, 2)), (1, 2, 2)),
          )
          / 2
      )
      assert couple(
          TensorProduct(
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(1, -1),
          ),
          ((1, 3), (1, 2)),
      ) == JzKetCoupled(
          2,
          -2,
          (S.Half, S.Half, 1),
          ((1, 3, Rational(3, 2)), (1, 2, 2)),
      )
      # j 1=1/2, j 2=1, j 3=1
      assert couple(
          TensorProduct(
              JzKet(S.Half, S.Half), JzKet(1, 1), JzKet(1, 1)
          ),
          ((1, 3), (1, 2)),
      ) == JzKetCoupled(
          Rational(5, 2),
          Rational(5, 2),
          (S.Half, 1, 1),
          ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, 1), JzKet(1, 0)
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          - 2
          * sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              S(5) / 2,
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, 1), JzKet(1, -1)
              ),
              ((1, 3), (1, 2)),
          )
          == -2
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, S.Half)),
          )
          / 3
          + sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, S.Half)),
          )
          / 6
          + sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          - 2
          * sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, 0), JzKet(1, 1)
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 5
          + sqrt(10)
          * JzKetCoupled(
              S(5) / 2,
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, 0), JzKet(1, 0)
              ),
              ((1, 3), (1, 2)),
          )
          == JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, S.Half)),
          )
          / 3
          - sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, S.Half)),
          )
          / 3
          + sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          + sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              S(5) / 2,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, 0), JzKet(1, -1)
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, S.Half)),
          )
          / 3
          - JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, S.Half)),
          )
          / 3
          + 2
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          - sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, -1), JzKet(1, 1)
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, S.Half)),
          )
          / 2
          + sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 5
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, -1), JzKet(1, 0)
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, S.Half)),
          )
          / 3
          + JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, S.Half)),
          )
          / 3
          + JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          + 4
          * sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half), JzKet(1, -1), JzKet(1, -1)
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          + sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 1),
                  JzKet(1, 1),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(6)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          - sqrt(30)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              S(5) / 2,
              Rational(3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 1),
                  JzKet(1, 0),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, S.Half)),
          )
          / 3
          + JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, S.Half)),
          )
          / 3
          - JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          - 4
          * sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              S(5) / 2,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 1),
                  JzKet(1, -1),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, S.Half)),
          )
          / 2
          - sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 5
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 0),
                  JzKet(1, 1),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, S.Half)),
          )
          / 3
          - JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, S.Half)),
          )
          / 3
          - 2
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          + sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(5)
          * JzKetCoupled(
              S(5) / 2,
              S.Half,
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 0),
                  JzKet(1, 0),
              ),
              ((1, 3), (1, 2)),
          )
          == JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, S.Half)),
          )
          / 3
          - sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, S.Half)),
          )
          / 3
          - sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          - sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, 0),
                  JzKet(1, -1),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 5
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, -1),
                  JzKet(1, 1),
              ),
              ((1, 3), (1, 2)),
          )
          == -2
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, S.Half)),
          )
          / 3
          + sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, S.Half)),
          )
          / 6
          - sqrt(2)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          + 2
          * sqrt(10)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-1, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(1, -1),
                  JzKet(1, 0),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 3, S.Half), (1, 2, Rational(3, 2))),
          )
          / 3
          + 2
          * sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(3, 2))),
          )
          / 15
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, 1, 1),
              ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert couple(
          TensorProduct(
              JzKet(S.Half, Rational(-1, 2)), JzKet(1, -1), JzKet(1, -1)
          ),
          ((1, 3), (1, 2)),
      ) == JzKetCoupled(
          S(5) / 2,
          Rational(-5, 2),
          (S.Half, 1, 1),
          ((1, 3, Rational(3, 2)), (1, 2, Rational(5, 2))),
      )
      # j1=1, 1, 1
      assert couple(
          TensorProduct(JzKet(1, 1), JzKet(1, 1), JzKet(1, 1)),
          ((1, 3), (1, 2)),
      ) == JzKetCoupled(3, 3, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
      assert (
          couple(
              TensorProduct(JzKet(1, 1), JzKet(1, 1), JzKet(1, 0)),
              ((1, 3), (1, 2)),
          )
          == sqrt(2)
          * JzKetCoupled(2, 2, (1, 1, 1), ((1, 3, 1), (1, 2, 2)))
          / 2
          - sqrt(6)
          * JzKetCoupled(2, 2, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + sqrt(3)
          * JzKetCoupled(3, 2, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 3
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 1), JzKet(1, 1), JzKet(1, -1)),
              ((1, 3), (1, 2)),
          )
          == sqrt(3)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 0), (1, 2, 1)))
          / 3
          - JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 30
          + JzKetCoupled(2, 1, (1, 1, 1), ((1, 3, 1), (1, 2, 2))) / 2
          - sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 1), JzKet(1, 0), JzKet(1, 1)),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(2, 2, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 3
          + sqrt(3)
          * JzKetCoupled(3, 2, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 3
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 1), JzKet(1, 0), JzKet(1, 0)),
              ((1, 3), (1, 2)),
          )
          == JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          - sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 10
          + JzKetCoupled(2, 1, (1, 1, 1), ((1, 3, 1), (1, 2, 2))) / 2
          + sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + 2
          * sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 1), JzKet(1, 0), JzKet(1, -1)),
              ((1, 3), (1, 2)),
          )
          == -sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 0)))
          / 6
          + sqrt(3)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 0), (1, 2, 1)))
          / 3
          - sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 15
          + sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 2)))
          / 3
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 10
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 1), JzKet(1, -1), JzKet(1, 1)),
              ((1, 3), (1, 2)),
          )
          == sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 5
          + sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 3
          + sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 1), JzKet(1, -1), JzKet(1, 0)),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 0)))
          / 6
          + JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 10
          + sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 2)))
          / 6
          + JzKetCoupled(2, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 2))) / 2
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 10
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 1), JzKet(1, -1), JzKet(1, -1)),
              ((1, 3), (1, 2)),
          )
          == sqrt(3)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 0), (1, 2, 1)))
          / 3
          + JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 30
          + JzKetCoupled(2, -1, (1, 1, 1), ((1, 3, 1), (1, 2, 2))) / 2
          + sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 0), JzKet(1, 1), JzKet(1, 1)),
              ((1, 3), (1, 2)),
          )
          == -sqrt(2)
          * JzKetCoupled(2, 2, (1, 1, 1), ((1, 3, 1), (1, 2, 2)))
          / 2
          - sqrt(6)
          * JzKetCoupled(2, 2, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + sqrt(3)
          * JzKetCoupled(3, 2, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 3
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 0), JzKet(1, 1), JzKet(1, 0)),
              ((1, 3), (1, 2)),
          )
          == -sqrt(3)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 0), (1, 2, 1)))
          / 3
          + sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 15
          - sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 3
          + 2
          * sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 0), JzKet(1, 1), JzKet(1, -1)),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 0)))
          / 6
          - JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 10
          + sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 2)))
          / 6
          - JzKetCoupled(2, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 2))) / 2
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 10
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 0), JzKet(1, 0), JzKet(1, 1)),
              ((1, 3), (1, 2)),
          )
          == -JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          - sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 10
          - JzKetCoupled(2, 1, (1, 1, 1), ((1, 3, 1), (1, 2, 2))) / 2
          + sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + 2
          * sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 0), JzKet(1, 0), JzKet(1, 0)),
              ((1, 3), (1, 2)),
          )
          == -sqrt(3)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 0), (1, 2, 1)))
          / 3
          - 2
          * sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 15
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 5
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 0), JzKet(1, 0), JzKet(1, -1)),
              ((1, 3), (1, 2)),
          )
          == -JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          - sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 10
          + JzKetCoupled(2, -1, (1, 1, 1), ((1, 3, 1), (1, 2, 2))) / 2
          - sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + 2
          * sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 0), JzKet(1, -1), JzKet(1, 1)),
              ((1, 3), (1, 2)),
          )
          == -sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 0)))
          / 6
          - JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 10
          - sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 2)))
          / 6
          + JzKetCoupled(2, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 2))) / 2
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 10
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 0), JzKet(1, -1), JzKet(1, 0)),
              ((1, 3), (1, 2)),
          )
          == -sqrt(3)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 0), (1, 2, 1)))
          / 3
          + sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 15
          + sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 3
          + 2
          * sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, 0), JzKet(1, -1), JzKet(1, -1)),
              ((1, 3), (1, 2)),
          )
          == sqrt(2)
          * JzKetCoupled(2, -2, (1, 1, 1), ((1, 3, 1), (1, 2, 2)))
          / 2
          + sqrt(6)
          * JzKetCoupled(2, -2, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + sqrt(3)
          * JzKetCoupled(3, -2, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 3
      )
      assert (
          couple(
              TensorProduct(JzKet(1, -1), JzKet(1, 1), JzKet(1, 1)),
              ((1, 3), (1, 2)),
          )
          == sqrt(3)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 0), (1, 2, 1)))
          / 3
          + JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 30
          - JzKetCoupled(2, 1, (1, 1, 1), ((1, 3, 1), (1, 2, 2))) / 2
          - sqrt(3)
          * JzKetCoupled(2, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + sqrt(15)
          * JzKetCoupled(3, 1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, -1), JzKet(1, 1), JzKet(1, 0)),
              ((1, 3), (1, 2)),
          )
          == -sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 0)))
          / 6
          + JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 10
          - sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 2)))
          / 6
          - JzKetCoupled(2, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 2))) / 2
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 10
      )
      assert (
          couple(
              TensorProduct(JzKet(1, -1), JzKet(1, 1), JzKet(1, -1)),
              ((1, 3), (1, 2)),
          )
          == sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 5
          - sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 3
          + sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, -1), JzKet(1, 0), JzKet(1, 1)),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(0, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 0)))
          / 6
          + sqrt(3)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 0), (1, 2, 1)))
          / 3
          - sqrt(15)
          * JzKetCoupled(1, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 15
          - sqrt(3)
          * JzKetCoupled(2, 0, (1, 1, 1), ((1, 3, 1), (1, 2, 2)))
          / 3
          + sqrt(10)
          * JzKetCoupled(3, 0, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 10
      )
      assert (
          couple(
              TensorProduct(JzKet(1, -1), JzKet(1, 0), JzKet(1, 0)),
              ((1, 3), (1, 2)),
          )
          == JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          - sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 10
          - JzKetCoupled(2, -1, (1, 1, 1), ((1, 3, 1), (1, 2, 2))) / 2
          - sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + 2
          * sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, -1), JzKet(1, 0), JzKet(1, -1)),
              ((1, 3), (1, 2)),
          )
          == -sqrt(6)
          * JzKetCoupled(2, -2, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 3
          + sqrt(3)
          * JzKetCoupled(3, -2, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 3
      )
      assert (
          couple(
              TensorProduct(JzKet(1, -1), JzKet(1, -1), JzKet(1, 1)),
              ((1, 3), (1, 2)),
          )
          == sqrt(3)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 0), (1, 2, 1)))
          / 3
          - JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 1), (1, 2, 1))) / 2
          + sqrt(15)
          * JzKetCoupled(1, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 1)))
          / 30
          - JzKetCoupled(2, -1, (1, 1, 1), ((1, 3, 1), (1, 2, 2))) / 2
          + sqrt(3)
          * JzKetCoupled(2, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + sqrt(15)
          * JzKetCoupled(3, -1, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 15
      )
      assert (
          couple(
              TensorProduct(JzKet(1, -1), JzKet(1, -1), JzKet(1, 0)),
              ((1, 3), (1, 2)),
          )
          == -sqrt(2)
          * JzKetCoupled(2, -2, (1, 1, 1), ((1, 3, 1), (1, 2, 2)))
          / 2
          + sqrt(6)
          * JzKetCoupled(2, -2, (1, 1, 1), ((1, 3, 2), (1, 2, 2)))
          / 6
          + sqrt(3)
          * JzKetCoupled(3, -2, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
          / 3
      )
      assert couple(
          TensorProduct(JzKet(1, -1), JzKet(1, -1), JzKet(1, -1)),
          ((1, 3), (1, 2)),
      ) == JzKetCoupled(3, -3, (1, 1, 1), ((1, 3, 2), (1, 2, 3)))
      # j1=1/2, j2=1/2, j3=3/2
      assert couple(
          TensorProduct(
              JzKet(S.Half, S.Half),
              JzKet(S.Half, S.Half),
              JzKet(Rational(3, 2), Rational(3, 2)),
          ),
          ((1, 3), (1, 2)),
      ) == JzKetCoupled(
          Rational(5, 2),
          S(5) / 2,
          (S.Half, S.Half, Rational(3, 2)),
          ((1, 3, 2), (1, 2, Rational(5, 2))),
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), S.Half),
              ),
              ((1, 3), (1, 2)),
          )
          == JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 2
          - sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 10
          + sqrt(15)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), Rational(-1, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 3
          - sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 5
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), Rational(-3, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(2)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 2
          + JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 2
          - sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 10
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              -S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), Rational(3, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == 2
          * sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 5
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), S.Half),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 6
          + 3
          * sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 10
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), Rational(-1, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 6
          + sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 3
          + sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 5
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              -S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, S.Half),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), Rational(-3, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 2
          + sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 10
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), Rational(3, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 2
          - sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 10
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), S.Half),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 6
          - sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 3
          - sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 5
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), Rational(-1, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 6
          - sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 6
          - 3
          * sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 10
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              -S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, S.Half),
                  JzKet(Rational(3, 2), Rational(-3, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == -2
          * sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 5
          + sqrt(5)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), Rational(3, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(2)
          * JzKetCoupled(
              S.Half,
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 2
          - JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 2
          + sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              S.Half,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 10
          + sqrt(10)
          * JzKetCoupled(
              Rational(5, 2),
              S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), S.Half),
              ),
              ((1, 3), (1, 2)),
          )
          == -sqrt(6)
          * JzKetCoupled(
              S.Half,
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, S.Half)),
          )
          / 6
          - sqrt(3)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 3
          + sqrt(5)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-1, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 5
          + sqrt(30)
          * JzKetCoupled(
              Rational(5, 2),
              -S(1) / 2,
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 10
      )
      assert (
          couple(
              TensorProduct(
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(S.Half, Rational(-1, 2)),
                  JzKet(Rational(3, 2), Rational(-1, 2)),
              ),
              ((1, 3), (1, 2)),
          )
          == -JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 1), (1, 2, Rational(3, 2))),
          )
          / 2
          + sqrt(15)
          * JzKetCoupled(
              Rational(3, 2),
              Rational(-3, 2),
              (S.Half, S.Half, Rational(3, 2)),
              ((1, 3, 2), (1, 2, Rational(3, 2))),
          )
          / 10
          + sqrt(15)
          * JzKetCoupled(
              Rational(5, 2),
              Rational(-3, 2),
              (S.Half, S.Half, S(3) / 2),
              ((1, 3, 2), (1, 2, Rational(5, 2))),
          )
          / 5
      )
      assert couple(
          TensorProduct(
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(S.Half, Rational(-1, 2)),
              JzKet(Rational(3, 2), Rational(-3, 2)),
          ),
          ((1, 3), (1, 2)),
      ) == JzKetCoupled(
          Rational(5, 2),
          -S(5) / 2,
          (S.Half, S.Half, Rational(3, 2)),
          ((1, 3, 2), (1, 2, Rational(5, 2))),
      )
  
imports: from sympy.concrete.summations import Sum
  from sympy.core.function import expand
  from sympy.core.numbers import I, Rational, pi
  from sympy.core.singleton import S
  from sympy.core.symbol import symbols
  from sympy.functions.elementary.exponential import exp
  from sympy.functions.elementary.miscellaneous import sqrt
  from sympy.functions.elementary.trigonometric import cos, sin
  from sympy.matrices.dense import Matrix
  from sympy.abc import alpha, beta, gamma, j, m
  from sympy.physics.quantum import hbar, represent, Commutator, InnerProduct
  from sympy.physics.quantum.qapply import qapply
  from sympy.physics.quantum.tensorproduct import TensorProduct
  from sympy.physics.quantum.cg import CG
  from sympy.physics.quantum.spin import Jx, Jy, Jz, Jplus, Jminus, J2, JxBra, JyBra, JzBra, JxKet, JyKet, JzKet, JxKetCoupled, JyKetCoupled, JzKetCoupled, couple, uncouple, Rotation, WignerD
  from sympy.testing.pytest import raises, slow
label_desc: Use image node_id 23 for calling a global function with example usage: test_couple_3_states_numerical() without return types