Using associated Legendre polynomials in PyStan model

I have recently started working with PyStan.
My problem is that I need to use associated legendre polynomails for a fit, however the implementation from scipy seems not to be compatible with PyStan.
Are there any preferable solutions to having to implement them all by hand ?

Best way is to call boost and add external cpp code until they’re added to stan. You can help add them or follow the issue Additional combinatorial functions and polynomials · Issue #2309 · stan-dev/math · GitHub