[libc-commits] [libc] [libc][NFC] Add sollya script to compute worst case range reduction. (PR #104803)

[via libc-commits]

https://github.com/lntue updated https://github.com/llvm/llvm-project/pull/104803

>From 6348e2c9ff2cb7ca55430f6b64e6cb07818feb03 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue at google.com>
Date: Mon, 19 Aug 2024 11:39:02 -0400
Subject: [PATCH 1/2] [libc] Add sollya script to compute wort case range
 reduction.

---
 libc/utils/mathtools/worst_case.sollya | 45 ++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)
 create mode 100644 libc/utils/mathtools/worst_case.sollya

diff --git a/libc/utils/mathtools/worst_case.sollya b/libc/utils/mathtools/worst_case.sollya
new file mode 100644
index 00000000000000..f420a2aa0cd520
--- /dev/null
+++ b/libc/utils/mathtools/worst_case.sollya
@@ -0,0 +1,45 @@
+// Implement WorstCase functions to compute the worst case for x mod C, with
+// the exponent of x ranges from emin to emax, and precision of x is p.
+// Adapted to Sollya from the Maple function in 
+//   J-M. Muller, "Elementary Functions", 3rd ed, Section 11.3.2.
+verbosity=0;
+procedure WorstCase(p, emin, emax, C, ndigits) {
+    epsilonmin = 12345.0;
+    Digits = ndigits;
+
+    powerofBoverC = 2^(emin - p) / C;
+    for e from emin - p + 1 to emax - p + 1 do {
+        powerofBoverC = 2 * powerofBoverC;
+        a = floor(powerofBoverC);
+        Plast = a;
+        r = round(1/(powerofBoverC - a), ndigits, RN);
+        a = floor(r);
+        Qlast = 1;
+        Q = a;
+        P = Plast * a + 1;
+        while (Q < 2^p - 1) do {
+            r = round(1/(r - a), ndigits, RN);
+            a = floor(r);
+            NewQ = Q * a + Qlast;
+            NewP = P * a + Plast;
+            Qlast = Q;
+            Plast = P;
+            Q = NewQ;
+            P = NewP;
+        };
+        epsilon = C * abs(Plast - Qlast * powerofBoverC);
+        if (epsilon < epsilonmin) then {
+            epsilonmin = epsilon;
+            numbermin = Qlast;
+            expmin = e;
+        };
+    };
+    display=decimal!;
+    print("numbermin : ", numbermin);
+    print("expmin    : ", expmin);
+    display=hexadecimal!;
+    print("Worst case: ", numbermin * 2^expmin);
+    display=decimal!;
+    ell = round(log2(epsilonmin), ndigits, RN);
+    print("numberofdigits : ", ell);
+};

>From 54b52d9200cd039caf4e48730d08b7ccef44dbcf Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue at google.com>
Date: Mon, 19 Aug 2024 15:09:18 -0400
Subject: [PATCH 2/2] Add examples.

---
 libc/utils/mathtools/worst_case.sollya | 35 ++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/libc/utils/mathtools/worst_case.sollya b/libc/utils/mathtools/worst_case.sollya
index f420a2aa0cd520..3a8d11b3da44d5 100644
--- a/libc/utils/mathtools/worst_case.sollya
+++ b/libc/utils/mathtools/worst_case.sollya
@@ -2,6 +2,41 @@
 // the exponent of x ranges from emin to emax, and precision of x is p.
 // Adapted to Sollya from the Maple function in 
 //   J-M. Muller, "Elementary Functions", 3rd ed, Section 11.3.2.
+//
+// Some examples:
+//
+// 1) Worst case for trig range reduction fast passes:
+//
+// Single precision
+// > WorstCase(24, -6, 32, pi/32, 128); 
+// numbermin :  10741887
+// expmin    :  7
+// Worst case:  0x1.47d0fep30
+// numberofdigits :  -32.888
+//
+// Double precision
+// > WorstCase(53, -8, 32, pi/128, 256);
+// numbermin :  6411027962775774
+// expmin    :  -53
+// Worst case:  0x1.6c6cbc45dc8dep-1
+// numberofdigits :  -66.4867
+//
+// 2) Worst case for exponential function range reduction:
+//
+// Single precision
+// > WorstCase(24, 1, 8, log(2), 128);
+// numbermin :  2907270
+// expmin    :  -22
+// Worst case:  0x1.62e43p-1
+// numberofdigits :  -28.9678
+//
+// Double precision
+// > WorstCase(53, 0, 11, log(2), 128);
+// numbermin :  7804143460206699
+// expmin    :  -51
+// Worst case:  0x1.bb9d3beb8c86bp1
+// numberofdigits :  -57.4931
+//
 verbosity=0;
 procedure WorstCase(p, emin, emax, C, ndigits) {
     epsilonmin = 12345.0;