From dc4122163b0bd2bdb045cedf9a980d58b53cdb3a Mon Sep 17 00:00:00 2001 From: Miloslav Metelka <miloslav.metelka@oracle.com> Date: Mon, 26 Mar 2018 16:17:13 +0200 Subject: [PATCH] Implemented Rf_gammafn and others. --- .../ffi/impl/common/JavaUpCallsRFFIImpl.java | 167 +- .../r/ffi/impl/nodes/FFIUpCallNode.java | 9 + .../r/ffi/impl/nodes/RandFunctionsNodes.java | 738 ++++++ .../r/ffi/impl/upcalls/StdUpCallsRFFI.java | 98 +- .../fficall/src/common/rffi_upcalls.h | 32 + .../fficall/src/common/rffi_upcallsindex.h | 380 +-- .../fficall/src/truffle_common/Rmath.c | 97 +- .../builtin/base/BaseBesselFunctions.java | 214 ++ .../builtin/base/BaseGammaFunctions.java | 491 +--- .../r/nodes/builtin/base/BasePackage.java | 7 + .../nodes/builtin/base/TrigExpFunctions.java | 25 +- .../truffle/r/nodes/builtin/base/Trunc.java | 7 +- .../truffle/r/nodes/builtin/InternalNode.java | 4 +- .../com/oracle/truffle/r/runtime/RError.java | 6 +- .../r/runtime/nmath/BesselFunctions.java | 2257 +++++++++++++++++ .../oracle/truffle/r/runtime/nmath/Beta.java | 51 + .../r/runtime/nmath/GammaFunctions.java | 509 +++- .../r/runtime/nmath/MathConstants.java | 27 +- .../oracle/truffle/r/runtime/nmath/RMath.java | 135 +- .../truffle/r/runtime/nmath/distr/Logis.java | 4 +- .../truffle/r/runtime/nmath/distr/Pnorm.java | 7 + .../packages/testrffi/testrffi/src/testrffi.c | 255 +- .../truffle/r/test/ExpectedTestOutput.test | 142 +- .../r/test/builtins/TestBuiltin_besselI.java | 15 +- .../r/test/builtins/TestBuiltin_besselJ.java | 20 +- .../r/test/builtins/TestBuiltin_besselK.java | 23 +- .../r/test/builtins/TestBuiltin_besselY.java | 24 +- .../gnu_r_ihaka_core.copyright.star.regex | 2 +- mx.fastr/copyrights/overrides | 4 +- 29 files changed, 4889 insertions(+), 861 deletions(-) create mode 100644 com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseBesselFunctions.java create mode 100644 com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/BesselFunctions.java create mode 100644 com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/Beta.java diff --git a/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/common/JavaUpCallsRFFIImpl.java b/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/common/JavaUpCallsRFFIImpl.java index c1c4043dd0..7fbdac2f8d 100644 --- a/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/common/JavaUpCallsRFFIImpl.java +++ b/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/common/JavaUpCallsRFFIImpl.java @@ -107,6 +107,7 @@ import com.oracle.truffle.r.runtime.ffi.UnsafeAdapter; import com.oracle.truffle.r.runtime.ffi.VectorRFFIWrapper; import com.oracle.truffle.r.runtime.gnur.SA_TYPE; import com.oracle.truffle.r.runtime.gnur.SEXPTYPE; +import com.oracle.truffle.r.runtime.nmath.RMath; import com.oracle.truffle.r.runtime.nodes.RNode; import com.oracle.truffle.r.runtime.nodes.RSyntaxNode; import com.oracle.truffle.r.runtime.rng.RRNG; @@ -1635,6 +1636,11 @@ public abstract class JavaUpCallsRFFIImpl implements UpCallsRFFI { throw implementedAsNode(); } + @Override + public void Rf_pnorm_both(double arg0, Object arg1, Object arg2, int arg3, int arg4) { + throw implementedAsNode(); + } + @Override public double Rf_dlnorm(double a, double b, double c, int d) { throw implementedAsNode(); @@ -1675,6 +1681,31 @@ public abstract class JavaUpCallsRFFIImpl implements UpCallsRFFI { throw implementedAsNode(); } + @Override + public double Rf_log1pmx(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_log1pexp(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_lgamma1p(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_logspace_add(double a, double b) { + throw implementedAsNode(); + } + + @Override + public double Rf_logspace_sub(double a, double b) { + throw implementedAsNode(); + } + @Override public double Rf_dbeta(double a, double b, double c, int d) { throw implementedAsNode(); @@ -2040,13 +2071,139 @@ public abstract class JavaUpCallsRFFIImpl implements UpCallsRFFI { throw implementedAsNode(); } + @Override + public double Rf_gammafn(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_lgammafn(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_lgammafn_sign(double a, Object b) { + throw implementedAsNode(); + } + + @Override + public void Rf_dpsifn(double a, int b, int c, int d, Object e, Object f, Object g) { + throw implementedAsNode(); + } + + @Override + public double Rf_psigamma(double a, double b) { + throw implementedAsNode(); + } + + @Override + public double Rf_digamma(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_trigamma(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_tetragamma(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_pentagamma(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_beta(double a, double b) { + throw implementedAsNode(); + } + + @Override + public double Rf_lbeta(double a, double b) { + throw implementedAsNode(); + } + + @Override + public double Rf_choose(double a, double b) { + throw implementedAsNode(); + } + + @Override + public double Rf_lchoose(double a, double b) { + throw implementedAsNode(); + } + + @Override + public double Rf_bessel_i(double a, double b, double c) { + throw implementedAsNode(); + } + + @Override + public double Rf_bessel_j(double a, double b) { + throw implementedAsNode(); + } + + @Override + public double Rf_bessel_k(double a, double b, double c) { + throw implementedAsNode(); + } + + @Override + public double Rf_bessel_y(double a, double b) { + throw implementedAsNode(); + } + + @Override + public double Rf_bessel_i_ex(double a, double b, double c, Object d) { + throw implementedAsNode(); + } + + @Override + public double Rf_bessel_j_ex(double a, double b, Object c) { + throw implementedAsNode(); + } + + @Override + public double Rf_bessel_k_ex(double a, double b, double c, Object d) { + throw implementedAsNode(); + } + + @Override + public double Rf_bessel_y_ex(double a, double b, Object c) { + throw implementedAsNode(); + } + + @Override + public double Rf_sign(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_fprec(double a, double b) { + throw implementedAsNode(); + } + @Override public double Rf_ftrunc(double a) { - if (a > 0) { - return Math.floor(a); - } else { - return Math.ceil(a); - } + return RMath.trunc(a); + } + + @Override + public double Rf_cospi(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_sinpi(double a) { + throw implementedAsNode(); + } + + @Override + public double Rf_tanpi(double a) { + throw implementedAsNode(); } @Override diff --git a/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/nodes/FFIUpCallNode.java b/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/nodes/FFIUpCallNode.java index 3b00756efd..b5726a116d 100644 --- a/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/nodes/FFIUpCallNode.java +++ b/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/nodes/FFIUpCallNode.java @@ -104,4 +104,13 @@ public abstract class FFIUpCallNode extends Node { return 6; } } + + public abstract static class Arg7 extends FFIUpCallNode { + public abstract Object executeObject(Object arg0, Object arg1, Object arg2, Object arg3, Object arg4, Object arg5, Object arg6); + + @Override + protected int numArgs() { + return 7; + } + } } diff --git a/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/nodes/RandFunctionsNodes.java b/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/nodes/RandFunctionsNodes.java index 50a546ea96..db0abf411a 100644 --- a/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/nodes/RandFunctionsNodes.java +++ b/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/nodes/RandFunctionsNodes.java @@ -26,8 +26,10 @@ import com.oracle.truffle.api.CompilerDirectives; import com.oracle.truffle.api.dsl.Cached; import com.oracle.truffle.api.dsl.Specialization; import com.oracle.truffle.api.interop.ForeignAccess; +import com.oracle.truffle.api.interop.InteropException; import com.oracle.truffle.api.interop.Message; import com.oracle.truffle.api.interop.TruffleObject; +import com.oracle.truffle.api.interop.UnknownIdentifierException; import com.oracle.truffle.api.interop.UnsupportedMessageException; import com.oracle.truffle.api.nodes.Node; import com.oracle.truffle.r.runtime.RInternalError; @@ -35,9 +37,17 @@ import com.oracle.truffle.r.runtime.data.RDataFactory; import com.oracle.truffle.r.runtime.data.RNull; import com.oracle.truffle.r.runtime.data.model.RAbstractDoubleVector; import com.oracle.truffle.r.runtime.data.model.RAbstractIntVector; +import com.oracle.truffle.r.runtime.data.nodes.GetReadonlyData; import com.oracle.truffle.r.runtime.data.nodes.VectorAccess; import com.oracle.truffle.r.runtime.data.nodes.VectorAccess.SequentialIterator; +import com.oracle.truffle.r.runtime.ffi.interop.UnsafeAdapter; import com.oracle.truffle.r.runtime.interop.ForeignArray2R; +import com.oracle.truffle.r.runtime.nmath.BesselFunctions; +import com.oracle.truffle.r.runtime.nmath.Beta; +import com.oracle.truffle.r.runtime.nmath.Choose; +import com.oracle.truffle.r.runtime.nmath.GammaFunctions; +import com.oracle.truffle.r.runtime.nmath.LBeta; +import com.oracle.truffle.r.runtime.nmath.MathConstants; import com.oracle.truffle.r.runtime.nmath.MathFunctions; import com.oracle.truffle.r.runtime.nmath.MathFunctions.Function2_1; import com.oracle.truffle.r.runtime.nmath.MathFunctions.Function2_2; @@ -45,12 +55,15 @@ import com.oracle.truffle.r.runtime.nmath.MathFunctions.Function3_1; import com.oracle.truffle.r.runtime.nmath.MathFunctions.Function3_2; import com.oracle.truffle.r.runtime.nmath.MathFunctions.Function4_1; import com.oracle.truffle.r.runtime.nmath.MathFunctions.Function4_2; +import com.oracle.truffle.r.runtime.nmath.RMath; import com.oracle.truffle.r.runtime.nmath.RandomFunctions.RandFunction1_Double; import com.oracle.truffle.r.runtime.nmath.RandomFunctions.RandFunction2_Double; import com.oracle.truffle.r.runtime.nmath.RandomFunctions.RandFunction3_Double; import com.oracle.truffle.r.runtime.nmath.RandomFunctions.RandomNumberProvider; import com.oracle.truffle.r.runtime.nmath.RMultinom; +import com.oracle.truffle.r.runtime.nmath.distr.Logis.PLogis; import com.oracle.truffle.r.runtime.nmath.distr.Rbinom; +import com.oracle.truffle.r.runtime.nmath.distr.Pnorm.PnormBoth; public final class RandFunctionsNodes { @@ -300,4 +313,729 @@ public final class RandFunctionsNodes { } + public abstract static class RfPnormBothNode extends FFIUpCallNode.Arg5 { + + @Child private Node cumIsPointerNode = Message.IS_POINTER.createNode(); + @Child private Node cumAsPointerNode; + @Child private Node cumRead; + @Child private Node cumWrite; + @Child private Node ccumIsPointerNode = Message.IS_POINTER.createNode(); + @Child private Node ccumAsPointerNode; + @Child private Node ccumWrite; + + @Specialization + protected Object evaluate(double x, Object cum, Object ccum, int lowerTail, int logP) { + // cum is R/W double* with size==1 and ccum is Writeonly double* with size==1 + double[] cumArr = new double[1]; + double[] ccumArr = new double[1]; + long cumPtr; + long ccumPtr; + TruffleObject cumTO = (TruffleObject) cum; + if (ForeignAccess.sendIsPointer(cumIsPointerNode, cumTO)) { + if (cumAsPointerNode == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + cumAsPointerNode = insert(Message.AS_POINTER.createNode()); + } + try { + cumPtr = ForeignAccess.sendAsPointer(cumAsPointerNode, cumTO); + cumArr[0] = UnsafeAdapter.UNSAFE.getDouble(cumPtr); + } catch (UnsupportedMessageException e) { + throw RInternalError.shouldNotReachHere("IS_POINTER message returned true, AS_POINTER should not fail"); + } + } else { + cumPtr = 0L; + if (cumRead == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + cumRead = Message.READ.createNode(); + cumWrite = Message.WRITE.createNode(); + } + try { + cumArr[0] = ((Number) ForeignAccess.sendRead(cumRead, cumTO, 0)).doubleValue(); + } catch (UnsupportedMessageException | UnknownIdentifierException e) { + throw RInternalError.shouldNotReachHere(e); + } + } + + TruffleObject ccumTO = (TruffleObject) ccum; + if (ForeignAccess.sendIsPointer(ccumIsPointerNode, ccumTO)) { + if (ccumAsPointerNode == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + ccumAsPointerNode = insert(Message.AS_POINTER.createNode()); + } + try { + ccumPtr = ForeignAccess.sendAsPointer(ccumAsPointerNode, ccumTO); + } catch (UnsupportedMessageException e) { + throw RInternalError.shouldNotReachHere("IS_POINTER message returned true, AS_POINTER should not fail"); + } + } else { + ccumPtr = 0L; + if (ccumWrite == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + ccumWrite = Message.WRITE.createNode(); + } + } + + PnormBoth.evaluate(x, cumArr, ccumArr, lowerTail != 0, logP != 0); + if (cumPtr != 0L) { + UnsafeAdapter.UNSAFE.putDouble(cumPtr, cumArr[0]); + } else { + try { + ForeignAccess.sendWrite(cumWrite, cumTO, 0, cumArr[0]); + } catch (InteropException e) { + throw RInternalError.shouldNotReachHere("WRITE message support expected"); + } + } + if (ccumPtr != 0L) { + UnsafeAdapter.UNSAFE.putDouble(ccumPtr, ccumArr[0]); + } else { + try { + ForeignAccess.sendWrite(ccumWrite, ccumTO, 0, ccumArr[0]); + } catch (InteropException e) { + throw RInternalError.shouldNotReachHere("WRITE message support expected"); + } + } + return RNull.instance; + } + + public static RfPnormBothNode create() { + return RandFunctionsNodesFactory.RfPnormBothNodeGen.create(); + } + } + + public abstract static class Log1pmxNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double a) { + return GammaFunctions.log1pmx(a); + } + + public static Log1pmxNode create() { + return RandFunctionsNodesFactory.Log1pmxNodeGen.create(); + } + } + + public abstract static class Log1pexpNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double a) { + return PLogis.log1pexp(a); + } + + public static Log1pexpNode create() { + return RandFunctionsNodesFactory.Log1pexpNodeGen.create(); + } + } + + public abstract static class Lgamma1pNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double a) { + return GammaFunctions.lgamma1p(a); + } + + public static Lgamma1pNode create() { + return RandFunctionsNodesFactory.Lgamma1pNodeGen.create(); + } + } + + public abstract static class LogspaceAddNode extends FFIUpCallNode.Arg2 { + + @Specialization + protected Object evaluate(double a, double b) { + return MathConstants.logspaceAdd(a, b); + } + + public static LogspaceAddNode create() { + return RandFunctionsNodesFactory.LogspaceAddNodeGen.create(); + } + } + + public abstract static class LogspaceSubNode extends FFIUpCallNode.Arg2 { + + @Specialization + protected Object evaluate(double a, double b) { + return MathConstants.logspaceSub(a, b); + } + + public static LogspaceSubNode create() { + return RandFunctionsNodesFactory.LogspaceSubNodeGen.create(); + } + } + + public abstract static class GammafnNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double a) { + return GammaFunctions.gammafn(a); + } + + public static GammafnNode create() { + return RandFunctionsNodesFactory.GammafnNodeGen.create(); + } + } + + public abstract static class LGammafnNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double a) { + return GammaFunctions.lgammafn(a); + } + + public static LGammafnNode create() { + return RandFunctionsNodesFactory.LGammafnNodeGen.create(); + } + } + + public abstract static class LGammafnSignNode extends FFIUpCallNode.Arg2 { + + @Child private Node sgnIsPointerNode = Message.IS_POINTER.createNode(); + @Child private Node sgnAsPointerNode; + @Child private Node sgnRead; + @Child private Node sgnWrite; + + @Specialization + protected Object evaluate(double a, Object sgn) { + int[] sgnArr = new int[1]; + long sgnPtr; + TruffleObject sgnTO = (TruffleObject) sgn; + if (ForeignAccess.sendIsPointer(sgnIsPointerNode, sgnTO)) { + if (sgnAsPointerNode == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + sgnAsPointerNode = insert(Message.AS_POINTER.createNode()); + } + try { + sgnPtr = ForeignAccess.sendAsPointer(sgnAsPointerNode, sgnTO); + sgnArr[0] = UnsafeAdapter.UNSAFE.getInt(sgnPtr); + } catch (UnsupportedMessageException e) { + throw RInternalError.shouldNotReachHere("IS_POINTER message returned true, AS_POINTER should not fail"); + } + } else { + sgnPtr = 0L; + if (sgnRead == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + sgnRead = Message.READ.createNode(); + sgnWrite = Message.WRITE.createNode(); + } + try { + sgnArr[0] = ((Number) ForeignAccess.sendRead(sgnRead, sgnTO, 0)).intValue(); + } catch (UnsupportedMessageException | UnknownIdentifierException e) { + throw RInternalError.shouldNotReachHere(e); + } + } + + double result = GammaFunctions.lgammafnSign(a, sgnArr); + if (sgnPtr != 0L) { + UnsafeAdapter.UNSAFE.putInt(sgnPtr, sgnArr[0]); + } else { + try { + ForeignAccess.sendWrite(sgnWrite, sgnTO, 0, sgnArr[0]); + } catch (InteropException e) { + throw RInternalError.shouldNotReachHere("WRITE message support expected"); + } + } + return result; + } + + public static LGammafnSignNode create() { + return RandFunctionsNodesFactory.LGammafnSignNodeGen.create(); + } + } + + public abstract static class DpsiFnNode extends FFIUpCallNode.Arg7 { + + @Child private Node ansIsPointerNode = Message.IS_POINTER.createNode(); + @Child private Node ansAsPointerNode; + @Child private Node ansRead; + @Child private Node ansWrite; + @Child private Node nzIsPointerNode = Message.IS_POINTER.createNode(); + @Child private Node nzAsPointerNode; + @Child private Node nzRead; + @Child private Node nzWrite; + @Child private Node ierrIsPointerNode = Message.IS_POINTER.createNode(); + @Child private Node ierrAsPointerNode; + @Child private Node ierrRead; + @Child private Node ierrWrite; + + @Specialization + protected Object evaluate(double x, int n, int kode, int m, Object ans, Object nz, Object ierr) { + // ans is R/W double* size==1 and nz is R/W int* size==1 + // ierr is R/W int* size==1 + double ansIn; + int nzIn; + int ierrIn; + TruffleObject ansTO = (TruffleObject) ans; + long ansPtr; + if (ForeignAccess.sendIsPointer(ansIsPointerNode, ansTO)) { + if (ansAsPointerNode == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + ansAsPointerNode = insert(Message.AS_POINTER.createNode()); + } + try { + ansPtr = ForeignAccess.sendAsPointer(ansAsPointerNode, ansTO); + ansIn = UnsafeAdapter.UNSAFE.getDouble(ansPtr); + } catch (UnsupportedMessageException e) { + throw RInternalError.shouldNotReachHere("IS_POINTER message returned true, AS_POINTER should not fail"); + } + } else { + ansPtr = 0L; + if (ansRead == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + ansRead = Message.READ.createNode(); + ansWrite = Message.WRITE.createNode(); + } + try { + ansIn = ((Number) ForeignAccess.sendRead(ansRead, ansTO, 0)).doubleValue(); + } catch (UnsupportedMessageException | UnknownIdentifierException e) { + throw RInternalError.shouldNotReachHere(e); + } + } + + TruffleObject nzTO = (TruffleObject) nz; + long nzPtr; + if (ForeignAccess.sendIsPointer(nzIsPointerNode, nzTO)) { + if (nzAsPointerNode == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + nzAsPointerNode = insert(Message.AS_POINTER.createNode()); + } + try { + nzPtr = ForeignAccess.sendAsPointer(nzAsPointerNode, nzTO); + nzIn = UnsafeAdapter.UNSAFE.getInt(nzPtr); + } catch (UnsupportedMessageException e) { + throw RInternalError.shouldNotReachHere("IS_POINTER message returned true, AS_POINTER should not fail"); + } + } else { + nzPtr = 0L; + if (nzRead == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + nzRead = Message.READ.createNode(); + nzWrite = Message.WRITE.createNode(); + } + try { + nzIn = ((Number) ForeignAccess.sendRead(nzRead, nzTO, 0)).intValue(); + } catch (UnsupportedMessageException | UnknownIdentifierException e) { + throw RInternalError.shouldNotReachHere(e); + } + } + + TruffleObject ierrTO = (TruffleObject) ierr; + long ierrPtr; + if (ForeignAccess.sendIsPointer(ierrIsPointerNode, ierrTO)) { + if (ierrAsPointerNode == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + ierrAsPointerNode = insert(Message.AS_POINTER.createNode()); + } + try { + ierrPtr = ForeignAccess.sendAsPointer(ierrAsPointerNode, ierrTO); + ierrIn = UnsafeAdapter.UNSAFE.getInt(ierrPtr); + } catch (UnsupportedMessageException e) { + throw RInternalError.shouldNotReachHere("IS_POINTER message returned true, AS_POINTER should not fail"); + } + } else { + ierrPtr = 0L; + if (ierrRead == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + ierrRead = Message.READ.createNode(); + ierrWrite = Message.WRITE.createNode(); + } + try { + ierrIn = ((Number) ForeignAccess.sendRead(ierrRead, ierrTO, 0)).intValue(); + } catch (UnsupportedMessageException | UnknownIdentifierException e) { + throw RInternalError.shouldNotReachHere(e); + } + } + + GammaFunctions.DpsiFnResult result = GammaFunctions.dpsifn(x, n, kode, ansIn, nzIn, ierrIn); + if (ansPtr != 0L) { + UnsafeAdapter.UNSAFE.putDouble(ansPtr, result.ans); + } else { + try { + ForeignAccess.sendWrite(ansWrite, ansTO, 0, result.ans); + } catch (InteropException e) { + throw RInternalError.shouldNotReachHere("WRITE message support expected"); + } + } + if (nzPtr != 0L) { + UnsafeAdapter.UNSAFE.putInt(nzPtr, result.nz); + } else { + try { + ForeignAccess.sendWrite(nzWrite, nzTO, 0, result.nz); + } catch (InteropException e) { + throw RInternalError.shouldNotReachHere("WRITE message support expected"); + } + } + if (ierrPtr != 0L) { + UnsafeAdapter.UNSAFE.putInt(ierrPtr, result.ierr); + } else { + try { + ForeignAccess.sendWrite(ierrWrite, ierrTO, 0, result.ierr); + } catch (InteropException e) { + throw RInternalError.shouldNotReachHere("WRITE message support expected"); + } + } + return RNull.instance; + } + + public static DpsiFnNode create() { + return RandFunctionsNodesFactory.DpsiFnNodeGen.create(); + } + } + + public abstract static class PsiGammaNode extends FFIUpCallNode.Arg2 { + + @Specialization + protected Object evaluate(double x, double deriv) { + return GammaFunctions.psigamma(x, deriv); + } + + public static PsiGammaNode create() { + return RandFunctionsNodesFactory.PsiGammaNodeGen.create(); + } + } + + public abstract static class DiGammaNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double x) { + return GammaFunctions.digamma(x); + } + + public static DiGammaNode create() { + return RandFunctionsNodesFactory.DiGammaNodeGen.create(); + } + } + + public abstract static class TriGammaNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double x) { + return GammaFunctions.trigamma(x); + } + + public static TriGammaNode create() { + return RandFunctionsNodesFactory.TriGammaNodeGen.create(); + } + } + + public abstract static class TetraGammaNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double x) { + return GammaFunctions.tetragamma(x); + } + + public static TetraGammaNode create() { + return RandFunctionsNodesFactory.TetraGammaNodeGen.create(); + } + } + + public abstract static class PentaGammaNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double x) { + return GammaFunctions.pentagamma(x); + } + + public static PentaGammaNode create() { + return RandFunctionsNodesFactory.PentaGammaNodeGen.create(); + } + } + + public abstract static class BetaNode extends FFIUpCallNode.Arg2 { + + @Specialization + protected Object evaluate(double a, double b) { + return Beta.beta(a, b); + } + + public static BetaNode create() { + return RandFunctionsNodesFactory.BetaNodeGen.create(); + } + } + + public abstract static class LBetaNode extends FFIUpCallNode.Arg2 { + + @Specialization + protected Object evaluate(double a, double b) { + return LBeta.lbeta(a, b); + } + + public static LBetaNode create() { + return RandFunctionsNodesFactory.LBetaNodeGen.create(); + } + } + + public abstract static class ChooseNode extends FFIUpCallNode.Arg2 { + + @Specialization + protected Object evaluate(double n, double k) { + return Choose.choose(n, k); + } + + public static ChooseNode create() { + return RandFunctionsNodesFactory.ChooseNodeGen.create(); + } + } + + public abstract static class LChooseNode extends FFIUpCallNode.Arg2 { + + @Specialization + protected Object evaluate(double n, double k) { + return Choose.lchoose(n, k); + } + + public static LChooseNode create() { + return RandFunctionsNodesFactory.LChooseNodeGen.create(); + } + } + + public abstract static class BesselINode extends FFIUpCallNode.Arg3 { + + @Specialization + protected Object evaluate(double a, double b, double c) { + return BesselFunctions.bessel_i(a, b, c); + } + + public static BesselINode create() { + return RandFunctionsNodesFactory.BesselINodeGen.create(); + } + } + + public abstract static class BesselJNode extends FFIUpCallNode.Arg2 { + + @Specialization + protected Object evaluate(double a, double b) { + return BesselFunctions.bessel_j(a, b); + } + + public static BesselJNode create() { + return RandFunctionsNodesFactory.BesselJNodeGen.create(); + } + } + + public abstract static class BesselKNode extends FFIUpCallNode.Arg3 { + + @Specialization + protected Object evaluate(double a, double b, double c) { + return BesselFunctions.bessel_k(a, b, c); + } + + public static BesselKNode create() { + return RandFunctionsNodesFactory.BesselKNodeGen.create(); + } + } + + public abstract static class BesselYNode extends FFIUpCallNode.Arg2 { + + @Specialization + protected Object evaluate(double a, double b) { + return BesselFunctions.bessel_y(a, b); + } + + public static BesselYNode create() { + return RandFunctionsNodesFactory.BesselYNodeGen.create(); + } + } + + public abstract static class BesselIExNode extends FFIUpCallNode.Arg4 { + + @Specialization + protected Object evaluate(final double a, final double b, final double c, Object d, + @Cached("create()") BesselExNode besselEx) { + return besselEx.execute(new BesselExCaller() { + @Override + public double call(double[] arr) { + return BesselFunctions.bessel_i_ex(a, b, c, arr); + } + + @Override + public int arrLen() { + return (int) Math.floor(b) + 1; + } + }, d); + } + + public static BesselIExNode create() { + return RandFunctionsNodesFactory.BesselIExNodeGen.create(); + } + } + + public abstract static class BesselJExNode extends FFIUpCallNode.Arg3 { + + @Specialization + protected Object evaluate(final double a, final double b, Object c, + @Cached("create()") BesselExNode besselEx) { + return besselEx.execute(new BesselExCaller() { + @Override + public double call(double[] arr) { + return BesselFunctions.bessel_j_ex(a, b, arr); + } + + @Override + public int arrLen() { + return (int) Math.floor(b) + 1; + } + }, c); + } + + public static BesselJExNode create() { + return RandFunctionsNodesFactory.BesselJExNodeGen.create(); + } + } + + public abstract static class BesselKExNode extends FFIUpCallNode.Arg4 { + + @Specialization + protected Object evaluate(double a, double b, double c, Object d, + @Cached("create()") BesselExNode besselEx) { + return besselEx.execute(new BesselExCaller() { + @Override + public double call(double[] arr) { + return BesselFunctions.bessel_k_ex(a, b, c, arr); + } + + @Override + public int arrLen() { + return (int) Math.floor(b) + 1; + } + }, d); + } + + public static BesselKExNode create() { + return RandFunctionsNodesFactory.BesselKExNodeGen.create(); + } + } + + public abstract static class BesselYExNode extends FFIUpCallNode.Arg3 { + + @Specialization + protected Object evaluate(double a, double b, Object c, + @Cached("create()") BesselExNode besselEx) { + return besselEx.execute(new BesselExCaller() { + @Override + public double call(double[] arr) { + return BesselFunctions.bessel_y_ex(a, b, arr); + } + + @Override + public int arrLen() { + return (int) Math.floor(b) + 1; + } + }, c); + } + + public static BesselYExNode create() { + return RandFunctionsNodesFactory.BesselYExNodeGen.create(); + } + } + + interface BesselExCaller { + + double call(double[] arr); + + int arrLen(); + + } + + abstract static class BesselExNode extends Node { + + @Child private Node bIsPointerNode = Message.IS_POINTER.createNode(); + @Child private Node bAsPointerNode; + @Child private ForeignArray2R bForeignArray2R; + + public abstract double execute(BesselExCaller caller, Object b); + + @Specialization + protected double besselEx(BesselExCaller caller, Object b, + @Cached("create()") GetReadonlyData.Double bReadonlyData) { + RAbstractDoubleVector bVec; + TruffleObject bTO = (TruffleObject) b; + if (ForeignAccess.sendIsPointer(bIsPointerNode, bTO)) { + if (bAsPointerNode == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + bAsPointerNode = insert(Message.AS_POINTER.createNode()); + } + long addr; + try { + addr = ForeignAccess.sendAsPointer(bAsPointerNode, bTO); + bVec = RDataFactory.createDoubleVectorFromNative(addr, caller.arrLen()); + } catch (UnsupportedMessageException e) { + throw RInternalError.shouldNotReachHere("IS_POINTER message returned true, AS_POINTER should not fail"); + } + } else { + if (bForeignArray2R == null) { + CompilerDirectives.transferToInterpreterAndInvalidate(); + bForeignArray2R = ForeignArray2R.create(); + } + bVec = (RAbstractDoubleVector) bForeignArray2R.convert(bTO); + } + + return caller.call(bReadonlyData.execute(bVec.materialize())); + } + + public static BesselExNode create() { + return RandFunctionsNodesFactory.BesselExNodeGen.create(); + } + + } + + public abstract static class SignNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double x) { + return RMath.sign(x); + } + + public static SignNode create() { + return RandFunctionsNodesFactory.SignNodeGen.create(); + } + } + + public abstract static class FPrecNode extends FFIUpCallNode.Arg2 { + + @Specialization + protected Object evaluate(double x, double digits) { + return RMath.fprec(x, digits); + } + + public static FPrecNode create() { + return RandFunctionsNodesFactory.FPrecNodeGen.create(); + } + } + + public abstract static class SinpiNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double x) { + return RMath.sinpi(x); + } + + public static SinpiNode create() { + return RandFunctionsNodesFactory.SinpiNodeGen.create(); + } + } + + public abstract static class CospiNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double x) { + return RMath.cospi(x); + } + + public static CospiNode create() { + return RandFunctionsNodesFactory.CospiNodeGen.create(); + } + } + + public abstract static class TanpiNode extends FFIUpCallNode.Arg1 { + + @Specialization + protected Object evaluate(double x) { + return RMath.tanpi(x); + } + + public static TanpiNode create() { + return RandFunctionsNodesFactory.TanpiNodeGen.create(); + } + } + } diff --git a/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/upcalls/StdUpCallsRFFI.java b/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/upcalls/StdUpCallsRFFI.java index 09df801945..dc0f358098 100644 --- a/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/upcalls/StdUpCallsRFFI.java +++ b/com.oracle.truffle.r.ffi.impl/src/com/oracle/truffle/r/ffi/impl/upcalls/StdUpCallsRFFI.java @@ -528,6 +528,9 @@ public interface StdUpCallsRFFI { @RFFIUpCallNode(value = RandFunctionsNodes.RandFunction2Node.class, functionClass = Rnorm.class) double Rf_rnorm(double a, double b); + @RFFIUpCallNode(value = RandFunctionsNodes.RfPnormBothNode.class) + void Rf_pnorm_both(double a, @RFFICpointer Object b, @RFFICpointer Object c, int d, int e); + @RFFIUpCallNode(value = RandFunctionsNodes.RandFunction3_1Node.class, functionClass = LogNormal.DLNorm.class) double Rf_dlnorm(double a, double b, double c, int d); @@ -552,6 +555,21 @@ public interface StdUpCallsRFFI { @RFFIUpCallNode(value = RandFunctionsNodes.RandFunction2Node.class, functionClass = RGamma.class) double Rf_rgamma(double a, double b); + @RFFIUpCallNode(RandFunctionsNodes.Log1pmxNode.class) + double Rf_log1pmx(double a); + + @RFFIUpCallNode(RandFunctionsNodes.Log1pexpNode.class) + double Rf_log1pexp(double a); + + @RFFIUpCallNode(RandFunctionsNodes.Lgamma1pNode.class) + double Rf_lgamma1p(double a); + + @RFFIUpCallNode(RandFunctionsNodes.LogspaceAddNode.class) + double Rf_logspace_add(double a, double b); + + @RFFIUpCallNode(RandFunctionsNodes.LogspaceSubNode.class) + double Rf_logspace_sub(double a, double b); + @RFFIUpCallNode(value = RandFunctionsNodes.RandFunction3_1Node.class, functionClass = DBeta.class) double Rf_dbeta(double a, double b, double c, int d); @@ -672,7 +690,6 @@ public interface StdUpCallsRFFI { @RFFIUpCallNode(value = RandFunctionsNodes.RandFunction2Node.class, functionClass = RNBinom.RNBinomMu.class) double Rf_rnbinom_mu(double a, double b); - @RFFIRunGC @RFFIUpCallNode(value = RandFunctionsNodes.RfRMultinomNode.class) void Rf_rmultinom(int a, @RFFICpointer Object b, int c, @RFFICpointer Object d); @@ -721,6 +738,7 @@ public interface StdUpCallsRFFI { @RFFIUpCallNode(value = RandFunctionsNodes.RandFunction4_2Node.class, functionClass = QNBeta.class) double Rf_qnbeta(double a, double b, double c, double d, int e, int f); + // Unable to find implementation of Rf_rnbeta // double Rf_rnbeta(double a, double b, double c); @RFFIUpCallNode(value = RandFunctionsNodes.RandFunction4_1Node.class, functionClass = Dnf.class) @@ -771,8 +789,86 @@ public interface StdUpCallsRFFI { @RFFIUpCallNode(value = RandFunctionsNodes.RandFunction1Node.class, functionClass = Signrank.RSignrank.class) double Rf_rsignrank(double a); + @RFFIUpCallNode(RandFunctionsNodes.GammafnNode.class) + double Rf_gammafn(double a); + + @RFFIUpCallNode(RandFunctionsNodes.LGammafnNode.class) + double Rf_lgammafn(double a); + + @RFFIUpCallNode(RandFunctionsNodes.LGammafnSignNode.class) + double Rf_lgammafn_sign(double a, @RFFICpointer Object b); + + @RFFIUpCallNode(RandFunctionsNodes.DpsiFnNode.class) + void Rf_dpsifn(double a, int b, int c, int d, @RFFICpointer Object e, @RFFICpointer Object f, @RFFICpointer Object g); + + @RFFIUpCallNode(RandFunctionsNodes.PsiGammaNode.class) + double Rf_psigamma(double a, double b); + + @RFFIUpCallNode(RandFunctionsNodes.DiGammaNode.class) + double Rf_digamma(double a); + + @RFFIUpCallNode(RandFunctionsNodes.TriGammaNode.class) + double Rf_trigamma(double a); + + @RFFIUpCallNode(RandFunctionsNodes.TetraGammaNode.class) + double Rf_tetragamma(double a); + + @RFFIUpCallNode(RandFunctionsNodes.PentaGammaNode.class) + double Rf_pentagamma(double a); + + @RFFIUpCallNode(RandFunctionsNodes.BetaNode.class) + double Rf_beta(double a, double b); + + @RFFIUpCallNode(RandFunctionsNodes.LBetaNode.class) + double Rf_lbeta(double a, double b); + + @RFFIUpCallNode(RandFunctionsNodes.ChooseNode.class) + double Rf_choose(double a, double b); + + @RFFIUpCallNode(RandFunctionsNodes.LChooseNode.class) + double Rf_lchoose(double a, double b); + + @RFFIUpCallNode(RandFunctionsNodes.BesselINode.class) + double Rf_bessel_i(double a, double b, double c); + + @RFFIUpCallNode(RandFunctionsNodes.BesselJNode.class) + double Rf_bessel_j(double a, double b); + + @RFFIUpCallNode(RandFunctionsNodes.BesselKNode.class) + double Rf_bessel_k(double a, double b, double c); + + @RFFIUpCallNode(RandFunctionsNodes.BesselYNode.class) + double Rf_bessel_y(double a, double b); + + @RFFIUpCallNode(RandFunctionsNodes.BesselIExNode.class) + double Rf_bessel_i_ex(double a, double b, double c, @RFFICpointer Object d); + + @RFFIUpCallNode(RandFunctionsNodes.BesselJExNode.class) + double Rf_bessel_j_ex(double a, double b, @RFFICpointer Object c); + + @RFFIUpCallNode(RandFunctionsNodes.BesselKExNode.class) + double Rf_bessel_k_ex(double a, double b, double c, @RFFICpointer Object d); + + @RFFIUpCallNode(RandFunctionsNodes.BesselYExNode.class) + double Rf_bessel_y_ex(double a, double b, @RFFICpointer Object c); + + @RFFIUpCallNode(RandFunctionsNodes.SignNode.class) + double Rf_sign(double a); + + @RFFIUpCallNode(RandFunctionsNodes.FPrecNode.class) + double Rf_fprec(double a, double b); + double Rf_ftrunc(double a); + @RFFIUpCallNode(RandFunctionsNodes.CospiNode.class) + double Rf_cospi(double a); + + @RFFIUpCallNode(RandFunctionsNodes.SinpiNode.class) + double Rf_sinpi(double a); + + @RFFIUpCallNode(RandFunctionsNodes.TanpiNode.class) + double Rf_tanpi(double a); + @RFFIUpCallNode(MiscNodes.NamesGetsNode.class) Object Rf_namesgets(Object vec, Object val); diff --git a/com.oracle.truffle.r.native/fficall/src/common/rffi_upcalls.h b/com.oracle.truffle.r.native/fficall/src/common/rffi_upcalls.h index b5c5f9feb3..a0dcd7923f 100644 --- a/com.oracle.truffle.r.native/fficall/src/common/rffi_upcalls.h +++ b/com.oracle.truffle.r.native/fficall/src/common/rffi_upcalls.h @@ -287,6 +287,7 @@ typedef double (*call_Rf_dnorm)(double a, double b, double c, int d); typedef double (*call_Rf_pnorm)(double a, double b, double c, int d, int e); typedef double (*call_Rf_qnorm)(double a, double b, double c, int d, int e); typedef double (*call_Rf_rnorm)(double a, double b); +typedef void (*call_Rf_pnorm_both)(double a, double* b, double* c, int d, int e); typedef double (*call_Rf_dlnorm)(double a, double b, double c, int d); typedef double (*call_Rf_plnorm)(double a, double b, double c, int d, int e); typedef double (*call_Rf_qlnorm)(double a, double b, double c, int d, int e); @@ -295,6 +296,11 @@ typedef double (*call_Rf_dgamma)(double a, double b, double c, int d); typedef double (*call_Rf_pgamma)(double a, double b, double c, int d, int e); typedef double (*call_Rf_qgamma)(double a, double b, double c, int d, int e); typedef double (*call_Rf_rgamma)(double a, double b); +typedef double (*call_Rf_log1pmx)(double a); +typedef double (*call_Rf_log1pexp)(double a); +typedef double (*call_Rf_lgamma1p)(double a); +typedef double (*call_Rf_logspace_add)(double a, double b); +typedef double (*call_Rf_logspace_sub)(double a, double b); typedef double (*call_Rf_dbeta)(double a, double b, double c, int d); typedef double (*call_Rf_pbeta)(double a, double b, double c, int d, int e); typedef double (*call_Rf_qbeta)(double a, double b, double c, int d, int e); @@ -367,7 +373,33 @@ typedef double (*call_Rf_dsignrank)(double a, double b, int c); typedef double (*call_Rf_psignrank)(double a, double b, int c, int d); typedef double (*call_Rf_qsignrank)(double a, double b, int c, int d); typedef double (*call_Rf_rsignrank)(double a); +typedef double (*call_Rf_gammafn)(double a); +typedef double (*call_Rf_lgammafn)(double a); +typedef double (*call_Rf_lgammafn_sign)(double a, int* b); +typedef void (*call_Rf_dpsifn)(double a, int b, int c, int d, double* e, int* f, int* g); +typedef double (*call_Rf_psigamma)(double a, double b); +typedef double (*call_Rf_digamma)(double a); +typedef double (*call_Rf_trigamma)(double a); +typedef double (*call_Rf_tetragamma)(double a); +typedef double (*call_Rf_pentagamma)(double a); +typedef double (*call_Rf_beta)(double a, double b); +typedef double (*call_Rf_lbeta)(double a, double b); +typedef double (*call_Rf_choose)(double a, double b); +typedef double (*call_Rf_lchoose)(double a, double b); +typedef double (*call_Rf_bessel_i)(double a, double b, double c); +typedef double (*call_Rf_bessel_j)(double a, double b); +typedef double (*call_Rf_bessel_k)(double a, double b, double c); +typedef double (*call_Rf_bessel_y)(double a, double b); +typedef double (*call_Rf_bessel_i_ex)(double a, double b, double c, double * d); +typedef double (*call_Rf_bessel_j_ex)(double a, double b, double * c); +typedef double (*call_Rf_bessel_k_ex)(double a, double b, double c, double * d); +typedef double (*call_Rf_bessel_y_ex)(double a, double b, double * c); +typedef double (*call_Rf_sign)(double a); +typedef double (*call_Rf_fprec)(double a, double b); typedef double (*call_Rf_ftrunc)(double a); +typedef double (*call_Rf_cospi)(double a); +typedef double (*call_Rf_sinpi)(double a); +typedef double (*call_Rf_tanpi)(double a); typedef SEXP (*call_Rf_match)(SEXP itable, SEXP ix, int nmatch); typedef Rboolean (*call_Rf_NonNullStringMatch)(SEXP s, SEXP t); typedef SEXP (*call_getvar)(); diff --git a/com.oracle.truffle.r.native/fficall/src/common/rffi_upcallsindex.h b/com.oracle.truffle.r.native/fficall/src/common/rffi_upcallsindex.h index 71610473d3..b785e1ccea 100644 --- a/com.oracle.truffle.r.native/fficall/src/common/rffi_upcallsindex.h +++ b/com.oracle.truffle.r.native/fficall/src/common/rffi_upcallsindex.h @@ -111,180 +111,212 @@ #define Rf_asInteger_x 105 #define Rf_asLogical_x 106 #define Rf_asReal_x 107 -#define Rf_classgets_x 108 -#define Rf_coerceVector_x 109 -#define Rf_cons_x 110 -#define Rf_copyListMatrix_x 111 -#define Rf_copyMatrix_x 112 -#define Rf_copyMostAttrib_x 113 -#define Rf_dbeta_x 114 -#define Rf_dbinom_x 115 -#define Rf_dcauchy_x 116 -#define Rf_dchisq_x 117 -#define Rf_defineVar_x 118 -#define Rf_dexp_x 119 -#define Rf_df_x 120 -#define Rf_dgamma_x 121 -#define Rf_dgeom_x 122 -#define Rf_dhyper_x 123 -#define Rf_dlnorm_x 124 -#define Rf_dlogis_x 125 -#define Rf_dnbeta_x 126 -#define Rf_dnbinom_x 127 -#define Rf_dnbinom_mu_x 128 -#define Rf_dnchisq_x 129 -#define Rf_dnf_x 130 -#define Rf_dnorm4_x 131 -#define Rf_dnt_x 132 -#define Rf_dpois_x 133 -#define Rf_dsignrank_x 134 -#define Rf_dt_x 135 -#define Rf_dunif_x 136 -#define Rf_duplicate_x 137 -#define Rf_dweibull_x 138 -#define Rf_dwilcox_x 139 -#define Rf_error_x 140 -#define Rf_errorcall_x 141 -#define Rf_eval_x 142 -#define Rf_findFun_x 143 -#define Rf_findVar_x 144 -#define Rf_findVarInFrame_x 145 -#define Rf_findVarInFrame3_x 146 -#define Rf_ftrunc_x 147 -#define Rf_getAttrib_x 148 -#define Rf_gsetVar_x 149 -#define Rf_inherits_x 150 -#define Rf_install_x 151 -#define Rf_installChar_x 152 -#define Rf_isNull_x 153 -#define Rf_isString_x 154 -#define Rf_lengthgets_x 155 -#define Rf_match_x 156 -#define Rf_mkCharLenCE_x 157 -#define Rf_namesgets_x 158 -#define Rf_ncols_x 159 -#define Rf_nrows_x 160 -#define Rf_pbeta_x 161 -#define Rf_pbinom_x 162 -#define Rf_pcauchy_x 163 -#define Rf_pchisq_x 164 -#define Rf_pexp_x 165 -#define Rf_pf_x 166 -#define Rf_pgamma_x 167 -#define Rf_pgeom_x 168 -#define Rf_phyper_x 169 -#define Rf_plnorm_x 170 -#define Rf_plogis_x 171 -#define Rf_pnbeta_x 172 -#define Rf_pnbinom_x 173 -#define Rf_pnbinom_mu_x 174 -#define Rf_pnchisq_x 175 -#define Rf_pnf_x 176 -#define Rf_pnorm5_x 177 -#define Rf_pnt_x 178 -#define Rf_ppois_x 179 -#define Rf_protect_x 180 -#define Rf_psignrank_x 181 -#define Rf_pt_x 182 -#define Rf_ptukey_x 183 -#define Rf_punif_x 184 -#define Rf_pweibull_x 185 -#define Rf_pwilcox_x 186 -#define Rf_qbeta_x 187 -#define Rf_qbinom_x 188 -#define Rf_qcauchy_x 189 -#define Rf_qchisq_x 190 -#define Rf_qexp_x 191 -#define Rf_qf_x 192 -#define Rf_qgamma_x 193 -#define Rf_qgeom_x 194 -#define Rf_qhyper_x 195 -#define Rf_qlnorm_x 196 -#define Rf_qlogis_x 197 -#define Rf_qnbeta_x 198 -#define Rf_qnbinom_x 199 -#define Rf_qnbinom_mu_x 200 -#define Rf_qnchisq_x 201 -#define Rf_qnf_x 202 -#define Rf_qnorm5_x 203 -#define Rf_qnt_x 204 -#define Rf_qpois_x 205 -#define Rf_qsignrank_x 206 -#define Rf_qt_x 207 -#define Rf_qtukey_x 208 -#define Rf_qunif_x 209 -#define Rf_qweibull_x 210 -#define Rf_qwilcox_x 211 -#define Rf_rbeta_x 212 -#define Rf_rbinom_x 213 -#define Rf_rcauchy_x 214 -#define Rf_rchisq_x 215 -#define Rf_rexp_x 216 -#define Rf_rf_x 217 -#define Rf_rgamma_x 218 -#define Rf_rgeom_x 219 -#define Rf_rhyper_x 220 -#define Rf_rlnorm_x 221 -#define Rf_rlogis_x 222 -#define Rf_rmultinom_x 223 -#define Rf_rnbinom_x 224 -#define Rf_rnbinom_mu_x 225 -#define Rf_rnchisq_x 226 -#define Rf_rnorm_x 227 -#define Rf_rpois_x 228 -#define Rf_rsignrank_x 229 -#define Rf_rt_x 230 -#define Rf_runif_x 231 -#define Rf_rweibull_x 232 -#define Rf_rwilcox_x 233 -#define Rf_setAttrib_x 234 -#define Rf_str2type_x 235 -#define Rf_unprotect_x 236 -#define Rf_unprotect_ptr_x 237 -#define Rf_warning_x 238 -#define Rf_warningcall_x 239 -#define Rprintf_x 240 -#define SETCAD4R_x 241 -#define SETCADDDR_x 242 -#define SETCADDR_x 243 -#define SETCADR_x 244 -#define SETCAR_x 245 -#define SETCDR_x 246 -#define SET_ATTRIB_x 247 -#define SET_BODY_x 248 -#define SET_CLOENV_x 249 -#define SET_FORMALS_x 250 -#define SET_NAMED_FASTR_x 251 -#define SET_OBJECT_x 252 -#define SET_RDEBUG_x 253 -#define SET_RSTEP_x 254 -#define SET_S4_OBJECT_x 255 -#define SET_STRING_ELT_x 256 -#define SET_SYMVALUE_x 257 -#define SET_TAG_x 258 -#define SET_TYPEOF_x 259 -#define SET_VECTOR_ELT_x 260 -#define STRING_ELT_x 261 -#define SYMVALUE_x 262 -#define TAG_x 263 -#define TYPEOF_x 264 -#define UNSET_S4_OBJECT_x 265 -#define VECTOR_ELT_x 266 -#define forceSymbols_x 267 -#define getCCallable_x 268 -#define getConnectionClassString_x 269 -#define getEmbeddingDLLInfo_x 270 -#define getOpenModeString_x 271 -#define getSummaryDescription_x 272 -#define isSeekable_x 273 -#define octsize_x 274 -#define registerCCallable_x 275 -#define registerRoutines_x 276 -#define restoreHandlerStacks_x 277 -#define setDotSymbolValues_x 278 -#define unif_rand_x 279 -#define useDynamicSymbols_x 280 +#define Rf_bessel_i_x 108 +#define Rf_bessel_i_ex_x 109 +#define Rf_bessel_j_x 110 +#define Rf_bessel_j_ex_x 111 +#define Rf_bessel_k_x 112 +#define Rf_bessel_k_ex_x 113 +#define Rf_bessel_y_x 114 +#define Rf_bessel_y_ex_x 115 +#define Rf_beta_x 116 +#define Rf_choose_x 117 +#define Rf_classgets_x 118 +#define Rf_coerceVector_x 119 +#define Rf_cons_x 120 +#define Rf_copyListMatrix_x 121 +#define Rf_copyMatrix_x 122 +#define Rf_copyMostAttrib_x 123 +#define Rf_cospi_x 124 +#define Rf_dbeta_x 125 +#define Rf_dbinom_x 126 +#define Rf_dcauchy_x 127 +#define Rf_dchisq_x 128 +#define Rf_defineVar_x 129 +#define Rf_dexp_x 130 +#define Rf_df_x 131 +#define Rf_dgamma_x 132 +#define Rf_dgeom_x 133 +#define Rf_dhyper_x 134 +#define Rf_digamma_x 135 +#define Rf_dlnorm_x 136 +#define Rf_dlogis_x 137 +#define Rf_dnbeta_x 138 +#define Rf_dnbinom_x 139 +#define Rf_dnbinom_mu_x 140 +#define Rf_dnchisq_x 141 +#define Rf_dnf_x 142 +#define Rf_dnorm4_x 143 +#define Rf_dnt_x 144 +#define Rf_dpois_x 145 +#define Rf_dpsifn_x 146 +#define Rf_dsignrank_x 147 +#define Rf_dt_x 148 +#define Rf_dunif_x 149 +#define Rf_duplicate_x 150 +#define Rf_dweibull_x 151 +#define Rf_dwilcox_x 152 +#define Rf_error_x 153 +#define Rf_errorcall_x 154 +#define Rf_eval_x 155 +#define Rf_findFun_x 156 +#define Rf_findVar_x 157 +#define Rf_findVarInFrame_x 158 +#define Rf_findVarInFrame3_x 159 +#define Rf_fprec_x 160 +#define Rf_ftrunc_x 161 +#define Rf_gammafn_x 162 +#define Rf_getAttrib_x 163 +#define Rf_gsetVar_x 164 +#define Rf_inherits_x 165 +#define Rf_install_x 166 +#define Rf_installChar_x 167 +#define Rf_isNull_x 168 +#define Rf_isString_x 169 +#define Rf_lbeta_x 170 +#define Rf_lchoose_x 171 +#define Rf_lengthgets_x 172 +#define Rf_lgamma1p_x 173 +#define Rf_lgammafn_x 174 +#define Rf_lgammafn_sign_x 175 +#define Rf_log1pexp_x 176 +#define Rf_log1pmx_x 177 +#define Rf_logspace_add_x 178 +#define Rf_logspace_sub_x 179 +#define Rf_match_x 180 +#define Rf_mkCharLenCE_x 181 +#define Rf_namesgets_x 182 +#define Rf_ncols_x 183 +#define Rf_nrows_x 184 +#define Rf_pbeta_x 185 +#define Rf_pbinom_x 186 +#define Rf_pcauchy_x 187 +#define Rf_pchisq_x 188 +#define Rf_pentagamma_x 189 +#define Rf_pexp_x 190 +#define Rf_pf_x 191 +#define Rf_pgamma_x 192 +#define Rf_pgeom_x 193 +#define Rf_phyper_x 194 +#define Rf_plnorm_x 195 +#define Rf_plogis_x 196 +#define Rf_pnbeta_x 197 +#define Rf_pnbinom_x 198 +#define Rf_pnbinom_mu_x 199 +#define Rf_pnchisq_x 200 +#define Rf_pnf_x 201 +#define Rf_pnorm5_x 202 +#define Rf_pnorm_both_x 203 +#define Rf_pnt_x 204 +#define Rf_ppois_x 205 +#define Rf_protect_x 206 +#define Rf_psigamma_x 207 +#define Rf_psignrank_x 208 +#define Rf_pt_x 209 +#define Rf_ptukey_x 210 +#define Rf_punif_x 211 +#define Rf_pweibull_x 212 +#define Rf_pwilcox_x 213 +#define Rf_qbeta_x 214 +#define Rf_qbinom_x 215 +#define Rf_qcauchy_x 216 +#define Rf_qchisq_x 217 +#define Rf_qexp_x 218 +#define Rf_qf_x 219 +#define Rf_qgamma_x 220 +#define Rf_qgeom_x 221 +#define Rf_qhyper_x 222 +#define Rf_qlnorm_x 223 +#define Rf_qlogis_x 224 +#define Rf_qnbeta_x 225 +#define Rf_qnbinom_x 226 +#define Rf_qnbinom_mu_x 227 +#define Rf_qnchisq_x 228 +#define Rf_qnf_x 229 +#define Rf_qnorm5_x 230 +#define Rf_qnt_x 231 +#define Rf_qpois_x 232 +#define Rf_qsignrank_x 233 +#define Rf_qt_x 234 +#define Rf_qtukey_x 235 +#define Rf_qunif_x 236 +#define Rf_qweibull_x 237 +#define Rf_qwilcox_x 238 +#define Rf_rbeta_x 239 +#define Rf_rbinom_x 240 +#define Rf_rcauchy_x 241 +#define Rf_rchisq_x 242 +#define Rf_rexp_x 243 +#define Rf_rf_x 244 +#define Rf_rgamma_x 245 +#define Rf_rgeom_x 246 +#define Rf_rhyper_x 247 +#define Rf_rlnorm_x 248 +#define Rf_rlogis_x 249 +#define Rf_rmultinom_x 250 +#define Rf_rnbinom_x 251 +#define Rf_rnbinom_mu_x 252 +#define Rf_rnchisq_x 253 +#define Rf_rnorm_x 254 +#define Rf_rpois_x 255 +#define Rf_rsignrank_x 256 +#define Rf_rt_x 257 +#define Rf_runif_x 258 +#define Rf_rweibull_x 259 +#define Rf_rwilcox_x 260 +#define Rf_setAttrib_x 261 +#define Rf_sign_x 262 +#define Rf_sinpi_x 263 +#define Rf_str2type_x 264 +#define Rf_tanpi_x 265 +#define Rf_tetragamma_x 266 +#define Rf_trigamma_x 267 +#define Rf_unprotect_x 268 +#define Rf_unprotect_ptr_x 269 +#define Rf_warning_x 270 +#define Rf_warningcall_x 271 +#define Rprintf_x 272 +#define SETCAD4R_x 273 +#define SETCADDDR_x 274 +#define SETCADDR_x 275 +#define SETCADR_x 276 +#define SETCAR_x 277 +#define SETCDR_x 278 +#define SET_ATTRIB_x 279 +#define SET_BODY_x 280 +#define SET_CLOENV_x 281 +#define SET_FORMALS_x 282 +#define SET_NAMED_FASTR_x 283 +#define SET_OBJECT_x 284 +#define SET_RDEBUG_x 285 +#define SET_RSTEP_x 286 +#define SET_S4_OBJECT_x 287 +#define SET_STRING_ELT_x 288 +#define SET_SYMVALUE_x 289 +#define SET_TAG_x 290 +#define SET_TYPEOF_x 291 +#define SET_VECTOR_ELT_x 292 +#define STRING_ELT_x 293 +#define SYMVALUE_x 294 +#define TAG_x 295 +#define TYPEOF_x 296 +#define UNSET_S4_OBJECT_x 297 +#define VECTOR_ELT_x 298 +#define forceSymbols_x 299 +#define getCCallable_x 300 +#define getConnectionClassString_x 301 +#define getEmbeddingDLLInfo_x 302 +#define getOpenModeString_x 303 +#define getSummaryDescription_x 304 +#define isSeekable_x 305 +#define octsize_x 306 +#define registerCCallable_x 307 +#define registerRoutines_x 308 +#define restoreHandlerStacks_x 309 +#define setDotSymbolValues_x 310 +#define unif_rand_x 311 +#define useDynamicSymbols_x 312 -#define UPCALLS_TABLE_SIZE 281 +#define UPCALLS_TABLE_SIZE 313 #endif // RFFI_UPCALLSINDEX_H diff --git a/com.oracle.truffle.r.native/fficall/src/truffle_common/Rmath.c b/com.oracle.truffle.r.native/fficall/src/truffle_common/Rmath.c index 7e763b29d3..a01625e10d 100644 --- a/com.oracle.truffle.r.native/fficall/src/truffle_common/Rmath.c +++ b/com.oracle.truffle.r.native/fficall/src/truffle_common/Rmath.c @@ -51,9 +51,8 @@ double Rf_rnorm(double a, double b) { return ((call_Rf_rnorm) callbacks[Rf_rnorm_x])(a, b); } -void Rf_pnorm_both(double a, double * b, double * c, int d, int e) { - unimplemented("Rf_rnorm"); -// return ((call_Rf_pnorm_both) callbacks[Rf_pnorm_both_x])(a, b, c, d, e); +void Rf_pnorm_both(double a, double* b, double* c, int d, int e) { + ((call_Rf_pnorm_both) callbacks[Rf_pnorm_both_x])(a, b, c, d, e); } double Rf_dunif(double a, double b, double c, int d) { @@ -89,28 +88,23 @@ double Rf_rgamma(double a, double b) { } double Rf_log1pmx(double a) { - unimplemented("Rf_log1pmx"); - return 0; + return ((call_Rf_log1pmx) callbacks[Rf_log1pmx_x])(a); } double Rf_log1pexp(double a) { - unimplemented("Rf_log1pexp"); - return 0; + return ((call_Rf_log1pexp) callbacks[Rf_log1pexp_x])(a); } double Rf_lgamma1p(double a) { - unimplemented("Rf_lgamma1p"); - return 0; + return ((call_Rf_lgamma1p) callbacks[Rf_lgamma1p_x])(a); } double Rf_logspace_add(double a, double b) { - unimplemented("Rf_logspace_add"); - return 0; + return ((call_Rf_logspace_add) callbacks[Rf_logspace_add_x])(a, b); } double Rf_logspace_sub(double a, double b) { - unimplemented("Rf_logspace_sub"); - return 0; + return ((call_Rf_logspace_sub) callbacks[Rf_logspace_sub_x])(a, b); } double Rf_dbeta(double a, double b, double c, int d) { @@ -455,107 +449,87 @@ double Rf_rsignrank(double a) { } double Rf_gammafn(double a) { - unimplemented("Rf_gammafn"); - return 0; + return ((call_Rf_gammafn) callbacks[Rf_gammafn_x])(a); } double Rf_lgammafn(double a) { - unimplemented("Rf_lgammafn"); - return 0; + return ((call_Rf_lgammafn) callbacks[Rf_lgammafn_x])(a); } double Rf_lgammafn_sign(double a, int* b) { - unimplemented("Rf_lgammafn_sign"); - return 0; + return ((call_Rf_lgammafn_sign) callbacks[Rf_lgammafn_sign_x])(a, b); } void Rf_dpsifn(double a, int b, int c, int d, double* e, int* f, int* g) { - unimplemented("Rf_dpsifn"); + return ((call_Rf_dpsifn) callbacks[Rf_dpsifn_x])(a, b, c, d, e, f, g); } double Rf_psigamma(double a, double b) { - unimplemented("Rf_psigamma"); - return 0; + return ((call_Rf_psigamma) callbacks[Rf_psigamma_x])(a, b); } double Rf_digamma(double a) { - unimplemented("Rf_digamma"); - return 0; + return ((call_Rf_digamma) callbacks[Rf_digamma_x])(a); } double Rf_trigamma(double a) { - unimplemented("Rf_trigamma"); - return 0; + return ((call_Rf_trigamma) callbacks[Rf_trigamma_x])(a); } double Rf_tetragamma(double a) { - unimplemented("Rf_tetragamma"); - return 0; + return ((call_Rf_tetragamma) callbacks[Rf_tetragamma_x])(a); } double Rf_pentagamma(double a) { - unimplemented("Rf_pentagamma"); - return 0; + return ((call_Rf_pentagamma) callbacks[Rf_pentagamma_x])(a); } double Rf_beta(double a, double b) { - unimplemented("Rf_beta"); - return 0; + return ((call_Rf_beta) callbacks[Rf_beta_x])(a, b); } double Rf_lbeta(double a, double b) { - unimplemented("Rf_lbeta"); - return 0; + return ((call_Rf_lbeta) callbacks[Rf_lbeta_x])(a, b); } double Rf_choose(double a, double b) { - unimplemented("Rf_choose"); - return 0; + return ((call_Rf_choose) callbacks[Rf_choose_x])(a, b); } double Rf_lchoose(double a, double b) { - unimplemented("Rf_lchoose"); - return 0; + return ((call_Rf_lchoose) callbacks[Rf_lchoose_x])(a, b); } double Rf_bessel_i(double a, double b, double c) { - unimplemented("Rf_bessel_i"); - return 0; + return ((call_Rf_bessel_i) callbacks[Rf_bessel_i_x])(a, b, c); } double Rf_bessel_j(double a, double b) { - unimplemented("Rf_bessel_j"); - return 0; + return ((call_Rf_bessel_j) callbacks[Rf_bessel_j_x])(a, b); } double Rf_bessel_k(double a, double b, double c) { - unimplemented("Rf_bessel_k"); - return 0; + return ((call_Rf_bessel_k) callbacks[Rf_bessel_k_x])(a, b, c); } double Rf_bessel_y(double a, double b) { - unimplemented("Rf_bessel_y"); - return 0; + return ((call_Rf_bessel_y) callbacks[Rf_bessel_y_x])(a, b); } double Rf_bessel_i_ex(double a, double b, double c, double * d) { - unimplemented("Rf_bessel_i_ex"); - return 0; + return ((call_Rf_bessel_i_ex) callbacks[Rf_bessel_i_ex_x])(a, b, c, d); } double Rf_bessel_j_ex(double a, double b, double * c) { - unimplemented("Rf_bessel_j_ex"); - return 0; + return ((call_Rf_bessel_j_ex) callbacks[Rf_bessel_j_ex_x])(a, b, c); } double Rf_bessel_k_ex(double a, double b, double c, double * d) { - unimplemented("Rf_bessel_k_ex"); - return 0; + return ((call_Rf_bessel_k_ex) callbacks[Rf_bessel_k_ex_x])(a, b, c, d); } double Rf_bessel_y_ex(double a, double b, double * c) { - unimplemented("Rf_bessel_y_ex"); - return 0; + return ((call_Rf_bessel_y_ex) callbacks[Rf_bessel_y_ex_x])(a, b, c); } int Rf_imax2(int x, int y) { @@ -575,13 +549,11 @@ double Rf_fmin2(double x, double y) { } double Rf_sign(double a) { - unimplemented("Rf_sign"); - return 0; + return ((call_Rf_sign) callbacks[Rf_sign_x])(a); } double Rf_fprec(double a, double b) { - unimplemented("Rf_fprec"); - return 0; + return ((call_Rf_fprec) callbacks[Rf_fprec_x])(a, b); } double Rf_fsign(double a, double b) { @@ -597,17 +569,14 @@ double Rf_ftrunc(double a) { } double Rf_cospi(double a) { - unimplemented("Rf_cospi"); - return 0; + return ((call_Rf_cospi) callbacks[Rf_cospi_x])(a); } double Rf_sinpi(double a) { - unimplemented("Rf_sinpi"); - return 0; + return ((call_Rf_sinpi) callbacks[Rf_sinpi_x])(a); } double Rf_tanpi(double a) { - unimplemented("Rf_tanpi"); - return 0; + return ((call_Rf_tanpi) callbacks[Rf_tanpi_x])(a); } diff --git a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseBesselFunctions.java b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseBesselFunctions.java new file mode 100644 index 0000000000..49a6fdc2e5 --- /dev/null +++ b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseBesselFunctions.java @@ -0,0 +1,214 @@ +/* + * Copyright (c) 2018, Oracle and/or its affiliates. All rights reserved. + * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. + * + * This code is free software; you can redistribute it and/or modify it + * under the terms of the GNU General Public License version 3 only, as + * published by the Free Software Foundation. + * + * This code is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * version 3 for more details (a copy is included in the LICENSE file that + * accompanied this code). + * + * You should have received a copy of the GNU General Public License version + * 3 along with this work; if not, write to the Free Software Foundation, + * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. + * + * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA + * or visit www.oracle.com if you need additional information or have any + * questions. + */ +package com.oracle.truffle.r.nodes.builtin.base; + +import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.numericValue; +import static com.oracle.truffle.r.runtime.builtins.RBehavior.PURE; +import static com.oracle.truffle.r.runtime.builtins.RBuiltinKind.INTERNAL; + +import com.oracle.truffle.api.dsl.Cached; +import com.oracle.truffle.api.dsl.Specialization; +import com.oracle.truffle.r.nodes.builtin.RBuiltinNode; +import com.oracle.truffle.r.runtime.RRuntime; +import com.oracle.truffle.r.runtime.builtins.RBuiltin; +import com.oracle.truffle.r.runtime.data.RDataFactory.VectorFactory; +import com.oracle.truffle.r.runtime.data.model.RAbstractDoubleVector; +import com.oracle.truffle.r.runtime.data.nodes.VectorAccess; +import com.oracle.truffle.r.runtime.data.nodes.VectorAccess.SequentialIterator; +import com.oracle.truffle.r.runtime.nmath.BesselFunctions; + +public class BaseBesselFunctions { + + @RBuiltin(name = "besselI", kind = INTERNAL, parameterNames = {"x", "nu", "expo"}, behavior = PURE) + public abstract static class BesselI extends RBuiltinNode.Arg3 { + + static { + Casts casts = new Casts(BesselI.class); + casts.arg(0).mustBe(numericValue()).asDoubleVector(); + casts.arg(1).mustBe(numericValue()).asDoubleVector(); + casts.arg(2).mustBe(numericValue()).asDoubleVector(); + } + + @Specialization(guards = {"xAccess.supports(x)", "nuAccess.supports(nu)", "expoAccess.supports(expo)"}) + protected RAbstractDoubleVector doFast(RAbstractDoubleVector x, RAbstractDoubleVector nu, RAbstractDoubleVector expo, + @Cached("x.access()") VectorAccess xAccess, + @Cached("nu.access()") VectorAccess nuAccess, + @Cached("expo.access()") VectorAccess expoAccess, + @Cached("create()") VectorFactory factory) { + SequentialIterator xIter = xAccess.access(x); + SequentialIterator nuIter = nuAccess.access(nu); + SequentialIterator expoIter = expoAccess.access(expo); + int xLen = xAccess.getLength(xIter); + int nuLen = nuAccess.getLength(nuIter); + int expoLen = expoAccess.getLength(expoIter); + int n = (xLen != 0 && nuLen != 0 && expoLen != 0) ? Math.max(Math.max(xLen, nuLen), expoLen) : 0; + double[] result = new double[n]; + for (int i = 0; i < n; i++) { + xAccess.nextWithWrap(xIter); + nuAccess.nextWithWrap(nuIter); + expoAccess.nextWithWrap(expoIter); + double xElem = xAccess.getDouble(xIter); + double nuElem = nuAccess.getDouble(nuIter); + double expoElem = expoAccess.getDouble(expoIter); + result[i] = (xAccess.na.check(xElem) || nuAccess.na.check(nuElem) || expoAccess.na.check(expoElem)) + ? RRuntime.DOUBLE_NA + : BesselFunctions.bessel_i(xElem, nuElem, expoElem); + } + return factory.createDoubleVector(result, xAccess.na.neverSeenNA() && nuAccess.na.neverSeenNA() && expoAccess.na.neverSeenNA()); + } + + @Specialization(replaces = "doFast") + protected RAbstractDoubleVector doGeneric(RAbstractDoubleVector x, RAbstractDoubleVector nu, RAbstractDoubleVector expo, + @Cached("create()") VectorFactory factory) { + return doFast(x, nu, expo, x.slowPathAccess(), nu.slowPathAccess(), expo.slowPathAccess(), factory); + } + + } + + @RBuiltin(name = "besselJ", kind = INTERNAL, parameterNames = {"x", "nu"}, behavior = PURE) + public abstract static class BesselJ extends RBuiltinNode.Arg2 { + + static { + Casts casts = new Casts(BesselJ.class); + casts.arg(0).mustBe(numericValue()).asDoubleVector(); + casts.arg(1).mustBe(numericValue()).asDoubleVector(); + } + + @Specialization(guards = {"xAccess.supports(x)", "nuAccess.supports(nu)"}) + protected RAbstractDoubleVector doFast(RAbstractDoubleVector x, RAbstractDoubleVector nu, + @Cached("x.access()") VectorAccess xAccess, + @Cached("nu.access()") VectorAccess nuAccess, + @Cached("create()") VectorFactory factory) { + SequentialIterator xIter = xAccess.access(x); + SequentialIterator nuIter = nuAccess.access(nu); + int xLen = xAccess.getLength(xIter); + int nuLen = nuAccess.getLength(nuIter); + int n = (xLen != 0 && nuLen != 0) ? Math.max(xLen, nuLen) : 0; + double[] result = new double[n]; + for (int i = 0; i < n; i++) { + xAccess.nextWithWrap(xIter); + nuAccess.nextWithWrap(nuIter); + double xElem = xAccess.getDouble(xIter); + double nuElem = nuAccess.getDouble(nuIter); + result[i] = (xAccess.na.check(xElem) || nuAccess.na.check(nuElem)) + ? RRuntime.DOUBLE_NA + : BesselFunctions.bessel_j(xElem, nuElem); + } + return factory.createDoubleVector(result, xAccess.na.neverSeenNA() && nuAccess.na.neverSeenNA()); + } + + @Specialization(replaces = "doFast") + protected RAbstractDoubleVector doGeneric(RAbstractDoubleVector x, RAbstractDoubleVector nu, + @Cached("create()") VectorFactory factory) { + return doFast(x, nu, x.slowPathAccess(), nu.slowPathAccess(), factory); + } + + } + + @RBuiltin(name = "besselK", kind = INTERNAL, parameterNames = {"x", "nu", "expo"}, behavior = PURE) + public abstract static class BesselK extends RBuiltinNode.Arg3 { + + static { + Casts casts = new Casts(BesselK.class); + casts.arg(0).mustBe(numericValue()).asDoubleVector(); + casts.arg(1).mustBe(numericValue()).asDoubleVector(); + casts.arg(2).mustBe(numericValue()).asDoubleVector(); + } + + @Specialization(guards = {"xAccess.supports(x)", "nuAccess.supports(nu)", "expoAccess.supports(expo)"}) + protected RAbstractDoubleVector doFast(RAbstractDoubleVector x, RAbstractDoubleVector nu, RAbstractDoubleVector expo, + @Cached("x.access()") VectorAccess xAccess, + @Cached("nu.access()") VectorAccess nuAccess, + @Cached("expo.access()") VectorAccess expoAccess, + @Cached("create()") VectorFactory factory) { + SequentialIterator xIter = xAccess.access(x); + SequentialIterator nuIter = nuAccess.access(nu); + SequentialIterator expoIter = expoAccess.access(expo); + int xLen = xAccess.getLength(xIter); + int nuLen = nuAccess.getLength(nuIter); + int expoLen = expoAccess.getLength(expoIter); + int n = (xLen != 0 && nuLen != 0 && expoLen != 0) ? Math.max(Math.max(xLen, nuLen), expoLen) : 0; + double[] result = new double[n]; + for (int i = 0; i < n; i++) { + xAccess.nextWithWrap(xIter); + nuAccess.nextWithWrap(nuIter); + expoAccess.nextWithWrap(expoIter); + double xElem = xAccess.getDouble(xIter); + double nuElem = nuAccess.getDouble(nuIter); + double expoElem = expoAccess.getDouble(expoIter); + result[i] = (xAccess.na.check(xElem) || nuAccess.na.check(nuElem) || expoAccess.na.check(expoElem)) + ? RRuntime.DOUBLE_NA + : BesselFunctions.bessel_k(xElem, nuElem, expoElem); + } + return factory.createDoubleVector(result, xAccess.na.neverSeenNA() && nuAccess.na.neverSeenNA() && expoAccess.na.neverSeenNA()); + } + + @Specialization(replaces = "doFast") + protected RAbstractDoubleVector doGeneric(RAbstractDoubleVector x, RAbstractDoubleVector nu, RAbstractDoubleVector expo, + @Cached("create()") VectorFactory factory) { + return doFast(x, nu, expo, x.slowPathAccess(), nu.slowPathAccess(), expo.slowPathAccess(), factory); + } + + } + + @RBuiltin(name = "besselY", kind = INTERNAL, parameterNames = {"x", "nu"}, behavior = PURE) + public abstract static class BesselY extends RBuiltinNode.Arg2 { + + static { + Casts casts = new Casts(BesselY.class); + casts.arg(0).mustBe(numericValue()).asDoubleVector(); + casts.arg(1).mustBe(numericValue()).asDoubleVector(); + } + + @Specialization(guards = {"xAccess.supports(x)", "nuAccess.supports(nu)"}) + protected RAbstractDoubleVector doFast(RAbstractDoubleVector x, RAbstractDoubleVector nu, + @Cached("x.access()") VectorAccess xAccess, + @Cached("nu.access()") VectorAccess nuAccess, + @Cached("create()") VectorFactory factory) { + SequentialIterator xIter = xAccess.access(x); + SequentialIterator nuIter = nuAccess.access(nu); + int xLen = xAccess.getLength(xIter); + int nuLen = nuAccess.getLength(nuIter); + int n = (xLen != 0 && nuLen != 0) ? Math.max(xLen, nuLen) : 0; + double[] result = new double[n]; + for (int i = 0; i < n; i++) { + xAccess.nextWithWrap(xIter); + nuAccess.nextWithWrap(nuIter); + double xElem = xAccess.getDouble(xIter); + double nuElem = nuAccess.getDouble(nuIter); + result[i] = (xAccess.na.check(xElem) || nuAccess.na.check(nuElem)) + ? RRuntime.DOUBLE_NA + : BesselFunctions.bessel_y(xElem, nuElem); + } + return factory.createDoubleVector(result, xAccess.na.neverSeenNA() && nuAccess.na.neverSeenNA()); + } + + @Specialization(replaces = "doFast") + protected RAbstractDoubleVector doGeneric(RAbstractDoubleVector x, RAbstractDoubleVector nu, + @Cached("create()") VectorFactory factory) { + return doFast(x, nu, x.slowPathAccess(), nu.slowPathAccess(), factory); + } + + } + +} diff --git a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseGammaFunctions.java b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseGammaFunctions.java index 34ef1ba276..02e4619a72 100644 --- a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseGammaFunctions.java +++ b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BaseGammaFunctions.java @@ -6,7 +6,7 @@ * Copyright (C) 1998 Ross Ihaka * Copyright (c) 1998--2012, The R Core Team * Copyright (c) 2004, The R Foundation - * Copyright (c) 2013, 2017, Oracle and/or its affiliates + * Copyright (c) 2013, 2018, Oracle and/or its affiliates * * All rights reserved. */ @@ -16,32 +16,25 @@ import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.complexValue import static com.oracle.truffle.r.nodes.builtin.CastBuilder.Predef.numericValue; import static com.oracle.truffle.r.runtime.RDispatch.MATH_GROUP_GENERIC; import static com.oracle.truffle.r.runtime.builtins.RBehavior.PURE; +import static com.oracle.truffle.r.runtime.builtins.RBuiltinKind.INTERNAL; import static com.oracle.truffle.r.runtime.builtins.RBuiltinKind.PRIMITIVE; -import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_EPSILON; -import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MANT_DIG; -import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MAX_EXP; -import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MIN_EXP; -import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_LOG10_2; -import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_PI; -import static com.oracle.truffle.r.runtime.nmath.RMath.fmax2; - -import com.oracle.truffle.api.CompilerDirectives; -import com.oracle.truffle.api.CompilerDirectives.CompilationFinal; + import com.oracle.truffle.api.dsl.Cached; import com.oracle.truffle.api.dsl.Specialization; -import com.oracle.truffle.api.nodes.Node; import com.oracle.truffle.r.nodes.attributes.UnaryCopyAttributesNode; import com.oracle.truffle.r.nodes.builtin.RBuiltinNode; -import com.oracle.truffle.r.nodes.builtin.base.BaseGammaFunctionsFactory.DpsiFnCalcNodeGen; +import com.oracle.truffle.r.nodes.primitive.BinaryMapFunctionNode; +import com.oracle.truffle.r.nodes.primitive.BinaryMapNode; import com.oracle.truffle.r.runtime.RError; import com.oracle.truffle.r.runtime.RInternalError; import com.oracle.truffle.r.runtime.RRuntime; +import com.oracle.truffle.r.runtime.RType; import com.oracle.truffle.r.runtime.builtins.RBuiltin; import com.oracle.truffle.r.runtime.data.RDataFactory; import com.oracle.truffle.r.runtime.data.RDoubleVector; import com.oracle.truffle.r.runtime.data.model.RAbstractDoubleVector; +import com.oracle.truffle.r.runtime.data.model.RAbstractVector; import com.oracle.truffle.r.runtime.nmath.GammaFunctions; -import com.oracle.truffle.r.runtime.nmath.RMath; import com.oracle.truffle.r.runtime.ops.na.NACheck; public class BaseGammaFunctions { @@ -107,7 +100,7 @@ public class BaseGammaFunctions { @Override protected double scalarFunction(double xv) { - throw RError.nyi(this, "trigamma"); + return GammaFunctions.trigamma(xv); } } @@ -127,444 +120,96 @@ public class BaseGammaFunctions { @RBuiltin(name = "digamma", kind = PRIMITIVE, parameterNames = {"x"}, dispatch = MATH_GROUP_GENERIC, behavior = PURE) public abstract static class DiGamma extends GammaBase { - @Child private DpsiFnCalc dpsiFnCalc = DpsiFnCalcNodeGen.create(); - static { casts(new Casts(DiGamma.class)); } @Override protected double scalarFunction(double xv) { - return -dpsiFnCalc.executeDouble(xv, 0, 1, 0); + return GammaFunctions.digamma(xv); } } - protected abstract static class DpsiFnCalc extends Node { - - // the following is transcribed from polygamma.c + @RBuiltin(name = "tetragamma", kind = PRIMITIVE, parameterNames = {"x"}, dispatch = MATH_GROUP_GENERIC, behavior = PURE) + public abstract static class TetraGamma extends GammaBase { - public abstract double executeDouble(double x, int n, int kode, double ans); + static { + casts(new Casts(TetraGamma.class)); + } - @Child private DpsiFnCalc dpsiFnCalc; + @Override + protected double scalarFunction(double xv) { + return GammaFunctions.tetragamma(xv); + } + } - @CompilationFinal(dimensions = 1) private static final double[] bvalues = new double[]{1.00000000000000000e+00, -5.00000000000000000e-01, 1.66666666666666667e-01, -3.33333333333333333e-02, - 2.38095238095238095e-02, -3.33333333333333333e-02, 7.57575757575757576e-02, -2.53113553113553114e-01, 1.16666666666666667e+00, -7.09215686274509804e+00, - 5.49711779448621554e+01, -5.29124242424242424e+02, 6.19212318840579710e+03, -8.65802531135531136e+04, 1.42551716666666667e+06, -2.72982310678160920e+07, - 6.01580873900642368e+08, -1.51163157670921569e+10, 4.29614643061166667e+11, -1.37116552050883328e+13, 4.88332318973593167e+14, -1.92965793419400681e+16}; + @RBuiltin(name = "pentagamma", kind = PRIMITIVE, parameterNames = {"x"}, dispatch = MATH_GROUP_GENERIC, behavior = PURE) + public abstract static class PentaGamma extends GammaBase { - private static final int n_max = 100; - // the following is actually a parameter in the original code, but it's always 1 and must be - // as the original code treats the "ans" value of type double as an array, which is legal - // only if a the first element of the array is accessed at all times - private static final int m = 1; + static { + casts(new Casts(PentaGamma.class)); + } - private double dpsiFnCalc(double x, int n, int kode, double ans) { - if (dpsiFnCalc == null) { - CompilerDirectives.transferToInterpreterAndInvalidate(); - dpsiFnCalc = insert(DpsiFnCalcNodeGen.create()); - } - return dpsiFnCalc.executeDouble(x, n, kode, ans); + @Override + protected double scalarFunction(double xv) { + return GammaFunctions.pentagamma(xv); } + } - // TODO: it's recursive - turn into AST recursion - @Specialization - double dpsifn(double xOld, int n, int kode, double ansOld) { - - double x = xOld; - double ans = ansOld; - - int mm; - int mx; - int nn; - int np; - int nx; - int fn; - double arg; - double den; - double elim; - double eps; - double fln; - double rln; - double r1m4; - double r1m5; - double s; - double slope; - double t; - double tk; - double tt; - double t1; - double t2; - double wdtol; - double xdmln; - double xdmy; - double xinc; - double xln = 0.0; - double xm; - double xmin; - double yint; - double[] trm = new double[23]; - double[] trmr = new double[n_max + 1]; - - // non-zero ierr always results in generating a NaN - // mVal.ierr = 0; - if (n < 0 || kode < 1 || kode > 2 || m < 1) { - return Double.NaN; - } - if (x <= 0.) { - /* - * use Abramowitz & Stegun 6.4.7 "Reflection Formula" psi(k, x) = (-1)^k psi(k, 1-x) - * - pi^{n+1} (d/dx)^n cot(x) - */ - if (x == RMath.round(x)) { - /* non-positive integer : +Inf or NaN depends on n */ - // for(j=0; j < m; j++) /* k = j + n : */ - // ans[j] = ((j+n) % 2) ? ML_POSINF : ML_NAN; - // m is always 1 - ans = (n % 2) != 0 ? Double.POSITIVE_INFINITY : Double.NaN; - return ans; - } - /* This could cancel badly */ - ans = dpsiFnCalc(1. - x, n, /* kode = */1, ans); - /* - * ans[j] == (-1)^(k+1) / gamma(k+1) * psi(k, 1 - x) for j = 0:(m-1) , k = n + j - */ - - /* Cheat for now: only work for m = 1, n in {0,1,2,3} : */ - if (m > 1 || n > 3) { /* doesn't happen for digamma() .. pentagamma() */ - /* not yet implemented */ - // non-zero ierr always results in generating a NaN - // mVal.ierr = 4; - return Double.NaN; - } - x *= M_PI; /* pi * x */ - if (n == 0) { - tt = Math.cos(x) / Math.sin(x); - } else if (n == 1) { - tt = -1 / Math.pow(Math.sin(x), 2); - } else if (n == 2) { - tt = 2 * Math.cos(x) / Math.pow(Math.sin(x), 3); - } else if (n == 3) { - tt = -2 * (2 * Math.pow(Math.cos(x), 2) + 1.) / Math.pow(Math.sin(x), 4); - } else { /* can not happen! */ - tt = RRuntime.DOUBLE_NA; - } - /* end cheat */ - - s = (n % 2) != 0 ? -1. : 1.; /* s = (-1)^n */ - /* - * t := pi^(n+1) * d_n(x) / gamma(n+1) , where d_n(x) := (d/dx)^n cot(x) - */ - t1 = t2 = s = 1.; - for (int k = 0, j = k - n; j < m; k++, j++, s = -s) { - /* k == n+j , s = (-1)^k */ - t1 *= M_PI; /* t1 == pi^(k+1) */ - if (k >= 2) { - t2 *= k; /* t2 == k! == gamma(k+1) */ - } - if (j >= 0) { /* by cheat above, tt === d_k(x) */ - // j must always be 0 - assert j == 0; - // ans[j] = s*(ans[j] + t1/t2 * tt); - ans = s * (ans + t1 / t2 * tt); - } - } - if (n == 0 && kode == 2) { /* unused from R, but "wrong": xln === 0 : */ - // ans[0] += xln; - ans += xln; - } - return ans; - } /* x <= 0 */ - - /* else : x > 0 */ - // nz not used - // mVal.nz = 0; - xln = Math.log(x); - if (kode == 1 /* && m == 1 */) { /* the R case --- for very large x: */ - double lrg = 1 / (2. * DBL_EPSILON); - if (n == 0 && x * xln > lrg) { - // ans[0] = -xln; - ans = -xln; - return ans; - } else if (n >= 1 && x > n * lrg) { - // ans[0] = exp(-n * xln)/n; /* == x^-n / n == 1/(n * x^n) */ - ans = Math.exp(-n * xln) / n; - return ans; - } - } - mm = m; - // nx = imin2(-Rf_i1mach(15), Rf_i1mach(16));/* = 1021 */ - nx = Math.min(-DBL_MIN_EXP, DBL_MAX_EXP); - assert (nx == 1021); - r1m5 = M_LOG10_2; // Rf_d1mach(5); - r1m4 = DBL_EPSILON * 0.5; // Rf_d1mach(4) * 0.5; - wdtol = fmax2(r1m4, 0.5e-18); /* 1.11e-16 */ - - /* elim = approximate exponential over and underflow limit */ - elim = 2.302 * (nx * r1m5 - 3.0); /* = 700.6174... */ - for (;;) { - nn = n + mm - 1; - fn = nn; - t = (fn + 1) * xln; - - /* overflow and underflow test for small and large x */ - - if (Math.abs(t) > elim) { - if (t <= 0.0) { - // nz not used - // mVal.nz = 0; - // non-zero ierr always results in generating a NaN - // mVal.ierr = 2; - return Double.NaN; - } - } else { - if (x < wdtol) { - // ans[0] = R_pow_di(x, -n-1); - ans = Math.pow(x, -n - 1); - if (mm != 1) { - // for(k = 1; k < mm ; k++) - // ans[k] = ans[k-1] / x; - assert mm < 2; - // int the original code, ans should not be accessed beyond the 0th - // index - } - if (n == 0 && kode == 2) { - // ans[0] += xln; - ans += xln; - } - return ans; - } + @RBuiltin(name = "psigamma", kind = INTERNAL, parameterNames = {"x", "deriv"}, behavior = PURE) + public abstract static class PsiGamma extends RBuiltinNode.Arg2 { - /* compute xmin and the number of terms of the series, fln+1 */ - - rln = r1m5 * DBL_MANT_DIG; // Rf_i1mach(14); - rln = Math.min(rln, 18.06); - fln = Math.max(rln, 3.0) - 3.0; - yint = 3.50 + 0.40 * fln; - slope = 0.21 + fln * (0.0006038 * fln + 0.008677); - xm = yint + slope * fn; - mx = (int) xm + 1; - xmin = mx; - if (n != 0) { - xm = -2.302 * rln - Math.min(0.0, xln); - arg = xm / n; - arg = Math.min(0.0, arg); - eps = Math.exp(arg); - xm = 1.0 - eps; - if (Math.abs(arg) < 1.0e-3) { - xm = -arg; - } - fln = x * xm / eps; - xm = xmin - x; - if (xm > 7.0 && fln < 15.0) { - break; - } - } - xdmy = x; - xdmln = xln; - xinc = 0.0; - if (x < xmin) { - nx = (int) x; - xinc = xmin - nx; - xdmy = x + xinc; - xdmln = Math.log(xdmy); - } + static { + Casts casts = new Casts(PsiGamma.class); + casts.arg(0).defaultError(RError.Message.NON_NUMERIC_MATH).mustBe(complexValue().not(), RError.Message.UNIMPLEMENTED_COMPLEX_FUN).mustBe(numericValue()).asDoubleVector(); + casts.arg(1).defaultError(RError.Message.NON_NUMERIC_MATH).mustBe(complexValue().not(), RError.Message.UNIMPLEMENTED_COMPLEX_FUN).mustBe(numericValue()).asDoubleVector(); + } - /* generate w(n+mm-1, x) by the asymptotic expansion */ + private final NACheck naCheck = NACheck.create(); - t = fn * xdmln; - t1 = xdmln + xdmln; - t2 = t + xdmln; - tk = Math.max(Math.abs(t), fmax2(Math.abs(t1), Math.abs(t2))); - if (tk <= elim) { /* for all but large x */ - return l10(t, tk, xdmy, xdmln, x, nn, nx, wdtol, fn, trm, trmr, xinc, mm, kode, ans); - } - } - // nz not used - // mVal.nz++; /* underflow */ - mm--; - // ans[mm] = 0.; - assert mm == 0; - ans = 0.; - if (mm == 0) { - return ans; - } - } /* end{for()} */ - nn = (int) fln + 1; - np = n + 1; - t1 = (n + 1) * xln; - t = Math.exp(-t1); - s = t; - den = x; - for (int i = 1; i <= nn; i++) { - den += 1.; - trm[i] = Math.pow(den, -np); - s += trm[i]; - } - // ans[0] = s; - ans = s; - if (n == 0 && kode == 2) { - // ans[0] = s + xln; - ans = s + xln; - } + @Specialization(guards = "binaryMapNode.isSupported(x, deriv)") + protected RAbstractDoubleVector psiGammaFast(RAbstractDoubleVector x, RAbstractDoubleVector deriv, + @Cached("createFastCached(x, deriv)") BinaryMapNode binaryMapNode) { + return (RAbstractDoubleVector) binaryMapNode.apply(x, deriv); + } - if (mm != 1) { /* generate higher derivatives, j > n */ - assert false; - // tol = wdtol / 5.0; - // for(j = 1; j < mm; j++) { - // t /= x; - // s = t; - // tols = t * tol; - // den = x; - // for(i=1; i <= nn; i++) { - // den += 1.; - // trm[i] /= den; - // s += trm[i]; - // if (trm[i] < tols) { - // break; - // } - // } - // ans[j] = s; - // } - } - return ans; + @Specialization(replaces = "psiGammaFast") + protected RAbstractDoubleVector psiGammaGeneric(RAbstractDoubleVector x, RAbstractDoubleVector deriv, + @Cached("createGeneric(x, deriv)") BinaryMapNode binaryMapNode) { + return (RAbstractDoubleVector) binaryMapNode.apply(x, deriv); + } + protected BinaryMapNode createFastCached(RAbstractDoubleVector x, RAbstractDoubleVector deriv) { + return createCached(x, deriv, false); } - private static double l10(double oldT, double oldTk, double xdmy, double xdmln, double x, double nn, double oldNx, double wdtol, double oldFn, double[] trm, double[] trmr, double xinc, - double mm, int kode, double ansOld) { - double t = oldT; - double tk = oldTk; - double nx = oldNx; - double fn = oldFn; - double ans = ansOld; - - double tss = Math.exp(-t); - double tt = 0.5 / xdmy; - double t1 = tt; - double tst = wdtol * tt; - if (nn != 0) { - t1 = tt + 1.0 / fn; - } - double rxsq = 1.0 / (xdmy * xdmy); - double ta = 0.5 * rxsq; - t = (fn + 1) * ta; - double s = t * bvalues[2]; - if (Math.abs(s) >= tst) { - tk = 2.0; - for (int k = 4; k <= 22; k++) { - t = t * ((tk + fn + 1) / (tk + 1.0)) * ((tk + fn) / (tk + 2.0)) * rxsq; - trm[k] = t * bvalues[k - 1]; - if (Math.abs(trm[k]) < tst) { - break; - } - s += trm[k]; - tk += 2.; - } - } - s = (s + t1) * tss; - if (xinc != 0.0) { - - /* backward recur from xdmy to x */ - - nx = (int) xinc; - double np = nn + 1; - if (nx > n_max) { - // nz not used - // mVal.nz = 0; - // non-zero ierr always results in generating a NaN - // mVal.ierr = 3; - return Double.NaN; - } else { - if (nn == 0) { - return l20(xdmln, xdmy, x, s, nx, kode, ans); - } - double xm = xinc - 1.0; - double fx = x + xm; - - /* this loop should not be changed. fx is accurate when x is small */ - for (int i = 1; i <= nx; i++) { - trmr[i] = Math.pow(fx, -np); - s += trmr[i]; - xm -= 1.; - fx = x + xm; - } - } - } - // ans[mm-1] = s; - assert (mm - 1) == 0; - ans = s; - if (fn == 0) { - return l30(xdmln, xdmy, x, s, kode, ans); - } + protected BinaryMapNode createGeneric(RAbstractDoubleVector x, RAbstractDoubleVector deriv) { + return createCached(x, deriv, true); + } - /* generate lower derivatives, j < n+mm-1 */ + protected BinaryMapNode createCached(RAbstractDoubleVector x, RAbstractDoubleVector deriv, boolean isGeneric) { + return BinaryMapNode.create(new PsiGammaFunction(), x, deriv, RType.Double, RType.Double, false, isGeneric); + } - for (int j = 2; j <= mm; j++) { - fn--; - tss *= xdmy; - t1 = tt; - if (fn != 0) { - t1 = tt + 1.0 / fn; - } - t = (fn + 1) * ta; - s = t * bvalues[2]; - if (Math.abs(s) >= tst) { - tk = 4 + fn; - for (int k = 4; k <= 22; k++) { - trm[k] = trm[k] * (fn + 1) / tk; - if (Math.abs(trm[k]) < tst) { - break; - } - s += trm[k]; - tk += 2.; - } - } - s = (s + t1) * tss; - if (xinc != 0.0) { - if (fn == 0) { - return l20(xdmln, xdmy, x, s, nx, kode, ans); - } - double xm = xinc - 1.0; - double fx = x + xm; - for (int i = 1; i <= nx; i++) { - trmr[i] = trmr[i] * fx; - s += trmr[i]; - xm -= 1.; - fx = x + xm; - } - } - // ans[mm - j] = s; - assert (mm - j) == 0; - ans = s; - if (fn == 0) { - return l30(xdmln, xdmy, x, s, kode, ans); - } - } - return ans; + } - } + static final class PsiGammaFunction extends BinaryMapFunctionNode { - private static double l20(double xdmln, double xdmy, double x, double oldS, double nx, int kode, double ans) { - double s = oldS; - for (int i = 1; i <= nx; i++) { - s += 1. / (x + (nx - i)); /* avoid disastrous cancellation, PR#13714 */ - } + @Override + public double applyDouble(double x, double deriv) { + return GammaFunctions.psigamma(x, deriv); + } - return l30(xdmln, xdmy, x, s, kode, ans); + @Override + public boolean mayFoldConstantTime(Class<? extends RAbstractVector> left, Class<? extends RAbstractVector> right) { + return false; } - private static double l30(double xdmln, double xdmy, double x, double s, int kode, double ansOld) { - double ans = ansOld; - if (kode != 2) { /* always */ - // ans[0] = s - xdmln; - ans = s - xdmln; - } else if (xdmy != x) { - double xq; - xq = xdmy / x; - // ans[0] = s - log(xq); - ans = s - Math.log(xq); - } - return ans; + @Override + public RAbstractVector tryFoldConstantTime(RAbstractVector left, int leftLength, RAbstractVector right, int rightLength) { + return null; } + } + } diff --git a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BasePackage.java b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BasePackage.java index 05b9d5d861..00c7379390 100644 --- a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BasePackage.java +++ b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/BasePackage.java @@ -268,6 +268,13 @@ public class BasePackage extends RBuiltinPackage { add(BaseGammaFunctions.Gamma.class, BaseGammaFunctionsFactory.GammaNodeGen::create); add(BaseGammaFunctions.Lgamma.class, BaseGammaFunctionsFactory.LgammaNodeGen::create); add(BaseGammaFunctions.TriGamma.class, BaseGammaFunctionsFactory.TriGammaNodeGen::create); + add(BaseGammaFunctions.TetraGamma.class, BaseGammaFunctionsFactory.TetraGammaNodeGen::create); + add(BaseGammaFunctions.PentaGamma.class, BaseGammaFunctionsFactory.PentaGammaNodeGen::create); + add(BaseGammaFunctions.PsiGamma.class, BaseGammaFunctionsFactory.PsiGammaNodeGen::create); + add(BaseBesselFunctions.BesselI.class, BaseBesselFunctionsFactory.BesselINodeGen::create); + add(BaseBesselFunctions.BesselJ.class, BaseBesselFunctionsFactory.BesselJNodeGen::create); + add(BaseBesselFunctions.BesselK.class, BaseBesselFunctionsFactory.BesselKNodeGen::create); + add(BaseBesselFunctions.BesselY.class, BaseBesselFunctionsFactory.BesselYNodeGen::create); add(ChooseFunctions.Choose.class, ChooseFunctionsFactory.ChooseNodeGen::create); add(ChooseFunctions.LChoose.class, ChooseFunctionsFactory.LChooseNodeGen::create); add(BetaFunctions.LBeta.class, BetaFunctionsFactory.LBetaNodeGen::create); diff --git a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/TrigExpFunctions.java b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/TrigExpFunctions.java index 28ccdc2d89..84e556a01d 100644 --- a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/TrigExpFunctions.java +++ b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/TrigExpFunctions.java @@ -48,6 +48,7 @@ import com.oracle.truffle.r.runtime.data.RDataFactory; import com.oracle.truffle.r.runtime.data.RDoubleVector; import com.oracle.truffle.r.runtime.data.model.RAbstractComplexVector; import com.oracle.truffle.r.runtime.data.model.RAbstractDoubleVector; +import com.oracle.truffle.r.runtime.nmath.RMath; import com.oracle.truffle.r.runtime.ops.BinaryArithmetic; import com.oracle.truffle.r.runtime.ops.BinaryArithmetic.Pow.CHypot; import com.oracle.truffle.r.runtime.ops.UnaryArithmetic; @@ -131,17 +132,7 @@ public class TrigExpFunctions { public static final class Sinpi extends UnaryArithmetic { @Override public double op(double op) { - double norm = op % 2d; - if (norm == 0d || norm == 1d || norm == -1d) { - return 0d; - } - if (norm == -1.5d || norm == 0.5d) { - return 1d; - } - if (norm == -0.5d || norm == 1.5d) { - return -1d; - } - return Math.sin(norm * Math.PI); + return RMath.sinpi(op); } @Override @@ -184,17 +175,7 @@ public class TrigExpFunctions { public static final class Cospi extends UnaryArithmetic { @Override public double op(double op) { - double norm = op % 2d; - if (norm == 0d) { - return 1d; - } - if (norm == -1d || norm == 1d) { - return -1d; - } - if (norm == -1.5d || norm == -0.5d || norm == 0.5d || norm == 1.5d) { - return 0d; - } - return Math.cos(norm * Math.PI); + return RMath.cospi(op); } @Override diff --git a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/Trunc.java b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/Trunc.java index c06470ea47..4ce3069d1d 100644 --- a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/Trunc.java +++ b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/Trunc.java @@ -27,6 +27,7 @@ import static com.oracle.truffle.r.runtime.builtins.RBehavior.PURE_ARITHMETIC; import static com.oracle.truffle.r.runtime.builtins.RBuiltinKind.PRIMITIVE; import com.oracle.truffle.r.runtime.builtins.RBuiltin; +import com.oracle.truffle.r.runtime.nmath.RMath; import com.oracle.truffle.r.runtime.ops.UnaryArithmetic; import com.oracle.truffle.r.runtime.ops.UnaryArithmeticFactory; @@ -37,10 +38,6 @@ public final class Trunc extends UnaryArithmetic { @Override public double op(double op) { - if (op > 0) { - return Math.floor(op); - } else { - return Math.ceil(op); - } + return RMath.trunc(op); } } diff --git a/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/InternalNode.java b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/InternalNode.java index c7d5eb0def..571ab65213 100644 --- a/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/InternalNode.java +++ b/com.oracle.truffle.r.nodes/src/com/oracle/truffle/r/nodes/builtin/InternalNode.java @@ -342,9 +342,9 @@ public abstract class InternalNode extends OperatorNode { private static final List<String> NOT_IMPLEMENTED = Arrays.asList( ".addTryHandlers", "interruptsSuspended", "restart", "max.col", "comment", "`comment<-`", "list2env", "tcrossprod", "beta", "dchisq", "pchisq", "qchisq", "dexp", "pexp", "qexp", "dgeom", "pgeom", "qgeom", "dpois", "ppois", "qpois", "dt", "pt", "qt", "dsignrank", "psignrank", - "qsignrank", "besselJ", "besselY", "psigamma", "dbeta", "pbeta", "qbeta", "dbinom", "pbinom", "qbinom", "dcauchy", "pcauchy", "qcauchy", "df", "pf", "qf", "dgamma", "pgamma", + "qsignrank", "dbeta", "pbeta", "qbeta", "dbinom", "pbinom", "qbinom", "dcauchy", "pcauchy", "qcauchy", "df", "pf", "qf", "dgamma", "pgamma", "qgamma", "dlnorm", "plnorm", "qlnorm", "dlogis", "plogis", "qlogis", "dnbinom", "pnbinom", "qnbinom", "dnorm", "pnorm", "qnorm", "dunif", "punif", "qunif", "dweibull", "pweibull", - "qweibull", "dnchisq", "pnchisq", "qnchisq", "dnt", "pnt", "qnt", "dwilcox", "pwilcox", "qwilcox", "besselI", "besselK", "dnbinom_mu", "pnbinom_mu", "qnbinom_mu", "dhyper", + "qweibull", "dnchisq", "pnchisq", "qnchisq", "dnt", "pnt", "qnt", "dwilcox", "pwilcox", "qwilcox", "dnbinom_mu", "pnbinom_mu", "qnbinom_mu", "dhyper", "phyper", "qhyper", "dnbeta", "pnbeta", "qnbeta", "dnf", "pnf", "qnf", "dtukey", "ptukey", "qtukey", "rchisq", "rexp", "rgeom", "rpois", "rt", "rsignrank", "rbeta", "rbinom", "rcauchy", "rf", "rgamma", "rlnorm", "rlogis", "rnbinom", "rnbinom_mu", "rnchisq", "rnorm", "runif", "rweibull", "rwilcox", "rhyper", "grepRaw", "regexec", "adist", "aregexec", "chartr", "strtrim", "eapply", "machine", "save", "dump", "prmatrix", "gcinfo", diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/RError.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/RError.java index 2f01fd7d28..4e358e1714 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/RError.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/RError.java @@ -963,7 +963,11 @@ public final class RError extends RuntimeException implements TruffleException { MUST_BE_COMPLEX_MATRIX("'%s' must be a complex matrix"), INVALID_FORMAL_ARG_LIST("invalid formal argument list for \"%s\""), SINGULAR_BACKSOLVE("singular matrix in 'backsolve'. First zero in diagonal [%d]"), - EOF_AFTER_BACKSLASH("\\ followed by EOF"); + EOF_AFTER_BACKSLASH("\\ followed by EOF"), + DERIV_OVER_N_MAX("deriv = %d > %d (= n_max)"), + BESSEL_ARG_RANGE("bessel_%s(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?"), + BESSEL_PRECISION_LOST("bessel_%s(%g,nu=%g): precision lost in result"), + BESSEL_NU_TOO_LARGE("bessel%s(x, nu): nu=%g too large for bessel_%s() algorithm"); public final String message; final boolean hasArgs; diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/BesselFunctions.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/BesselFunctions.java new file mode 100644 index 0000000000..9423d8c77a --- /dev/null +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/BesselFunctions.java @@ -0,0 +1,2257 @@ +/* + * This material is distributed under the GNU General Public License + * Version 2. You may review the terms of this license at + * http://www.gnu.org/licenses/gpl-2.0.html + * + * Copyright (C) 1998-2014 Ross Ihaka + * Copyright (c) 2002-3, The R Core Team + * Copyright (c) 2018, Oracle and/or its affiliates + * + * All rights reserved. + */ +package com.oracle.truffle.r.runtime.nmath; + +import static com.oracle.truffle.r.runtime.nmath.Arithmetic.pow; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_EPSILON; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MAX; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MIN; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_1_PI; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_PI; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_PI_2; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_SQRT_2dPI; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.ML_NAN; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.ML_NEGINF; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.ML_POSINF; +import static com.oracle.truffle.r.runtime.nmath.RMath.cospi; +import static com.oracle.truffle.r.runtime.nmath.RMath.fmax2; +import static com.oracle.truffle.r.runtime.nmath.RMath.sinpi; +import static com.oracle.truffle.r.runtime.nmath.RMath.trunc; +import static com.oracle.truffle.r.runtime.nmath.TOMS708.fabs; + +import com.oracle.truffle.r.runtime.RError; + +// Checkstyle: stop +// @formatter:off +public class BesselFunctions { + + // GNUR from bessel.h + + private static final int nsig_BESS = 16; + private static final double ensig_BESS = 1e16; + private static final double rtnsig_BESS = 1e-4; + private static final double enmten_BESS = 8.9e-308; + private static final double enten_BESS = 1e308; + + private static final double exparg_BESS = 709.; + private static final double xlrg_BESS_IJ = 1e5; + private static final double xlrg_BESS_Y = 1e8; + private static final double thresh_BESS_Y = 16.; + + private static final double xmax_BESS_K = 705.342; /* maximal x for UNscaled answer */ + + /* sqrt(DBL_MIN) = 1.491668e-154 */ + private static final double sqxmin_BESS_K = 1.49e-154; + + /* + * x < eps_sinc <==> sin(x)/x == 1 (particularly "==>"); Linux (around 2001-02) gives + * 2.14946906753213e-08 Solaris 2.5.1 gives 2.14911933289084e-08 + */ + private static final double M_eps_sinc = 2.149e-8; + + // GNUR from bessel_i.c + + /* .Internal(besselI(*)) : */ + public static double bessel_i(double x, double alpha, double expo) { + /* NaNs propagated correctly */ + if (Double.isNaN(x) || Double.isNaN(alpha)) + return x + alpha; + if (x < 0) { + RMathError.error(RMathError.MLError.RANGE, "bessel_i"); + return ML_NAN; + } + int ize = (int) expo; + double na = Math.floor(alpha); + if (alpha < 0) { + /* + * Using Abramowitz & Stegun 9.6.2 & 9.6.6 this may not be quite optimal (CPU and + * accuracy wise) + */ + return (bessel_i(x, -alpha, expo) + + ((alpha == na) ? /* sin(pi * alpha) = 0 */ 0 : bessel_k(x, -alpha, expo) * + ((ize == 1) ? 2. : 2. * Math.exp(-2. * x)) / M_PI * sinpi(-alpha))); + } + int nb = 1 + (int) na;/* nb-1 <= alpha < nb */ + alpha -= (double) (nb - 1); + double[] bi = new double[nb]; + int ncalc = I_bessel(x, alpha, nb, ize, bi); + if (ncalc != nb) {/* error input */ + if (ncalc < 0) + RMathError.warning(RError.Message.BESSEL_ARG_RANGE, "i", x, ncalc, nb, alpha); + else + RMathError.warning(RError.Message.BESSEL_PRECISION_LOST, "i", x, alpha + (double) nb - 1); + } + return bi[nb - 1]; + } + + /* + * modified version of bessel_i that accepts a work array instead of allocating one. + */ + public static double bessel_i_ex(double x, double alpha, double expo, double[] bi) { + /* NaNs propagated correctly */ + if (Double.isNaN(x) || Double.isNaN(alpha)) + return x + alpha; + if (x < 0) { + RMathError.error(RMathError.MLError.RANGE, "bessel_i"); + return ML_NAN; + } + int ize = (int) expo; + double na = Math.floor(alpha); + if (alpha < 0) { + /* + * Using Abramowitz & Stegun 9.6.2 & 9.6.6 this may not be quite optimal (CPU and + * accuracy wise) + */ + return (bessel_i_ex(x, -alpha, expo, bi) + + ((alpha == na) ? 0 : bessel_k_ex(x, -alpha, expo, bi) * + ((ize == 1) ? 2. : 2. * Math.exp(-2. * x)) / M_PI * sinpi(-alpha))); + } + int nb = 1 + (int) na;/* nb-1 <= alpha < nb */ + alpha -= (double) (nb - 1); + int ncalc = I_bessel(x, alpha, nb, ize, bi); + if (ncalc != nb) {/* error input */ + if (ncalc < 0) + RMathError.warning(RError.Message.BESSEL_ARG_RANGE, "i", x, ncalc, nb, alpha); + else + RMathError.warning(RError.Message.BESSEL_PRECISION_LOST, "i", x, alpha + (double) nb - 1); + } + x = bi[nb - 1]; + return x; + } + + /* + * ------------------------------------------------------------------- + * + * This routine calculates Bessel functions I_(N+ALPHA) (X) for non-negative argument X, and + * non-negative order N+ALPHA, with or without exponential scaling. + * + * + * Explanation of variables in the calling sequence + * + * X - Non-negative argument for which I's or exponentially scaled I's (I*EXP(-X)) are to be + * calculated. If I's are to be calculated, X must be less than exparg_BESS (IZE=1) or + * xlrg_BESS_IJ (IZE=2), (see bessel.h). ALPHA - Fractional part of order for which I's or + * exponentially scaled I's (I*EXP(-X)) are to be calculated. 0 <= ALPHA < 1.0. NB - Number of + * functions to be calculated, NB > 0. The first function calculated is of order ALPHA, and the + * last is of order (NB - 1 + ALPHA). IZE - Type. IZE = 1 if unscaled I's are to be calculated, + * = 2 if exponentially scaled I's are to be calculated. BI - Output vector of length NB. If the + * routine terminates normally (NCALC=NB), the vector BI contains the functions I(ALPHA,X) + * through I(NB-1+ALPHA,X), or the corresponding exponentially scaled functions. NCALC - Output + * variable indicating possible errors. Before using the vector BI, the user should check that + * NCALC=NB, i.e., all orders have been calculated to the desired accuracy. See error returns + * below. + ******************************************************************* + ******************************************************************* + * + * + * + * Error returns + * + * In case of an error, NCALC != NB, and not all I's are calculated to the desired accuracy. + * + * NCALC < 0: An argument is out of range. For example, NB <= 0, IZE is not 1 or 2, or IZE=1 and + * ABS(X) >= EXPARG_BESS. In this case, the BI-vector is not calculated, and NCALC is set to + * MIN0(NB,0)-1 so that NCALC != NB. + * + * NB > NCALC > 0: Not all requested function values could be calculated accurately. This + * usually occurs because NB is much larger than ABS(X). In this case, BI[N] is calculated to + * the desired accuracy for N <= NCALC, but precision is lost for NCALC < N <= NB. If BI[N] does + * not vanish for N > NCALC (because it is too small to be represented), and BI[N]/BI[NCALC] = + * 10**(-K), then only the first NSIG-K significant figures of BI[N] can be trusted. + * + * + * Intrinsic functions required are: + * + * DBLE, EXP, gamma_cody, INT, MAX, MIN, REAL, SQRT + * + * + * Acknowledgement + * + * This program is based on a program written by David J. Sookne (2) that computes values of the + * Bessel functions J or I of float argument and long order. Modifications include the + * restriction of the computation to the I Bessel function of non-negative float argument, the + * extension of the computation to arbitrary positive order, the inclusion of optional + * exponential scaling, and the elimination of most underflow. An earlier version was published + * in (3). + * + * References: "A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne, D. J., + * Math. Comp. 26, 1972, pp 941-947. + * + * "Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS Jour. of Res. B. + * 77B, 1973, pp 125-132. + * + * "ALGORITHM 597, Sequence of Modified Bessel Functions of the First Kind," Cody, W. J., Trans. + * Math. Soft., 1983, pp. 242-245. + * + * Latest modification: May 30, 1989 + * + * Modified by: W. J. Cody and L. Stoltz Applied Mathematics Division Argonne National + * Laboratory Argonne, IL 60439 + */ + + /*------------------------------------------------------------------- + Mathematical constants + -------------------------------------------------------------------*/ + private static final double const__ = 1.585; + + private static int I_bessel(double x, double alpha, int nb, int ize, double[] bi) { + /* Local variables */ + int nend, intx, nbmx, k, l, n, nstart; + double pold, test, p, em, en, empal, emp2al, halfx, + aa, bb, cc, psave, plast, tover, psavel, sum, nu, twonu; + + /* Parameter adjustments */ + nu = alpha; + twonu = nu + nu; + + /*------------------------------------------------------------------- + Check for X, NB, OR IZE out of range. + ------------------------------------------------------------------- */ + int ncalc; + if (nb > 0 && x >= 0. && (0. <= nu && nu < 1.) && + (1 <= ize && ize <= 2)) { + + ncalc = nb; + if (ize == 1 && x > exparg_BESS) { + for (k = 1; k <= nb; k++) + bi[k - 1] = ML_POSINF; /* the limit *is* = Inf */ + return ncalc; + } + if (ize == 2 && x > xlrg_BESS_IJ) { + for (k = 1; k <= nb; k++) + bi[k - 1] = 0.; /* The limit exp(-x) * I_nu(x) --> 0 : */ + return ncalc; + } + intx = (int) (x);/* fine, since *x <= xlrg_BESS_IJ <<< LONG_MAX */ + if (x >= rtnsig_BESS) { /* "non-small" x ( >= 1e-4 ) */ + /* + * ------------------------------------------------------------------- Initialize + * the forward sweep, the P-sequence of Olver + * ------------------------------------------------------------------- + */ + nbmx = nb - intx; + n = intx + 1; + en = (double) (n + n) + twonu; + plast = 1.; + p = en / x; + /* + * ------------------------------------------------ Calculate general significance + * test ------------------------------------------------ + */ + test = ensig_BESS + ensig_BESS; + if (intx << 1 > nsig_BESS * 5) { + test = Math.sqrt(test * p); + } else { + test /= Arithmetic.powDi(const__, intx); + } + L120W: while (true) { + if (nbmx >= 3) { + /* + * -------------------------------------------------- Calculate P-sequence + * until N = NB-1 Check for possible overflow. + * ------------------------------------------------ + */ + tover = enten_BESS / ensig_BESS; + nstart = intx + 2; + nend = nb - 1; + for (k = nstart; k <= nend; ++k) { + n = k; + en += 2.; + pold = plast; + plast = p; + p = en * plast / x + pold; + if (p > tover) { + /* + * ------------------------------------------------ To avoid + * overflow, divide P-sequence by TOVER. Calculate P-sequence until + * ABS(P) > 1. ---------------------------------------------- + */ + tover = enten_BESS; + p /= tover; + plast /= tover; + psave = p; + psavel = plast; + nstart = n + 1; + do { + ++n; + en += 2.; + pold = plast; + plast = p; + p = en * plast / x + pold; + } while (p <= 1.); + + bb = en / x; + /* + * ------------------------------------------------ Calculate + * backward test, and find NCALC, the highest N such that the test + * is passed. ------------------------------------------------ + */ + test = pold * plast / ensig_BESS; + test *= .5 - .5 / (bb * bb); + p = plast * tover; + --n; + en -= 2.; + nend = Math.min(nb, n); + L90W: while (true) { + for (l = nstart; l <= nend; ++l) { + ncalc = l; + pold = psavel; + psavel = psave; + psave = en * psavel / x + pold; + if (psave * psavel > test) { + break L90W; + } + } + ncalc = nend + 1; + break L90W; + } + L90: --(ncalc); + break L120W; + } + } + n = nend; + en = (double) (n + n) + twonu; + /*--------------------------------------------------- + Calculate special significance test for NBMX > 2. + --------------------------------------------------- */ + test = fmax2(test, Math.sqrt(plast * ensig_BESS) * Math.sqrt(p + p)); + } + /* + * -------------------------------------------------------- Calculate P-sequence + * until significance test passed. + * -------------------------------------------------------- + */ + do { + ++n; + en += 2.; + pold = plast; + plast = p; + p = en * plast / x + pold; + } while (p < test); + break; // of extra while (true) + } + + L120: + /* + * ------------------------------------------------------------------- Initialize + * the backward recursion and the normalization sum. + * ------------------------------------------------------------------- + */ + ++n; + en += 2.; + bb = 0.; + aa = 1. / p; + em = (double) n - 1.; + empal = em + nu; + emp2al = em - 1. + twonu; + sum = aa * empal * emp2al / em; + nend = n - nb; + L230W: while (true) { + L220W: while (true) { + if (nend < 0) { + /* + * ----------------------------------------------------- N < NB, so + * store BI[N] and set higher orders to 0.. + * ----------------------------------------------------- + */ + bi[n - 1] = aa; + nend = -nend; + for (l = 1; l <= nend; ++l) { + bi[n] = 0.; + } + } else { + if (nend > 0) { + /* + * ----------------------------------------------------- Recur + * backward via difference equation, calculating (but not storing) + * BI[N], until N = NB. + * --------------------------------------------------- + */ + + for (l = 1; l <= nend; ++l) { + --n; + en -= 2.; + cc = bb; + bb = aa; + /* + * for x ~= 1500, sum would overflow to 'inf' here, and the + * final bi[] /= sum would give 0 wrongly; RE-normalize (aa, + * sum) here -- no need to undo + */ + if (nend > 100 && aa > 1e200) { + /* multiply by 2^-900 = 1.18e-271 */ + cc = Math.scalb(cc, -900); + bb = Math.scalb(bb, -900); + sum = Math.scalb(sum, -900); + } + aa = en * bb / x + cc; + em -= 1.; + emp2al -= 1.; + if (n == 1) { + break; + } + if (n == 2) { + emp2al = 1.; + } + empal -= 1.; + sum = (sum + aa * empal) * emp2al / em; + } + } + /* + * --------------------------------------------------- Store BI[NB] + * --------------------------------------------------- + */ + bi[n - 1] = aa; + if (nb <= 1) { + sum = sum + sum + aa; + break L230W; + } + /* + * ------------------------------------------------- Calculate and Store + * BI[NB-1] ------------------------------------------------- + */ + --n; + en -= 2.; + bi[n - 1] = en * aa / x + bb; + if (n == 1) { + break L220W; + } + em -= 1.; + if (n == 2) + emp2al = 1.; + else + emp2al -= 1.; + + empal -= 1.; + sum = (sum + bi[n - 1] * empal) * emp2al / em; + } + nend = n - 2; + if (nend > 0) { + /* + * -------------------------------------------- Calculate via difference + * equation and store BI[N], until N = 2. + * ------------------------------------------ + */ + for (l = 1; l <= nend; ++l) { + --n; + en -= 2.; + bi[n - 1] = en * bi[n] / x + bi[n + 1]; + em -= 1.; + if (n == 2) + emp2al = 1.; + else + emp2al -= 1.; + empal -= 1.; + sum = (sum + bi[n - 1] * empal) * emp2al / em; + } + } + /* + * ---------------------------------------------- Calculate BI[1] + * -------------------------------------------- + */ + bi[0] = 2. * empal * bi[1] / x + bi[2]; + break; // of extra while (true) + } + L220: sum = sum + sum + bi[0]; + break; // of extra while (true) + } + + L230: + /* + * --------------------------------------------------------- Normalize. Divide all + * BI[N] by sum. --------------------------------------------------------- + */ + if (nu != 0.) + sum *= (gammaCody(1. + nu) * pow(x * .5, -nu)); + if (ize == 1) + sum *= Math.exp(-(x)); + aa = enmten_BESS; + if (sum > 1.) + aa *= sum; + for (n = 1; n <= nb; ++n) { + if (bi[n - 1] < aa) + bi[n - 1] = 0.; + else + bi[n - 1] /= sum; + } + return ncalc; + } else { /* small x < 1e-4 */ + /* + * ----------------------------------------------------------- Two-term ascending + * series for small X. ----------------------------------------------------------- + */ + aa = 1.; + empal = 1. + nu; + /* No need to check for underflow */ + halfx = .5 * x; + if (nu != 0.) + aa = pow(halfx, nu) / gammaCody(empal); + if (ize == 2) + aa *= Math.exp(-(x)); + bb = halfx * halfx; + bi[0] = aa + aa * bb / empal; + if (x != 0. && bi[0] == 0.) + ncalc = 0; + if (nb > 1) { + if (x == 0.) { + for (n = 2; n <= nb; ++n) + bi[n - 1] = 0.; + } else { + /* + * ------------------------------------------------- Calculate higher-order + * functions. ------------------------------------------------- + */ + cc = halfx; + tover = (enmten_BESS + enmten_BESS) / x; + if (bb != 0.) + tover = enmten_BESS / bb; + for (n = 2; n <= nb; ++n) { + aa /= empal; + empal += 1.; + aa *= cc; + if (aa <= tover * empal) + bi[n - 1] = aa = 0.; + else + bi[n - 1] = aa + aa * bb / empal; + if (bi[n - 1] == 0. && ncalc > n) + ncalc = n - 1; + } + } + } + } + } else { /* argument out of range */ + ncalc = Math.min(nb, 0) - 1; + } + return ncalc; + } + + // GNUR from bessel_j.c + + // unused now from R + public static double bessel_j(double x, double alpha) { + /* NaNs propagated correctly */ + if (Double.isNaN(x) || Double.isNaN(alpha)) + return x + alpha; + if (x < 0) { + RMathError.error(RMathError.MLError.RANGE, "bessel_j"); + return ML_NAN; + } + double na = Math.floor(alpha); + if (alpha < 0) { + /* + * Using Abramowitz & Stegun 9.1.2 this may not be quite optimal (CPU and accuracy wise) + */ + return (((alpha - na == 0.5) ? 0 : bessel_j(x, -alpha) * cospi(alpha)) + + ((alpha == na) ? 0 : bessel_y(x, -alpha) * sinpi(alpha))); + } else if (alpha > 1e7) { + RMathError.warning(RError.Message.BESSEL_NU_TOO_LARGE, "J", alpha, "j"); + return ML_NAN; + } + int nb = 1 + (int) na; /* nb-1 <= alpha < nb */ + alpha -= (double) (nb - 1); // ==> alpha' in [0, 1) + double[] bj = new double[nb]; + int ncalc = J_bessel(x, alpha, nb, bj); + if (ncalc != nb) {/* error input */ + if (ncalc < 0) + RMathError.warning(RError.Message.BESSEL_ARG_RANGE, "j", x, ncalc, nb, alpha); + else + RMathError.warning(RError.Message.BESSEL_PRECISION_LOST, "j", x, alpha + (double) nb - 1); + } + x = bj[nb - 1]; + return x; + } + + /* + * Called from R: modified version of bessel_j(), accepting a work array instead of allocating + * one. + */ + public static double bessel_j_ex(double x, double alpha, double[] bj) { + /* NaNs propagated correctly */ + if (Double.isNaN(x) || Double.isNaN(alpha)) + return x + alpha; + if (x < 0) { + RMathError.error(RMathError.MLError.RANGE, "bessel_j"); + return ML_NAN; + } + double na = Math.floor(alpha); + if (alpha < 0) { + /* + * Using Abramowitz & Stegun 9.1.2 this may not be quite optimal (CPU and accuracy wise) + */ + return (((alpha - na == 0.5) ? 0 : bessel_j_ex(x, -alpha, bj) * cospi(alpha)) + + ((alpha == na) ? 0 : bessel_y_ex(x, -alpha, bj) * sinpi(alpha))); + } else if (alpha > 1e7) { + RMathError.warning(RError.Message.BESSEL_NU_TOO_LARGE, "J", alpha, "j"); + return ML_NAN; + } + int nb = 1 + (int) na; /* nb-1 <= alpha < nb */ + alpha -= (double) (nb - 1); // ==> alpha' in [0, 1) + int ncalc = J_bessel(x, alpha, nb, bj); + if (ncalc != nb) {/* error input */ + if (ncalc < 0) + RMathError.warning(RError.Message.BESSEL_ARG_RANGE, "j", x, ncalc, nb, alpha); + else + RMathError.warning(RError.Message.BESSEL_PRECISION_LOST, "j", x, alpha + (double) nb - 1); + } + x = bj[nb - 1]; + return x; + } + + /* + * Calculates Bessel functions J_{n+alpha} (x) for non-negative argument x, and non-negative + * order n+alpha, n = 0,1,..,nb-1. + * + * Explanation of variables in the calling sequence. + * + * X - Non-negative argument for which J's are to be calculated. ALPHA - Fractional part of + * order for which J's are to be calculated. 0 <= ALPHA < 1. NB - Number of functions to be + * calculated, NB >= 1. The first function calculated is of order ALPHA, and the last is of + * order (NB - 1 + ALPHA). B - Output vector of length NB. If RJBESL terminates normally + * (NCALC=NB), the vector B contains the functions J/ALPHA/(X) through J/NB-1+ALPHA/(X). NCALC - + * Output variable indicating possible errors. Before using the vector B, the user should check + * that NCALC=NB, i.e., all orders have been calculated to the desired accuracy. See the + * following + **************************************************************** + * + * + * Error return codes + * + * In case of an error, NCALC != NB, and not all J's are calculated to the desired accuracy. + * + * NCALC < 0: An argument is out of range. For example, NBES <= 0, ALPHA < 0 or > 1, or X is too + * large. In this case, b[1] is set to zero, the remainder of the B-vector is not calculated, + * and NCALC is set to MIN(NB,0)-1 so that NCALC != NB. + * + * NB > NCALC > 0: Not all requested function values could be calculated accurately. This + * usually occurs because NB is much larger than ABS(X). In this case, b[N] is calculated to the + * desired accuracy for N <= NCALC, but precision is lost for NCALC < N <= NB. If b[N] does not + * vanish for N > NCALC (because it is too small to be represented), and b[N]/b[NCALC] = + * 10^(-K), then only the first NSIG - K significant figures of b[N] can be trusted. + * + * + * Acknowledgement + * + * This program is based on a program written by David J. Sookne (2) that computes values of the + * Bessel functions J or I of float argument and long order. Modifications include the + * restriction of the computation to the J Bessel function of non-negative float argument, the + * extension of the computation to arbitrary positive order, and the elimination of most + * underflow. + * + * References: + * + * Olver, F.W.J., and Sookne, D.J. (1972) "A Note on Backward Recurrence Algorithms"; Math. + * Comp. 26, 941-947. + * + * Sookne, D.J. (1973) "Bessel Functions of Real Argument and Integer Order"; NBS Jour. of Res. + * B. 77B, 125-132. + * + * Latest modification: March 19, 1990 + * + * Author: W. J. Cody Applied Mathematics Division Argonne National Laboratory Argonne, IL 60439 + ******************************************************************* + */ + + /* + * --------------------------------------------------------------------- Mathematical constants + * + * PI2 = 2 / PI TWOPI1 = first few significant digits of 2 * PI TWOPI2 = (2*PI - TWOPI1) to + * working precision, i.e., TWOPI1 + TWOPI2 = 2 * PI to extra precision. + * --------------------------------------------------------------------- + */ + private static final double pi2 = .636619772367581343075535; + private static final double twopi1 = 6.28125; + private static final double twopi2 = .001935307179586476925286767; + + /*--------------------------------------------------------------------- + * Factorial(N) + *--------------------------------------------------------------------- */ + private static final double[] fact = new double[]{1., 1., 2., 6., 24., 120., 720., 5040., 40320., + 362880., 3628800., 39916800., 479001600., 6227020800., 87178291200., + 1.307674368e12, 2.0922789888e13, 3.55687428096e14, 6.402373705728e15, + 1.21645100408832e17, 2.43290200817664e18, 5.109094217170944e19, + 1.12400072777760768e21, 2.585201673888497664e22, + 6.2044840173323943936e23}; + + private static int J_bessel(double x, double alpha, int nb, double[] b) { + /* Local variables */ + int nend, intx, nbmx, i, j, k, l, m, n, nstart; + + double nu, twonu, capp, capq, pold, vcos, test, vsin; + double p, s, t, z, alpem, halfx, aa, bb, cc, psave, plast; + double tover, t1, alp2em, em, en, xc, xk, xm, psavel, gnu, xin, sum; + + /* Parameter adjustment */ + nu = alpha; + twonu = nu + nu; + + int ncalc; + /*------------------------------------------------------------------- + Check for out of range arguments. + -------------------------------------------------------------------*/ + if (nb > 0 && x >= 0. && 0. <= nu && nu < 1.) { + + ncalc = nb; + if (x > xlrg_BESS_IJ) { + RMathError.error(RMathError.MLError.RANGE, "J_bessel"); + /* + * indeed, the limit is 0, but the cutoff happens too early + */ + for (i = 1; i <= nb; i++) + b[i - 1] = 0.; /* was ML_POSINF (really nonsense) */ + return ncalc; + } + intx = (int) (x); + /* Initialize result array to zero. */ + for (i = 1; i <= nb; ++i) + b[i - 1] = 0.; + + /* + * =================================================================== Branch into 3 + * cases : 1) use 2-term ascending series for small X 2) use asymptotic form for large X + * when NB is not too large 3) use recursion otherwise + * =================================================================== + */ + + if (x < rtnsig_BESS) { + /* + * --------------------------------------------------------------- Two-term + * ascending series for small X. + * --------------------------------------------------------------- + */ + alpem = 1. + nu; + + halfx = (x > enmten_BESS) ? .5 * x : 0.; + aa = (nu != 0.) ? pow(halfx, nu) / (nu * gammaCody(nu)) : 1.; + bb = (x + 1. > 1.) ? -halfx * halfx : 0.; + b[0] = aa + aa * bb / alpem; + if (x != 0. && b[0] == 0.) + ncalc = 0; + + if (nb != 1) { + if (x <= 0.) { + for (n = 2; n <= nb; ++n) + b[n - 1] = 0.; + } else { + /* + * ---------------------------------------------- Calculate higher order + * functions. ---------------------------------------------- + */ + if (bb == 0.) + tover = (enmten_BESS + enmten_BESS) / x; + else + tover = enmten_BESS / bb; + cc = halfx; + for (n = 2; n <= nb; ++n) { + aa /= alpem; + alpem += 1.; + aa *= cc; + if (aa <= tover * alpem) + aa = 0.; + + b[n - 1] = aa + aa * bb / alpem; + if (b[n - 1] == 0. && ncalc > n) + ncalc = n - 1; + } + } + } + } else if (x > 25. && nb <= intx + 1) { + /* + * ------------------------------------------------------------ Asymptotic series + * for X > 25 (and not too large nb) + * ------------------------------------------------------------ + */ + xc = Math.sqrt(pi2 / x); + xin = 1 / (64 * x * x); + if (x >= 130.) + m = 4; + else if (x >= 35.) + m = 8; + else + m = 11; + xm = 4. * (double) m; + /* + * ------------------------------------------------ Argument reduction for SIN and + * COS routines. ------------------------------------------------ + */ + t = trunc(x / (twopi1 + twopi2) + .5); + z = (x - t * twopi1) - t * twopi2 - (nu + .5) / pi2; + vsin = Math.sin(z); + vcos = Math.cos(z); + gnu = twonu; + for (i = 1; i <= 2; ++i) { + s = (xm - 1. - gnu) * (xm - 1. + gnu) * xin * .5; + t = (gnu - (xm - 3.)) * (gnu + (xm - 3.)); + t1 = (gnu - (xm + 1.)) * (gnu + (xm + 1.)); + k = m + m; + capp = s * t / fact[k]; + capq = s * t1 / fact[k + 1]; + xk = xm; + for (; k >= 4; k -= 2) {/* k + 2(j-2) == 2m */ + xk -= 4.; + s = (xk - 1. - gnu) * (xk - 1. + gnu); + t1 = t; + t = (gnu - (xk - 3.)) * (gnu + (xk - 3.)); + capp = (capp + 1. / fact[k - 2]) * s * t * xin; + capq = (capq + 1. / fact[k - 1]) * s * t1 * xin; + + } + capp += 1.; + capq = (capq + 1.) * (gnu * gnu - 1.) * (.125 / x); + b[i - 1] = xc * (capp * vcos - capq * vsin); + if (nb == 1) + return ncalc; + + /* vsin <--> vcos */ t = vsin; + vsin = -vcos; + vcos = t; + gnu += 2.; + } + /* + * ----------------------------------------------- If NB > 2, compute J(X,ORDER+I) + * for I = 2, NB-1 ----------------------------------------------- + */ + if (nb > 2) + for (gnu = twonu + 2., j = 3; j <= nb; j++, gnu += 2.) + b[j - 1] = gnu * b[j - 2] / x - b[j - 3]; + } else { + /* + * rtnsig_BESS <= x && ( x <= 25 || intx+1 < *nb ) : + * -------------------------------------------------------- Use recurrence to + * generate results. First initialize the calculation of P*S. + * -------------------------------------------------------- + */ + nbmx = nb - intx; + n = intx + 1; + en = (double) (n + n) + twonu; + plast = 1.; + p = en / x; + /* + * --------------------------------------------------- Calculate general + * significance test. --------------------------------------------------- + */ + test = ensig_BESS + ensig_BESS; + L190W: while (true) { + if (nbmx >= 3) { + /* + * ------------------------------------------------------------ Calculate + * P*S until N = NB-1. Check for possible overflow. + * ---------------------------------------------------------- + */ + tover = enten_BESS / ensig_BESS; + nstart = intx + 2; + nend = nb - 1; + en = (double) (nstart + nstart) - 2. + twonu; + for (k = nstart; k <= nend; ++k) { + n = k; + en += 2.; + pold = plast; + plast = p; + p = en * plast / x - pold; + if (p > tover) { + /* + * ------------------------------------------- To avoid overflow, + * divide P*S by TOVER. Calculate P*S until ABS(P) > 1. + * ------------------------------------------- + */ + tover = enten_BESS; + p /= tover; + plast /= tover; + psave = p; + psavel = plast; + nstart = n + 1; + do { + ++n; + en += 2.; + pold = plast; + plast = p; + p = en * plast / x - pold; + } while (p <= 1.); + + bb = en / x; + /* + * ----------------------------------------------- Calculate + * backward test and find NCALC, the highest N such that the test is + * passed. ----------------------------------------------- + */ + test = pold * plast * (.5 - .5 / (bb * bb)); + test /= ensig_BESS; + p = plast * tover; + --n; + en -= 2.; + nend = Math.min(nb, n); + for (l = nstart; l <= nend; ++l) { + pold = psavel; + psavel = psave; + psave = en * psavel / x - pold; + if (psave * psavel > test) { + ncalc = l - 1; + break L190W; + } + } + ncalc = nend; + break L190W; + } + } + n = nend; + en = (double) (n + n) + twonu; + /* + * ----------------------------------------------------- Calculate special + * significance test for NBMX > 2. + * ----------------------------------------------------- + */ + test = fmax2(test, Math.sqrt(plast * ensig_BESS) * Math.sqrt(p + p)); + } + /* + * ------------------------------------------------ Calculate P*S until + * significance test passes. + */ + do { + ++n; + en += 2.; + pold = plast; + plast = p; + p = en * plast / x - pold; + } while (p < test); + break; // of extra while (true) + } + + L190: + /*--------------------------------------------------------------- + Initialize the backward recursion and the normalization sum. + --------------------------------------------------------------- */ + ++n; + en += 2.; + bb = 0.; + aa = 1. / p; + m = n / 2; + em = (double) m; + m = (n << 1) - (m << 2);/* + * = 2 n - 4 (n/2) = 0 for even, 2 for odd n + */ + if (m == 0) + sum = 0.; + else { + alpem = em - 1. + nu; + alp2em = em + em + nu; + sum = aa * alpem * alp2em / em; + } + nend = n - nb; + /* if (nend > 0) */ + /* + * -------------------------------------------------------- Recur backward via + * difference equation, calculating (but not storing) b[N], until N = NB. + * -------------------------------------------------------- + */ + for (l = 1; l <= nend; ++l) { + --n; + en -= 2.; + cc = bb; + bb = aa; + aa = en * bb / x - cc; + m = (m != 0) ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */ + if (m != 0) { + em -= 1.; + alp2em = em + em + nu; + if (n == 1) + break; + + alpem = em - 1. + nu; + if (alpem == 0.) + alpem = 1.; + sum = (sum + aa * alp2em) * alpem / em; + } + } + /*-------------------------------------------------- + Store b[NB]. + --------------------------------------------------*/ + b[n - 1] = aa; + L250W: while (true) { + L240W: while (true) { + if (nend >= 0) { + if (nb <= 1) { + if (nu + 1. == 1.) + alp2em = 1.; + else + alp2em = nu; + sum += b[0] * alp2em; + break L250W; + } else {/*-- nb >= 2 : --------------------------- + Calculate and store b[NB-1]. + ----------------------------------------*/ + --n; + en -= 2.; + b[n - 1] = en * aa / x - bb; + if (n == 1) + break L240W; + + m = (m != 0) ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */ + if (m != 0) { + em -= 1.; + alp2em = em + em + nu; + alpem = em - 1. + nu; + if (alpem == 0.) + alpem = 1.; + sum = (sum + b[n - 1] * alp2em) * alpem / em; + } + } + } + + /* if (n - 2 != 0) */ + /* + * -------------------------------------------------------- Calculate via + * difference equation and store b[N], until N = 2. + * -------------------------------------------------------- + */ + for (n = n - 1; n >= 2; n--) { + en -= 2.; + b[n - 1] = en * b[n] / x - b[n + 1]; + m = (m != 0) ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */ + if (m != 0) { + em -= 1.; + alp2em = em + em + nu; + alpem = em - 1. + nu; + if (alpem == 0.) + alpem = 1.; + sum = (sum + b[n - 1] * alp2em) * alpem / em; + } + } + /* + * --------------------------------------- Calculate b[1]. + * ----------------------------------------- + */ + b[0] = 2. * (nu + 1.) * b[1] / x - b[2]; + break; // of extra while (true) + } + + L240: em -= 1.; + alp2em = em + em + nu; + if (alp2em == 0.) + alp2em = 1.; + sum += b[0] * alp2em; + break; // of extra while (true) + } + + L250: + /* + * --------------------------------------------------- Normalize. Divide all b[N] by + * sum. --------------------------------------------------- + */ + /* if (nu + 1. != 1.) poor test */ + if (fabs(nu) > 1e-15) + sum *= (gammaCody(nu) * pow(.5 * x, -nu)); + + aa = enmten_BESS; + if (sum > 1.) + aa *= sum; + for (n = 1; n <= nb; ++n) { + if (fabs(b[n - 1]) < aa) + b[n - 1] = 0.; + else + b[n - 1] /= sum; + } + } + } else { + /* Error return -- X, NB, or ALPHA is out of range : */ + b[0] = 0.; + ncalc = Math.min(nb, 0) - 1; + } + return ncalc; + } + + // GNUR from bessel_k.c + + public static double bessel_k(double x, double alpha, double expo) { + /* NaNs propagated correctly */ + if (Double.isNaN(x) || Double.isNaN(alpha)) + return x + alpha; + if (x < 0) { + RMathError.error(RMathError.MLError.RANGE, "bessel_k"); + return ML_NAN; + } + int ize = (int) expo; + if (alpha < 0) + alpha = -alpha; + int nb = 1 + (int) Math.floor(alpha);/* nb-1 <= |alpha| < nb */ + alpha -= (double) (nb - 1); + double[] bk = new double[nb]; + int ncalc = K_bessel(x, alpha, nb, ize, bk); + if (ncalc != nb) {/* error input */ + if (ncalc < 0) + RMathError.warning(RError.Message.BESSEL_ARG_RANGE, "k", x, ncalc, nb, alpha); + else + RMathError.warning(RError.Message.BESSEL_PRECISION_LOST, "k", x, alpha + (double) nb - 1); + } + x = bk[nb - 1]; + return x; + } + + /* + * modified version of bessel_k that accepts a work array instead of allocating one. + */ + public static double bessel_k_ex(double x, double alpha, double expo, double[] bk) { + /* NaNs propagated correctly */ + if (Double.isNaN(x) || Double.isNaN(alpha)) + return x + alpha; + if (x < 0) { + RMathError.error(RMathError.MLError.RANGE, "bessel_k"); + return ML_NAN; + } + int ize = (int) expo; + if (alpha < 0) + alpha = -alpha; + int nb = 1 + (int) Math.floor(alpha);/* nb-1 <= |alpha| < nb */ + alpha -= (double) (nb - 1); + int ncalc = K_bessel(x, alpha, nb, ize, bk); + if (ncalc != nb) {/* error input */ + if (ncalc < 0) + RMathError.warning(RError.Message.BESSEL_ARG_RANGE, "k", x, ncalc, nb, alpha); + else + RMathError.warning(RError.Message.BESSEL_PRECISION_LOST, "k", x, alpha + (double) nb - 1); + } + x = bk[nb - 1]; + return x; + } + + /*------------------------------------------------------------------- + + This routine calculates modified Bessel functions + of the third kind, K_(N+ALPHA) (X), for non-negative + argument X, and non-negative order N+ALPHA, with or without + exponential scaling. + + Explanation of variables in the calling sequence + + X - Non-negative argument for which + K's or exponentially scaled K's (K*EXP(X)) + are to be calculated. If K's are to be calculated, + X must not be greater than XMAX_BESS_K. + ALPHA - Fractional part of order for which + K's or exponentially scaled K's (K*EXP(X)) are + to be calculated. 0 <= ALPHA < 1.0. + NB - Number of functions to be calculated, NB > 0. + The first function calculated is of order ALPHA, and the + last is of order (NB - 1 + ALPHA). + IZE - Type. IZE = 1 if unscaled K's are to be calculated, + = 2 if exponentially scaled K's are to be calculated. + BK - Output vector of length NB. If the + routine terminates normally (NCALC=NB), the vector BK + contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X), + or the corresponding exponentially scaled functions. + If (0 < NCALC < NB), BK(I) contains correct function + values for I <= NCALC, and contains the ratios + K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. + NCALC - Output variable indicating possible errors. + Before using the vector BK, the user should check that + NCALC=NB, i.e., all orders have been calculated to + the desired accuracy. See error returns below. + + + ******************************************************************* + + Error returns + + In case of an error, NCALC != NB, and not all K's are + calculated to the desired accuracy. + + NCALC < -1: An argument is out of range. For example, + NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= XMAX_BESS_K. + In this case, the B-vector is not calculated, + and NCALC is set to MIN0(NB,0)-2 so that NCALC != NB. + NCALC = -1: Either K(ALPHA,X) >= XINF or + K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) >= XINF. In this case, + the B-vector is not calculated. Note that again + NCALC != NB. + + 0 < NCALC < NB: Not all requested function values could + be calculated accurately. BK(I) contains correct function + values for I <= NCALC, and contains the ratios + K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. + + + Intrinsic functions required are: + + ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT + + + Acknowledgement + + This program is based on a program written by J. B. Campbell + (2) that computes values of the Bessel functions K of float + argument and float order. Modifications include the addition + of non-scaled functions, parameterization of machine + dependencies, and the use of more accurate approximations + for SINH and SIN. + + References: "On Temme's Algorithm for the Modified Bessel + Functions of the Third Kind," Campbell, J. B., + TOMS 6(4), Dec. 1980, pp. 581-586. + + "A FORTRAN IV Subroutine for the Modified Bessel + Functions of the Third Kind of Real Order and Real + Argument," Campbell, J. B., Report NRC/ERB-925, + National Research Council, Canada. + + Latest modification: May 30, 1989 + + Modified by: W. J. Cody and L. Stoltz + Applied Mathematics Division + Argonne National Laboratory + Argonne, IL 60439 + + ------------------------------------------------------------------- + */ + /*--------------------------------------------------------------------- + * Mathematical constants + * A = LOG(2) - Euler's constant + * D = SQRT(2/PI) + ---------------------------------------------------------------------*/ + private static final double a = .11593151565841244881; + + /*--------------------------------------------------------------------- + P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA + Euler's constant + Coefficients converted from hex to decimal and modified + by W. J. Cody, 2/26/82 */ + private static final double[] p = new double[]{.805629875690432845, 20.4045500205365151, + 157.705605106676174, 536.671116469207504, 900.382759291288778, + 730.923886650660393, 229.299301509425145, .822467033424113231}; + private static final double[] q = new double[]{29.4601986247850434, 277.577868510221208, + 1206.70325591027438, 2762.91444159791519, 3443.74050506564618, + 2210.63190113378647, 572.267338359892221}; + /* R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) */ + private static final double[] r = new double[]{-.48672575865218401848, 13.079485869097804016, + -101.96490580880537526, 347.65409106507813131, + 3.495898124521934782e-4}; + private static final double[] s = new double[]{-25.579105509976461286, 212.57260432226544008, + -610.69018684944109624, 422.69668805777760407}; + /* T - Approximation for SINH(Y)/Y */ + private static final double[] t = new double[]{1.6125990452916363814e-10, + 2.5051878502858255354e-8, 2.7557319615147964774e-6, + 1.9841269840928373686e-4, .0083333333333334751799, + .16666666666666666446}; + /*---------------------------------------------------------------------*/ + private static final double[] estm = new double[]{52.0583, 5.7607, 2.7782, 14.4303, 185.3004, 9.3715}; + private static final double[] estf = new double[]{41.8341, 7.1075, 6.4306, 42.511, 1.35633, 84.5096, 20.}; + + private static int K_bessel(double x, double alpha, int nb, int ize, double[] bk) { + /* Local variables */ + int iend, i, j, k, m, ii, mplus1; + double x2by4, twox, c, blpha, ratio, wminf; + double d1, d2, d3, f0, f1, f2, p0, q0, t1, t2, twonu; + double dm, ex, bk1, bk2, nu; + + ii = 0; /* -Wall */ + + ex = x; + nu = alpha; + int ncalc = Math.min(nb, 0) - 2; + if (nb > 0 && (0. <= nu && nu < 1.) && (1 <= ize && ize <= 2)) { + if (ex <= 0 || (ize == 1 && ex > xmax_BESS_K)) { + if (ex <= 0) { + if (ex < 0) { + RMathError.error(RMathError.MLError.RANGE, "K_bessel"); + } + for (i = 0; i < nb; i++) + bk[i] = ML_POSINF; + } else /* would only have underflow */ + for (i = 0; i < nb; i++) + bk[i] = 0.; + ncalc = nb; + return ncalc; + } + k = 0; + if (nu < sqxmin_BESS_K) { + nu = 0.; + } else if (nu > .5) { + k = 1; + nu -= 1.; + } + twonu = nu + nu; + iend = nb + k - 1; + c = nu * nu; + d3 = -c; + L420W: while (true) { + if (ex <= 1.) { + /* + * ------------------------------------------------------------ Calculation of + * P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA + * ------------------------------------------------------------ + */ + d1 = 0.; + d2 = p[0]; + t1 = 1.; + t2 = q[0]; + for (i = 2; i <= 7; i += 2) { + d1 = c * d1 + p[i - 1]; + d2 = c * d2 + p[i]; + t1 = c * t1 + q[i - 1]; + t2 = c * t2 + q[i]; + } + d1 = nu * d1; + t1 = nu * t1; + f1 = Math.log(ex); + f0 = a + nu * (p[7] - nu * (d1 + d2) / (t1 + t2)) - f1; + q0 = Math.exp(-nu * (a - nu * (p[7] + nu * (d1 - d2) / (t1 - t2)) - f1)); + f1 = nu * f0; + p0 = Math.exp(f1); + /* + * ----------------------------------------------------------- Calculation of F0 + * = ----------------------------------------------------------- + */ + d1 = r[4]; + t1 = 1.; + for (i = 0; i < 4; ++i) { + d1 = c * d1 + r[i]; + t1 = c * t1 + s[i]; + } + /* + * d2 := sinh(f1)/ nu = sinh(f1)/(f1/f0) = f0 * sinh(f1)/f1 + */ + if (fabs(f1) <= .5) { + f1 *= f1; + d2 = 0.; + for (i = 0; i < 6; ++i) { + d2 = f1 * d2 + t[i]; + } + d2 = f0 + f0 * f1 * d2; + } else { + d2 = Math.sinh(f1) / nu; + } + f0 = d2 - nu * d1 / (t1 * p0); + if (ex <= 1e-10) { + /* + * --------------------------------------------------------- X <= 1.0E-10 + * Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X) + * --------------------------------------------------------- + */ + bk[0] = f0 + ex * f0; + if (ize == 1) { + bk[0] -= ex * bk[0]; + } + ratio = p0 / f0; + c = ex * DBL_MAX; + if (k != 0) { + /* + * --------------------------------------------------- Calculation of + * K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X), ALPHA >= 1/2 + * --------------------------------------------------- + */ + ncalc = -1; + if (bk[0] >= c / ratio) { + return ncalc; + } + bk[0] = ratio * bk[0] / ex; + twonu += 2.; + ratio = twonu; + } + ncalc = 1; + if (nb == 1) + return ncalc; + + /* + * ----------------------------------------------------- Calculate + * K(ALPHA+L,X)/K(ALPHA+L-1,X), L = 1, 2, ... , NB-1 + * ----------------------------------------------------- + */ + ncalc = -1; + for (i = 1; i < nb; ++i) { + if (ratio >= c) + return ncalc; + + bk[i] = ratio / ex; + twonu += 2.; + ratio = twonu; + } + ncalc = 1; + break L420W; + } else { + /* + * ------------------------------------------------------ 10^-10 < X <= 1.0 + * ------------------------------------------------------ + */ + c = 1.; + x2by4 = ex * ex / 4.; + p0 = .5 * p0; + q0 = .5 * q0; + d1 = -1.; + d2 = 0.; + bk1 = 0.; + bk2 = 0.; + f1 = f0; + f2 = p0; + do { + d1 += 2.; + d2 += 1.; + d3 = d1 + d3; + c = x2by4 * c / d2; + f0 = (d2 * f0 + p0 + q0) / d3; + p0 /= d2 - nu; + q0 /= d2 + nu; + t1 = c * f0; + t2 = c * (p0 - d2 * f0); + bk1 += t1; + bk2 += t2; + } while (fabs(t1 / (f1 + bk1)) > DBL_EPSILON || + fabs(t2 / (f2 + bk2)) > DBL_EPSILON); + bk1 = f1 + bk1; + bk2 = 2. * (f2 + bk2) / ex; + if (ize == 2) { + d1 = Math.exp(ex); + bk1 *= d1; + bk2 *= d1; + } + wminf = estf[0] * ex + estf[1]; + } + } else if (DBL_EPSILON * ex > 1.) { + /* + * ------------------------------------------------- X > 1./EPS + * ------------------------------------------------- + */ + ncalc = nb; + bk1 = 1. / (M_SQRT_2dPI * Math.sqrt(ex)); + for (i = 0; i < nb; ++i) + bk[i] = bk1; + return ncalc; + + } else { + /* + * ------------------------------------------------------- X > 1.0 + * ------------------------------------------------------- + */ + twox = ex + ex; + blpha = 0.; + ratio = 0.; + if (ex <= 4.) { + /* + * ---------------------------------------------------------- Calculation of + * K(ALPHA+1,X)/K(ALPHA,X), 1.0 <= X <= 4.0 + * ---------------------------------------------------------- + */ + d2 = trunc(estm[0] / ex + estm[1]); + m = (int) d2; + d1 = d2 + d2; + d2 -= .5; + d2 *= d2; + for (i = 2; i <= m; ++i) { + d1 -= 2.; + d2 -= d1; + ratio = (d3 + d2) / (twox + d1 - ratio); + } + /* + * ----------------------------------------------------------- Calculation + * of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward recurrence and K(ALPHA,X) + * from the wronskian + * ----------------------------------------------------------- + */ + d2 = trunc(estm[2] * ex + estm[3]); + m = (int) d2; + c = fabs(nu); + d3 = c + c; + d1 = d3 - 1.; + f1 = DBL_MIN; + f0 = (2. * (c + d2) / ex + .5 * ex / (c + d2 + 1.)) * DBL_MIN; + for (i = 3; i <= m; ++i) { + d2 -= 1.; + f2 = (d3 + d2 + d2) * f0; + blpha = (1. + d1 / d2) * (f2 + blpha); + f2 = f2 / ex + f1; + f1 = f0; + f0 = f2; + } + f1 = (d3 + 2.) * f0 / ex + f1; + d1 = 0.; + t1 = 1.; + for (i = 1; i <= 7; ++i) { + d1 = c * d1 + p[i - 1]; + t1 = c * t1 + q[i - 1]; + } + p0 = Math.exp(c * (a + c * (p[7] - c * d1 / t1) - Math.log(ex))) / ex; + f2 = (c + .5 - ratio) * f1 / ex; + bk1 = p0 + (d3 * f0 - f2 + f0 + blpha) / (f2 + f1 + f0) * p0; + if (ize == 1) { + bk1 *= Math.exp(-ex); + } + wminf = estf[2] * ex + estf[3]; + } else { + /* + * --------------------------------------------------------- Calculation of + * K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by backward recurrence, for X > + * 4.0 ---------------------------------------------------------- + */ + dm = trunc(estm[4] / ex + estm[5]); + m = (int) dm; + d2 = dm - .5; + d2 *= d2; + d1 = dm + dm; + for (i = 2; i <= m; ++i) { + dm -= 1.; + d1 -= 2.; + d2 -= d1; + ratio = (d3 + d2) / (twox + d1 - ratio); + blpha = (ratio + ratio * blpha) / dm; + } + bk1 = 1. / ((M_SQRT_2dPI + M_SQRT_2dPI * blpha) * Math.sqrt(ex)); + if (ize == 1) + bk1 *= Math.exp(-ex); + wminf = estf[4] * (ex - fabs(ex - estf[6])) + estf[5]; + } + /* + * --------------------------------------------------------- Calculation of + * K(ALPHA+1,X) from K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X) + * --------------------------------------------------------- + */ + bk2 = bk1 + bk1 * (nu + .5 - ratio) / ex; + } + /*-------------------------------------------------------------------- + Calculation of 'NCALC', K(ALPHA+I,X), I = 0, 1, ... , NCALC-1, + & K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1 + -------------------------------------------------------------------*/ + ncalc = nb; + bk[0] = bk1; + if (iend == 0) + return ncalc; + + j = 1 - k; + if (j >= 0) + bk[j] = bk2; + + if (iend == 1) + return ncalc; + + m = Math.min((int) (wminf - nu), iend); + for (i = 2; i <= m; ++i) { + t1 = bk1; + bk1 = bk2; + twonu += 2.; + if (ex < 1.) { + if (bk1 >= DBL_MAX / twonu * ex) + break; + } else { + if (bk1 / ex >= DBL_MAX / twonu) + break; + } + bk2 = twonu / ex * bk1 + t1; + ii = i; + ++j; + if (j >= 0) { + bk[j] = bk2; + } + } + + m = ii; + if (m == iend) { + return ncalc; + } + ratio = bk2 / bk1; + mplus1 = m + 1; + ncalc = -1; + for (i = mplus1; i <= iend; ++i) { + twonu += 2.; + ratio = twonu / ex + 1. / ratio; + ++j; + if (j >= 1) { + bk[j] = ratio; + } else { + if (bk2 >= DBL_MAX / ratio) + return ncalc; + + bk2 *= ratio; + } + } + ncalc = Math.max(1, mplus1 - k); + if (ncalc == 1) + bk[0] = bk2; + if (nb == 1) + return ncalc; + break; // of extra while (true) + } + + L420: for (i = ncalc; i < nb; ++i) { /* i == *ncalc */ + bk[i] *= bk[i - 1]; + (ncalc)++; + } + } + return ncalc; + } + + // GNUR from bessel_y.c + + // unused now from R + public static double bessel_y(double x, double alpha) { + /* NaNs propagated correctly */ + if (Double.isNaN(x) || Double.isNaN(alpha)) + return x + alpha; + if (x < 0) { + RMathError.error(RMathError.MLError.RANGE, "bessel_y"); + return ML_NAN; + } + double na = Math.floor(alpha); + if (alpha < 0) { + /* + * Using Abramowitz & Stegun 9.1.2 this may not be quite optimal (CPU and accuracy wise) + */ + return (((alpha - na == 0.5) ? 0 : bessel_y(x, -alpha) * cospi(alpha)) - + ((alpha == na) ? 0 : bessel_j(x, -alpha) * sinpi(alpha))); + } else if (alpha > 1e7) { + RMathError.warning(RError.Message.BESSEL_NU_TOO_LARGE, "Y", alpha, "y"); + return ML_NAN; + } + int nb = 1 + (int) na;/* nb-1 <= alpha < nb */ + alpha -= (double) (nb - 1); + double[] by = new double[nb]; + int ncalc = Y_bessel(x, alpha, nb, by); + if (ncalc != nb) {/* error input */ + if (ncalc == -1) { + return ML_POSINF; + } else if (ncalc < -1) + RMathError.warning(RError.Message.BESSEL_ARG_RANGE, "y", x, ncalc, nb, alpha); + else /* ncalc >= 0 */ + RMathError.warning(RError.Message.BESSEL_PRECISION_LOST, "y", x, alpha + (double) nb - 1); + } + x = by[nb - 1]; + return x; + } + + /* + * Called from R: modified version of bessel_y(), accepting a work array instead of allocating + * one. + */ + public static double bessel_y_ex(double x, double alpha, double[] by) { + /* NaNs propagated correctly */ + if (Double.isNaN(x) || Double.isNaN(alpha)) + return x + alpha; + if (x < 0) { + RMathError.error(RMathError.MLError.RANGE, "bessel_y"); + return ML_NAN; + } + double na = Math.floor(alpha); + if (alpha < 0) { + /* + * Using Abramowitz & Stegun 9.1.2 this may not be quite optimal (CPU and accuracy wise) + */ + return (((alpha - na == 0.5) ? 0 : bessel_y_ex(x, -alpha, by) * cospi(alpha)) - + ((alpha == na) ? 0 : bessel_j_ex(x, -alpha, by) * sinpi(alpha))); + } else if (alpha > 1e7) { + RMathError.warning(RError.Message.BESSEL_NU_TOO_LARGE, "Y", alpha, "y"); + return ML_NAN; + } + int nb = 1 + (int) na;/* nb-1 <= alpha < nb */ + alpha -= (double) (nb - 1); + int ncalc = Y_bessel(x, alpha, nb, by); + if (ncalc != nb) {/* error input */ + if (ncalc == -1) + return ML_POSINF; + else if (ncalc < -1) + RMathError.warning(RError.Message.BESSEL_ARG_RANGE, "y", x, ncalc, nb, alpha); + else /* ncalc >= 0 */ + RMathError.warning(RError.Message.BESSEL_PRECISION_LOST, "y", x, alpha + (double) nb - 1); + } + x = by[nb - 1]; + return x; + } + + /* + * ---------------------------------------------------------------------- + * + * This routine calculates Bessel functions Y_(N+ALPHA) (X) v for non-negative argument X, and + * non-negative order N+ALPHA. + * + * + * Explanation of variables in the calling sequence + * + * X - Non-negative argument for which Y's are to be calculated. ALPHA - Fractional part of + * order for which Y's are to be calculated. 0 <= ALPHA < 1.0. NB - Number of functions to be + * calculated, NB > 0. The first function calculated is of order ALPHA, and the last is of order + * (NB - 1 + ALPHA). BY - Output vector of length NB. If the routine terminates normally + * (NCALC=NB), the vector BY contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X), If (0 < + * NCALC < NB), BY(I) contains correct function values for I <= NCALC, and contains the ratios + * Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array. NCALC - Output variable indicating + * possible errors. Before using the vector BY, the user should check that NCALC=NB, i.e., all + * orders have been calculated to the desired accuracy. See error returns below. + ******************************************************************* + * + * + * + * Error returns + * + * In case of an error, NCALC != NB, and not all Y's are calculated to the desired accuracy. + * + * NCALC < -1: An argument is out of range. For example, NB <= 0, IZE is not 1 or 2, or IZE=1 + * and ABS(X) >= XMAX. In this case, BY[0] = 0.0, the remainder of the BY-vector is not + * calculated, and NCALC is set to MIN0(NB,0)-2 so that NCALC != NB. NCALC = -1: Y(ALPHA,X) >= + * XINF. The requested function values are set to 0.0. 1 < NCALC < NB: Not all requested + * function values could be calculated accurately. BY(I) contains correct function values for I + * <= NCALC, and and the remaining NB-NCALC array elements contain 0.0. + * + * + * Intrinsic functions required are: + * + * DBLE, EXP, INT, MAX, MIN, REAL, SQRT + * + * + * Acknowledgement + * + * This program draws heavily on Temme's Algol program for Y(a,x) and Y(a+1,x) and on Campbell's + * programs for Y_nu(x). Temme's scheme is used for x < THRESH, and Campbell's scheme is used in + * the asymptotic region. Segments of code from both sources have been translated into Fortran + * 77, merged, and heavily modified. Modifications include parameterization of machine + * dependencies, use of a new approximation for ln(gamma(x)), and built-in protection against + * over/underflow. + * + * References: "Bessel functions J_nu(x) and Y_nu(x) of float order and float argument," + * Campbell, J. B., Comp. Phy. Comm. 18, 1979, pp. 133-142. + * + * "On the numerical evaluation of the ordinary Bessel function of the second kind," Temme, N. + * M., J. Comput. Phys. 21, 1976, pp. 343-350. + * + * Latest modification: March 19, 1990 + * + * Modified by: W. J. Cody Applied Mathematics Division Argonne National Laboratory Argonne, IL + * 60439 ---------------------------------------------------------------------- + */ + + /* + * ---------------------------------------------------------------------- Mathematical constants + * FIVPI = 5*PI PIM5 = 5*PI - 15 + * ---------------------------------------------------------------------- + */ + private static final double fivpi = 15.707963267948966192; + private static final double pim5 = .70796326794896619231; + + /*---------------------------------------------------------------------- + Coefficients for Chebyshev polynomial expansion of + 1/gamma(1-x), abs(x) <= .5 + ----------------------------------------------------------------------*/ + private static final double[] ch = new double[]{-6.7735241822398840964e-24, + -6.1455180116049879894e-23, 2.9017595056104745456e-21, + 1.3639417919073099464e-19, 2.3826220476859635824e-18, + -9.0642907957550702534e-18, -1.4943667065169001769e-15, + -3.3919078305362211264e-14, -1.7023776642512729175e-13, + 9.1609750938768647911e-12, 2.4230957900482704055e-10, + 1.7451364971382984243e-9, -3.3126119768180852711e-8, + -8.6592079961391259661e-7, -4.9717367041957398581e-6, + 7.6309597585908126618e-5, .0012719271366545622927, + .0017063050710955562222, -.07685284084478667369, + -.28387654227602353814, .92187029365045265648}; + + private static int Y_bessel(double x, double alpha, int nb, double[] by) { + /* Local variables */ + int i, k, na; + + double alfa, div, ddiv, even, gamma, term, cosmu, sinmu, + b, c, d, e, f, g, h, p, q, r, s, d1, d2, q0, pa, pa1, qa, qa1, + en, en1, nu, ex, ya, ya1, twobyx, den, odd, aye, dmu, x2, xna; + + en1 = ya = ya1 = 0; /* -Wall */ + + ex = x; + nu = alpha; + int ncalc; + if (nb > 0 && 0. <= nu && nu < 1.) { + if (ex < DBL_MIN || ex > xlrg_BESS_Y) { + /* + * Warning is not really appropriate, give proper limit: ML_ERROR(ME_RANGE, + * "Y_bessel"); + */ + ncalc = nb; + if (ex > xlrg_BESS_Y) + by[0] = 0.; /* was ML_POSINF */ + else if (ex < DBL_MIN) + by[0] = ML_NEGINF; + for (i = 0; i < nb; i++) + by[i] = by[0]; + return ncalc; + } + xna = trunc(nu + .5); + na = (int) xna; + if (na == 1) {/* <==> .5 <= *alpha < 1 <==> -5. <= nu < 0 */ + nu -= xna; + } + if (nu == -.5) { + p = M_SQRT_2dPI / Math.sqrt(ex); + ya = p * Math.sin(ex); + ya1 = -p * Math.cos(ex); + } else if (ex < 3.) { + /* + * ------------------------------------------------------------- Use Temme's scheme + * for small X ------------------------------------------------------------- + */ + b = ex * .5; + d = -Math.log(b); + f = nu * d; + e = pow(b, -nu); + if (fabs(nu) < M_eps_sinc) + c = M_1_PI; + else + c = nu / sinpi(nu); + + /* + * ------------------------------------------------------------ Computation of + * sinh(f)/f ------------------------------------------------------------ + */ + if (fabs(f) < 1.) { + x2 = f * f; + en = 19.; + s = 1.; + for (i = 1; i <= 9; ++i) { + s = s * x2 / en / (en - 1.) + 1.; + en -= 2.; + } + } else { + s = (e - 1. / e) * .5 / f; + } + /* + * -------------------------------------------------------- Computation of + * 1/gamma(1-a) using Chebyshev polynomials + */ + x2 = nu * nu * 8.; + aye = ch[0]; + even = 0.; + alfa = ch[1]; + odd = 0.; + for (i = 3; i <= 19; i += 2) { + even = -(aye + aye + even); + aye = -even * x2 - aye + ch[i - 1]; + odd = -(alfa + alfa + odd); + alfa = -odd * x2 - alfa + ch[i]; + } + even = (even * .5 + aye) * x2 - aye + ch[20]; + odd = (odd + alfa) * 2.; + gamma = odd * nu + even; + /* + * End of computation of 1/gamma(1-a) + * ----------------------------------------------------------- + */ + g = e * gamma; + e = (e + 1. / e) * .5; + f = 2. * c * (odd * e + even * s * d); + e = nu * nu; + p = g * c; + q = M_1_PI / g; + c = nu * M_PI_2; + if (fabs(c) < M_eps_sinc) + r = 1.; + else + r = sinpi(nu / 2) / c; + + r = M_PI * c * r * r; + c = 1.; + d = -b * b; + h = 0.; + ya = f + r * q; + ya1 = p; + en = 1.; + + while (fabs(g / (1. + fabs(ya))) + + fabs(h / (1. + fabs(ya1))) > DBL_EPSILON) { + f = (f * en + p + q) / (en * en - e); + c *= (d / en); + p /= en - nu; + q /= en + nu; + g = c * (f + r * q); + h = c * p - en * g; + ya += g; + ya1 += h; + en += 1.; + } + ya = -ya; + ya1 = -ya1 / b; + } else if (ex < thresh_BESS_Y) { + /* + * -------------------------------------------------------------- Use Temme's scheme + * for moderate X : 3 <= x < 16 + * -------------------------------------------------------------- + */ + c = (.5 - nu) * (.5 + nu); + b = ex + ex; + e = ex * M_1_PI * cospi(nu) / DBL_EPSILON; + e *= e; + p = 1.; + q = -ex; + r = 1. + ex * ex; + s = r; + en = 2.; + while (r * en * en < e) { + en1 = en + 1.; + d = (en - 1. + c / en) / s; + p = (en + en - p * d) / en1; + q = (-b + q * d) / en1; + s = p * p + q * q; + r *= s; + en = en1; + } + f = p / s; + p = f; + g = -q / s; + q = g; + L220: while (true) { + en -= 1.; + if (en > 0.) { + r = en1 * (2. - p) - 2.; + s = b + en1 * q; + d = (en - 1. + c / en) / (r * r + s * s); + p = d * r; + q = d * s; + e = f + 1.; + f = p * e - g * q; + g = q * e + p * g; + en1 = en; + continue; // goto L220; + } + break; // of extra while (true) + } + f = 1. + f; + d = f * f + g * g; + pa = f / d; + qa = -g / d; + d = nu + .5 - p; + q += ex; + pa1 = (pa * q - qa * d) / ex; + qa1 = (qa * q + pa * d) / ex; + b = ex - M_PI_2 * (nu + .5); + c = Math.cos(b); + s = Math.sin(b); + d = M_SQRT_2dPI / Math.sqrt(ex); + ya = d * (pa * s + qa * c); + ya1 = d * (qa1 * s - pa1 * c); + } else { /* x > thresh_BESS_Y */ + /* + * ---------------------------------------------------------- Use Campbell's + * asymptotic scheme. ---------------------------------------------------------- + */ + na = 0; + d1 = trunc(ex / fivpi); + i = (int) d1; + dmu = ex - 15. * d1 - d1 * pim5 - (alpha + .5) * M_PI_2; + if (i - (i / 2 << 1) == 0) { + cosmu = Math.cos(dmu); + sinmu = Math.sin(dmu); + } else { + cosmu = -Math.cos(dmu); + sinmu = -Math.sin(dmu); + } + ddiv = 8. * ex; + dmu = alpha; + den = Math.sqrt(ex); + for (k = 1; k <= 2; ++k) { + p = cosmu; + cosmu = sinmu; + sinmu = -p; + d1 = (2. * dmu - 1.) * (2. * dmu + 1.); + d2 = 0.; + div = ddiv; + p = 0.; + q = 0.; + q0 = d1 / div; + term = q0; + for (i = 2; i <= 20; ++i) { + d2 += 8.; + d1 -= d2; + div += ddiv; + term = -term * d1 / div; + p += term; + d2 += 8.; + d1 -= d2; + div += ddiv; + term *= (d1 / div); + q += term; + if (fabs(term) <= DBL_EPSILON) { + break; + } + } + p += 1.; + q += q0; + if (k == 1) + ya = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den; + else + ya1 = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den; + dmu += 1.; + } + } + if (na == 1) { + h = 2. * (nu + 1.) / ex; + if (h > 1.) { + if (fabs(ya1) > DBL_MAX / h) { + h = 0.; + ya = 0.; + } + } + h = h * ya1 - ya; + ya = ya1; + ya1 = h; + } + + /* + * --------------------------------------------------------------- Now have first one or + * two Y's --------------------------------------------------------------- + */ + by[0] = ya; + ncalc = 1; + L450W: while (true) { + if (nb > 1) { + by[1] = ya1; + if (ya1 != 0.) { + aye = 1. + alpha; + twobyx = 2. / ex; + ncalc = 2; + for (i = 2; i < nb; ++i) { + if (twobyx < 1.) { + if (fabs(by[i - 1]) * twobyx >= DBL_MAX / aye) + break L450W; + } else { + if (fabs(by[i - 1]) >= DBL_MAX / aye / twobyx) + break L450W; + } + by[i] = twobyx * aye * by[i - 1] - by[i - 2]; + aye += 1.; + ++(ncalc); + } + } + } + break; // of extra while (true) + } + L450: for (i = ncalc; i < nb; ++i) + by[i] = ML_NEGINF;/* was 0 */ + + } else { + by[0] = 0.; + ncalc = Math.min(nb, 0) - 1; + } + return ncalc; + } + + // GNUR from gamma_cody.c + + /* + * ---------------------------------------------------------------------- Mathematical constants + * ---------------------------------------------------------------------- + */ + private static final double sqrtpi = .9189385332046727417803297; /* == ??? */ + + private static double xbig = 171.624; + /* ML_POSINF == const double xinf = 1.79e308; */ + /* DBL_EPSILON = const double eps = 2.22e-16; */ + /* DBL_MIN == const double xminin = 2.23e-308; */ + + /*---------------------------------------------------------------------- + Numerator and denominator coefficients for rational minimax + approximation over (1,2). + ----------------------------------------------------------------------*/ + private static double[] pGC = new double[]{ + -1.71618513886549492533811, + 24.7656508055759199108314, -379.804256470945635097577, + 629.331155312818442661052, 866.966202790413211295064, + -31451.2729688483675254357, -36144.4134186911729807069, + 66456.1438202405440627855}; + private static double[] qGC = new double[]{ + -30.8402300119738975254353, + 315.350626979604161529144, -1015.15636749021914166146, + -3107.77167157231109440444, 22538.1184209801510330112, + 4755.84627752788110767815, -134659.959864969306392456, + -115132.259675553483497211}; + /*---------------------------------------------------------------------- + Coefficients for minimax approximation over (12, INF). + ----------------------------------------------------------------------*/ + private static double cGC[] = new double[]{ + -.001910444077728, 8.4171387781295e-4, + -5.952379913043012e-4, 7.93650793500350248e-4, + -.002777777777777681622553, .08333333333333333331554247, + .0057083835261}; + + private static final double gammaCody(double x) { + + /* + * ---------------------------------------------------------------------- + * + * This routine calculates the GAMMA function for a float argument X. Computation is based + * on an algorithm outlined in reference [1]. The program uses rational functions that + * approximate the GAMMA function to at least 20 significant decimal digits. Coefficients + * for the approximation over the interval (1,2) are unpublished. Those for the + * approximation for X >= 12 are from reference [2]. The accuracy achieved depends on the + * arithmetic system, the compiler, the intrinsic functions, and proper selection of the + * machine-dependent constants. + ******************************************************************* + * + * + * Error returns + * + * The program returns the value XINF for singularities or when overflow would occur. The + * computation is believed to be free of underflow and overflow. + * + * Intrinsic functions required are: + * + * INT, DBLE, EXP, LOG, REAL, SIN + * + * + * References: [1] "An Overview of Software Development for Special Functions", W. J. Cody, + * Lecture Notes in Mathematics, 506, Numerical Analysis Dundee, 1975, G. A. Watson (ed.), + * Springer Verlag, Berlin, 1976. + * + * [2] Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968. + * + * Latest modification: October 12, 1989 + * + * Authors: W. J. Cody and L. Stoltz Applied Mathematics Division Argonne National + * Laboratory Argonne, IL 60439 + * ---------------------------------------------------------------------- + */ + + /* + * ******************************************************************* + * + * Explanation of machine-dependent constants + * + * beta - radix for the floating-point representation maxexp - the smallest positive power + * of beta that overflows XBIG - the largest argument for which GAMMA(X) is representable in + * the machine, i.e., the solution to the equation GAMMA(XBIG) = beta**maxexp XINF - the + * largest machine representable floating-point number; approximately beta**maxexp EPS - the + * smallest positive floating-point number such that 1.0+EPS > 1.0 XMININ - the smallest + * positive floating-point number such that 1/XMININ is machine representable + * + * Approximate values for some important machines are: + * + * beta maxexp XBIG + * + * CRAY-1 (S.P.) 2 8191 966.961 Cyber 180/855 under NOS (S.P.) 2 1070 177.803 IEEE (IBM/XT, + * SUN, etc.) (S.P.) 2 128 35.040 IEEE (IBM/XT, SUN, etc.) (D.P.) 2 1024 171.624 IBM 3033 + * (D.P.) 16 63 57.574 VAX D-Format (D.P.) 2 127 34.844 VAX G-Format (D.P.) 2 1023 171.489 + * + * XINF EPS XMININ + * + * CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 Cyber 180/855 under NOS (S.P.) 1.26E+322 + * 3.55E-15 3.14E-294 IEEE (IBM/XT, SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 IEEE + * (IBM/XT, SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 IBM 3033 (D.P.) 7.23D+75 2.22D-16 + * 1.39D-76 VAX D-Format (D.P.) 1.70D+38 1.39D-17 5.88D-39 VAX G-Format (D.P.) 8.98D+307 + * 1.11D-16 1.12D-308 + ******************************************************************* + * + * + * ---------------------------------------------------------------------- Machine dependent + * parameters ---------------------------------------------------------------------- + */ + + /* Local variables */ + int i, n; + double fact, xden, xnum, y, z, yi, res, sum, ysq; + + boolean parity = false; + fact = 1.; + n = 0; + y = x; + if (y <= 0.) { + /* + * ------------------------------------------------------------- Argument is negative + * ------------------------------------------------------------- + */ + y = -x; + yi = trunc(y); + res = y - yi; + if (res != 0.) { + if (yi != trunc(yi * .5) * 2.) + parity = true; + fact = -M_PI / sinpi(res); + y += 1.; + } else { + return (ML_POSINF); + } + } + /* + * ----------------------------------------------------------------- Argument is positive + * ----------------------------------------------------------------- + */ + if (y < DBL_EPSILON) { + /* + * -------------------------------------------------------------- Argument < EPS + * -------------------------------------------------------------- + */ + if (y >= DBL_MIN) { + res = 1. / y; + } else { + return (ML_POSINF); + } + } else if (y < 12.) { + yi = y; + if (y < 1.) { + /* + * --------------------------------------------------------- EPS < argument < 1 + * --------------------------------------------------------- + */ + z = y; + y += 1.; + } else { + /* + * ----------------------------------------------------------- 1 <= argument < 12, + * reduce argument if necessary + * ----------------------------------------------------------- + */ + n = (int) y - 1; + y -= (double) n; + z = y - 1.; + } + /* + * --------------------------------------------------------- Evaluate approximation for + * 1. < argument < 2. --------------------------------------------------------- + */ + xnum = 0.; + xden = 1.; + for (i = 0; i < 8; ++i) { + xnum = (xnum + pGC[i]) * z; + xden = xden * z + qGC[i]; + } + res = xnum / xden + 1.; + if (yi < y) { + /* + * -------------------------------------------------------- Adjust result for case + * 0. < argument < 1. -------------------------------------------------------- + */ + res /= yi; + } else if (yi > y) { + /* + * ---------------------------------------------------------- Adjust result for case + * 2. < argument < 12. ---------------------------------------------------------- + */ + for (i = 0; i < n; ++i) { + res *= y; + y += 1.; + } + } + } else { + /* + * ------------------------------------------------------------- Evaluate for argument + * >= 12., ------------------------------------------------------------- + */ + if (y <= xbig) { + ysq = y * y; + sum = cGC[6]; + for (i = 0; i < 6; ++i) { + sum = sum / ysq + cGC[i]; + } + sum = sum / y - y + sqrtpi; + sum += (y - .5) * Math.log(y); + res = Math.exp(sum); + } else { + return (ML_POSINF); + } + } + /* + * ---------------------------------------------------------------------- Final adjustments + * and return ---------------------------------------------------------------------- + */ + if (parity) + res = -res; + if (fact != 1.) + res = fact / res; + return res; + } + +} +// Checkstyle: resume +// @formatter:on diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/Beta.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/Beta.java new file mode 100644 index 0000000000..ba0a5b30de --- /dev/null +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/Beta.java @@ -0,0 +1,51 @@ +/* + * This material is distributed under the GNU General Public License + * Version 2. You may review the terms of this license at + * http://www.gnu.org/licenses/gpl-2.0.html + * + * Copyright (C) 1998 Ross Ihaka + * Copyright (c) 2000-2014, The R Core Team + * Copyright (c) 2018, Oracle and/or its affiliates + * + * All rights reserved. + */ +package com.oracle.truffle.r.runtime.nmath; + +import com.oracle.truffle.r.runtime.RRuntime; +import static com.oracle.truffle.r.runtime.nmath.GammaFunctions.gammafn; +import static com.oracle.truffle.r.runtime.nmath.LBeta.lbeta; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.ML_POSINF; + +public class Beta { + + private static final double xmin = GammaFunctions.gfn_xmin; + private static final double xmax = GammaFunctions.gfn_xmax; + private static final double lnsml = -708.39641853226412; + + public static double beta(double a, double b) { + if (Double.isNaN(a) || Double.isNaN(b)) { + return a + b; + } + if (a < 0 || b < 0) + return RMathError.defaultError(); // ML_ERR_return_NAN + else if (a == 0 || b == 0) + return ML_POSINF; + else if (!RRuntime.isFinite(a) || !RRuntime.isFinite(b)) + return 0; + + if (a + b < xmax) {/* ~= 171.61 for IEEE */ + // return gammafn(a) * gammafn(b) / gammafn(a+b); + /* + * All the terms are positive, and all can be large for large or small arguments. They + * are never much less than one. gammafn(x) can still overflow for x ~ 1e-308, but the + * result would too. + */ + return (1 / gammafn(a + b)) * gammafn(a) * gammafn(b); + } else { + double val = lbeta(a, b); + // underflow to 0 is not harmful per se; exp(-999) also gives no warning + return Math.exp(val); + } + } + +} diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/GammaFunctions.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/GammaFunctions.java index 91284883c1..4c6a90ec0f 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/GammaFunctions.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/GammaFunctions.java @@ -8,7 +8,7 @@ * Copyright (c) 1998--2014, The R Core Team * Copyright (c) 2002--2010, The R Foundation * Copyright (C) 2005--2006, Morten Welinder - * Copyright (c) 2014, 2017, Oracle and/or its affiliates + * Copyright (c) 2014, 2018, Oracle and/or its affiliates * * based on AS 91 (C) 1979 Royal Statistical Society * and on AS 111 (C) 1977 Royal Statistical Society @@ -27,11 +27,20 @@ import static com.oracle.truffle.r.runtime.nmath.DPQ.rdtclog; import static com.oracle.truffle.r.runtime.nmath.DPQ.rdtqiv; import static com.oracle.truffle.r.runtime.nmath.DPQ.rlog1exp; import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_EPSILON; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MANT_DIG; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MAX_EXP; import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MIN; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MIN_EXP; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_LOG10_2; import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_LN2; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_PI; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.ML_NAN; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.ML_POSINF; +import static com.oracle.truffle.r.runtime.nmath.RMath.fmax2; import com.oracle.truffle.api.CompilerDirectives.CompilationFinal; import com.oracle.truffle.api.CompilerDirectives.TruffleBoundary; +import com.oracle.truffle.r.runtime.RError; import com.oracle.truffle.r.runtime.RError.Message; import com.oracle.truffle.r.runtime.RRuntime; import com.oracle.truffle.r.runtime.Utils; @@ -121,8 +130,8 @@ public abstract class GammaFunctions { * non-trivial, see ./gammalims.c xsml = exp(.01)*DBL_MIN dxrel = sqrt(DBL_EPSILON) = 2 ^ -26 */ private static final int ngam = 22; - private static final double gfn_xmin = -170.5674972726612; - private static final double gfn_xmax = 171.61447887182298; + static final double gfn_xmin = -170.5674972726612; + static final double gfn_xmax = 171.61447887182298; private static final double gfn_xsml = 2.2474362225598545e-308; private static final double gfn_dxrel = 1.490116119384765696e-8; @@ -258,12 +267,12 @@ public abstract class GammaFunctions { private static final double M_LN_SQRT_PId2 = 0.225791352644727432363097614947; // convert C's int* sgn to a 1-element int[] - static double lgammafnSign(double x, int[] sgn) { + public static double lgammafnSign(double x, int[] sgn) { double ans; double y; double sinpiy; - if (sgn[0] != 0) { + if (sgn != null) { sgn[0] = 1; } @@ -271,10 +280,8 @@ public abstract class GammaFunctions { return x; } - if (x < 0 && RMath.fmod(Math.floor(-x), 2.) == 0) { - if (sgn[0] != 0) { - sgn[0] = 1; - } + if (sgn != null && x < 0 && RMath.fmod(Math.floor(-x), 2.) == 0) { + sgn[0] = -1; } if (x <= 0 && x == (long) x) { /* Negative integer argument */ @@ -672,7 +679,7 @@ public abstract class GammaFunctions { private static final double minLog1Value = -0.79149064; /* Accurate calculation of log(1+x)-x, particularly for small x. */ - private static double log1pmx(double x) { + public static double log1pmx(double x) { if (x > 1 || x < minLog1Value) { return RMath.log1p(x) - x; } else { /* @@ -712,7 +719,7 @@ public abstract class GammaFunctions { private static final double tol_logcf = 1e-14; /* Compute log(gamma(a+1)) accurately also for small a (0 < a < 0.5). */ - private static double lgamma1p(double a) { + public static double lgamma1p(double a) { double lgam; int i; @@ -1184,4 +1191,484 @@ public abstract class GammaFunctions { pr = DPois.dpoisRaw(shape - 1, x / scale, giveLog); return giveLog ? pr - Math.log(scale) : pr / scale; } + + // the following is transcribed from polygamma.c + + @CompilationFinal(dimensions = 1) private static final double[] bvalues = new double[]{1.00000000000000000e+00, -5.00000000000000000e-01, 1.66666666666666667e-01, -3.33333333333333333e-02, + 2.38095238095238095e-02, -3.33333333333333333e-02, 7.57575757575757576e-02, -2.53113553113553114e-01, 1.16666666666666667e+00, -7.09215686274509804e+00, + 5.49711779448621554e+01, -5.29124242424242424e+02, 6.19212318840579710e+03, -8.65802531135531136e+04, 1.42551716666666667e+06, -2.72982310678160920e+07, + 6.01580873900642368e+08, -1.51163157670921569e+10, 4.29614643061166667e+11, -1.37116552050883328e+13, 4.88332318973593167e+14, -1.92965793419400681e+16}; + + private static final int n_max = 100; + // the following is actually a parameter in the original code, but it's always 1 and must be + // as the original code treats the "ans" value of type double as an array, which is legal + // only if a the first element of the array is accessed at all times + private static final int m = 1; + + public static final class DpsiFnResult { + + public double ans; + public int nz; + public int ierr; + + DpsiFnResult(double ans) { + this.ans = ans; + } + + DpsiFnResult(double ans, int nz, int ierr) { + this.ans = ans; + this.nz = nz; + this.ierr = ierr; + } + } + + public static DpsiFnResult dpsifn(double x, int n, int kode, double ans) { + DpsiFnResult result = new DpsiFnResult(ans); + dpsifn(x, n, kode, result); + return result; + } + + public static DpsiFnResult dpsifn(double x, int n, int kode, double ans, int nz, int ierr) { + DpsiFnResult result = new DpsiFnResult(ans, nz, ierr); + dpsifn(x, n, kode, result); + return result; + } + + @TruffleBoundary + private static void dpsifn(double x, int n, int kode, DpsiFnResult result) { + int mm; + int mx; + int nn; + int np; + int nx; + int fn; + double arg; + double den; + double elim; + double eps; + double fln; + double rln; + double r1m4; + double r1m5; + double s; + double slope; + double t; + double tk; + double tt; + double t1; + double t2; + double wdtol; + double xdmln; + double xdmy; + double xinc; + double xln = 0.0; + double xm; + double xmin; + double yint; + double[] trm = new double[23]; + double[] trmr = new double[n_max + 1]; + + // non-zero ierr always results in generating a NaN + result.ierr = 0; + if (n < 0 || kode < 1 || kode > 2 || m < 1) { + result.ans = Double.NaN; + return; + } + if (x <= 0.) { + /* + * use Abramowitz & Stegun 6.4.7 "Reflection Formula" psi(k, x) = (-1)^k psi(k, 1-x) - + * pi^{n+1} (d/dx)^n cot(x) + */ + if (x == RMath.round(x)) { + /* non-positive integer : +Inf or NaN depends on n */ + // for(j=0; j < m; j++) /* k = j + n : */ + // ans[j] = ((j+n) % 2) ? ML_POSINF : ML_NAN; + // m is always 1 + result.ans = (n % 2) != 0 ? ML_POSINF : ML_NAN; + return; + } + /* This could cancel badly */ + dpsifn(1. - x, n, /* kode = */1, result); + /* + * ans[j] == (-1)^(k+1) / gamma(k+1) * psi(k, 1 - x) for j = 0:(m-1) , k = n + j + */ + + /* Cheat for now: only work for m = 1, n in {0,1,2,3} : */ + if (m > 1 || n > 3) { /* doesn't happen for digamma() .. pentagamma() */ + /* not yet implemented */ + // non-zero ierr always results in generating a NaN + result.ierr = 4; + result.ans = Double.NaN; + return; + } + x *= M_PI; /* pi * x */ + if (n == 0) { + tt = Math.cos(x) / Math.sin(x); + } else if (n == 1) { + tt = -1 / Math.pow(Math.sin(x), 2); + } else if (n == 2) { + tt = 2 * Math.cos(x) / Math.pow(Math.sin(x), 3); + } else if (n == 3) { + tt = -2 * (2 * Math.pow(Math.cos(x), 2) + 1.) / Math.pow(Math.sin(x), 4); + } else { /* can not happen! */ + tt = RRuntime.DOUBLE_NA; + } + /* end cheat */ + + s = (n % 2) != 0 ? -1. : 1.; /* s = (-1)^n */ + /* + * t := pi^(n+1) * d_n(x) / gamma(n+1) , where d_n(x) := (d/dx)^n cot(x) + */ + t1 = t2 = s = 1.; + for (int k = 0, j = k - n; j < m; k++, j++, s = -s) { + /* k == n+j , s = (-1)^k */ + t1 *= M_PI; /* t1 == pi^(k+1) */ + if (k >= 2) { + t2 *= k; /* t2 == k! == gamma(k+1) */ + } + if (j >= 0) { /* by cheat above, tt === d_k(x) */ + // j must always be 0 + assert j == 0; + // ans[j] = s*(ans[j] + t1/t2 * tt); + result.ans = s * (result.ans + t1 / t2 * tt); + } + } + if (n == 0 && kode == 2) { /* unused from R, but "wrong": xln === 0 : */ + // ans[0] += xln; + result.ans += xln; + } + return; + } /* x <= 0 */ + + /* else : x > 0 */ + // nz not used + result.nz = 0; + xln = Math.log(x); + if (kode == 1 /* && m == 1 */) { /* the R case --- for very large x: */ + double lrg = 1 / (2. * DBL_EPSILON); + if (n == 0 && x * xln > lrg) { + // ans[0] = -xln; + result.ans = -xln; + return; + } else if (n >= 1 && x > n * lrg) { + // ans[0] = exp(-n * xln)/n; /* == x^-n / n == 1/(n * x^n) */ + result.ans = Math.exp(-n * xln) / n; + return; + } + } + mm = m; + // nx = imin2(-Rf_i1mach(15), Rf_i1mach(16));/* = 1021 */ + nx = Math.min(-DBL_MIN_EXP, DBL_MAX_EXP); + assert (nx == 1021); + r1m5 = M_LOG10_2; // Rf_d1mach(5); + r1m4 = DBL_EPSILON * 0.5; // Rf_d1mach(4) * 0.5; + wdtol = fmax2(r1m4, 0.5e-18); /* 1.11e-16 */ + + /* elim = approximate exponential over and underflow limit */ + elim = 2.302 * (nx * r1m5 - 3.0); /* = 700.6174... */ + for (;;) { + nn = n + mm - 1; + fn = nn; + t = (fn + 1) * xln; + + /* overflow and underflow test for small and large x */ + + if (Math.abs(t) > elim) { + if (t <= 0.0) { + // nz not used + result.nz = 0; + // non-zero ierr always results in generating a NaN + result.ierr = 2; + result.ans = Double.NaN; + return; + } + } else { + if (x < wdtol) { + // ans[0] = R_pow_di(x, -n-1); + result.ans = Math.pow(x, -n - 1); + if (mm != 1) { + // for(k = 1; k < mm ; k++) + // ans[k] = ans[k-1] / x; + assert mm < 2; + // int the original code, ans should not be accessed beyond the 0th + // index + } + if (n == 0 && kode == 2) { + // ans[0] += xln; + result.ans += xln; + } + return; + } + + /* compute xmin and the number of terms of the series, fln+1 */ + + rln = r1m5 * DBL_MANT_DIG; // Rf_i1mach(14); + rln = Math.min(rln, 18.06); + fln = Math.max(rln, 3.0) - 3.0; + yint = 3.50 + 0.40 * fln; + slope = 0.21 + fln * (0.0006038 * fln + 0.008677); + xm = yint + slope * fn; + mx = (int) xm + 1; + xmin = mx; + if (n != 0) { + xm = -2.302 * rln - Math.min(0.0, xln); + arg = xm / n; + arg = Math.min(0.0, arg); + eps = Math.exp(arg); + xm = 1.0 - eps; + if (Math.abs(arg) < 1.0e-3) { + xm = -arg; + } + fln = x * xm / eps; + xm = xmin - x; + if (xm > 7.0 && fln < 15.0) { + break; + } + } + xdmy = x; + xdmln = xln; + xinc = 0.0; + if (x < xmin) { + nx = (int) x; + xinc = xmin - nx; + xdmy = x + xinc; + xdmln = Math.log(xdmy); + } + + /* generate w(n+mm-1, x) by the asymptotic expansion */ + + t = fn * xdmln; + t1 = xdmln + xdmln; + t2 = t + xdmln; + tk = Math.max(Math.abs(t), fmax2(Math.abs(t1), Math.abs(t2))); + if (tk <= elim) { /* for all but large x */ + l10(t, tk, xdmy, xdmln, x, nn, nx, wdtol, fn, trm, trmr, xinc, mm, kode, result); + return; + } + } + // nz not used + result.nz++; /* underflow */ + mm--; + // ans[mm] = 0.; + assert mm == 0; + result.ans = 0.; + if (mm == 0) { + return; + } + } /* end{for()} */ + nn = (int) fln + 1; + np = n + 1; + t1 = (n + 1) * xln; + t = Math.exp(-t1); + s = t; + den = x; + for (int i = 1; i <= nn; i++) { + den += 1.; + trm[i] = Math.pow(den, -np); + s += trm[i]; + } + // ans[0] = s; + result.ans = s; + if (n == 0 && kode == 2) { + // ans[0] = s + xln; + result.ans = s + xln; + } + + if (mm != 1) { /* generate higher derivatives, j > n */ + assert false; + // tol = wdtol / 5.0; + // for(j = 1; j < mm; j++) { + // t /= x; + // s = t; + // tols = t * tol; + // den = x; + // for(i=1; i <= nn; i++) { + // den += 1.; + // trm[i] /= den; + // s += trm[i]; + // if (trm[i] < tols) { + // break; + // } + // } + // ans[j] = s; + // } + } + } + + private static void l10(double oldT, double oldTk, double xdmy, double xdmln, double x, double nn, double oldNx, double wdtol, double oldFn, double[] trm, double[] trmr, double xinc, + double mm, int kode, DpsiFnResult result) { + double t = oldT; + double tk = oldTk; + double nx = oldNx; + double fn = oldFn; + + double tss = Math.exp(-t); + double tt = 0.5 / xdmy; + double t1 = tt; + double tst = wdtol * tt; + if (nn != 0) { + t1 = tt + 1.0 / fn; + } + double rxsq = 1.0 / (xdmy * xdmy); + double ta = 0.5 * rxsq; + t = (fn + 1) * ta; + double s = t * bvalues[2]; + if (Math.abs(s) >= tst) { + tk = 2.0; + for (int k = 4; k <= 22; k++) { + t = t * ((tk + fn + 1) / (tk + 1.0)) * ((tk + fn) / (tk + 2.0)) * rxsq; + trm[k] = t * bvalues[k - 1]; + if (Math.abs(trm[k]) < tst) { + break; + } + s += trm[k]; + tk += 2.; + } + } + s = (s + t1) * tss; + if (xinc != 0.0) { + + /* backward recur from xdmy to x */ + + nx = (int) xinc; + double np = nn + 1; + if (nx > n_max) { + // nz not used + result.nz = 0; + // non-zero ierr always results in generating a NaN + result.ierr = 3; + result.ans = Double.NaN; + return; + } else { + if (nn == 0) { + l20(xdmln, xdmy, x, s, nx, kode, result); + return; + } + double xm = xinc - 1.0; + double fx = x + xm; + + /* this loop should not be changed. fx is accurate when x is small */ + for (int i = 1; i <= nx; i++) { + trmr[i] = Math.pow(fx, -np); + s += trmr[i]; + xm -= 1.; + fx = x + xm; + } + } + } + // ans[mm-1] = s; + assert (mm - 1) == 0; + result.ans = s; + if (fn == 0) { + l30(xdmln, xdmy, x, s, kode, result); + return; + } + + /* generate lower derivatives, j < n+mm-1 */ + + for (int j = 2; j <= mm; j++) { + fn--; + tss *= xdmy; + t1 = tt; + if (fn != 0) { + t1 = tt + 1.0 / fn; + } + t = (fn + 1) * ta; + s = t * bvalues[2]; + if (Math.abs(s) >= tst) { + tk = 4 + fn; + for (int k = 4; k <= 22; k++) { + trm[k] = trm[k] * (fn + 1) / tk; + if (Math.abs(trm[k]) < tst) { + break; + } + s += trm[k]; + tk += 2.; + } + } + s = (s + t1) * tss; + if (xinc != 0.0) { + if (fn == 0) { + l20(xdmln, xdmy, x, s, nx, kode, result); + return; + } + double xm = xinc - 1.0; + double fx = x + xm; + for (int i = 1; i <= nx; i++) { + trmr[i] = trmr[i] * fx; + s += trmr[i]; + xm -= 1.; + fx = x + xm; + } + } + // ans[mm - j] = s; + assert (mm - j) == 0; + result.ans = s; + if (fn == 0) { + l30(xdmln, xdmy, x, s, kode, result); + return; + } + } + + } + + private static void l20(double xdmln, double xdmy, double x, double oldS, double nx, int kode, DpsiFnResult result) { + double s = oldS; + for (int i = 1; i <= nx; i++) { + s += 1. / (x + (nx - i)); /* avoid disastrous cancellation, PR#13714 */ + } + + l30(xdmln, xdmy, x, s, kode, result); + } + + private static void l30(double xdmln, double xdmy, double x, double s, int kode, DpsiFnResult result) { + if (kode != 2) { /* always */ + // ans[0] = s - xdmln; + result.ans = s - xdmln; + } else if (xdmy != x) { + double xq; + xq = xdmy / x; + // ans[0] = s - log(xq); + result.ans = s - Math.log(xq); + } + } + + public static double psigamma(double x, double deriv) { + // polygamma.c + /* n-th derivative of psi(x); e.g., psigamma(x,0) == digamma(x) */ + if (Double.isNaN(x)) { + return x; + } + deriv = RMath.forceint(deriv); + int n = (int) deriv; + if (n > n_max) { + RMathError.warning(RError.Message.DERIV_OVER_N_MAX, n, n_max); + return ML_NAN; + } + DpsiFnResult result = dpsifn(x, n, 1, 0); + // ML_TREAT_psigam(ierr); + /* Now, ans == A := (-1)^(n+1) / gamma(n+1) * psi(n, x) */ + result.ans = -result.ans; /* = (-1)^(0+1) * gamma(0+1) * A */ + for (int k = 1; k <= n; k++) { + result.ans *= (-k); /* = (-1)^(k+1) * gamma(k+1) * A */ + } + return result.ans; /* = psi(n, x) */ + } + + public static double digamma(double x) { + return -dpsifn(x, 0, 1, 0).ans; + } + + public static double trigamma(double x) { + return dpsifn(x, 1, 1, 1).ans; + } + + public static double tetragamma(double x) { + return -2.0 * dpsifn(x, 2, 1, 1).ans; + } + + public static double pentagamma(double x) { + return 6.0 * dpsifn(x, 3, 1, 1).ans; + } + } diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/MathConstants.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/MathConstants.java index 543ce27dc9..4b772eb36e 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/MathConstants.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/MathConstants.java @@ -6,7 +6,7 @@ * Copyright (C) 1998 Ross Ihaka * Copyright (c) 1998--2012, The R Core Team * Copyright (c) 2004, The R Foundation - * Copyright (c) 2013, 2017, Oracle and/or its affiliates + * Copyright (c) 2013, 2018, Oracle and/or its affiliates * * All rights reserved. */ @@ -86,8 +86,17 @@ public final class MathConstants { public static final double DBL_EPSILON = Math.ulp(1.0); + /* Max.expon. of 10 (=308.2547) */ + public static final int MAX10E = (int) (DBL_MAX_EXP * M_LOG10_2); + + public static final long LONG_MAX = Long.MAX_VALUE; + public static final double ML_NAN = Double.NaN; + public static final double ML_NEGINF = Double.NEGATIVE_INFINITY; + + public static final double ML_POSINF = Double.POSITIVE_INFINITY; + // Different to Double.MIN_VALUE! public static final double DBL_MIN = Double.MIN_NORMAL; @@ -106,4 +115,20 @@ public final class MathConstants { public static double logspaceAdd(double logx, double logy) { return Math.max(logx, logy) + Math.log1p(Math.exp(-Math.abs(logx - logy))); } + + /** + * Compute the log of a difference from logs of terms, i.e., + * + * log (exp (logx) - exp (logy)) + * + * without causing overflows and without throwing away large handfuls of accuracy. + */ + public static double logspaceSub(double logx, double logy) { + return logx + log1Exp(logy - logx); + } + + private static double log1Exp(double x) { + return (x) > -M_LN2 ? Math.log(-Math.expm1(x)) : Math.log1p(-Math.exp(x)); + } + } diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/RMath.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/RMath.java index 5db1157cc3..6ffed815d4 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/RMath.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/RMath.java @@ -6,13 +6,17 @@ * Copyright (C) 1998 Ross Ihaka * Copyright (c) 1998--2012, The R Core Team * Copyright (c) 2004, The R Foundation - * Copyright (c) 2013, 2017, Oracle and/or its affiliates + * Copyright (c) 2013, 2018, Oracle and/or its affiliates * * All rights reserved. */ package com.oracle.truffle.r.runtime.nmath; +import static com.oracle.truffle.r.runtime.nmath.Arithmetic.powDi; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_EPSILON; import static com.oracle.truffle.r.runtime.nmath.MathConstants.DBL_MIN; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.LONG_MAX; +import static com.oracle.truffle.r.runtime.nmath.MathConstants.MAX10E; import static com.oracle.truffle.r.runtime.nmath.MathConstants.M_LN_SQRT_2PI; import static com.oracle.truffle.r.runtime.nmath.TOMS708.fabs; @@ -65,6 +69,107 @@ public final class RMath { return forceint(x); } + public static double trunc(double x) { + if (x > 0) { + return Math.floor(x); + } else { + return Math.ceil(x); + } + } + + // GNUR from fprec.c + + private static final int MAX_DIGITS = 22; + + public static double fprec(double x, double digits) { + if (Double.isNaN(x) || Double.isNaN(digits)) + return x + digits; + if (!RRuntime.isFinite(x)) { + return x; + } + if (!RRuntime.isFinite(digits)) { + if (digits > 0.0) { + return x; + } else { + digits = 1.0; + } + } + if (x == 0) { + return x; + } + int dig = (int) round(digits); + if (dig > MAX_DIGITS) { + return x; + } else if (dig < 1) { + dig = 1; + } + + double sgn = 1.0; + if (x < 0.0) { + sgn = -sgn; + x = -x; + } + double l10 = Math.log10(x); + int e10 = (int) (dig - 1 - Math.floor(l10)); + if (fabs(l10) < MAX10E - 2) { + double p10 = 1.0; + if (e10 > MAX10E) { /* numbers less than 10^(dig-1) * 1e-308 */ + p10 = powDi(10., e10 - MAX10E); + e10 = MAX10E; + } + if (e10 > 0) { /* + * Try always to have pow >= 1 and so exactly representable + */ + double pow10 = powDi(10., e10); + return (sgn * (rint((x * pow10) * p10) / pow10) / p10); + } else { + double pow10 = powDi(10., -e10); + return (sgn * (rint((x / pow10)) * pow10)); + } + } else { /* -- LARGE or small -- */ + boolean do_round = (MAX10E - l10 >= powDi(10., -dig)); + int e2 = dig + ((e10 > 0) ? 1 : -1) * MAX_DIGITS; + double p10 = powDi(10., e2); + x *= p10; + double P10 = powDi(10., e10 - e2); + x *= P10; + /*-- p10 * P10 = 10 ^ e10 */ + if (do_round) { + x += 0.5; + } + x = Math.floor(x) / p10; + return (sgn * x / P10); + } + } + + // GNUR from fround.c + + public static double rint(double x) { + double tmp, sgn = 1.0; + long ltmp; + + if (x != x) { /* NaN */ + return x; + } + + if (x < 0.0) { + x = -x; + sgn = -1.0; + } + + if (x < (double) LONG_MAX) { + ltmp = (long) (x + 0.5); + /* implement round to even */ + if (fabs(x + 0.5 - ltmp) < 10 * DBL_EPSILON && (ltmp % 2 == 1)) + ltmp--; + tmp = ltmp; + } else { + /* ignore round to even: too small a point to bother */ + tmp = Math.floor(x + 0.5); + } + return sgn * tmp; + } + public static double fsign(double x, double y) { if (Double.isNaN(x) || Double.isNaN(y)) { return x + y; @@ -85,6 +190,34 @@ public final class RMath { return tmp - Math.floor(tmp / b) * b; } + public static double sinpi(double x) { + double norm = x % 2d; + if (norm == 0d || norm == 1d || norm == -1d) { + return 0d; + } + if (norm == -1.5d || norm == 0.5d) { + return 1d; + } + if (norm == -0.5d || norm == 1.5d) { + return -1d; + } + return Math.sin(norm * Math.PI); + } + + public static double cospi(double x) { + double norm = x % 2d; + if (norm == 0d) { + return 1d; + } + if (norm == -1d || norm == 1d) { + return -1d; + } + if (norm == -1.5d || norm == -0.5d || norm == 0.5d || norm == 1.5d) { + return 0d; + } + return Math.cos(norm * Math.PI); + } + public static double tanpi(double x) { if (Double.isNaN(x)) { return x; diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Logis.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Logis.java index 3bb60b17a6..c1ed6b93d2 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Logis.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Logis.java @@ -5,7 +5,7 @@ * * Copyright (C) 1998 Ross Ihaka * Copyright (c) 1998--2008, The R Core Team - * Copyright (c) 2016, 2017, Oracle and/or its affiliates + * Copyright (c) 2016, 2018, Oracle and/or its affiliates * * All rights reserved. */ @@ -112,7 +112,7 @@ public final class Logis { } } - private static double log1pexp(double x) { + public static double log1pexp(double x) { if (x <= 18.) { return RMath.log1p(Math.exp(x)); } diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Pnorm.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Pnorm.java index c8ae38cd2b..b75a57edee 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Pnorm.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Pnorm.java @@ -78,6 +78,13 @@ public final class Pnorm implements Function3_2 { public static final class PnormBoth { + public static void evaluate(double x, double[] cum, double[] ccum, boolean lowerTail, boolean logP) { + PnormBoth pnormBoth = new PnormBoth(cum[0]); + pnormBoth.pnormBoth(x, lowerTail, logP); + cum[0] = pnormBoth.cum; + ccum[0] = pnormBoth.ccum; + } + private double cum; private double ccum; diff --git a/com.oracle.truffle.r.test.native/packages/testrffi/testrffi/src/testrffi.c b/com.oracle.truffle.r.test.native/packages/testrffi/testrffi/src/testrffi.c index b880edfeff..63a67dd670 100644 --- a/com.oracle.truffle.r.test.native/packages/testrffi/testrffi/src/testrffi.c +++ b/com.oracle.truffle.r.test.native/packages/testrffi/testrffi/src/testrffi.c @@ -590,106 +590,165 @@ SEXP test_RfEvalWithPromiseInPairList() { return result; } +//typedef double Rf_cospi(double a); +//typedef double Rf_sinpi(double a); +//typedef double Rf_tanpi(double a); SEXP test_RfRandomFunctions() { - SEXP v; - PROTECT(v = allocVector(REALSXP, 95)); + SEXP v, vNames; + int vLen = 128; + PROTECT(v = allocVector(VECSXP, vLen)); + PROTECT(vNames = allocVector(STRSXP, vLen)); int n = 0; - REAL(v)[n++] = Rf_dunif(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_qunif(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_punif(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_runif(0, 1); - REAL(v)[n++] = Rf_dchisq(0, 1, FALSE); - REAL(v)[n++] = Rf_pchisq(0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qchisq(0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rchisq(1); - REAL(v)[n++] = Rf_dnchisq(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pnchisq(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qnchisq(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rnchisq(0, 1); - REAL(v)[n++] = Rf_dnorm4(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pnorm5(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qnorm5(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rnorm(1, 0); - REAL(v)[n++] = Rf_dlnorm(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_plnorm(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qlnorm(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rlnorm(1, 0); - REAL(v)[n++] = Rf_dgamma(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pgamma(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qgamma(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rgamma(1, 0); - REAL(v)[n++] = Rf_dbeta(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pbeta(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qbeta(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rbeta(1, 0); - REAL(v)[n++] = Rf_df(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pf(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qf(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rf(1, 0); - REAL(v)[n++] = Rf_dt(1, 0, FALSE); - REAL(v)[n++] = Rf_pt(1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_qt(1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_rt(1); - REAL(v)[n++] = Rf_dbinom(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pbinom(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qbinom(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rbinom(1, 0); - REAL(v)[n++] = Rf_dcauchy(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pcauchy(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qcauchy(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rcauchy(1, 0); - REAL(v)[n++] = Rf_dexp(1, 0, FALSE); - REAL(v)[n++] = Rf_pexp(1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_qexp(1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_rexp(1); - REAL(v)[n++] = Rf_dgeom(1, 0, FALSE); - REAL(v)[n++] = Rf_pgeom(1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_qgeom(1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_rgeom(1); - REAL(v)[n++] = Rf_dhyper(1, 0, 1, 0, FALSE); - REAL(v)[n++] = Rf_phyper(1, 0, 1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_qhyper(1, 0, 1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_rhyper(1, 0, 1); - REAL(v)[n++] = Rf_dnbinom(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pnbinom(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qnbinom(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rnbinom(1, 0); - REAL(v)[n++] = Rf_dnbinom_mu(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pnbinom_mu(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qnbinom_mu(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rnbinom_mu(1, 0); - REAL(v)[n++] = Rf_dpois(1, 0, FALSE); - REAL(v)[n++] = Rf_ppois(1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_qpois(1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_rpois(1); - REAL(v)[n++] = Rf_dweibull(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pweibull(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qweibull(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rweibull(1, 0); - REAL(v)[n++] = Rf_dlogis(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_plogis(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qlogis(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rlogis(1, 0); - REAL(v)[n++] = Rf_dnbeta(1, 0, 1, 0, FALSE); - REAL(v)[n++] = Rf_pnbeta(1, 0, 1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_qnbeta(1, 0, 1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_dnf(1, 0, 1, 0, FALSE); - REAL(v)[n++] = Rf_pnf(1, 0, 1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_qnf(1, 0, 1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_dnt(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pnt(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qnt(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_ptukey(1, 0, 1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_qtukey(1, 0, 1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_dwilcox(1, 0, 1, FALSE); - REAL(v)[n++] = Rf_pwilcox(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_qwilcox(1, 0, 1, TRUE, FALSE); - REAL(v)[n++] = Rf_rwilcox(1, 0); - REAL(v)[n++] = Rf_dsignrank(1, 0, FALSE); - REAL(v)[n++] = Rf_psignrank(1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_qsignrank(1, 0, TRUE, FALSE); - REAL(v)[n++] = Rf_rsignrank(1); - UNPROTECT(1); + SET_STRING_ELT(vNames, n, mkChar("Rf_dunif")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dunif(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qunif")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qunif(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_punif")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_punif(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_runif")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_runif(0, 1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dchisq")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dchisq(0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pchisq")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pchisq(0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qchisq")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qchisq(0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rchisq")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rchisq(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dnchisq")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dnchisq(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pnchisq")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pnchisq(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qnchisq")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qnchisq(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rnchisq")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rnchisq(0, 1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dnorm4")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dnorm4(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pnorm5")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pnorm5(1, 0, 1, TRUE, FALSE))); + + double cum = 1.0; + double ccum = 0.0; + Rf_pnorm_both(2, &cum, &ccum, TRUE, FALSE); + SET_STRING_ELT(vNames, n, mkChar("Rf_pnorm_both-arg-cum")); SET_VECTOR_ELT(v, n++, ScalarReal(cum)); // Index 15 + SET_STRING_ELT(vNames, n, mkChar("Rf_pnorm_both-arg-ccum")); SET_VECTOR_ELT(v, n++, ScalarReal(ccum)); + + SET_STRING_ELT(vNames, n, mkChar("Rf_qnorm5")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qnorm5(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rnorm")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rnorm(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dlnorm")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dlnorm(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_plnorm")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_plnorm(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qlnorm")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qlnorm(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rlnorm")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rlnorm(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dgamma")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dgamma(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pgamma")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pgamma(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qgamma")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qgamma(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rgamma")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rgamma(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_log1pmx")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_log1pmx(42))); + SET_STRING_ELT(vNames, n, mkChar("Rf_log1pexp")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_log1pexp(6))); + SET_STRING_ELT(vNames, n, mkChar("Rf_lgamma1p")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_lgamma1p(42))); + SET_STRING_ELT(vNames, n, mkChar("Rf_logspace_add")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_logspace_add(42, 6))); + SET_STRING_ELT(vNames, n, mkChar("Rf_logspace_sub")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_logspace_sub(42, 6))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dbeta")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dbeta(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pbeta")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pbeta(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qbeta")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qbeta(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rbeta")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rbeta(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_df")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_df(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pf")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pf(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qf")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qf(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rf")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rf(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dt")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dt(1, 0, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pt")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pt(1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qt")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qt(1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rt")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rt(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dbinom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dbinom(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pbinom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pbinom(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qbinom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qbinom(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rbinom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rbinom(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dcauchy")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dcauchy(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pcauchy")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pcauchy(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qcauchy")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qcauchy(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rcauchy")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rcauchy(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dexp")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dexp(1, 0, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pexp")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pexp(1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qexp")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qexp(1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rexp")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rexp(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dgeom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dgeom(1, 0, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pgeom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pgeom(1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qgeom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qgeom(1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rgeom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rgeom(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dhyper")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dhyper(1, 0, 1, 0, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_phyper")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_phyper(1, 0, 1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qhyper")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qhyper(1, 0, 1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rhyper")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rhyper(1, 0, 1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dnbinom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dnbinom(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pnbinom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pnbinom(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qnbinom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qnbinom(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rnbinom")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rnbinom(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dnbinom_mu")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dnbinom_mu(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pnbinom_mu")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pnbinom_mu(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qnbinom_mu")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qnbinom_mu(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rnbinom_mu")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rnbinom_mu(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dpois")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dpois(1, 0, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_ppois")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_ppois(1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qpois")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qpois(1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rpois")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rpois(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dweibull")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dweibull(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pweibull")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pweibull(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qweibull")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qweibull(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rweibull")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rweibull(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dlogis")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dlogis(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_plogis")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_plogis(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qlogis")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qlogis(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rlogis")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rlogis(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dnbeta")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dnbeta(1, 0, 1, 0, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pnbeta")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pnbeta(1, 0, 1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qnbeta")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qnbeta(1, 0, 1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dnf")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dnf(1, 0, 1, 0, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pnf")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pnf(1, 0, 1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qnf")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qnf(1, 0, 1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dnt")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dnt(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pnt")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pnt(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qnt")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qnt(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_ptukey")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_ptukey(1, 0, 1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qtukey")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qtukey(1, 0, 1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dwilcox")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dwilcox(1, 0, 1, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pwilcox")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pwilcox(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qwilcox")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qwilcox(1, 0, 1, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rwilcox")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rwilcox(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_dsignrank")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_dsignrank(1, 0, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_psignrank")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_psignrank(1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_qsignrank")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_qsignrank(1, 0, TRUE, FALSE))); + SET_STRING_ELT(vNames, n, mkChar("Rf_rsignrank")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_rsignrank(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_gammafn")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_gammafn(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_lgammafn")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_lgammafn(1))); + + int signRet = 0; + SET_STRING_ELT(vNames, n, mkChar("Rf_lgammafn_sign")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_lgammafn_sign(1, &signRet))); + SET_STRING_ELT(vNames, n, mkChar("Rf_lgammafn_sign-arg-signRet")); SET_VECTOR_ELT(v, n++, ScalarReal(signRet)); + + double ans = 0.0; + int nz = 0; + int ierr = 0; + Rf_dpsifn(/*x*/1, /*n>=0*/2, /*kode>=1&&<=2*/1, /*m==1*/ 1, &ans, &nz, &ierr); + SET_STRING_ELT(vNames, n, mkChar("Rf_dpsinfn-arg-ans")); SET_VECTOR_ELT(v, n++, ScalarReal(ans)); + SET_STRING_ELT(vNames, n, mkChar("Rf_dpsinfn-arg-nz")); SET_VECTOR_ELT(v, n++, ScalarReal(nz)); + SET_STRING_ELT(vNames, n, mkChar("Rf_dpsinfn-arg-ierr")); SET_VECTOR_ELT(v, n++, ScalarReal(ierr)); + + SET_STRING_ELT(vNames, n, mkChar("Rf_psigamma")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_psigamma(1, 0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_digamma")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_digamma(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_trigamma")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_trigamma(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_tetragamma")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_tetragamma(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_pentagamma")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_pentagamma(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_beta")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_beta(1, 1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_lbeta")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_lbeta(1, 1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_choose")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_choose(1, 1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_lchoose")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_lchoose(1, 1))); + + double nb[3] = { 0.0, 0.0, 0.0 }; // work-array of size floor(second-param) + 1 + SET_STRING_ELT(vNames, n, mkChar("Rf_bessel_i")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_bessel_i(42.0, 2.0, 1.0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_bessel_i_ex")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_bessel_i_ex(1.0, 2.0, 1.0, nb))); + SET_STRING_ELT(vNames, n, mkChar("Rf_bessel_j")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_bessel_j(42.0, 2.0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_bessel_j_ex")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_bessel_j_ex(1.0, 2.0, nb))); + SET_STRING_ELT(vNames, n, mkChar("Rf_bessel_k")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_bessel_k(42.0, 2.0, 1.0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_bessel_k_ex")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_bessel_k_ex(1.0, 2.0, 1.0, nb))); + SET_STRING_ELT(vNames, n, mkChar("Rf_bessel_y")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_bessel_y(42.0, 2.0))); + SET_STRING_ELT(vNames, n, mkChar("Rf_bessel_y_ex")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_bessel_y_ex(1.0, 2.0, nb))); + +// SET_STRING_ELT(vNames, n, mkChar("Rf_cospi")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_cospi(1))); +// SET_STRING_ELT(vNames, n, mkChar("Rf_sinpi")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_sinpi(1))); +// SET_STRING_ELT(vNames, n, mkChar("Rf_tanpi")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_tanpi(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_sign")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_sign(1))); + SET_STRING_ELT(vNames, n, mkChar("Rf_fprec")); SET_VECTOR_ELT(v, n++, ScalarReal(Rf_fprec(1,2))); + setAttrib(v, R_NamesSymbol, vNames); +// printf("n=%d, vLen=%d\n", n, vLen); + UNPROTECT(2); return v; } diff --git a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test index aa4e0e0934..a589460550 100644 --- a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test +++ b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test @@ -11189,23 +11189,51 @@ character(0) #.Internal(bcVersion()) [1] 10 -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI1#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI1# #argv <- list(FALSE, FALSE, 1); .Internal(besselI(argv[[1]], argv[[2]], argv[[3]])) [1] 1 -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI2#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI1# +#besselI(1,c(NA,1)) +[1] NA 0.5651591 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI1# +#besselI(c(1,2),1) +[1] 0.5651591 1.5906369 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI1# +#besselI(c(1,2,3),c(1,2)) +[1] 0.5651591 0.6889484 3.9533702 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI1# +#besselI(c(1,2,3),c(1,2),c(3,4,5,6,7,8)) +[1] 0.20791042 0.09323903 0.19682671 0.04993878 0.21526929 0.11178255 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI1# +#besselI(c(1,2,3),c(1,2),c(3,4,NA,6,7,8)) +[1] 0.20791042 0.09323903 NA 0.04993878 0.21526929 0.11178255 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI1# +#besselI(c(1,2,3),c(NA,2)) +[1] NA 0.6889484 NA + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI1# +#besselI(c(1,NA),1) +[1] 0.5651591 NA + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselI.testbesselI2# #argv <- list(logical(0), logical(0), 1); .Internal(besselI(argv[[1]], argv[[2]], argv[[3]])) numeric(0) -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ1#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ1# #argv <- list(logical(0), logical(0)); .Internal(besselJ(argv[[1]], argv[[2]])) numeric(0) -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ2#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ2# #argv <- list(FALSE, FALSE); .Internal(besselJ(argv[[1]], argv[[2]])) [1] 1 -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ3#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ3# #argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 2.5); .Internal(besselJ(argv[[1]], argv[[2]])) [1] 4.724426e-17 2.672539e-16 1.511816e-15 8.552124e-15 4.837812e-14 [6] 2.736680e-13 1.548100e-12 8.757375e-12 4.953919e-11 2.802360e-10 @@ -11215,15 +11243,45 @@ numeric(0) [26] -8.858053e-02 -9.352409e-02 -4.969565e-02 4.984926e-02 -2.597979e-03 [31] 3.880718e-03 -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK1#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ3# +#besselJ(1,c(NA,1)) +[1] NA 0.4400506 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ3# +#besselJ(c(1,2),1) +[1] 0.4400506 0.5767248 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ3# +#besselJ(c(1,2,3),c(1,2)) +[1] 0.4400506 0.3528340 0.3390590 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ3# +#besselJ(c(1,2,3),c(1,2),c(3,4,5,6,7,8)) +Error in besselJ(c(1, 2, 3), c(1, 2), c(3, 4, 5, 6, 7, 8)) : + unused argument (c(3, 4, 5, 6, 7, 8)) + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ3# +#besselJ(c(1,2,3),c(1,2),c(3,4,NA,6,7,8)) +Error in besselJ(c(1, 2, 3), c(1, 2), c(3, 4, NA, 6, 7, 8)) : + unused argument (c(3, 4, NA, 6, 7, 8)) + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ3# +#besselJ(c(1,2,3),c(NA,2)) +[1] NA 0.352834 NA + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselJ.testbesselJ3# +#besselJ(c(1,NA),1) +[1] 0.4400506 NA + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK1# #argv <- list(FALSE, FALSE, 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]])) [1] Inf -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK2#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK2# #argv <- list(logical(0), logical(0), 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]])) numeric(0) -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK3#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK3# #argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 3, 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]])) [1] 9.223372e+18 1.152922e+18 1.441152e+17 1.801440e+16 2.251800e+15 [6] 2.814750e+14 3.518437e+13 4.398047e+12 5.497558e+11 6.871947e+10 @@ -11233,7 +11291,7 @@ numeric(0) [26] 3.209956e-15 2.688919e-29 2.948133e-57 5.271814e-113 2.445443e-224 [31] 0.000000e+00 -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK4#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK4# #argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 3.5, 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]])) [1] 2.219478e+22 1.961760e+21 1.733967e+20 1.532625e+19 1.354662e+18 [6] 1.197363e+17 1.058330e+16 9.354401e+14 8.268201e+13 7.308126e+12 @@ -11243,7 +11301,35 @@ numeric(0) [26] 3.374310e-15 2.757500e-29 2.985649e-57 5.305318e-113 2.453209e-224 [31] 0.000000e+00 -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY1#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK4# +#besselK(1,c(NA,1)) +[1] NA 0.6019072 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK4# +#besselK(c(1,2),1) +[1] 0.6019072 0.1398659 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK4# +#besselK(c(1,2,3),c(1,2)) +[1] 0.60190723 0.25375975 0.04015643 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK4# +#besselK(c(1,2,3),c(1,2),c(3,4,5,6,7,8)) +[1] 1.6361535 1.8750451 0.8065635 4.4167701 1.0334768 1.2354706 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK4# +#besselK(c(1,2,3),c(1,2),c(3,4,NA,6,7,8)) +[1] 1.636153 1.875045 NA 4.416770 1.033477 1.235471 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK4# +#besselK(c(1,2,3),c(NA,2)) +[1] NA 0.2537598 NA + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselK.testbesselK4# +#besselK(c(1,NA),1) +[1] 0.6019072 NA + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY1# #argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 20.5); .Internal(besselY(argv[[1]], argv[[2]])) [1] -6.747747e+146 -4.550341e+140 -3.068520e+134 -2.069255e+128 -1.395401e+122 [6] -9.409884e+115 -6.345551e+109 -4.279120e+103 -2.885623e+97 -1.945918e+91 @@ -11253,7 +11339,7 @@ numeric(0) [26] -7.076470e-02 2.381079e-02 4.744158e-02 -3.516090e-02 3.336562e-02 [31] -2.491015e-02 -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY2#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY2# #argv <- list(2, c(3, 8.94, 14.88, 20.82, 26.76, 32.7, 38.64, 44.58, 50.52, 56.46, 62.4, 68.34, 74.28, 80.22, 86.16, 92.1, 98.04, 103.98, 109.92, 115.86, 121.8, 127.74, 133.68, 139.62, 145.56, 151.5, 157.44, 163.38, 169.32, 175.26, 181.2, 187.14, 193.08, 199.02, 204.96, 210.9, 216.84, 222.78, 228.72, 234.66, 240.6, 246.54, 252.48, 258.42, 264.36, 270.3, 276.24, 282.18, 288.12, 294.06, 300)); .Internal(besselY(argv[[1]], argv[[2]])) [1] -1.127784e+00 -1.282008e+04 -2.165098e+10 -4.733004e+17 -6.084569e+25 [6] -3.046226e+34 -4.601250e+43 -1.761838e+53 -1.506980e+63 -2.615568e+73 @@ -11268,7 +11354,7 @@ numeric(0) [51] -Inf There were 22 warnings (use warnings() to see them) -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY3#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY3# #argv <- list(c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10), -0.2); .Internal(besselY(argv[[1]], argv[[2]])) [1] -Inf -1.129937637 -0.690945975 -0.435238869 -0.251890636 [6] -0.108032918 0.010318976 0.110293104 0.195945764 0.269765825 @@ -11292,10 +11378,40 @@ There were 22 warnings (use warnings() to see them) [96] 0.102417825 0.077752074 0.052569412 0.027122694 0.001664807 [101] -0.023553761 -##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY4#Ignored.Unimplemented# +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY4# #argv <- list(logical(0), logical(0)); .Internal(besselY(argv[[1]], argv[[2]])) numeric(0) +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY4# +#besselY(1,c(NA,1)) +[1] NA -0.7812128 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY4# +#besselY(c(1,2),1) +[1] -0.7812128 -0.1070324 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY4# +#besselY(c(1,2,3),c(1,2)) +[1] -0.7812128 -0.6174081 0.3246744 + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY4# +#besselY(c(1,2,3),c(1,2),c(3,4,5,6,7,8)) +Error in besselY(c(1, 2, 3), c(1, 2), c(3, 4, 5, 6, 7, 8)) : + unused argument (c(3, 4, 5, 6, 7, 8)) + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY4# +#besselY(c(1,2,3),c(1,2),c(3,4,NA,6,7,8)) +Error in besselY(c(1, 2, 3), c(1, 2), c(3, 4, NA, 6, 7, 8)) : + unused argument (c(3, 4, NA, 6, 7, 8)) + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY4# +#besselY(c(1,2,3),c(NA,2)) +[1] NA -0.6174081 NA + +##com.oracle.truffle.r.test.builtins.TestBuiltin_besselY.testbesselY4# +#besselY(c(1,NA),1) +[1] -0.7812128 NA + ##com.oracle.truffle.r.test.builtins.TestBuiltin_beta.testbeta1#Ignored.Unimplemented# #argv <- list(FALSE, FALSE); .Internal(beta(argv[[1]], argv[[2]])) [1] Inf diff --git a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselI.java b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselI.java index 78c36a510b..c26fa879f5 100644 --- a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselI.java +++ b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselI.java @@ -4,7 +4,7 @@ * http://www.gnu.org/licenses/gpl-2.0.html * * Copyright (c) 2014, Purdue University - * Copyright (c) 2014, 2017, Oracle and/or its affiliates + * Copyright (c) 2014, 2018, Oracle and/or its affiliates * * All rights reserved. */ @@ -19,13 +19,18 @@ public class TestBuiltin_besselI extends TestBase { @Test public void testbesselI1() { - // FIXME RInternalError: not implemented: .Internal besselI - assertEval(Ignored.Unimplemented, "argv <- list(FALSE, FALSE, 1); .Internal(besselI(argv[[1]], argv[[2]], argv[[3]]))"); + assertEval("besselI(1,c(NA,1))"); + assertEval("besselI(c(1,2),1)"); + assertEval("besselI(c(1,2,3),c(1,2))"); + assertEval("besselI(c(1,2,3),c(1,2),c(3,4,5,6,7,8))"); + assertEval("besselI(c(1,NA),1)"); + assertEval("besselI(c(1,2,3),c(NA,2))"); + assertEval("besselI(c(1,2,3),c(1,2),c(3,4,NA,6,7,8))"); + assertEval("argv <- list(FALSE, FALSE, 1); .Internal(besselI(argv[[1]], argv[[2]], argv[[3]]))"); } @Test public void testbesselI2() { - // FIXME RInternalError: not implemented: .Internal besselI - assertEval(Ignored.Unimplemented, "argv <- list(logical(0), logical(0), 1); .Internal(besselI(argv[[1]], argv[[2]], argv[[3]]))"); + assertEval("argv <- list(logical(0), logical(0), 1); .Internal(besselI(argv[[1]], argv[[2]], argv[[3]]))"); } } diff --git a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselJ.java b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselJ.java index 9906603888..4af22569f9 100644 --- a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselJ.java +++ b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselJ.java @@ -4,7 +4,7 @@ * http://www.gnu.org/licenses/gpl-2.0.html * * Copyright (c) 2014, Purdue University - * Copyright (c) 2014, 2017, Oracle and/or its affiliates + * Copyright (c) 2014, 2018, Oracle and/or its affiliates * * All rights reserved. */ @@ -19,20 +19,24 @@ public class TestBuiltin_besselJ extends TestBase { @Test public void testbesselJ1() { - // FIXME RInternalError: not implemented: .Internal besselJ - assertEval(Ignored.Unimplemented, "argv <- list(logical(0), logical(0)); .Internal(besselJ(argv[[1]], argv[[2]]))"); + assertEval("argv <- list(logical(0), logical(0)); .Internal(besselJ(argv[[1]], argv[[2]]))"); } @Test public void testbesselJ2() { - // FIXME RInternalError: not implemented: .Internal besselJ - assertEval(Ignored.Unimplemented, "argv <- list(FALSE, FALSE); .Internal(besselJ(argv[[1]], argv[[2]]))"); + assertEval("argv <- list(FALSE, FALSE); .Internal(besselJ(argv[[1]], argv[[2]]))"); } @Test public void testbesselJ3() { - // FIXME RInternalError: not implemented: .Internal besselJ - assertEval(Ignored.Unimplemented, - "argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 2.5); .Internal(besselJ(argv[[1]], argv[[2]]))"); + assertEval("argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 2.5); .Internal(besselJ(argv[[1]], argv[[2]]))"); + assertEval("besselJ(1,c(NA,1))"); + assertEval("besselJ(c(1,2),1)"); + assertEval("besselJ(c(1,2,3),c(1,2))"); + assertEval("besselJ(c(1,2,3),c(1,2),c(3,4,5,6,7,8))"); + assertEval("besselJ(c(1,NA),1)"); + assertEval("besselJ(c(1,2,3),c(NA,2))"); + assertEval("besselJ(c(1,2,3),c(1,2),c(3,4,NA,6,7,8))"); + } } diff --git a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselK.java b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselK.java index 100011e365..36ed394ce3 100644 --- a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselK.java +++ b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselK.java @@ -4,7 +4,7 @@ * http://www.gnu.org/licenses/gpl-2.0.html * * Copyright (c) 2014, Purdue University - * Copyright (c) 2014, 2017, Oracle and/or its affiliates + * Copyright (c) 2014, 2018, Oracle and/or its affiliates * * All rights reserved. */ @@ -19,27 +19,28 @@ public class TestBuiltin_besselK extends TestBase { @Test public void testbesselK1() { - // FIXME RInternalError: not implemented: .Internal besselK - assertEval(Ignored.Unimplemented, "argv <- list(FALSE, FALSE, 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]]))"); + assertEval("argv <- list(FALSE, FALSE, 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]]))"); } @Test public void testbesselK2() { - // FIXME RInternalError: not implemented: .Internal besselK - assertEval(Ignored.Unimplemented, "argv <- list(logical(0), logical(0), 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]]))"); + assertEval("argv <- list(logical(0), logical(0), 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]]))"); } @Test public void testbesselK3() { - // FIXME RInternalError: not implemented: .Internal besselK - assertEval(Ignored.Unimplemented, - "argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 3, 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]]))"); + assertEval("argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 3, 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]]))"); } @Test public void testbesselK4() { - // FIXME RInternalError: not implemented: .Internal besselK - assertEval(Ignored.Unimplemented, - "argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 3.5, 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]]))"); + assertEval("argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 3.5, 1); .Internal(besselK(argv[[1]], argv[[2]], argv[[3]]))"); + assertEval("besselK(1,c(NA,1))"); + assertEval("besselK(c(1,2),1)"); + assertEval("besselK(c(1,2,3),c(1,2))"); + assertEval("besselK(c(1,2,3),c(1,2),c(3,4,5,6,7,8))"); + assertEval("besselK(c(1,NA),1)"); + assertEval("besselK(c(1,2,3),c(NA,2))"); + assertEval("besselK(c(1,2,3),c(1,2),c(3,4,NA,6,7,8))"); } } diff --git a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselY.java b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselY.java index 4d751a26ec..91bcb680bc 100644 --- a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselY.java +++ b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/builtins/TestBuiltin_besselY.java @@ -4,7 +4,7 @@ * http://www.gnu.org/licenses/gpl-2.0.html * * Copyright (c) 2014, Purdue University - * Copyright (c) 2014, 2017, Oracle and/or its affiliates + * Copyright (c) 2014, 2018, Oracle and/or its affiliates * * All rights reserved. */ @@ -19,28 +19,28 @@ public class TestBuiltin_besselY extends TestBase { @Test public void testbesselY1() { - // FIXME RInternalError: not implemented: .Internal besselY - assertEval(Ignored.Unimplemented, - "argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 20.5); .Internal(besselY(argv[[1]], argv[[2]]))"); + assertEval("argv <- list(c(9.5367431640625e-07, 1.9073486328125e-06, 3.814697265625e-06, 7.62939453125e-06, 1.52587890625e-05, 3.0517578125e-05, 6.103515625e-05, 0.0001220703125, 0.000244140625, 0.00048828125, 0.0009765625, 0.001953125, 0.00390625, 0.0078125, 0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024), 20.5); .Internal(besselY(argv[[1]], argv[[2]]))"); } @Test public void testbesselY2() { - // FIXME RInternalError: not implemented: .Internal besselY - assertEval(Ignored.Unimplemented, - "argv <- list(2, c(3, 8.94, 14.88, 20.82, 26.76, 32.7, 38.64, 44.58, 50.52, 56.46, 62.4, 68.34, 74.28, 80.22, 86.16, 92.1, 98.04, 103.98, 109.92, 115.86, 121.8, 127.74, 133.68, 139.62, 145.56, 151.5, 157.44, 163.38, 169.32, 175.26, 181.2, 187.14, 193.08, 199.02, 204.96, 210.9, 216.84, 222.78, 228.72, 234.66, 240.6, 246.54, 252.48, 258.42, 264.36, 270.3, 276.24, 282.18, 288.12, 294.06, 300)); .Internal(besselY(argv[[1]], argv[[2]]))"); + assertEval("argv <- list(2, c(3, 8.94, 14.88, 20.82, 26.76, 32.7, 38.64, 44.58, 50.52, 56.46, 62.4, 68.34, 74.28, 80.22, 86.16, 92.1, 98.04, 103.98, 109.92, 115.86, 121.8, 127.74, 133.68, 139.62, 145.56, 151.5, 157.44, 163.38, 169.32, 175.26, 181.2, 187.14, 193.08, 199.02, 204.96, 210.9, 216.84, 222.78, 228.72, 234.66, 240.6, 246.54, 252.48, 258.42, 264.36, 270.3, 276.24, 282.18, 288.12, 294.06, 300)); .Internal(besselY(argv[[1]], argv[[2]]))"); } @Test public void testbesselY3() { - // FIXME RInternalError: not implemented: .Internal besselY - assertEval(Ignored.Unimplemented, - "argv <- list(c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10), -0.2); .Internal(besselY(argv[[1]], argv[[2]]))"); + assertEval("argv <- list(c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10), -0.2); .Internal(besselY(argv[[1]], argv[[2]]))"); } @Test public void testbesselY4() { - // FIXME RInternalError: not implemented: .Internal besselY - assertEval(Ignored.Unimplemented, "argv <- list(logical(0), logical(0)); .Internal(besselY(argv[[1]], argv[[2]]))"); + assertEval("argv <- list(logical(0), logical(0)); .Internal(besselY(argv[[1]], argv[[2]]))"); + assertEval("besselY(1,c(NA,1))"); + assertEval("besselY(c(1,2),1)"); + assertEval("besselY(c(1,2,3),c(1,2))"); + assertEval("besselY(c(1,2,3),c(1,2),c(3,4,5,6,7,8))"); + assertEval("besselY(c(1,NA),1)"); + assertEval("besselY(c(1,2,3),c(NA,2))"); + assertEval("besselY(c(1,2,3),c(1,2),c(3,4,NA,6,7,8))"); } } diff --git a/mx.fastr/copyrights/gnu_r_ihaka_core.copyright.star.regex b/mx.fastr/copyrights/gnu_r_ihaka_core.copyright.star.regex index 8ede89dd35..5671a32698 100644 --- a/mx.fastr/copyrights/gnu_r_ihaka_core.copyright.star.regex +++ b/mx.fastr/copyrights/gnu_r_ihaka_core.copyright.star.regex @@ -1 +1 @@ -/\*\n \* This material is distributed under the GNU General Public License\n \* Version 2. You may review the terms of this license at\n \* http://www.gnu.org/licenses/gpl-2.0.html\n \*\n \* Copyright \(C\) 1998 Ross Ihaka\n \* Copyright \(c\) (?:[1-2][09][0-9][0-9]--?)?(?:[1-2][09])?[0-9][0-9], The R Core Team\n \* Copyright \(c\) (?:(20[0-9][0-9]), )?(20[0-9][0-9]), Oracle and/or its affiliates\n \*\n \* All rights reserved.\n \*/\n.* +/\*\n \* This material is distributed under the GNU General Public License\n \* Version 2. You may review the terms of this license at\n \* http://www.gnu.org/licenses/gpl-2.0.html\n \*\n \* Copyright \(C\) (?:[1-2][09][0-9][0-9]-?)?(?:[1-2][09][0-9][0-9]) Ross Ihaka\n \* Copyright \(c\) (?:[1-2][09][0-9][0-9]--?)?(?:[1-2][09])?[0-9]?[0-9], The R Core Team\n \* Copyright \(c\) (?:(20[0-9][0-9]), )?(20[0-9][0-9]), Oracle and/or its affiliates\n \*\n \* All rights reserved.\n \*/\n.* diff --git a/mx.fastr/copyrights/overrides b/mx.fastr/copyrights/overrides index 33d09ed2aa..da411e3341 100644 --- a/mx.fastr/copyrights/overrides +++ b/mx.fastr/copyrights/overrides @@ -208,6 +208,8 @@ com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/ffi/DLL.java,gnu_r com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/gnur/SA_TYPE.java,gnu_r.copyright com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/gnur/SEXPTYPE.java,gnu_r.copyright com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/Arithmetic.java,gnu_r_gentleman_ihaka.copyright +com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/BesselFunctions.java,gnu_r_ihaka_core.copyright +com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/Beta.java,gnu_r_ihaka_core.copyright com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/Choose.java,gnu_r.copyright com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Cauchy.java,gnu_r_ihaka_core.copyright com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/nmath/distr/Chisq.java,gnu_r_ihaka_core.copyright @@ -969,4 +971,4 @@ com.oracle.truffle.r.pkgs/rJava/src/java/RJavaImport.java,rjava.copyright com.oracle.truffle.r.pkgs/rJava/src/java/RJavaTools.java,rjava_romain.copyright com.oracle.truffle.r.pkgs/rJava/src/java/RectangularArrayBuilder.java,rjava.copyright com.oracle.truffle.r.pkgs/rJava/src/java/RectangularArrayExamples.java,rjava.copyright -com.oracle.truffle.r.pkgs/rJava/src/java/RectangularArraySummary.java,rjava.copyright \ No newline at end of file +com.oracle.truffle.r.pkgs/rJava/src/java/RectangularArraySummary.java,rjava.copyright -- GitLab