Skip to content

Hypergeometric functions for combinations of p and q other than (p = q = 0, p = 1 and q = 0, p = 0 and q = 1, p = 2 and q = 0, p = q = 1, p = 0 and q = 2, p = 2 and q = 1, p = 3 and q = 2) don't match Mathematica's results #99

New issue
New issue
Closed
#100
Closed
Hypergeometric functions for combinations of p and q other than (p = q = 0, p = 1 and q = 0, p = 0 and q = 1, p = 2 and q = 0, p = q = 1, p = 0 and q = 2, p = 2 and q = 1, p = 3 and q = 2) don't match Mathematica's results#99
#100

Description

@avachon100510
avachon100510
opened on Jul 20, 2026
Image Image

As you can see from the two images, there's a mismatch. Using the Julia interpreter,

julia> using HypergeometricFunctions

julia> pFq((2.0, 3.0, 5.0, 7.0, 11.0), (2.0, 3.0, 5.0, 7.0), 2)
-57.81767556925831

It clearly doesn't match Mathematica's output:

In[50]:= HypergeometricPFQ[{2, 3, 5, 7, 11}, {2, 3, 5, 7}, 2]
Out[50]:= -1

I suspect it's a bug.
Thanks in advance

Activity

  1. TSGut commented on Jul 20, 2026

    @TSGut
    TSGut
    on Jul 20, 2026
    Contributor

    It's a cancellation bug typical for hypergeometric type functions. I will prep a PR.

  2. avachon100510 commented on Jul 21, 2026

    @avachon100510
    avachon100510
    on Jul 21, 2026
    Author

    Also, I've been trying to implement the Weniger transformation method in Ultra Fractal, but there's that black wedge that I can't seem to get rid of...

    Image

    Can you check the code if I did something wrong?

           ComplexArray workA = new ComplexArray(p)
       ComplexArray workB = new ComplexArray(q)

       int i = 0
       while i < p
         workA.m_Elements[i] = A.m_Elements[i]
         i = i + 1
       endwhile
       i = 0
       while i < q
         workB.m_Elements[i] = B.m_Elements[i]
         i = i + 1
       endwhile

       CancelCommonParameters(workA, p, workB, q)

       ComplexArray absalphaarray = new ComplexArray(p + 1)
       absalphaarray.m_Elements[0] = 0.0
       int n = 1
       while n < absalphaarray.getArrayLength()
         absalphaarray.m_Elements[n] = cabs(workA.m_Elements[n - 1])
         n = n + 1
       endwhile
       n = 1
       absaalphaarrayprod = 1
       while n < absalphaarray.getArrayLength()
         absaalphaarrayprod = absaalphaarrayprod*absalphaarray.m_Elements[n]
         n = n + 1
       endwhile
       complex gamma = HypPFQGetDynamicGamma(p, q, A, B, z)
       complex zeta = recip(z)
       int r = trunc(max(p, q + 1))
       ComplexArray bigN = new ComplexArray(r + 4)
       ComplexArray nextN = new ComplexArray(r + 4)
       bigN.m_Elements[0] = 0
       int k = 1
       while k < bigN.getArrayLength()
         bigN.m_Elements[k] = 0.0
         nextN.m_Elements[k] = 0.0
         k = k + 1
       endwhile
       ComplexArray bigD = new ComplexArray(r + 4)
       ComplexArray nextD = new ComplexArray(r + 4)
       bigD.m_Elements[0] = 0
       k = 0
       while k < bigD.getArrayLength()
         bigD.m_Elements[k] = 0.0
         nextD.m_Elements[k] = 0.0
         k = k + 1
       endwhile
       ComplexArray bigR = new ComplexArray(r + 4)
       ComplexArray nextR = new ComplexArray(r + 4)
       bigR.m_Elements[0] = 0
       k = 0
       while k < bigR.getArrayLength()
         bigR.m_Elements[k] = 0.0
         nextR.m_Elements[k] = 0.0
         k = k + 1
       endwhile
       complex prodBeta = gmullimit(workB, q)
       complex prodAlpha = gmullimit(workA, p)
       complex initialVal = (prodBeta*zeta)/prodAlpha/(gamma - 1)
       bigN.m_Elements[r + 3] = initialVal
       bigD.m_Elements[r + 3] = initialVal
       bigR.m_Elements[r + 3] = bigN.m_Elements[r + 3]/bigD.m_Elements[r + 3]
       complex err = 1
       int j = 1
       while j <= p
         err = err*(absalphaarray.m_Elements[j] + 1)
         j = j + 1
       endwhile
       ComplexArray P_hat = new ComplexArray(r + 4)
       k = 0
       while k < P_hat.getArrayLength()
         P_hat.m_Elements[k] = 0.0
         k = k + 1
       endwhile
       complex t = 1
       j = 1
       while j <= p
         t = t*(workA.m_Elements[j - 1] + 1)
         j = j + 1
       endwhile
       P_hat.m_Elements[1] = t
       ComplexArray bigQ = new ComplexArray(r + 3)
       t = 2
       j = 1
       while j <= q
         t = t*(workB.m_Elements[j - 1] + 1)
         j = j + 1
       endwhile
       bigQ.m_Elements[1] = t
       complex P_hat_R = gamma + 2
       complex QR = 1
       k = 0
       while k <= r + 2 || (k < kmax && cabs(bigR.m_Elements[r + 3] - bigR.m_Elements[r + 2]) > max(cabs(bigR.m_Elements[r + 3]), cabs(bigR.m_Elements[r + 2]))*8*2.22e-16)
         j = 1
         while j <= r + 2
           bigN.m_Elements[j] = bigN.m_Elements[j + 1]
           bigD.m_Elements[j] = bigD.m_Elements[j + 1]
           bigR.m_Elements[j] = bigR.m_Elements[j + 1]
           j = j + 1
         endwhile
         complex t1 = 0.0
         j = 0
         while j <= r
           t1 = t1 + bigQ.m_Elements[j + 1]*(gamma + 2*k - 2*j - 1)*bigN.m_Elements[r - j + 2]
           j = j + 1
         endwhile
         if k <= r
           j = 0
           while j <= k
             complex t2 = 1.0
             int i = 1
             while i <= q
               t2 = t2*(workB.m_Elements[i - 1] + j + 1)
               i = i + 1
             endwhile
             t2 = t2*(j + 2)
             t1 = t1 + nCr(k, j)*(-1)^(k - j)*t2
             j = j + 1
           endwhile
         endif
         t2 = 0.0
         t2 = t2 + P_hat.m_Elements[1]*(gamma + k - 1)*bigN.m_Elements[r + 2]
         j = 1
         while j <= r + 1
           t2 = t2 + P_hat.m_Elements[j + 1]*(bigN.m_Elements[r - j + 3] + (gamma + k - j - 1)*bigN.m_Elements[r - j + 2])
           j = j + 1
         endwhile
         bigN.m_Elements[r + 3] = zeta*t1 - t2
         t1 = 0.0
         j = 0
         while j <= r
           t1 = t1 + bigQ.m_Elements[j + 1]*(gamma + 2*k - 2*j - 1)*bigD.m_Elements[r - j + 2]
           j = j + 1
         endwhile
         t2 = 0.0
         t2 = t2 + P_hat.m_Elements[1]*(gamma + k - 1)*bigD.m_Elements[r + 2]
         j = 1
         while j <= r + 1
           t2 = t2 + P_hat.m_Elements[j + 1]*(bigD.m_Elements[r - j + 3] + (gamma + k - j - 1)*bigD.m_Elements[r - j + 2])
           j = j + 1
         endwhile
         bigD.m_Elements[r + 3] = zeta*t1 - t2
         bigR.m_Elements[r + 3] = bigN.m_Elements[r + 3]/bigD.m_Elements[r + 3]
         if inf_norm(P_hat.m_Elements[1]) < 2.22e-16*err
           return bigR.m_Elements[r + 3]
         endif

         bigN.m_Elements[r + 3] = bigN.m_Elements[r + 3]/P_hat.m_Elements[1]
         bigD.m_Elements[r + 3] = bigD.m_Elements[r + 3]/P_hat.m_Elements[1]
         k = k + 1
         err = 1
         j = 1
         while j <= p
           err = err*(absalphaarray.m_Elements[j] + 1)
           j = j + 1
         endwhile
         if k <= r + 1
           t2 = 1
           j = 1
           while j <= p
             t2 = t2*(workA.m_Elements[j - 1] + k + 1)
             j = j + 1
           endwhile
           t2 = t2/P_hat_R
           P_hat_R = P_hat_R*(((gamma + 2*k + 1)*(gamma + 2*k + 2))/(gamma + k + 1))
           t1 = k*((gamma + 2*k)*t2 - P_hat.m_Elements[1])
           j = 2
           while j <= k
             complex s = ((k - j + 1)*((gamma + 2*k - j + 1)*t1 + k*P_hat.m_Elements[j - 1]) - k*P_hat.m_Elements[j])/j
             P_hat.m_Elements[j - 1] = t2
             t2 = t1
             t1 = s
             j = j + 1
           endwhile
           P_hat.m_Elements[k] = t2
           P_hat.m_Elements[k + 1] = t1
           t2 = k + 2
           j = 1
           while j <= q
             t2 = t2*(workB.m_Elements[j - 1] + k + 1)
             j = j + 1
           endwhile
           QR = QR*(((gamma + 2*k - 2)*(gamma + 2*k - 1))/(gamma + k - 1))
           t2 = t2/QR
           t1 = k*((gamma + 2*k - 1)*t2 - bigQ.m_Elements[1])
           j = 2
           while j <= trunc(min(k, r))
             s = ((k - j + 1)*((gamma + 2*k - j)*t1 + k*bigQ.m_Elements[j - 1]) - k*bigQ.m_Elements[j])/j
             bigQ.m_Elements[j - 1] = t2
             t2 = t1
             t1 = s
             j = j + 1
           endwhile
           bigQ.m_Elements[trunc(min(k, r))] = t2
           bigQ.m_Elements[trunc(min(k, r)) + 1] = t1
         else
           t2 = gamma + k
           j = 1
           while j <= p
             t2 = t2*(workA.m_Elements[j - 1] + k + 1)
             j = j + 1
           endwhile
           t2 = t2/P_hat_R
           P_hat_R = P_hat_R*(((gamma + 2*k + 1)*(gamma + 2*k + 2))/((gamma + 2*k - r - 1)*(gamma + 2*k - r)))
           t1 = k*((gamma + 2*k)*t2 - (gamma + 2*k - 1 - r - 2)*P_hat.m_Elements[1])
           j = 2
           while j <= r + 1
             complex s = ((k - j + 1)*((gamma + 2*k - j + 1)*t1 - (r - j + 3)*k*P_hat.m_Elements[j - 1]) - (gamma + 2*k - j - r - 2)*k*P_hat.m_Elements[j])/j
             P_hat.m_Elements[j - 1] = t2
             t2 = t1
             t1 = s
             j = j + 1
           endwhile
           P_hat.m_Elements[r + 1] = t2
           P_hat.m_Elements[r + 2] = t1
           t2 = k + 2
           j = 1
           while j <= q
             t2 = t2*(workB.m_Elements[j - 1] + k + 1)
             j = j + 1
           endwhile
           QR = QR*(((gamma + 2*k - 2)*(gamma + 2*k - 1))/((gamma + 2*k - r - 3)*(gamma + 2*k - r - 2)))
           t2 = t2/QR
           t1 = k*((gamma + 2*k - 1)*t2 - (gamma + 2*k - 1 - r - 2)*bigQ.m_Elements[1])
           j = 2
           while j <= r
             s = ((k - j + 1)*((gamma + 2*k - j)*t1 - (r - j + 2)*k*bigQ.m_Elements[j - 1]) - (gamma + 2*k - j - r - 2)*k*bigQ.m_Elements[j])/j
             bigQ.m_Elements[j - 1] = t2
             t2 = t1
             t1 = s
             j = j + 1
           endwhile
           bigQ.m_Elements[r] = t2
           bigQ.m_Elements[r + 1] = t1
         endif
       endwhile
       if isfinite(bigR.m_Elements[r + 3])
         return bigR.m_Elements[r + 3]
       else
         return bigR.m_Elements[r + 2]
       endif

  3. avachon100510 commented on Jul 22, 2026

    @avachon100510
    avachon100510
    on Jul 22, 2026
    Author

    What causes this NaN region?

  4. avachon100510 commented on Jul 22, 2026

    @avachon100510
    avachon100510
    on Jul 22, 2026
    Author

    I figured it out, turns out it was floating-point overflow, so I added this:

    
         complex max_val = 1.0
         int check_idx = 1
         while check_idx <= r + 3
           complex cur_n = bigN.m_Elements[check_idx]
           complex cur_d = bigD.m_Elements[check_idx]
           if cabs(cur_n) > cabs(max_val)
             max_val = cur_n
           endif
           if cabs(cur_d) > cabs(max_val)
             max_val = cur_d
           endif
           check_idx = check_idx + 1
         endwhile

         float norm_mag = inf_norm(max_val)
         if norm_mag > 1e20 || norm_mag < 1e-20
           int norm_idx = 1
           while norm_idx <= r + 3
             bigN.m_Elements[norm_idx] = bigN.m_Elements[norm_idx]/max_val
             bigD.m_Elements[norm_idx] = bigD.m_Elements[norm_idx]/max_val
             norm_idx = norm_idx + 1
           endwhile
         endif

    It worked.

  5. MikaelSlevinsky
    closed this as completedin #100on Jul 28, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions