Repository: fujiisoup/py3nj

Issue:
Unable to compute some values
```python
py3nj.wigner3j(4, 0, 0, 0, 0, 0)
```
gives an error
```python
IndexError: index 4 is out of bounds for axis 0 with size 1
```


Discussion:

Patch:
diff --git a/doc/conf.py b/doc/conf.py
index e2b803b..bd93ae3 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -29,25 +29,26 @@ print("py3nj: %s, %s" % (py3nj.__version__, py3nj.__file__))
 # -- General configuration ------------------------------------------------
 
 # If your documentation needs a minimal Sphinx version, state it here.
-#needs_sphinx = '1.0'
+# needs_sphinx = '1.0'
 
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
 extensions = [
-    'sphinx.ext.autodoc',
-    'sphinx.ext.autosummary',
-    'sphinx.ext.intersphinx',
-    'sphinx.ext.extlinks',
-    'sphinx.ext.mathjax',
-    'numpydoc',
-    'IPython.sphinxext.ipython_directive',
-    'IPython.sphinxext.ipython_console_highlighting',
+    "sphinx.ext.autodoc",
+    "sphinx.ext.autosummary",
+    "sphinx.ext.intersphinx",
+    "sphinx.ext.extlinks",
+    "sphinx.ext.mathjax",
+    "numpydoc",
+    "IPython.sphinxext.ipython_directive",
+    "IPython.sphinxext.ipython_console_highlighting",
 ]
 
-extlinks = {'issue': ('https://github.com/fujiisoup/py3nj/issues/%s', 'GH'),
-            'pull': ('https://github.com/fujiisoup/py3nj/pull/%s', 'PR'),
-            }
+extlinks = {
+    "issue": ("https://github.com/fujiisoup/py3nj/issues/%s", "GH"),
+    "pull": ("https://github.com/fujiisoup/py3nj/pull/%s", "PR"),
+}
 
 autosummary_generate = True
 
@@ -55,20 +56,20 @@ numpydoc_class_members_toctree = True
 numpydoc_show_class_members = False
 
 # Add any paths that contain templates here, relative to this directory.
-templates_path = ['_templates']
+templates_path = ["_templates"]
 
 # The suffix of source filenames.
-source_suffix = '.rst'
+source_suffix = ".rst"
 
 # The encoding of source files.
-#source_encoding = 'utf-8-sig'
+# source_encoding = 'utf-8-sig'
 
 # The master toctree document.
-master_doc = 'index'
+master_doc = "index"
 
 # General information about the project.
-project = 'py3nj'
-copyright = '2018-%s, py3nj Developers' % datetime.datetime.now().year
+project = "py3nj"
+copyright = "2018-%s, py3nj Developers" % datetime.datetime.now().year
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
@@ -81,41 +82,41 @@ release = py3nj.__version__
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
-#language = None
+# language = None
 
 # There are two options for replacing |today|: either, you set today to some
 # non-false value, then it is used:
-#today = ''
+# today = ''
 # Else, today_fmt is used as the format for a strftime call.
-today_fmt = '%Y-%m-%d'
+today_fmt = "%Y-%m-%d"
 
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
-exclude_patterns = ['_build']
+exclude_patterns = ["_build"]
 
 # The reST default role (used for this markup: `text`) to use for all
 # documents.
-#default_role = None
+# default_role = None
 
 # If true, '()' will be appended to :func: etc. cross-reference text.
-#add_function_parentheses = True
+# add_function_parentheses = True
 
 # If true, the current module name will be prepended to all description
 # unit titles (such as .. function::).
-#add_module_names = True
+# add_module_names = True
 
 # If true, sectionauthor and moduleauthor directives will be shown in the
 # output. They are ignored by default.
-#show_authors = False
+# show_authors = False
 
 # The name of the Pygments (syntax highlighting) style to use.
-pygments_style = 'sphinx'
+pygments_style = "sphinx"
 
 # A list of ignored prefixes for module index sorting.
-#modindex_common_prefix = []
+# modindex_common_prefix = []
 
 # If true, keep warnings as "system message" paragraphs in the built documents.
-#keep_warnings = False
+# keep_warnings = False
 
 
 # -- Options for HTML output ----------------------------------------------
@@ -125,11 +126,12 @@ pygments_style = 'sphinx'
 
 # on_rtd is whether we are on readthedocs.org, this line of code grabbed from
 # docs.readthedocs.org
-on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
+on_rtd = os.environ.get("READTHEDOCS", None) == "True"
 
 if not on_rtd:  # only import and set the theme if we're building docs locally
     import sphinx_rtd_theme
-    html_theme = 'sphinx_rtd_theme'
+
+    html_theme = "sphinx_rtd_theme"
     html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 
 # otherwise, readthedocs.org uses their theme by default, so no need to specify it
@@ -137,44 +139,45 @@ if not on_rtd:  # only import and set the theme if we're building docs locally
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
 # documentation.
-#html_theme_options = {}
+# html_theme_options = {}
 
 # Add any paths that contain custom themes here, relative to this directory.
-#html_theme_path = []
+# html_theme_path = []
 
 # The name for this set of Sphinx documents.  If None, it defaults to
 # "<project> v<release> documentation".
-#html_title = None
+# html_title = None
 
 # A shorter title for the navigation bar.  Default is the same as html_title.
-#html_short_title = None
+# html_short_title = None
 
 # The name of an image file (relative to this directory) to place at the top
 # of the sidebar.
-#html_logo = None
+# html_logo = None
 
 # The name of an image file (within the static path) to use as favicon of the
 # docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
 # pixels large.
-#html_favicon = None
+# html_favicon = None
 
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
+html_static_path = ["_static"]
 
 # Sometimes the savefig directory doesn't exist and needs to be created
 # https://github.com/ipython/ipython/issues/8733
 # becomes obsolete when we can pin ipython>=5.2; see doc/environment.yml
-ipython_savefig_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)),
-                                   '_build','html','_static')
+ipython_savefig_dir = os.path.join(
+    os.path.dirname(os.path.abspath(__file__)), "_build", "html", "_static"
+)
 if not os.path.exists(ipython_savefig_dir):
     os.makedirs(ipython_savefig_dir)
 
 # Add any extra paths that contain custom files (such as robots.txt or
 # .htaccess) here, relative to this directory. These files are copied
 # directly to the root of the documentation.
-#html_extra_path = []
+# html_extra_path = []
 
 # If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
 # using the given strftime format.
@@ -182,98 +185,92 @@ html_last_updated_fmt = today_fmt
 
 # If true, SmartyPants will be used to convert quotes and dashes to
 # typographically correct entities.
-#html_use_smartypants = True
+# html_use_smartypants = True
 
 # Custom sidebar templates, maps document names to template names.
-#html_sidebars = {}
+# html_sidebars = {}
 
 # Additional templates that should be rendered to pages, maps page names to
 # template names.
-#html_additional_pages = {}
+# html_additional_pages = {}
 
 # If false, no module index is generated.
-#html_domain_indices = True
+# html_domain_indices = True
 
 # If false, no index is generated.
 # html_use_index = True
 
 # If true, the index is split into individual pages for each letter.
-#html_split_index = False
+# html_split_index = False
 
 # If true, links to the reST sources are added to the pages.
-#html_show_sourcelink = True
+# html_show_sourcelink = True
 
 # If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
-#html_show_sphinx = True
+# html_show_sphinx = True
 
 # If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
-#html_show_copyright = True
+# html_show_copyright = True
 
 # If true, an OpenSearch description file will be output, and all pages will
 # contain a <link> tag referring to it.  The value of this option must be the
 # base URL from which the finished HTML is served.
-#html_use_opensearch = ''
+# html_use_opensearch = ''
 
 # This is the file name suffix for HTML files (e.g. ".xhtml").
-#html_file_suffix = None
+# html_file_suffix = None
 
 # Output file base name for HTML help builder.
-htmlhelp_basename = 'py3njdoc'
+htmlhelp_basename = "py3njdoc"
 
 
 # -- Options for LaTeX output ---------------------------------------------
 
 latex_elements = {
-# The paper size ('letterpaper' or 'a4paper').
-#'papersize': 'letterpaper',
-
-# The font size ('10pt', '11pt' or '12pt').
-#'pointsize': '10pt',
-
-# Additional stuff for the LaTeX preamble.
-#'preamble': '',
+    # The paper size ('letterpaper' or 'a4paper').
+    #'papersize': 'letterpaper',
+    # The font size ('10pt', '11pt' or '12pt').
+    #'pointsize': '10pt',
+    # Additional stuff for the LaTeX preamble.
+    #'preamble': '',
 }
 
 # Grouping the document tree into LaTeX files. List of tuples
 # (source start file, target name, title,
 #  author, documentclass [howto, manual, or own class]).
 latex_documents = [
-  ('index', 'py3nj.tex', 'py3nj Documentation',
-   'xr-scipy Developers', 'manual'),
+    ("index", "py3nj.tex", "py3nj Documentation", "xr-scipy Developers", "manual"),
 ]
 
 # The name of an image file (relative to this directory) to place at the top of
 # the title page.
-#latex_logo = None
+# latex_logo = None
 
 # For "manual" documents, if this is true, then toplevel headings are parts,
 # not chapters.
-#latex_use_parts = False
+# latex_use_parts = False
 
 # If true, show page references after internal links.
-#latex_show_pagerefs = False
+# latex_show_pagerefs = False
 
 # If true, show URL addresses after external links.
-#latex_show_urls = False
+# latex_show_urls = False
 
 # Documents to append as an appendix to all manuals.
-#latex_appendices = []
+# latex_appendices = []
 
 # If false, no module index is generated.
-#latex_domain_indices = True
+# latex_domain_indices = True
 
 
 # -- Options for manual page output ---------------------------------------
 
 # One entry per manual page. List of tuples
 # (source start file, name, description, authors, manual section).
-man_pages = [
-    ('index', 'py3nj', 'py3nj Documentation',
-     ['xr-scipy Developers'], 1)
-]
+man_pages = [("index", "py3nj", "py3nj Documentation", ["xr-scipy Developers"], 1)]
 
 # If true, show URL addresses after external links.
-#man_show_urls = False
+# man_show_urls = False
 
 
 # -- Options for Texinfo output -------------------------------------------
@@ -282,27 +279,33 @@ man_pages = [
 # (source start file, target name, title, author,
 #  dir menu entry, description, category)
 texinfo_documents = [
-  ('index', 'py3nj', 'py3nj Documentation',
-   'py3nj Developers', 'py3nj', 'Wigner symbols for python',
-   'Miscellaneous'),
+    (
+        "index",
+        "py3nj",
+        "py3nj Documentation",
+        "py3nj Developers",
+        "py3nj",
+        "Wigner symbols for python",
+        "Miscellaneous",
+    ),
 ]
 
 # Documents to append as an appendix to all manuals.
-#texinfo_appendices = []
+# texinfo_appendices = []
 
 # If false, no module index is generated.
-#texinfo_domain_indices = True
+# texinfo_domain_indices = True
 
 # How to display URL addresses: 'footnote', 'no', or 'inline'.
-#texinfo_show_urls = 'footnote'
+# texinfo_show_urls = 'footnote'
 
 # If true, do not generate a @detailmenu in the "Top" node's menu.
-#texinfo_no_detailmenu = False
+# texinfo_no_detailmenu = False
 
 
 # Example configuration for intersphinx: refer to the Python standard library.
 intersphinx_mapping = {
-    'python': ('https://docs.python.org/3/', None),
-    'pandas': ('https://pandas.pydata.org/pandas-docs/stable/', None),
-    'xarray': ('http://xarray.pydata.org/', None),
+    "python": ("https://docs.python.org/3/", None),
+    "pandas": ("https://pandas.pydata.org/pandas-docs/stable/", None),
+    "xarray": ("http://xarray.pydata.org/", None),
 }
diff --git a/py3nj/__init__.py b/py3nj/__init__.py
index b630a1d..f29907c 100644
--- a/py3nj/__init__.py
+++ b/py3nj/__init__.py
@@ -1,11 +1,12 @@
 # https://github.com/fujiisoup/py3nj/issues/6
 import os
 import sys
-extra_dll_dir = os.path.join(os.path.dirname(__file__), '.libs')
 
-if sys.platform == 'win32' and os.path.isdir(extra_dll_dir):
-    os.environ.setdefault('PATH', '')
-    os.environ['PATH'] += os.pathsep + extra_dll_dir
+extra_dll_dir = os.path.join(os.path.dirname(__file__), ".libs")
+
+if sys.platform == "win32" and os.path.isdir(extra_dll_dir):
+    os.environ.setdefault("PATH", "")
+    os.environ["PATH"] += os.pathsep + extra_dll_dir
 
 from ._version import __version__
 
diff --git a/py3nj/_version.py b/py3nj/_version.py
index df9144c..b3f4756 100644
--- a/py3nj/_version.py
+++ b/py3nj/_version.py
@@ -1,1 +1,1 @@
-__version__ = '0.1.1'
+__version__ = "0.1.2"
diff --git a/py3nj/wigner.py b/py3nj/wigner.py
index 01fd012..fa2f3e5 100644
--- a/py3nj/wigner.py
+++ b/py3nj/wigner.py
@@ -4,8 +4,8 @@ from py3nj import _wigner
 
 def int_broadcast(*args):
     args = np.broadcast_arrays(*args)
-    if any(a.dtype.kind not in 'iu' for a in args):
-        raise ValueError('Argument should be integer.')
+    if any(a.dtype.kind not in "iu" for a in args):
+        raise ValueError("Argument should be integer.")
     return args
 
 
@@ -31,11 +31,15 @@ def clebsch_gordan(two_j1, two_j2, two_j3, two_m1, two_m2, two_m3):
         arguments.
     """
     two_j1, two_j2, two_j3, two_m1, two_m2, two_m3 = int_broadcast(
-        two_j1, two_j2, two_j3, two_m1, two_m2, two_m3)
+        two_j1, two_j2, two_j3, two_m1, two_m2, two_m3
+    )
 
     phase = (two_j1 - two_j2 + two_m3) % 4 - 1
-    return -phase * np.sqrt(two_j3 + 1) * wigner3j(two_j1, two_j2, two_j3,
-                                                  two_m1, two_m2, -two_m3)
+    return (
+        -phase
+        * np.sqrt(two_j3 + 1)
+        * wigner3j(two_j1, two_j2, two_j3, two_m1, two_m2, -two_m3)
+    )
 
 
 def wigner3j(two_l1, two_l2, two_l3, two_m1, two_m2, two_m3):
@@ -61,10 +65,11 @@ def wigner3j(two_l1, two_l2, two_l3, two_m1, two_m2, two_m3):
         The value of 3J symbol with the same shape of the arguments.
     """
     two_l1, two_l2, two_l3, two_m1, two_m2, two_m3 = int_broadcast(
-        two_l1, two_l2, two_l3, two_m1, two_m2, two_m3)
+        two_l1, two_l2, two_l3, two_m1, two_m2, two_m3
+    )
 
     if (two_l1 < 0).any():
-        raise ValueError('Some of l values are negative')
+        raise ValueError("Some of l values are negative")
 
     l, thrcof = drc3jj(two_l2, two_l3, two_m2, two_m3)
 
@@ -74,8 +79,11 @@ def wigner3j(two_l1, two_l2, two_l3, two_m1, two_m2, two_m3):
     two_l1 = two_l1.ravel()
     valid = (two_l1 < l1max) * ((two_m1 + two_m2 + two_m3).ravel() == 0)
 
