Inverse Harmonic and Inverse Digamma Functions - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: HP Prime Software Library (/forum-15.html) +--- Thread: Inverse Harmonic and Inverse Digamma Functions (/thread-15524.html) Inverse Harmonic and Inverse Digamma Functions - Albert Chan - 09-01-2020 09:06 PM I don't actually have HP Prime, below code tested with XCas invPsi_guess(k) based on inverse harmonic function f3, optimized for small k = e^x For small x, Ψ(x) ≈ -1/x. Flip the order, we have x ≈ -1/Ψ(x) Guess is improved by Halley's method, 2 times Code: ```invPsi_guess(k) := k - k/(24.*k*k+1.) + 0.5;  // k = e^x invPsi(x) := {  local g, f0, f1;  g := when(x < -1.8, -1./x, invPsi_guess(e^x));  g -= (f0 := Psi(g)-x) * (f1 := Psi(g,1)) / (f1*f1 - 0.5*f0*Psi(g,2));  g -= (f0 := Psi(g)-x) * (f1 := Psi(g,1)) / (f1*f1 - 0.5*f0*Psi(g,2)); }; invH(h) := invPsi(h-euler_gamma) - 1;``` Test functions for round-trip ... XCas> map(x -> Psi(invPsi(x)), range(-5,5)) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → [-5.0, -4.0, -3.0, -2.0, -1.0, 3.19839640883e-17, 1.0, 2.0, 3.0, 4.0] XCas> map(k -> invPsi(Psi(10.^k)), range(-5,5)) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → [1e-05, 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, 100.0, 1000.0, 10000.0] XCas> map(n -> invH(harmonic(n)), range(1,10)) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0] Update: replaced ifte by when, to avoided the ifte bug. RE: Inverse Harmonic and Inverse Digamma Functions - Albert Chan - 09-02-2020 11:03 AM I try to plot invPsi(x), but failed. Why ? XCas> plot(invPsi(x), x=-1 .. 1) "Ifte: Unable to check test Error: Bad Argument Value" RE: Inverse Harmonic and Inverse Digamma Functions - Albert Chan - 09-02-2020 07:48 PM Compare these expressions, which should be equivalent: XCas> subst(x>0, x=1) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → true XCas> subst(ifte(x>0, true, false), x=1) ﻿ ﻿ ﻿ ﻿ ﻿ → "Ifte: Unable to check test Error: Bad Argument Value" ifte() tried to evaluate the expression, before substitution. However, "x>0", cannot be evaluated. XCas> subst(ifte(x-x, true, false), x=1) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → false XCas> subst(0/x, x=0) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 0 Again, expressions were evaluated first, before substitution. This is branchless invPsi(x), to avoid above ifte() issues. This version is able to do "plot(invPsi(x), x=-1 .. 1)" Code: ```invPsi(x) := {  local g, f0, f1;  g := [x >= -1.8, x < -1.8] * [invPsi_guess(e^x), -1/(x+not(x))];  g -= (f0 := Psi(g)-x) * (f1 := Psi(g,1)) / (f1*f1 - 0.5*f0*Psi(g,2));  g -= (f0 := Psi(g)-x) * (f1 := Psi(g,1)) / (f1*f1 - 0.5*f0*Psi(g,2)); };``` Both "then" and "else" branch were evaluated, which is not ideal. Also, -1/x guess replaced by -1/(x+not(x)), to patch against divide-by-zero. XCas> subst(invPsi(x), x = 0) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 1.46163214497 XCas> subst(invPsi(x), x = -euler_gamma) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 1.0 RE: Inverse Harmonic and Inverse Digamma Functions - Albert Chan - 09-03-2020 02:05 PM Something unexpected when I try to test invPsi() accuracy. Note: Psi(1) = -euler_gamma = -0.57721566490153286 ... XCas> 1 - invPsi(-0.57721566490153286) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 5.55111512313e-16 XCas> 1 - invPsi(-euler_gamma) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 1.713074127e-13 ??? Above were done on XCas 1.5.0-63 (win32), with 12-digits display float format. Raising to 17-digits, expressions now agreed, error = 5.5511151231258e-16 Rounding constant with changing float display format is likely a XCas bug. This bug seems to affect only euler_gamma. Example, setting float display format to 6 digits: XCas> 3.1415926535897931 - float(pi) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 0.0 XCas> 0.57721566490153286 - float(euler_gamma) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 0.470199e-6 ??? RE: Inverse Harmonic and Inverse Digamma Functions - Albert Chan - 10-30-2020 09:57 PM (09-02-2020 11:03 AM)Albert Chan Wrote:  I try to plot invPsi(x), but failed. Why ? XCas> plot(invPsi(x), x=-1 .. 1) "Ifte: Unable to check test Error: Bad Argument Value" I found the fix ! Replace Ifte by when (piecewise also work)