1# -*- julia -*-
2#
3# Generate tgamma128.h, containing polynomials and constants used by
4# tgamma128.c.
5#
6# Copyright (c) 2006-2023, Arm Limited.
7# SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
8
9# This Julia program depends on the 'Remez' and 'SpecialFunctions'
10# library packages. To install them, run this at the interactive Julia
11# prompt:
12#
13#   import Pkg; Pkg.add(["Remez", "SpecialFunctions"])
14#
15# Tested on Julia 1.4.1 (Ubuntu 20.04) and 1.9.0 (22.04).
16
17import Printf
18import Remez
19import SpecialFunctions
20
21# Round a BigFloat to 128-bit long double and format it as a C99 hex
22# float literal.
23function quadhex(x)
24    sign = " "
25    if x < 0
26        sign = "-"
27        x = -x
28    end
29
30    exponent = BigInt(floor(log2(x)))
31    exponent = max(exponent, -16382)
32    @assert(exponent <= 16383) # else overflow
33
34    x /= BigFloat(2)^exponent
35    @assert(1 <= x < 2)
36    x *= BigFloat(2)^112
37    mantissa = BigInt(round(x))
38
39    mantstr = string(mantissa, base=16, pad=29)
40    return Printf.@sprintf("%s0x%s.%sp%+dL", sign, mantstr[1], mantstr[2:end],
41                           exponent)
42end
43
44# Round a BigFloat to 128-bit long double and return it still as a
45# BigFloat.
46function quadval(x, round=0)
47    sign = +1
48    if x.sign < 0
49        sign = -1
50        x = -x
51    end
52
53    exponent = BigInt(floor(log2(x)))
54    exponent = max(exponent, -16382)
55    @assert(exponent <= 16383) # else overflow
56
57    x /= BigFloat(2)^exponent
58    @assert(1 <= x < 2)
59    x *= BigFloat(2)^112
60    if round < 0
61        mantissa = floor(x)
62    elseif round > 0
63        mantissa = ceil(x)
64    else
65        mantissa = round(x)
66    end
67
68    return sign * mantissa * BigFloat(2)^(exponent - 112)
69end
70
71# Output an array of BigFloats as a C array declaration.
72function dumparray(a, name)
73    println("static const long double ", name, "[] = {")
74    for x in N
75        println("    ", quadhex(x), ",")
76    end
77    println("};")
78end
79
80print("/*
81 * Polynomial coefficients and other constants for tgamma128.c.
82 *
83 * Copyright (c) 2006,2009,2023 Arm Limited.
84 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
85 */
86")
87
88Base.MPFR.setprecision(512)
89
90e = exp(BigFloat(1))
91
92print("
93/* The largest positive value for which 128-bit tgamma does not overflow. */
94")
95lo = BigFloat("1000")
96hi = BigFloat("2000")
97while true
98    global lo
99    global hi
100    global max_x
101
102    mid = (lo + hi) / 2
103    if mid == lo || mid == hi
104        max_x = mid
105        break
106    end
107    if SpecialFunctions.logabsgamma(mid)[1] < 16384 * log(BigFloat(2))
108        lo = mid
109    else
110        hi = mid
111    end
112end
113max_x = quadval(max_x, -1)
114println("static const long double max_x = ", quadhex(max_x), ";")
115
116print("
117/* Coefficients of the polynomial used in the tgamma_large() subroutine */
118")
119N, D, E, X = Remez.ratfn_minimax(
120    x -> x==0 ? sqrt(BigFloat(2)*pi/e) :
121                exp(SpecialFunctions.logabsgamma(1/x)[1] +
122                    (1/x-0.5)*(1+log(x))),
123    (0, 1/BigFloat(8)),
124    24, 0,
125    (x, y) -> 1/y
126)
127dumparray(N, "coeffs_large")
128
129print("
130/* Coefficients of the polynomial used in the tgamma_tiny() subroutine */
131")
132N, D, E, X = Remez.ratfn_minimax(
133    x -> x==0 ? 1 : 1/(x*SpecialFunctions.gamma(x)),
134    (0, 1/BigFloat(32)),
135    13, 0,
136)
137dumparray(N, "coeffs_tiny")
138
139print("
140/* The location within the interval [1,2] where gamma has a minimum.
141 * Specified as the sum of two 128-bit values, for extra precision. */
142")
143lo = BigFloat("1.4")
144hi = BigFloat("1.5")
145while true
146    global lo
147    global hi
148    global min_x
149
150    mid = (lo + hi) / 2
151    if mid == lo || mid == hi
152        min_x = mid
153        break
154    end
155    if SpecialFunctions.digamma(mid) < 0
156        lo = mid
157    else
158        hi = mid
159    end
160end
161min_x_hi = quadval(min_x, -1)
162println("static const long double min_x_hi = ", quadhex(min_x_hi), ";")
163println("static const long double min_x_lo = ", quadhex(min_x - min_x_hi), ";")
164
165print("
166/* The actual minimum value that gamma takes at that location.
167 * Again specified as the sum of two 128-bit values. */
168")
169min_y = SpecialFunctions.gamma(min_x)
170min_y_hi = quadval(min_y, -1)
171println("static const long double min_y_hi = ", quadhex(min_y_hi), ";")
172println("static const long double min_y_lo = ", quadhex(min_y - min_y_hi), ";")
173
174function taylor_bodge(x)
175    # Taylor series generated by Wolfram Alpha for (gamma(min_x+x)-min_y)/x^2.
176    # Used in the Remez calls below for x values very near the origin, to avoid
177    # significance loss problems when trying to compute it directly via that
178    # formula (even in MPFR's extra precision).
179    return BigFloat("0.428486815855585429730209907810650582960483696962660010556335457558784421896667728014324097132413696263704801646004585959298743677879606168187061990204432200")+x*(-BigFloat("0.130704158939785761928008749242671025181542078105370084716141350308119418619652583986015464395882363802104154017741656168641240436089858504560718773026275797")+x*(BigFloat("0.160890753325112844190519489594363387594505844658437718135952967735294789599989664428071656484587979507034160383271974554122934842441540146372016567834062876")+x*(-BigFloat("0.092277030213334350126864106458600575084335085690780082222880945224248438672595248111704471182201673989215223667543694847795410779036800385804729955729659506"))))
180end
181
182print("
183/* Coefficients of the polynomial used in the tgamma_central() subroutine
184 * for computing gamma on the interval [1,min_x] */
185")
186N, D, E, X = Remez.ratfn_minimax(
187    x -> x < BigFloat(0x1p-64) ? taylor_bodge(-x) :
188        (SpecialFunctions.gamma(min_x - x) - min_y) / (x*x),
189    (0, min_x - 1),
190    31, 0,
191    (x, y) -> x^2,
192)
193dumparray(N, "coeffs_central_neg")
194
195print("
196/* Coefficients of the polynomial used in the tgamma_central() subroutine
197 * for computing gamma on the interval [min_x,2] */
198")
199N, D, E, X = Remez.ratfn_minimax(
200    x -> x < BigFloat(0x1p-64) ? taylor_bodge(x) :
201        (SpecialFunctions.gamma(min_x + x) - min_y) / (x*x),
202    (0, 2 - min_x),
203    28, 0,
204    (x, y) -> x^2,
205)
206dumparray(N, "coeffs_central_pos")
207
208print("
209/* 128-bit float value of pi, used by the sin_pi_x_over_pi subroutine
210 */
211")
212println("static const long double pi = ", quadhex(BigFloat(pi)), ";")
213