-    if shape == ():  # scale case
+    # temporary set the invalid l1 to zero
+    two_l1 = np.where(valid, two_l1, 0)
+    if shape == ():  # scalar case
         return np.where(valid, thrcof[two_l1], 0.0)[0]
+    thrcof = thrcof.reshape(-1, l1max)
     thrcof = np.where(valid, thrcof[np.arange(len(two_l1)), two_l1], 0.0)
     return thrcof.reshape(shape)
 
@@ -103,10 +111,11 @@ def wigner6j(two_l1, two_l2, two_l3, two_l4, two_l5, two_l6):
         The value of 6J symbol with the same shape of the arguments.
     """
     two_l1, two_l2, two_l3, two_l4, two_l5, two_l6 = int_broadcast(
-        two_l1, two_l2, two_l3, two_l4, two_l5, two_l6)
+        two_l1, two_l2, two_l3, two_l4, two_l5, two_l6
+    )
 
     if (two_l1 < 0).any():
-        raise ValueError('Some of l values are negative')
+        raise ValueError("Some of l values are negative")
 
     l, sixcof = drc6j(two_l2, two_l3, two_l4, two_l5, two_l6)
 
@@ -116,14 +125,16 @@ def wigner6j(two_l1, two_l2, two_l3, two_l4, two_l5, two_l6):
     two_l1 = two_l1.ravel()
     valid = two_l1 < l1max
 
+    # temporary set the invalid l1 to the first value
+    two_l1 = np.where(valid, two_l1, 0)
     if shape == ():  # scale case
         return np.where(valid, sixcof[two_l1], 0.0)[0]
+    sixcof = sixcof.reshape(-1, l1max)
     sixcof = np.where(valid, sixcof[np.arange(len(two_l1)), two_l1], 0.0)
     return sixcof.reshape(shape)
 
 
-def wigner9j(two_l1, two_l2, two_l3, two_l4, two_l5, two_l6,
-             two_l7, two_l8, two_l9):
+def wigner9j(two_l1, two_l2, two_l3, two_l4, two_l5, two_l6, two_l7, two_l8, two_l9):
     """
     Calculate wigner 9j symbol
     (L1 L2 L3)
@@ -157,19 +168,32 @@ def wigner9j(two_l1, two_l2, two_l3, two_l4, two_l5, two_l6,
     threej: array
         The value of 9J symbol with the same shape of the arguments.
     """
-    (two_la, two_lb, two_lc, two_ld, two_le, two_lf,
-     two_lg, two_lh, two_lj) = int_broadcast(
-        two_l1, two_l2, two_l3, two_l4, two_l5, two_l6, two_l7, two_l8, two_l9)
+    (
+        two_la,
+        two_lb,
+        two_lc,
+        two_ld,
+        two_le,
+        two_lf,
+        two_lg,
+        two_lh,
+        two_lj,
+    ) = int_broadcast(
+        two_l1, two_l2, two_l3, two_l4, two_l5, two_l6, two_l7, two_l8, two_l9
+    )
 
     # calculat summation directory from 6-j symbol
-    x, sixj = drc6j(np.stack([two_lf, two_lb, two_la]),
-                    np.stack([two_lb, two_lf, two_lj]),
-                    np.stack([two_lc, two_le, two_lg]),
-                    np.stack([two_la, two_ld, two_lh]),
-                    np.stack([two_lj, two_lh, two_ld]))
+    x, sixj = drc6j(
+        np.stack([two_lf, two_lb, two_la]),
+        np.stack([two_lb, two_lf, two_lj]),
+        np.stack([two_lc, two_le, two_lg]),
+        np.stack([two_la, two_ld, two_lh]),
+        np.stack([two_lj, two_lh, two_ld]),
+    )
     neg_phase = (2 * x) % 4 - 1
     return np.sum(-neg_phase * (x + 1) * sixj[0] * sixj[1] * sixj[2], axis=-1)
 
+
 def drc3jj(two_l2, two_l3, two_m2, two_m3):
     """
     Calculate Wigner's 3j symbol
@@ -193,42 +217,44 @@ def drc3jj(two_l2, two_l3, two_m2, two_m3):
     threej: array, shape (..., n)
         The value of 3J symbol
     """
-    two_l2, two_l3, two_m2, two_m3 = int_broadcast(
-        two_l2, two_l3, two_m2, two_m3)
+    two_l2, two_l3, two_m2, two_m3 = int_broadcast(two_l2, two_l3, two_m2, two_m3)
 
     if (two_l2 < 0).any() or (two_l3 < 0).any():
-        raise ValueError('Some of l values are negative.')
+        raise ValueError("Some of l values are negative.")
     if (two_l2 < np.abs(two_m2)).any() or (two_l3 < np.abs(two_m3)).any():
-        raise ValueError('Some of m values are larger than l.')
+        raise ValueError("Some of m values are larger than l.")
 
     shape = two_l2.shape
     l1max = int(np.max(two_l2 + two_l3) + 1)
 
     thrcof, ier = _wigner.drc3jj_vec(
-        two_l2=two_l2.ravel(), two_l3=two_l3.ravel(),
-        two_m2=two_m2.ravel(), two_m3=two_m3.ravel(),
-        nvec=two_l2.size, ndim=l1max)
-    return np.arange(l1max), thrcof.reshape(shape + (l1max, ))
+        two_l2=two_l2.ravel(),
+        two_l3=two_l3.ravel(),
+        two_m2=two_m2.ravel(),
+        two_m3=two_m3.ravel(),
+        nvec=two_l2.size,
+        ndim=l1max,
+    )
+    return np.arange(l1max), thrcof.reshape(shape + (l1max,))
 
 
 def _drc3jj(two_l2, two_l3, two_m2, two_m3):
     """ scalar version of drc3jj """
-    l1max = (two_l2 + two_l3)
+    l1max = two_l2 + two_l3
     l1min = max(np.abs(two_l2 - two_l3), np.abs(two_m2 + two_m3))
     ndim = int((l1max - l1min + 2) / 2)
 
-    l1min, l1max, thrcof, ier = _wigner.drc3jj_int(
-        two_l2, two_l3, two_m2, two_m3, ndim)
+    l1min, l1max, thrcof, ier = _wigner.drc3jj_int(two_l2, two_l3, two_m2, two_m3, ndim)
     if ier == 1:
-        raise ValueError('Either L2.LT.ABS(M2) or L3.LT.ABS(M3).')
+        raise ValueError("Either L2.LT.ABS(M2) or L3.LT.ABS(M3).")
     elif ier == 2:
-        raise ValueError('Either L2+ABS(M2) or L3+ABS(M3) non-integer.')
+        raise ValueError("Either L2+ABS(M2) or L3+ABS(M3) non-integer.")
     elif ier == 3:
-        raise ValueError('L1MAX-L1MIN not an integer.')
+        raise ValueError("L1MAX-L1MIN not an integer.")
     elif ier == 4:
-        raise ValueError('L1MAX less than L1MIN.')
+        raise ValueError("L1MAX less than L1MIN.")
     elif ier == 5:
-        raise ValueError('NDIM less than L1MAX-L1MIN+1.')
+        raise ValueError("NDIM less than L1MAX-L1MIN+1.")
 
     return np.arange(l1min, l1min + ndim * 2, 2), thrcof
 
@@ -258,21 +284,31 @@ def drc6j(two_l2, two_l3, two_l4, two_l5, two_l6):
         The value of 3J symbol
     """
     two_l2, two_l3, two_l4, two_l5, two_l6 = int_broadcast(
-        two_l2, two_l3, two_l4, two_l5, two_l6)
-
-    if ((two_l2 < 0).any() or (two_l3 < 0).any() or (two_l4 < 0).any() or
-            (two_l4 < 0).any() or (two_l5 < 0).any()):
-        raise ValueError('Some of l values are negative')
+        two_l2, two_l3, two_l4, two_l5, two_l6
+    )
+
+    if (
+        (two_l2 < 0).any()
+        or (two_l3 < 0).any()
+        or (two_l4 < 0).any()
+        or (two_l4 < 0).any()
+        or (two_l5 < 0).any()
+    ):
+        raise ValueError("Some of l values are negative")
 
     shape = two_l2.shape
     l1max = int(np.max(two_l2 + two_l3) + 1)
 
     sixcof, ier = _wigner.drc6j_vec(
-        two_l2=two_l2.ravel(), two_l3=two_l3.ravel(),
-        two_l4=two_l4.ravel(), two_l5=two_l5.ravel(),
+        two_l2=two_l2.ravel(),
+        two_l3=two_l3.ravel(),
+        two_l4=two_l4.ravel(),
+        two_l5=two_l5.ravel(),
         two_l6=two_l6.ravel(),
-        nvec=two_l2.size, ndim=l1max)
-    return np.arange(l1max), sixcof.reshape(shape + (l1max, ))
+        nvec=two_l2.size,
+        ndim=l1max,
+    )
+    return np.arange(l1max), sixcof.reshape(shape + (l1max,))
 
 
 def _drc6j(two_l2, two_l3, two_l4, two_l5, two_l6):
@@ -282,16 +318,20 @@ def _drc6j(two_l2, two_l3, two_l4, two_l5, two_l6):
     ndim = int((two_l1max - two_l1min) / 2 + 1)
 
     l1min, l1max, sixcof, ier = _wigner.drc6j_int(
-        two_l2, two_l3, two_l4, two_l5, two_l6, ndim)
+        two_l2, two_l3, two_l4, two_l5, two_l6, ndim
+    )
 
     if ier == 1:
-        raise ValueError('L2+L3+L5+L6 and L2+L4+L6 must be integers')
+        raise ValueError("L2+L3+L5+L6 and L2+L4+L6 must be integers")
     elif ier == 2:
-        raise ValueError('ABS(L2-L4).LE.L6.LE.L2+L4 must be satisfied.')
+        raise ValueError("ABS(L2-L4).LE.L6.LE.L2+L4 must be satisfied.")
     elif ier == 3:
-        raise ValueError('ABS(L4-L5).LE.L3.LE.L4+L5 must be satisfied.')
+        raise ValueError("ABS(L4-L5).LE.L3.LE.L4+L5 must be satisfied.")
     elif ier == 4:
-        raise ValueError('L1MAX-L1MIN must be a non-negative integer, '
-                         'where L1MAX=MIN(L2+L3,L5+L6) and'  'L1MIN=MAX(ABS(L2-L3),ABS(L5-L6)).')
+        raise ValueError(
+            "L1MAX-L1MIN must be a non-negative integer, "
+            "where L1MAX=MIN(L2+L3,L5+L6) and"
+            "L1MIN=MAX(ABS(L2-L3),ABS(L5-L6))."
+        )
 
     return np.arange(l1min, l1min + ndim * 2, 2), sixcof
diff --git a/setup.py b/setup.py
index 672b42d..428c219 100644
--- a/setup.py
+++ b/setup.py
@@ -15,31 +15,40 @@ else:
     raise RuntimeError("Unable to find version string in %s." % (VERSIONFILE,))
 
 # fortran extension
-ext = Extension(name='py3nj._wigner',
-                sources=['fortran/_wigner.pyf', 'fortran/drc.f90',
-                         'fortran/xermsg.f', 'fortran/d1mach.f',
-                         'fortran/drc3jj.f',
-                         'fortran/drc6j.f'],
-                extra_compile_args=['-fopenmp'])
+ext = Extension(
+    name="py3nj._wigner",
+    sources=[
+        "fortran/_wigner.pyf",
+        "fortran/drc.f90",
+        "fortran/xermsg.f",
+        "fortran/d1mach.f",
+        "fortran/drc3jj.f",
+        "fortran/drc6j.f",
+    ],
+    extra_compile_args=["-fopenmp"],
+)
 
 # module
-setup(name='py3nj',
-      version=verstr,
-      author="Keisuke Fujii",
-      author_email="fujiisoup@gmail.com",
-      description=("numpy compatible wigner 3n-J symbols"),
-      license="BSD 3-clause",
-      keywords="atomic physics, quantum physics",
-      url="http://github.com/fujiisoup/py3nj",
-      packages=["py3nj", ],
-      package_dir={'py3nj': 'py3nj'},
-      py_modules=['py3nj.__init__'],
-      ext_modules=[ext],
-      classifiers=['License :: OSI Approved :: BSD License',
-                   'Natural Language :: English',
-                   'Operating System :: MacOS :: MacOS X',
-                   'Operating System :: Microsoft :: Windows',
-                   'Operating System :: POSIX :: Linux',
-                   'Programming Language :: Python :: 3.6',
-                   'Topic :: Scientific/Engineering :: Physics']
-      )
+setup(
+    name="py3nj",
+    version=verstr,
+    author="Keisuke Fujii",
+    author_email="fujiisoup@gmail.com",
+    description=("numpy compatible wigner 3n-J symbols"),
+    license="BSD 3-clause",
+    keywords="atomic physics, quantum physics",
+    url="http://github.com/fujiisoup/py3nj",
+    packages=["py3nj",],
+    package_dir={"py3nj": "py3nj"},
+    py_modules=["py3nj.__init__"],
+    ext_modules=[ext],
+    classifiers=[
+        "License :: OSI Approved :: BSD License",
+        "Natural Language :: English",
+        "Operating System :: MacOS :: MacOS X",
+        "Operating System :: Microsoft :: Windows",
+        "Operating System :: POSIX :: Linux",
+        "Programming Language :: Python :: 3.6",
+        "Topic :: Scientific/Engineering :: Physics",
+    ],
+)


Test patch:
diff --git a/tests/test_3j.py b/tests/test_3j.py
index 521abfa..2348583 100644
--- a/tests/test_3j.py
+++ b/tests/test_3j.py
@@ -6,28 +6,31 @@ from py3nj import wigner, wigner3j
 rng = np.random.RandomState(0)
 # precalculated 3j symbols. Input, result
 THREE_J = (
-    ((1, 1, 0, 0, 0, 0), -np.sqrt(1.0/3)),
-    ((0, 1, 1, 0, 0, 0), -np.sqrt(1.0/3)),
-    ((2, 2, 0, 0, 0, 0),  np.sqrt(1.0/5)),
-    ((2, 1, 1, 0, 0, 0),  np.sqrt(2.0/15)),
-    ((2, 1, 1, 0, 0, 0), np.sqrt(2.0/15)),
-    ((0, 1, 1, 0,-1, 1), np.sqrt(1.0/3)),
-    ((0.5, 0.5, 1.0, -0.5, 0.5, 0), np.sqrt(1.0/6.0)),
-    ((0.5, 0.5, 1.0,  0.5,-0.5, 0), np.sqrt(1.0/6.0)),
-    ((0.5, 0.5, 0.0,  0.5,-0.5, 0), np.sqrt(1.0/2.0)),
-    ((0.5, 0.5, 0.0, -0.5, 0.5, 0), -np.sqrt(1.0/2.0)),
-    ((1, 2, 3, 0, 2, -2), -np.sqrt(1.0/21)),
-    ((1, 2, 3, 1, 2, -3), np.sqrt(1.0/7)),
+    ((1, 1, 0, 0, 0, 0), -np.sqrt(1.0 / 3)),
+    ((0, 1, 1, 0, 0, 0), -np.sqrt(1.0 / 3)),
+    ((2, 2, 0, 0, 0, 0), np.sqrt(1.0 / 5)),
+    ((2, 1, 1, 0, 0, 0), np.sqrt(2.0 / 15)),
+    ((2, 1, 1, 0, 0, 0), np.sqrt(2.0 / 15)),
+    ((0, 1, 1, 0, -1, 1), np.sqrt(1.0 / 3)),
+    ((0.5, 0.5, 1.0, -0.5, 0.5, 0), np.sqrt(1.0 / 6.0)),
+    ((0.5, 0.5, 1.0, 0.5, -0.5, 0), np.sqrt(1.0 / 6.0)),
+    ((0.5, 0.5, 0.0, 0.5, -0.5, 0), np.sqrt(1.0 / 2.0)),
+    ((0.5, 0.5, 0.0, -0.5, 0.5, 0), -np.sqrt(1.0 / 2.0)),
+    ((1, 2, 3, 0, 2, -2), -np.sqrt(1.0 / 21)),
+    ((1, 2, 3, 1, 2, -3), np.sqrt(1.0 / 7)),
+    ((4, 0, 0, 0, 0, 0), 0),
 )
 
 
 def test_drc3jj():
     for three_j, result in THREE_J:
         three_j = (np.array(three_j) * 2).astype(int)
-        l, thrcof = wigner._drc3jj(three_j[1], three_j[2], three_j[4],
-                                   three_j[5])
-        assert three_j[0] in l
-        assert np.allclose(thrcof[l == three_j[0]], result)
+        l, thrcof = wigner._drc3jj(three_j[1], three_j[2], three_j[4], three_j[5])
+
+        if three_j[0] in l:
+            assert np.allclose(thrcof[l == three_j[0]], result)
+        else:
+            assert np.allclose(thrcof, 1)
 
 
 def test_invalid():
@@ -38,7 +41,7 @@ def test_invalid():
         l, thrcof = wigner._drc3jj(1, 3, 4, 0)
 
 
-@pytest.mark.parametrize('half_integer', [1, 0])
+@pytest.mark.parametrize("half_integer", [1, 0])
 def test_drc3jj_vec(half_integer):
     n = 1000
     l2 = rng.randint(0, 10, size=n) * 2 + half_integer
@@ -47,9 +50,9 @@ def test_drc3jj_vec(half_integer):
     m3 = np.zeros(n, dtype=int)
     for i in range(n):
         if l2[i] > 0:
-            m2[i] = rng.randint(-l2[i], l2[i]+1)
+            m2[i] = rng.randint(-l2[i], l2[i] + 1)
         if l3[i] > 0:
-            m3[i] = rng.randint(-l3[i], l3[i]+1)
+            m3[i] = rng.randint(-l3[i], l3[i] + 1)
 
     l, actual = wigner.drc3jj(l2, l3, m2, m3)
     for _ in range(2):
@@ -76,7 +79,7 @@ def test_drc3jj_vec_err():
         _, actual = wigner.drc3jj(l2, l3, m2, m3)
 
 
-@pytest.mark.parametrize('half_integer', [1, 0])
+@pytest.mark.parametrize("half_integer", [1, 0])
 def test_wigner3j(half_integer):
     n = 10000
     n2 = 10
@@ -86,9 +89,9 @@ def test_wigner3j(half_integer):
     m3 = np.zeros(n, dtype=int)
     for i in range(n):
         if l2[i] > 0:
-            m2[i] = rng.randint(-l2[i], l2[i]+1)
+            m2[i] = rng.randint(-l2[i], l2[i] + 1)
         if l3[i] > 0:
-            m3[i] = rng.randint(-l3[i], l3[i]+1)
+            m3[i] = rng.randint(-l3[i], l3[i] + 1)
 
     l, expected_all = wigner.drc3jj(l2, l3, m2, m3)
     for _ in range(n2):
@@ -106,6 +109,15 @@ def test_wigner3j_value():
         three_j = (np.array(three_j) * 2).astype(int)
         actual = wigner3j(*three_j)
         assert np.allclose(actual, expected)
+        # broadcast the first argument
+        actual = wigner3j(np.array([three_j[0]])[:, np.newaxis], *three_j[1:])
+        assert np.allclose(actual, expected)
+        # further broadcasting
+        actual = wigner3j(
+            np.array([three_j[0]])[:, np.newaxis, np.newaxis], *three_j[1:]
+        )
+        assert np.allclose(actual, expected)
+
     # test vector
     three_j = (np.array([thr for thr, _ in THREE_J]).T * 2).astype(int)
     expected = np.array([value for _, value in THREE_J]).T
@@ -119,10 +131,10 @@ def test0d():
         l3 = rng.randint(0, 40) * 2
         m2 = 0
         if l2 > 0:
-            m2 = rng.randint(-l2, l2+1)
+            m2 = rng.randint(-l2, l2 + 1)
         m3 = 0
         if l3 > 0:
-            m3 = rng.randint(-l3, l3+1)
+            m3 = rng.randint(-l3, l3 + 1)
         l, actual = wigner.drc3jj(l2, l3, m2, m3)
         assert actual.ndim == 1
 
@@ -147,9 +159,9 @@ def test2d():
         for i1 in range(n[1]):
             i = (i0, i1)
             if l2[i] > 0:
-                m2[i] = rng.randint(-l2[i], l2[i]+1)
+                m2[i] = rng.randint(-l2[i], l2[i] + 1)
             if l3[i] > 0:
-                m3[i] = rng.randint(-l3[i], l3[i]+1)
+                m3[i] = rng.randint(-l3[i], l3[i] + 1)
 
     l, actual = wigner.drc3jj(l2, l3, m2, m3)
     for i0 in range(n[0]):
diff --git a/tests/test_6j.py b/tests/test_6j.py
index 0874118..b6f81b5 100644
--- a/tests/test_6j.py
+++ b/tests/test_6j.py
@@ -7,11 +7,11 @@ rng = np.random.RandomState(0)
 # precalculated 3j symbols. Input, result
 SIX_J = (
     ((0, 0, 0, 0, 0, 0), 1.0),
-    ((2, 2, 2, 2, 2, 2), -3.0/70),
-    ((4, 2, 2, 2, 2, 2), 2.0/35),
-    ((4, 3, 2, 2, 2, 2), -(1.0/7) * np.sqrt(1.0/2)),
-    ((0.5, 0.5, 1, 0.5, 0.5, 1), 1.0/6),
-    ((0.5, 1.5, 1, 0.5, 0.5, 1), -(1.0/3)),
+    ((2, 2, 2, 2, 2, 2), -3.0 / 70),
+    ((4, 2, 2, 2, 2, 2), 2.0 / 35),
+    ((4, 3, 2, 2, 2, 2), -(1.0 / 7) * np.sqrt(1.0 / 2)),
+    ((0.5, 0.5, 1, 0.5, 0.5, 1), 1.0 / 6),
+    ((0.5, 1.5, 1, 0.5, 0.5, 1), -(1.0 / 3)),
 )
 
 
@@ -76,3 +76,13 @@ def test_wigner6j_value():
     expected = np.array([value for _, value in SIX_J]).T
     actual = wigner6j(*six_j.astype(int))
     assert np.allclose(actual, expected)
+
+
+def test_wigner6j_GH13():
+    """
+    Test for https://github.com/fujiisoup/py3nj/issues/13
+    """
+    N = np.arange(10)
+    J = np.arange(11)[:, np.newaxis]
+    # should not fail
+    wigner6j(0, N * 2, N * 2, J * 2, 2, 2)
\ No newline at end of file
diff --git a/tests/test_9j.py b/tests/test_9j.py
index 78ba8bc..1befb2a 100644
--- a/tests/test_9j.py
+++ b/tests/test_9j.py
@@ -6,17 +6,11 @@ from py3nj import wigner9j
 rng = np.random.RandomState(0)
 # precalculated 3j symbols. Input, result
 NINE_J = (
-    ((1, 1, 1,
-      1, 1, 1,
-      1, 1, 1), 0.0),
+    ((1, 1, 1, 1, 1, 1, 1, 1, 1), 0.0),
     #
-    ((1, 2, 1,
-      1, 1, 1,
-      1, 1, 1), 1.0 / 18),
+    ((1, 2, 1, 1, 1, 1, 1, 1, 1), 1.0 / 18),
     #
-    ((1, 2, 1,
-      2, 3, 1,
-      1, 2, 1), 1.0 / 5.0 * np.sqrt(1.0/30)),
+    ((1, 2, 1, 2, 3, 1, 1, 2, 1), 1.0 / 5.0 * np.sqrt(1.0 / 30)),
 )
 
 
@@ -35,7 +29,7 @@ def test_wigner9j():
 
 
 def test_issue2():
-  # https://github.com/fujiisoup/py3nj/issues/2
-  actual = wigner9j(0, 0, 0, 2 * 1, 0, 2 * 1, 2 * 1, 0, 2 * 1)
-  expected = 1.0 / 3.0
-  assert np.allclose(actual, expected)
\ No newline at end of file
+    # https://github.com/fujiisoup/py3nj/issues/2
+    actual = wigner9j(0, 0, 0, 2 * 1, 0, 2 * 1, 2 * 1, 0, 2 * 1)
+    expected = 1.0 / 3.0
+    assert np.allclose(actual, expected)
diff --git a/tests/test_cg.py b/tests/test_cg.py
index 0bc2a68..9dc7cf5 100644
--- a/tests/test_cg.py
+++ b/tests/test_cg.py
@@ -11,21 +11,14 @@ rng = np.random.RandomState(0)
 
 CG = (
     # j1 = 1/2, j2 = 1/2
-    (( 1, 1, 2,
-       1, 1, 2), 1),
-    (( 1, 1, 2,
-      -1,-1,-2), 1),
-    (( 1, 1, 2,
-       1,-1, 0), np.sqrt(0.5)),
-    (( 1, 1, 2,
-      -1, 1, 0), np.sqrt(0.5)),
-    (( 1, 1, 0,
-       1,-1, 0), np.sqrt(0.5)),
-    (( 1, 1, 0,
-      -1, 1, 0), -np.sqrt(0.5)),
+    ((1, 1, 2, 1, 1, 2), 1),
+    ((1, 1, 2, -1, -1, -2), 1),
+    ((1, 1, 2, 1, -1, 0), np.sqrt(0.5)),
+    ((1, 1, 2, -1, 1, 0), np.sqrt(0.5)),
+    ((1, 1, 0, 1, -1, 0), np.sqrt(0.5)),
+    ((1, 1, 0, -1, 1, 0), -np.sqrt(0.5)),
     # j1 = 1, j2 = 1/2
-    (( 2, 2, 4,
-       2, 0, 2), np.sqrt(1.0/2)),
+    ((2, 2, 4, 2, 0, 2), np.sqrt(1.0 / 2)),
 )